Noise handling demo

The noise framework employed by the simulator might also be useful on its own, as it allows to generate the exact same noises as in the simulator, and provides access to the exact noise definitions and PSDs.

This notebook shows the various noise types used in the simulator, demonstrates how to generate noise samples independent from the simulator, and how to obtain the PSDs for any noise type. It shows how to extract the default noise definitions from the simulator and plots those as examples.

import matplotlib_inline

%matplotlib inline
from matplotlib import pyplot as plt

plt.rcParams["figure.constrained_layout.use"] = True
plt.rcParams["figure.figsize"] = [9, 6]
plt.rcParams["figure.dpi"] = 50
plt.rcParams["font.size"] = 20
matplotlib_inline.backend_inline.set_matplotlib_formats("jpeg", dpi=50)
# reduces notebook file size

import numpy as np

from lisainstrument import Instrument, MosaID, SatID
from lisainstrument.noisy import (
    NoiseDefBase,
    NoiseDefCumSum,
    NoiseDefPowerLaw,
    NoiseDefRed,
    NoiseDefViolet,
    NoiseDefInfraRedStat,
    NoiseDefFalloff6Stat,
    NoiseDefFalloff8Stat,
    NoiseDefWhite,
    NoiseDefPink,
    psd_twosided_from_asd_onesided,
    NoiseDefModulation,
    NoiseDefLaser,
    NoiseDefBacklink,
    NoiseDefReciprocalBacklink,
    NoiseDefTestMass,
    LaserNoiseShape,
    TestMassNoiseShape,
)
from lisainstrument.noisy.estimate_psd import estimate_psd_twosided, PSDEstimate
from lisainstrument.noisy import NoiseGen

def plotnoise(noisedef, fmin, nseg=25, nseeds=10, conf_prob=0.95, detrend="constant"):
    """Generate noise samples, compute actual and expected PSD, plot both"""
    nperseg = int(noisedef.f_sample / fmin)
    size = nperseg * nseg
    plt.close()
    for seed in np.arange(nseeds):
        ngen = NoiseGen(noisedef, seed)
        noise = ngen.generate(size)
        est = estimate_psd_twosided(
            noise, noisedef.f_sample, nseg=nseg, conf_prob=conf_prob, detrend=detrend
        )
        lbl = "Welch estimate" if seed == 0 else None
        plt.loglog(est.freq, est.psd_twosided, "k-", alpha=1 / nseeds, label=lbl)

    f = np.exp(np.linspace(np.log(est.freq[0]), np.log(est.freq[-1]), 400))
    psdex = noisedef.discrete_psd_twosided(f)
    plt.loglog(f, psdex, "r-", label="Expected")

    plt.xlabel(r"$f\,[\mathrm{Hz}]$")
    plt.ylabel(r"PSD $[U/\mathrm{Hz}]$")
    plt.title(f"{noisedef.name} ({noisedef.__class__.__name__})")
    plt.legend()

White noise

fsamp = 4.0
asd1s = 10.0
ndef = NoiseDefWhite(f_sample=fsamp, asd_onesided_const=asd1s, name="White")
plotnoise(ndef, fsamp * 1e-3, nseeds=100)
print(ndef)

Noise type NoiseDefWhite
name = White
f_sample [Hz] = 4.0
asd_onesided_const [1/sqrt(Hz)] = 10.0 

../_images/931df25473db89d8f09115e70848e813093ad853bbc9f43df0ccea425e12873f.jpg

wnoise = NoiseGen(ndef, 13).generate(100000)
whist, wedge = np.histogram(wnoise, bins=100, density=True)
plt.close()
plt.stairs(whist, edges=wedge, color="k", label="Sample distribution")
x = np.linspace(wedge[0], wedge[-1], 500)
plt.plot(x, ndef.sample_distribution(x), "g-", label="expected")
plt.legend();

../_images/0e2cfb2a57eb1e09beddfb348219f23ed91d0ece6d241dc1466ec641b95e879d.jpg

Red noise

ndef = NoiseDefRed(f_sample=fsamp, asd_onesided_at_1hz=asd1s, f_min_hz=1e-2, name="Red")
plotnoise(ndef, 1e-4)
plt.axvline(x=ndef.f_min_hz, label=f"f_min_hz = {ndef.f_min_hz}")
plt.loglog(
    [1],
    [psd_twosided_from_asd_onesided(ndef.asd_onesided_at_1hz)],
    "xg",
    ms=10,
    label="Fiducial amplitude",
)
f = np.linspace(1e-3, 2, 200)
plt.loglog(
    f,
    psd_twosided_from_asd_onesided(ndef.asd_onesided_at_1hz) * f**-2,
    "k:",
    label="$f^{-2}$",
)
plt.legend()
print(ndef)

