Instrumental noise#

The ability to add noise to your detector timestreams is supported through the function add_noise_to_observations() and the low-level versions add_noise(), add_white_noise(), and add_one_over_f_noise().

White noise#

Here is a short example that shows how to add white noise to timelines:

import litebird_sim as lbs
import numpy as np

# Create a simulation lasting 100 seconds
sim = lbs.Simulation(
    base_path='./output',
    start_time=0,
    duration_s=100,
    random_seed=12345,
)

# Create a detector object
det = lbs.DetectorInfo(
  net_ukrts=100,
  sampling_rate_hz=10
)

obs = sim.create_observations(detectors=[det])

# Here we add white noise using the detector
# noise parameters from the `det` object.
# We use the random number generator provided
# by `sim`, which is initialized with the
# seed we passed to the Simulation constructor
# to ensure repeatability.
lbs.noise.add_noise_to_observations(obs, 'white', dets_random=sim.dets_random)

for i in range(10):
    print(f"{obs[0].tod[0][i]:.5e}")

-1.37982e-04
3.65642e-04
2.47778e-04
1.78779e-04
-5.03410e-05
4.21404e-04
5.90033e-04
5.07347e-04
-9.98478e-05
-5.19765e-05

Note that we pass sim.dets_random as the detector-level random number generators to use. This is a list member variable that is initialized by the constructor of the class Simulation, and it is safe to be used with multiple MPI processes as it ensures that each detector has its own random number generator with a different seed. You can also pass another list of random number generators, as long as each has the normal method. More information on the generation of random numbers can be found in Random numbers.

To add white noise using a custom white noise sigma, in µK, we can call the low level function directly:

import litebird_sim as lbs

sim = lbs.Simulation(
    base_path='./output',
    start_time=0,
    duration_s=100,
    random_seed=12345,
)

det = lbs.DetectorInfo(
  net_ukrts=100,
  sampling_rate_hz=10,
)

obs = sim.create_observations(detectors=[det])

custom_sigma_uk = 1234
lbs.noise.add_white_noise(obs[0].tod[0], custom_sigma_uk, random=sim.dets_random[0])

1/f Models and Engines#

The framework supports also the generatio of 1/f noise. Here you can choose the computational engine (how it is calculated) and the physical model (the shape of the power spectrum).

Engines#

The engine is selected via the engine parameter:

  1. “fft” (Default): Generates noise in the Fourier domain.

  2. “ducc”: Uses time-domain infinite impulse response (IIR) filtering provided by the ducc0 library. It only supports the “keshner” model.

Models#

The physical shape of the Power Spectral Density (PSD) is selected via the model parameter:

  1. “toast” (Default): The classic power-law ratio, also implemented in https://github.com/hpc4cmb/toast/blob/372fa7642bbe61a5f01d239e707c04b80ad4bf46/src/toast/tod/sim_noise.py#L74. The PSD is proportional to:

    \[P(f) \propto \frac{f^\alpha + f_{knee}^\alpha}{f^\alpha + f_{min}^\alpha}\]

  2. “keshner”: Corresponds to a sum of relaxation processes. This is the native model of the ducc engine. The PSD is proportional to:

    \[P(f) \propto \left( \frac{f^2 + f_{knee}^2}{f^2 + f_{min}^2} \right)^{\alpha/2}\]

This call allows to add 1/f noise:

import litebird_sim as lbs

sim = lbs.Simulation(
    base_path='./output',
    start_time=0,
    duration_s=100,
    random_seed=12345,
)

det = lbs.DetectorInfo(
  net_ukrts=100,
  sampling_rate_hz=10,
  alpha=1,
  fknee_mhz=10
)

obs = sim.create_observations(detectors=[det])

# Here we add 1/f noise using the detector noise
# parameters from the detector object
lbs.noise.add_noise_to_observations(obs, 'one_over_f', dets_random=sim.dets_random)

Again, to generate noise with custom parameters, we can either use the low-level function or edit the Observation object to contain the desired noise parameters.

import litebird_sim as lbs
import numpy as np

sim = lbs.Simulation(
    base_path='./output',
    start_time=0,
    duration_s=100,
    random_seed=12345,
)

det = lbs.DetectorInfo(
  net_ukrts=100,
  sampling_rate_hz=10,
  alpha=1,
  fknee_mhz=10,
  fmin_hz=0.001,
)

obs = sim.create_observations(detectors=[det])

custom_sigma_uk = 1234
custom_fknee_mhz = 12.34
custom_alpha = 1.234
custom_fmin_hz = 0.0123

# Option 1: we call the low-level function directly
lbs.noise.add_one_over_f_noise(
    obs[0].tod[0],
    custom_fknee_mhz,
    custom_fmin_hz,
    custom_alpha,
    custom_sigma_uk,
    obs[0].sampling_rate_hz,
    sim.dets_random[0],
)

# Option 2: we change the values in `obs`
obs[0].fknee_mhz[0] = custom_fknee_mhz
obs[0].fmin_hz[0] = custom_fmin_hz
obs[0].alpha[0] = custom_alpha
obs[0].net_ukrts[0] = (
    custom_sigma_uk / np.sqrt(obs[0].sampling_rate_hz)
)

