White noise NET in Radio-astronomy and Cosmology

Create a white noise map and compare with power spectrum expected from the NET
cosmology
python
healpy
Published

April 24, 2020

import healpy as hp
import numpy as np
%matplotlib inline
import matplotlib.pyplot as plt
import astropy.units as u

Noise-Equivalent-Temperature, it is a measure of sensitivity of a detector, in cosmology, it is often quoted in \(\mu K \sqrt(s)\), i.e. it is the sensitivity per unit time and can be divided by the integration time to get the actual standard deviation of the white noise of the instrument.

For example let’s consider a white noise NET of \(200 \mu K \sqrt(s)\)

it means that if you integrate for 100 seconds for each pixel, the standard deviation will be \(20 \mu K\).

net = 200 * u.Unit("uK * sqrt(s)")
net

\(200 \; \mathrm{\mu K\,s^{1/2}}\)

integration_time_per_pixel = 100 * u.s
standard_deviation = net / np.sqrt(integration_time_per_pixel)

Create a white noise map

Now that we have an estimate of the standard deviation per pixel, we can use numpy to create a map of gaussian white noise.

nside = 128
npix = hp.nside2npix(nside)
m = np.random.normal(scale = standard_deviation.value, size=npix) * standard_deviation.unit
hp.mollview(m, unit=m.unit, title="White noise map")

Power spectrum

Finally we can compute the angular power spectrum with anafast, i.e. the power as a function of the angular scales, from low \(\ell\) values for large angular scales, to high \(\ell\) values for small angular scales.

At low \(\ell\) there is not much statistics and the power spectrum is biased, but if we exclude lower ells, we can have an estimate of the white noise \(C_\ell\) coefficients. We can then compare with the theoretical power computed as:

\[ C_\ell = \Omega_{pix}\sigma^2 \]

Where: \(\Omega_{pix}\) is the pixel are in square-radians and \(\sigma^2\) is the white noise standard variance.

cl = hp.anafast(m)
cl[100:].mean()
0.02516797566530168
pixel_area = hp.nside2pixarea(nside)
white_noise_cl = standard_deviation.value**2 * pixel_area
white_noise_cl
0.025566346464760685
plt.figure(figsize=(6,4))
plt.loglog(cl, label="Map power spectrum", alpha=.7)
plt.hlines(white_noise_cl, 0, len(cl), label="White noise level")
plt.xlabel("$\ell$")
plt.ylabel("$C_\ell [\mu K ^ 2]$");

Masking

In case we are removing some pixels from a map, for example to mask out a strong signal (e.g. the Milky Way), our estimate of the power spectrum on the partial sky is lower. However we assume that the properties of the noise will be the same also in the masked region. At first order, for simple masks, we can just correct for the amplitude by dividing the power spectrum by the sky fraction.

m.value[len(m)//2-30000:len(m)//2+30000] = hp.UNSEEN
hp.mollview(m, unit=m.unit, title="White noise map")

cl_masked = hp.anafast(m)
plt.figure(figsize=(6,4))
plt.loglog(cl, label="Map power spectrum", alpha=.7)
plt.loglog(cl_masked, label="Map power spectrum (Masked)", alpha=.7)
plt.hlines(white_noise_cl, 0, len(cl), label="White noise level")
plt.xlabel("$\ell$")
plt.ylabel("$C_\ell [\mu K ^ 2]$")
plt.legend();

sky_fraction = hp.mask_good(m).sum() / len(m)
print(sky_fraction)
0.69482421875
plt.figure(figsize=(6,4))
plt.loglog(cl, label="Map power spectrum", alpha=.7)
plt.loglog(cl_masked / sky_fraction, label="Map power spectrum (Masked) - corrected", alpha=.7)
plt.hlines(white_noise_cl, 0, len(cl), label="White noise level")
plt.xlabel("$\ell$")
plt.ylabel("$C_\ell [\mu K ^ 2]$")
plt.legend();