Noise type NoiseDefRed
name = Red
f_sample [Hz] = 4.0
asd_onesided_at_1hz [1/sqrt(Hz)] = 10.0 
f_min_hz [Hz] = 0.01 

../_images/604bd69861fa9fbf18d44b9f3bb3e92cbad4b6709e9162354007e058c970e2c9.jpg

Infrared noise

Note there are two “infrared” noises. The original one is called NoiseDefInfrared and should not be used, except for comparison with old results. It is not stationary and the PSD has a low frequency divergence. Instead one should employ the new variant called NoiseDefInfraredStat which is shown below and does exhibit nonstationary behavior.

ndef = NoiseDefInfraRedStat(
    f_sample=fsamp, asd_onesided_at_1hz=asd1s, f_min_hz=1e-2, name="Infrared"
)
plotnoise(ndef, 1e-4)
plt.axvline(x=ndef.f_min_hz, label=f"f_min_hz = {ndef.f_min_hz}")
plt.loglog(
    [1],
    [psd_twosided_from_asd_onesided(ndef.asd_onesided_at_1hz)],
    "xg",
    ms=10,
    label="Fiducial amplitude",
)
f = np.linspace(1e-3, 2, 200)
plt.loglog(
    f,
    psd_twosided_from_asd_onesided(ndef.asd_onesided_at_1hz) * f**-4,
    "k:",
    label="$f^{-4}$",
)
plt.legend()
print(ndef)

Noise type NoiseDefInfraRedStat
name = Infrared
f_sample [Hz] = 4.0
asd_onesided_at_1hz [1/sqrt(Hz)] = 10.0 
f_min_hz [Hz] = 0.01 

../_images/68b7d8339b4c47f2cc5a8bf9c78fbcb69980af01adb8bc0eeba5dded5fcd3095.jpg

Violet noise

ndef = NoiseDefViolet(f_sample=fsamp, asd_onesided_at_1hz=asd1s, name="Violet")
plotnoise(ndef, 1e-3, nseeds=20)
plt.loglog(
    [1],
    [psd_twosided_from_asd_onesided(ndef.asd_onesided_at_1hz)],
    "xg",
    ms=10,
    label="Fiducial amplitude",
)
f = np.linspace(1e-3, 2, 200)
plt.loglog(
    f,
    psd_twosided_from_asd_onesided(ndef.asd_onesided_at_1hz) * f**2,
    "k:",
    label="$f^{2}$",
)
plt.legend()
print(ndef)

Noise type NoiseDefViolet
name = Violet
f_sample [Hz] = 4.0
asd_onesided_at_1hz [1/sqrt(Hz)] = 10.0 

../_images/1acbe91a6ef3d60372de229d43baf4ea8e22cb7f0c05ae8f399eec45df7b2921.jpg

Pink noise

ndef = NoiseDefPink(
    f_sample=fsamp, asd_onesided_at_1hz=asd1s, f_min_hz=1e-2, name="Pink"
)
plotnoise(ndef, 2e-4, nseeds=5)
plt.axvline(x=ndef.f_min_hz, label=f"f_min_hz = {ndef.f_min_hz}")
f = np.linspace(1e-3, 2, 200)
plt.loglog(
    f,
    psd_twosided_from_asd_onesided(ndef.asd_onesided_at_1hz) * f**-1,
    "k:",
    label="$f^{-1}$",
)
plt.loglog(
    [1],
    [psd_twosided_from_asd_onesided(ndef.asd_onesided_at_1hz)],
    "xg",
    ms=10,
    label="Fiducial amplitude",
)
plt.legend()
print(ndef)

Noise type NoiseDefPink
name = Pink
f_sample [Hz] = 4.0
asd_onesided_at_1hz [1/sqrt(Hz)] = 10.0 
f_min_hz [Hz] = 0.01 

../_images/ed3837836b882ff438ce7da56cb26a09852e748a94e37b153d11e59a74a1392e.jpg

Simulator default noise parameters

For the examples below, we show the instrumental noise types for simulations set up with default parameters of the Instrument class. Extracting the noise definitions for any Instrument instance works the same.

instru = Instrument()
defaults = instru.get_noise_defs()

