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:

“fft” (Default): Generates noise in the Fourier domain.
“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:

“toast” (Default): The classic power-law ratio, also implemented in hpc4cmb/toast. The PSD is proportional to:

\[P(f) \propto \frac{f^\alpha + f_{knee}^\alpha}{f^\alpha + f_{min}^\alpha}\]
“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.

Correlated noise#

Real detector arrays are often affected by noise that is partially correlated across detectors — for example through a common thermal bath, shared readout electronics, or optical leakage. The function add_correlated_noise() (and its high-level wrapper add_noise_to_observations()) supports two models for inter-detector correlations.

Common-mode model#

Each detector \(i\) belonging to group \(g\) receives:

\[n_i(t) = \sqrt{\rho_i}\,c_g(t) + \sqrt{1 - \rho_i}\,u_i(t)\]

where \(c_g(t)\) is a shared noise stream for the whole group and \(u_i(t)\) is a detector-unique stream. The parameter \(\rho_i \in [0, 1]\) controls the fraction of variance contributed by the common mode: \(\rho_i = 0\) gives fully independent detectors while \(\rho_i = 1\) makes all detectors in a group identical.

Detectors are assigned to groups via the group_by key (name of a per-detector attribute, e.g. "wafer") or an explicit integer label array passed as groups.

import litebird_sim as lbs
import numpy as np

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

# Two detectors with identical noise parameters
dets = [
    lbs.DetectorInfo(name="det_A", net_ukrts=50, sampling_rate_hz=10,
                     alpha=1.0, fknee_mhz=10, fmin_hz=0.001),
    lbs.DetectorInfo(name="det_B", net_ukrts=50, sampling_rate_hz=10,
                     alpha=1.0, fknee_mhz=10, fmin_hz=0.001),
]

obs = sim.create_observations(detectors=dets)

# Both detectors share the same common-mode stream (rho=0.5 means
# half of the variance is common, half is independent).
lbs.noise.add_noise_to_observations(
    obs,
    noise_type='correlated',
    dets_random=sim.dets_random,
    correlation={
        "groups": [0, 0],   # both detectors in group 0
        "rho": 0.5,
    },
)

print("det_A sample:", f"{obs[0].tod[0][0]:.3e}")
print("det_B sample:", f"{obs[0].tod[1][0]:.3e}")

det_A sample: ...
det_B sample: ...

When detectors should be grouped by a named attribute (e.g. focal-plane wafer), pass group_by="wafer" instead of an explicit array. All detectors that share the same value of the attribute will receive the same common-mode stream.

Cholesky model#

For a richer correlation structure, you can specify a full \(n_{det} \times n_{det}\) correlation matrix \(\mathbf{R}\). The function performs a Cholesky decomposition \(\mathbf{R} = \mathbf{L}\mathbf{L}^T\), generates \(n_{det}\) independent unit-variance noise streams \(z_j(t)\) (one per detector, each with the correct 1/f PSD), and mixes them:

\[n_i(t) = \sigma_i \sum_j L_{ij}\,z_j(t)\]

This allows arbitrary positive-semi-definite correlation structures, including block-diagonal layouts and continuously varying off-diagonal elements.

import litebird_sim as lbs
import numpy as np

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

dets = [
    lbs.DetectorInfo(name=f"det_{i}", net_ukrts=50, sampling_rate_hz=10,
                     alpha=1.0, fknee_mhz=10, fmin_hz=0.001)
    for i in range(3)
]

obs = sim.create_observations(detectors=dets)

# Define a correlation matrix: strong correlation between det_0/det_1,
# det_2 is weakly correlated with the others.
R = np.array([
    [1.0, 0.9, 0.1],
    [0.9, 1.0, 0.1],
    [0.1, 0.1, 1.0],
])

lbs.noise.add_noise_to_observations(
    obs,
    noise_type='correlated',
    dets_random=sim.dets_random,
    correlation={"corr_matrix": R},
)

print("det_0 sample:", f"{obs[0].tod[0][0]:.3e}")
print("det_1 sample:", f"{obs[0].tod[1][0]:.3e}")
print("det_2 sample:", f"{obs[0].tod[2][0]:.3e}")

det_0 sample: ...
det_1 sample: ...
det_2 sample: ...

Note

The matrix \(\mathbf{R}\) must be symmetric and positive-semi-definite (PSD). If the Cholesky decomposition fails because the matrix is numerically singular, add a small diagonal regularisation before passing it:

R += 1e-10 * np.eye(n)

The diagonal of \(\mathbf{R}\) should be 1 (unit variance); per-detector scaling is handled automatically by the NET values stored in the observation.

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_correlated_noise(tod, sampling_rate_hz, net_ukrts, fknee_mhz, fmin_hz, alpha, dets_random, groups=None, rho=0.25, scale=1.0, engine='fft', model='toast', common_mode_type='one_over_f', corr_matrix=None)#