lbs.noise.add_noise_to_observations(obs, 'one_over_f', dets_random=sim.dets_random)

Warning

It’s crucial to grasp the distinction between the noise level in a timestream and the noise level in a map. While the latter is dependent on the former, the conversion is influenced by several factors. This understanding will empower you in your data analysis tasks.

A common mistake is to use the mission time divided by the number of pixels in the map in a call to func:.add_white_noise. This is wrong, as the noise level per pixel depends on the overall integration time, which is always less than the mission time because of cosmic ray loss, repointing maneuvers, etc. These effects reduce the number of samples in the timeline that can be used to estimate the map, but they do not affect the noise of the timeline.

Methods of the Simulation class#

The class Simulation provides the function Simulation.add_noise() which adds noise to the timelines. All the details of the noise are provided in the class observation and the interface is simplified.

import litebird_sim as lbs
from astropy.time import Time
import numpy as np

start_time = 0
time_span_s = 1000.0
sampling_hz = 10.0
nside = 128

sim = lbs.Simulation(
    start_time=start_time,
    duration_s=time_span_s,
    random_seed=12345,
)

# We pick a simple scanning strategy where the spin axis is aligned
# with the Sun-Earth axis, and the spacecraft spins once every minute
sim.set_scanning_strategy(
    lbs.SpinningScanningStrategy(
        spin_sun_angle_rad=np.deg2rad(0),
        precession_rate_hz=0,
        spin_rate_hz=1 / 60,
        start_time=start_time,
    ),
    delta_time_s=5.0,
 )

# We simulate an instrument whose boresight is perpendicular to
# the spin axis.
sim.set_instrument(
    lbs.InstrumentInfo(
        boresight_rotangle_rad=0.0,
        spin_boresight_angle_rad=np.deg2rad(90),
        spin_rotangle_rad=np.deg2rad(75),
    )
)

# A simple detector looking along the boresight direction
det = lbs.DetectorInfo(
    name="Boresight_detector",
    sampling_rate_hz=sampling_hz,
    bandcenter_ghz=100.0,
    net_ukrts=50.0,
)

sim.create_observations(detectors=det)

sim.add_noise(noise_type='one_over_f')

for i in range(5):
    print(f"{sim.observations[0].tod[0][i]:.5e}")

-6.90763e-05
1.82736e-04
1.23804e-04
8.93039e-05
-2.52559e-05

API reference#

litebird_sim.noise.add_noise(tod, noise_type, sampling_rate_hz, net_ukrts, fknee_mhz, fmin_hz, alpha, dets_random, scale=1.0, engine='fft', model='toast')#

Adds noise (white or 1/f) to a TOD array for a specific detector.

This function handles the correct broadcasting if net_ukrts, fknee_mhz, etc., are arrays (indicating multiple detectors) while the tod is processed one detector at a time.

Parameters:

  • tod (ndarray) – The Time-Ordered Data array of shape (n_detectors, n_samples).

  • noise_type (str) – The type of noise to add: "white" or "one_over_f".

  • sampling_rate_hz (float) – Sampling rate in Hz.

  • net_ukrts (float or array-like) – NET in uK*sqrt(s). Can be a scalar or an array of length n_detectors.

  • fknee_mhz (float or array-like) – Knee frequency in mHz. Can be a scalar or an array of length n_detectors.

  • fmin_hz (float or array-like) – Minimum frequency in Hz. Can be a scalar or an array of length n_detectors.

  • alpha (float or array-like) – Spectral slope. Can be a scalar or an array of length n_detectors.

  • dets_random (list of numpy.random.Generator) – List of random number generators (one per detector).

  • scale (float, optional) – A multiplicative scaling factor applied to the NET. Defaults to 1.0.

  • engine (str, optional) – Computation engine ("fft" or "ducc"). Defaults to "fft".

  • model (str, optional) – Physical noise model ("toast" or "keshner"). Defaults to "toast".

litebird_sim.noise.add_noise_to_observations(observations, noise_type, user_seed=None, dets_random=None, scale=1.0, component='tod', engine='fft', model='toast')#

Adds instrumental noise (white or 1/f) to a list of Observations.

This is the high-level interface for noise simulation. It iterates over detectors and observations, handling random number generator initialization and parameter broadcasting. It modifies the observations in-place.