Using default set of locking beatnote frequencies might cause interferometric beatnote frequencies to fall outside the requirement range of 5..25 MHz

Modulation noise

ndef = defaults.noise_def_modulation(mosa=MosaID.MOSA_12)
plotnoise(ndef, 1e-4, nseeds=5)
plt.axvline(x=ndef.f_min_hz, label=f"f_min_hz = {ndef.f_min_hz}")
plt.loglog(
    [1],
    [psd_twosided_from_asd_onesided(ndef.asd_onesided_at_1hz)],
    "xg",
    ms=10,
    label="Fiducial amplitude",
)
f = np.linspace(1e-3, 2, 200)
plt.loglog(
    f,
    psd_twosided_from_asd_onesided(ndef.asd_onesided_at_1hz) * (f ** (2 / 3)),
    "k:",
    label="$f^{2/3}$",
)
plt.legend()
print(ndef)

Noise type NoiseDefModulation
name = modulation_12
f_sample [Hz] = 16.0
asd_onesided_at_1hz [1/sqrt(Hz)] = 5.2e-14 
f_min_hz [Hz] = 5e-05 

../_images/bade142a369b0b7cc9f08935f7e407e85e21d8da25ff30246acbef6edbf00d03.jpg

Laser

Note the laser noise has several PSD shape options. The default in the simulator is LaserNoiseShape.WHITE_PLUS_INFRARED_STAT. It used to be LaserNoiseShape.WHITE_PLUS_INFRARED which is still available but should not be used in production. The difference is that the new variant is stationary, in contrast to the old one.

ndef = defaults.noise_def_laser(MosaID.MOSA_12)
plotnoise(ndef, 5e-6, nseeds=5)
plt.axvline(x=ndef.f_min_hz, label=f"f_min_hz = {ndef.f_min_hz}")
plt.axvline(x=ndef.f_knee_hz, label=f"f_knee_hz = {ndef.f_knee_hz}")
plt.loglog(
    [1],
    [psd_twosided_from_asd_onesided(ndef.asd_onesided_const)],
    "xg",
    ms=10,
    label="Fiducial amplitude",
)
plt.legend(loc="upper right")
print(ndef)

Noise type NoiseDefLaser
name = laser_12
f_sample [Hz] = 16.0
shape = white+infrared_stationary
asd_onesided_const [Hz/sqrt(Hz)] = 30.0 
f_min_hz [Hz] = 5e-05 
f_knee_hz [Hz] = 0.002 

../_images/8760c6c9b2faaf6eed6e5133d05f03138079cdc87ce02d8b5cb31b8368d10b75.jpg

Backlink

Backlink noise is split into reciprocal and non-reciprocal noise.

ndef = defaults.noise_def_backlink(MosaID.MOSA_12)
plotnoise(ndef, 1e-4)
plt.axvline(x=ndef.f_min_hz, label=f"f_min_hz = {ndef.f_min_hz}")
plt.axvline(x=ndef.f_knee_hz, label=f"f_knee_hz = {ndef.f_knee_hz}")
plt.loglog(
    [1],
    [psd_twosided_from_asd_onesided(ndef.asd_onesided_at_1hz)],
    "xg",
    ms=10,
    label="Fiducial amplitude",
)
f = np.linspace(1e-3, 2, 200)
plt.loglog(
    f,
    psd_twosided_from_asd_onesided(ndef.asd_onesided_at_1hz) * (f**2),
    "k:",
    label="$f^2$",
)
plt.legend()
print(ndef)

Noise type NoiseDefBacklink
name = backlink_12
f_sample [Hz] = 16.0
displacement_asd_onesided [m/sqrt(Hz)] = 3e-12 
f_min_hz [Hz] = 5e-05 
f_knee_hz [Hz] = 0.002 

../_images/ce4b1b0b10dfe1bf32106c6342dd0010fd45cb3eadccfb1994135a1dd14bf1d9.jpg

ndef = defaults.noise_def_reciprocal_backlink(SatID.SAT_1)
plotnoise(ndef, 1e-4)
plt.loglog(
    [1],
    [psd_twosided_from_asd_onesided(ndef.asd_onesided_at_1hz)],
    "xg",
    ms=10,
    label="Fiducial amplitude",
)
f = np.linspace(1e-3, 2, 200)
plt.loglog(
    f,
    psd_twosided_from_asd_onesided(ndef.asd_onesided_at_1hz) * (f**2),
    "k:",
    label="$f^2$",
)
plt.legend()
print(ndef)

