# White noise NET in Radio-astronomy and Cosmology

Create a white noise map and compare with power spectrum expected from the NET

```
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
```

```
integration_time_per_pixel = 100 * u.s
```

```
standard_deviation = net / np.sqrt(integration_time_per_pixel)
```

```
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()
```

```
pixel_area = hp.nside2pixarea(nside)
```

```
white_noise_cl = standard_deviation.value**2 * pixel_area
```

```
white_noise_cl
```

```
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)
```

```
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();
```