Add correlated noise to a TOD.

Two correlation models are supported, selected by which parameter is provided:

Cholesky model (corr_matrix is given)

A full \(n_{det} \times n_{det}\) correlation matrix \(\mathbf{R}\) is factored as \(\mathbf{R} = \mathbf{L} \mathbf{L}^T\) via Cholesky decomposition. Then \(n_{det}\) independent unit-variance noise streams \(z_j(t)\) are generated (each with the PSD of detector j) and mixed:

\[n_i(t) = \sigma_i \sum_j L_{ij}\,z_j(t)\]

This supports arbitrary positive-semi-definite correlation structures.

Common-mode model (groups is given, default)

Each detector i in group g receives:

\[n_i(t) = \sqrt{\rho_i}\,c_g(t) + \sqrt{1 - \rho_i}\,u_i(t)\]

where \(c_g(t)\) is a shared stream for the whole group and \(u_i(t)\) is a detector-unique stream. This is a rank-1 approximation within each group.

Parameters:

tod (ndarray, shape (n_detectors, n_samples)) – Time-ordered data array, modified in-place.
sampling_rate_hz (float) – Sampling rate in Hz.
net_ukrts (float or array-like) – Noise Equivalent Temperature in μK√s. Scalar or per-detector array.
fknee_mhz (float or array-like) – Knee frequency in mHz. Scalar or per-detector array.
fmin_hz (float or array-like) – Minimum frequency in Hz. Scalar or per-detector array.
alpha (float or array-like) – Spectral slope. Scalar or per-detector array.
dets_random (list of numpy.random.Generator) – One RNG per detector.
groups (array-like of int or None, optional) – Integer group labels, length n_detectors. Used only by the common-mode model. Required when corr_matrix is None.
rho (float or array-like, optional) – Fraction of variance in the common mode. Must be in [0, 1]. Used only by the common-mode model. Defaults to 0.25.
scale (float, optional) – Multiplicative factor applied to NET before unit conversion. Defaults to 1.0.
engine (str, optional) – Noise-generation engine ("fft" or "ducc"). Defaults to "fft".
model (str, optional) – PSD model ("toast" or "keshner"). Defaults to "toast".
common_mode_type (str, optional) – PSD shape for the common-mode stream: "one_over_f" (default) or "white". Used only by the common-mode model.
corr_matrix (array-like of shape (n_detectors, n_detectors) or None, optional) – Symmetric positive-semi-definite correlation matrix. When provided, the Cholesky model is used and groups / rho / common_mode_type are ignored. The diagonal should be 1 (unit variance before per-detector sigma scaling).

Raises:

ValueError – If corr_matrix is provided but is not square, not the right size, or not positive-semi-definite.
ValueError – If corr_matrix is None and groups is also None.
ValueError – If common_mode_type is not recognised (common-mode path).

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', correlation=None)#

Adds noise to a TOD array for all detectors.

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", "one_over_f", or "correlated".
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".
correlation (dict or None, optional) –
Required when noise_type="correlated". Supported keys:
- "corr_matrix" (ndarray): full \((n_{det}, n_{det})\) correlation matrix. When present, the Cholesky model is used and "groups" / "rho" / "common_mode_type" are ignored.
- "groups" (array-like): integer group labels (one per detector). Required when "corr_matrix" is absent.
- "rho" (float or array-like): fraction of variance in the common mode, in [0, 1]. Defaults to 0.25.
- "common_mode_type" (str): PSD shape of the common-mode stream: "one_over_f" (default) or "white".

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', correlation=None)#

Adds instrumental noise (white, 1/f, or correlated) 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": 1/f noise characterised by knee frequency and spectral index.
- "correlated": Common-mode + detector-unique 1/f noise. Requires the correlation argument.
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}\)
correlation (dict or None, optional) –
Required when noise_type="correlated". Supported keys:
- "group_by" (str or None): name of a per-detector attribute on each Observation (e.g. "wafer"). Detectors with the same value share a common-mode stream. None (default) puts all detectors in one group.
- "groups" (array-like): explicit integer group-label array. Takes precedence over "group_by" when present.
- "corr_matrix" (ndarray): full \((n_{det}, n_{det})\) correlation matrix. When present, the Cholesky model is used and "group_by" / "groups" / "rho" / "common_mode_type" are ignored.
- "group_by" (str or None): name of a per-detector attribute on each Observation (e.g. "wafer"). Detectors with the same value share a common-mode stream. None (default) puts all detectors in one group. Used only when "corr_matrix" is absent.
- "groups" (array-like): explicit integer group-label array. Takes precedence over "group_by" when present.
- "rho" (float or array-like): fraction of variance in the common mode. Must be in [0, 1]. Defaults to 0.25.
- "common_mode_type" (str): PSD shape of the common-mode stream: "one_over_f" (default) or "white".

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