Noise type NoiseDefReciprocalBacklink
name = reciprocal_backlink_sc_1
f_sample [Hz] = 16.0
displacement_asd_onesided [m/sqrt(Hz)] = 1.7e-09 

../_images/4a85b43b34eedb16be3ca53a04b3b4c33aeecdf48ee875e8cdb41ef87c414cdf.jpg

Testmass

This noise type has several options for the PSD shape. Two of them, TestMassNoiseShape.ORIGINAL and TestMassNoiseShape.LOWFREQ_RELAX, should not be used in production as they exhibit nonstationary behavior. Instead, one should employ the replacements TestMassNoiseShape.ORIGINAL_STAT and TestMassNoiseShape.LOWFREQ_RELAX_STAT, which do not have this problem. Below we show the simulator default TestMassNoiseShape.ORIGINAL_STAT.

ndef = defaults.noise_def_testmass(MosaID.MOSA_12)
plotnoise(ndef, 5e-6, nseeds=5)
plt.axvline(x=ndef.f_min_hz, label=f"f_min_hz = {ndef.f_min_hz}")
plt.axvline(x=ndef.f_knee_hz, label=f"f_knee_hz = {ndef.f_knee_hz}")
plt.loglog(
    [1],
    [
        psd_twosided_from_asd_onesided(
            ndef.accel_asd_onesided / (2 * np.pi * ndef.f_break_hz**2)
        )
    ],
    "xg",
    ms=10,
    label="Fiducial amplitude",
)
plt.legend(loc="lower right")
print(ndef)

Noise type NoiseDefTestMass
name = testmass_12
f_sample [Hz] = 16.0
shape = original-stationary 
accel_asd_onesided [(m/s^2)/sqrt(Hz)] = 2.4e-15 
f_min_hz [Hz] = 5e-05 
f_knee_hz [Hz] = 0.0004 
f_break_hz [Hz] = 0.008 
f_relax_hz [Hz] = 8e-05 

../_images/aeb98cd0591c163f5deeba3d93b8ec8a7e1ac224aff7b3a0433ffc5cdb71bd9d.jpg

Clock noise

ndef = defaults.noise_def_clock(SatID.SAT_1)
plotnoise(ndef, 5e-6, nseeds=5)
plt.axvline(x=ndef.f_min_hz, label=f"f_min_hz = {ndef.f_min_hz}")
plt.loglog(
    [1],
    [psd_twosided_from_asd_onesided(ndef.asd_onesided_at_1hz)],
    "xg",
    ms=10,
    label="Fiducial amplitude",
)
f = np.linspace(1e-3, 2, 200)
plt.legend(loc="upper right")
print(ndef)

Noise type NoiseDefClock
name = clock_1
f_sample [Hz] = 16.0
asd_onesided_at_1hz [1/sqrt(Hz)] = 6.32e-14 
f_min_hz [Hz] = 1e-05 

../_images/6cefa96c407d22b5185f01c8c23a6831d8929f9f0ab1a9fbcd69f1a650cb1c21.jpg

Ranging noise

ndef = defaults.noise_def_ranging(MosaID.MOSA_12)
plotnoise(ndef, 5e-5, nseeds=5)
plt.loglog(
    [1],
    [psd_twosided_from_asd_onesided(ndef.asd_onesided_const)],
    "xg",
    ms=10,
    label="Fiducial amplitude",
)
f = np.linspace(1e-3, 2, 200)
plt.legend()
print(ndef)

Noise type NoiseDefRanging
name = ranging_12
f_sample [Hz] = 16.0
asd_onesided_const [s/sqrt(Hz)] = 3e-09 

../_images/5e3834cc7c4cd985493b61988b586b97eafbe5adb1921b6d3d64f474245cc062.jpg

OMS SCI Carrier

ndef = defaults.noise_def_oms_sci_carrier(MosaID.MOSA_12)
plotnoise(ndef, 5e-5, nseeds=5)
plt.loglog(
    [1],
    [psd_twosided_from_asd_onesided(ndef.asd_onesided_at_1hz)],
    "xg",
    ms=10,
    label="Fiducial amplitude",
)
plt.legend()
print(ndef)

Noise type NoiseDefOMS
name = oms_sci_carrier_12
f_sample [Hz] = 16.0
displacement_asd_onesided [m/sqrt(Hz)] = 6.35e-12 
f_min_hz [Hz] = 5e-05 
f_knee_hz [Hz] = 0.002 