Parameters:

  • observations (Observation or list of Observation) – The observation(s) to which noise will be added. Can be a single Observation instance or a list of them.

  • noise_type (str) –

    The type of noise to generate. Options are:

    • "white": Uncorrelated Gaussian noise based on NET.

    • "one_over_f": Correlated noise characterized by knee frequency and spectral index.

  • user_seed (int, optional) – A master integer seed used to initialize the random number generators if dets_random is not provided. If None and dets_random is also None, the generators will be initialized unpredictably (usually from the OS entropy source).

  • dets_random (list of numpy.random.Generator, optional) – A list of pre-initialized random number generators, one per detector. If provided, user_seed is ignored. This allows for precise control over the RNG state for reproducibility across different calls.

  • scale (float, optional) – A multiplicative scaling factor applied to the noise level (NET). Defaults to 1.0. Useful for simulating different noise realizations or scaling noise down for debugging.

  • component (str, optional) – The name of the attribute in the Observation objects where the noise should be added. Defaults to "tod". Can be changed (e.g., to "noise_tod") to store noise separately from the signal.

  • engine (str, optional) –

    The computational backend for 1/f noise generation. Defaults to "fft".

    • "fft": Uses Fourier synthesis. Supports all models.

    • "ducc": Uses the ducc0 library’s IIR filter. Supports only the "keshner" model.

  • model (str, optional) –

    The physical model for the 1/f noise Power Spectral Density (PSD). Defaults to "toast".

    • "toast": \(P(f) \propto (f^\alpha + f_{knee}^\alpha) / (f^\alpha + f_{min}^\alpha)\)

    • "keshner": \(P(f) \propto ((f^2 + f_{knee}^2) / (f^2 + f_{min}^2))^{\alpha/2}\)

Raises:

  • ValueError – If noise_type is not recognized.

  • ValueError – If the ducc engine is requested with the toast model.

litebird_sim.noise.add_one_over_f_noise(data, fknee_mhz: float, fmin_hz: float, alpha: float, sigma: float, sampling_rate_hz: float, random, engine: str = 'fft', model: str = 'toast')#

Adds 1/f noise to the data array using a specific engine and physical model.

This function supports multiple noise generation engines (FFT-based synthesis and DUCC0’s time-domain filtering) and multiple physical models for the noise power spectral density (PSD).

Parameters:

  • data (1-D numpy array) – The input time-ordered data (TOD) array (modified in-place).

  • fknee_mhz (float) – The knee frequency in mHz.

  • fmin_hz (float) – The minimum frequency in Hz below which the spectrum flattens.

  • alpha (float) – The spectral slope (e.g., 1.0 for pink noise).

  • sigma (float) – The white noise standard deviation (RMS) per sample.

  • sampling_rate_hz (float) – The sampling rate of the data in Hz.

  • random (numpy.random.Generator) – The random number generator instance.

  • engine (str, optional) –

    The computational method used to generate the noise. Defaults to "fft".

    • "fft": Generates noise in the Fourier domain. Supports all model types.

    • "ducc": Uses ducc0.misc.OofaNoise (time-domain filtering). Supports only the "keshner" model. Very efficient for long streams.

  • model (str, optional) –

    The physical model for the Power Spectral Density (PSD). Defaults to "toast".

    • "toast": The classic power-law ratio model. \(P(f) \propto (f^\alpha + f_{knee}^\alpha) / (f^\alpha + f_{min}^\alpha)\)

    • "keshner": The model implemented by DUCC0 (sum of relaxation processes). \(P(f) \propto ((f^2 + f_{knee}^2) / (f^2 + f_{min}^2))^{\alpha/2}\)

Raises:

ValueError – If the engine is unknown or if the selected engine does not support the requested model.

litebird_sim.noise.add_white_noise(data, sigma: float, random)#

Adds white noise with the given sigma to the array data.

Parameters:

  • data (1-D numpy array) – The input data array (modified in-place).

  • sigma (float) – The white noise level per sample. Be sure not to include cosmic ray loss, repointing maneuvers, etc., as these affect the integration time but not the white noise per sample.

  • random (numpy.random.Generator) – A random number generator that implements the normal method. This is typically obtained from the RNGHierarchy of the Simulation class.

litebird_sim.noise.apply_transfer_function(ft, freqs, fknee_mhz, fmin_hz, alpha, sigma, model_id)#

Applies the selected transfer function model to the Fourier coefficients.

Parameters:

  • ft (array-like) – The Fourier transform of the white noise (modified in-place).

  • freqs (array-like) – The frequency bins corresponding to ft.

  • fknee_mhz (float) – The knee frequency in mHz.

  • fmin_hz (float) – The minimum frequency in Hz.

  • alpha (float) – The spectral index (slope).

  • sigma (float) – The white noise level (standard deviation).

  • model_id (int) – The identifier for the model to use: 0 = ‘toast’ 1 = ‘keshner’

litebird_sim.noise.nearest_pow2(data)#

Returns the next largest power of 2 that will encompass the full data set.

Parameters:

data (1-D numpy array) – The input data array.

litebird_sim.noise.rescale_noise(net_ukrts: float, sampling_rate_hz: float, scale: float)#

Converts NET [uK*sqrt(s)] to sigma per sample [K].

Parameters:

  • net_ukrts (float or array-like) – Noise Equivalent Temperature in micro-Kelvin sqrt(seconds).

  • sampling_rate_hz (float) – The sampling rate in Hz.

  • scale (float) – A multiplicative scaling factor applied to the NET.

Returns:

The standard deviation (sigma) of the white noise per sample in Kelvin.

Return type:

float or array-like