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 nosie 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,
    NoiseDefInfraRed,
    NoiseDefFalloff6,
    NoiseDefFalloff8,
    NoiseDefWhite,
    NoiseDefPink,
    psd_twosided_from_asd_onesided,
    NoiseDefModulation,
    NoiseDefLaser,
    NoiseDefBacklink,
    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, ndef.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

ndef = NoiseDefInfraRed(
    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 NoiseDefInfraRed
name = Infrared
f_sample [Hz] = 4.0
asd_onesided_at_1hz [1/sqrt(Hz)] = 10.0 
f_min_hz [Hz] = 0.01 

../_images/924abc00ac22fb370a8a1b31a17a0c669c1500abdc55a266cd50a47f858e0cda.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 use the default parameters of the simulator for the instrumental noise types.

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

ndef = defaults.noise_def_laser(MosaID.MOSA_12)
plotnoise(ndef, 5e-5, 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()
print(ndef)

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

../_images/18fe5f3fbffdd34c6754d80c1a939eaf8b8fe04ad8e551d70792078f2782ba28.jpg

Backlink

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

Testmass

This noise definition is problematic because the PSD model has a low-frequency divergence for the “lowfreq-relax” option. Reproducing the model PSD for generated noise was unsucessful, even when tweaking the Welch method parameters. The default parameters of the simulator do not use the “lowfreq-relax” option. The small deviation in the plot below from the expected value converges when computing longer noise sequences.

ndef = defaults.noise_def_testmass(MosaID.MOSA_12)
plotnoise(ndef, 5e-5, 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()
print(ndef)

Noise type NoiseDefTestMass
name = testmass_12
f_sample [Hz] = 16.0
shape = original 
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/659a7d0f8503714c54626adf2a894c21ca67c953f71f231c8f0dc418abc52e53.jpg

Clock noise

ndef = defaults.noise_def_clock(SatID.SAT_1)
plotnoise(ndef, 5e-5, 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()
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/a4624f7942e3ed04d9917b185fddfa5a2179a452c78614147d2823ec15646f14.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 