Modify the normalization of the noise in DC2.1rc6 reported by JM and aurelien

[mjtueros]

Hello all,
I strongly believe but, needless to say, happy to be wrong that there is a factor 1./ sqrt(2) too much in the generation of the galactic noise in Grandlib.
More precisely
in galaxy.py (grand/grand/sim/noise/) at line 186 I think it should be
v_complex = np.abs(amp * size_out / np.sqrt(2)) * np.exp(1j * phase)

rather than
v_complex = np.abs(amp * size_out / 2) * np.exp(1j * phase)

This would explain the low level of Galactic noise in the latest DC2 release.

I am going to make some plots to backup my claim, but maybe people will agree with me, and I won't have too 😀

[mjtueros]

See DC2 Analysis on May27 and June13

[luckyjim]

Hello,

I've reviewed in detail the generation of noise from a PSD and written tutorials here.
https://github.com/luckyjim/mogwai/tree/main/doc

I then performed the code review of the galaxy module, starting at the end of the v_amplitude_infile calculation. In fact, it's done in just a few lines of code. Here they are.

   # SL
   nb_freq = len(freqs_mhz)
   freq_res = freqs_mhz[1] - freqs_mhz[0]
   v_amplitude_infile = v_amplitude_infile * np.sqrt(freq_res)
   v_amplitude = np.zeros((nb_freq, 3))
   v_amplitude[:, 0] = interpol_at_new_x(gala_freq[:, 0], v_amplitude_infile[:, 0], freqs_mhz)
   v_amplitude[:, 1] = interpol_at_new_x(gala_freq[:, 0], v_amplitude_infile[:, 1], freqs_mhz)
   v_amplitude[:, 2] = interpol_at_new_x(gala_freq[:, 0], v_amplitude_infile[:, 2], freqs_mhz)

   v_amplitude = v_amplitude.T
   rng   = np.random.default_rng(seed)     
   amp   = rng.normal(loc=0, scale=v_amplitude[np.newaxis,...], size=(nb_ant, 3, nb_freq))
   phase = 2 * np.pi * rng.random(size=(nb_ant, 3, nb_freq))
   v_complex = np.abs(amp * size_out / sqrt(2)) * np.exp(1j * phase)

   return v_complex

I assumed that v_amplitude_infile is the ASD and has a correct value.

With the modification mentioned by Matias, the normalization is correct but unclear; there are two terms in the code:

• sqrt(freq_res) at the top
• size_out / sqrt(2) at the bottom

Indeed, there are two normalizations to take into account:

1. SciPy "backward" FFT convention
   Since we are working directly in Fourier space and will ultimately perform an inverse FFT (irfft()) using SciPy's backward convention, there is a factor of sqrt(nb_sample_in_time_series=n_s) to take into account. See here
   https://github.com/luckyjim/mogwai/blob/main/doc/fft_and_white_noise.ipynb

2. Unit constraint with ASD (Amplitude Spectrum Density)
   We multiply an FFT in volts with an ASD in volts/sqrt(Hz). We must normalize with a sqrt(Hz) to obtain volts. The value of this frequency is the Nyquist frequency, see here

https://github.com/luckyjim/mogwai/blob/main/doc/spectrum_with_Welch_method.ipynb

So the "global" normalization is sqrt(f_nyquist*n_s).

The normalization applied is
N = np.sqrt(freq_res) * size_out / sqrt(2)
Since
freq_res = f_s/n_s
we find
N = sqrt(f_nyquist * n_s)

The two normalization terms compensate each other without really understanding the value of each of them...

Then, for noise generation, there is an error in the random selection of the modulus, which is biased by a factor of sqrt(2/pi) ~ 0.8.

Indeed, we want the modulus of the noise FFT to be identical to the ASD, and the code performs a random selection centered on ... 0.

amp = rng.normal(loc=0, scale= ...

then an absolute value of the whole, probably to get rid of negative values. We actually have a folded normal distribution:

https://en.wikipedia.org/wiki/Folded_normal_distribution

(from which I took the factor sqrt(2/pi))

If we want to perform a random draw on the modulus in addition to the phase, we must perform a draw centered on the ASD, like

amp = rng.normal(loc=asd, scale=asd/1000, ....

with small variance relative to ASD.
But this is not necessary unless we want to "noise" the ASD before generating the noise. For a draw based on the phase, this is already sufficient; see here.

https://github.com/luckyjim/mogwai/blob/main/doc/noise_wf.ipynb

Finally, the galactic_noise() function rereads all the files related to the galactic signal model and performs the ASD unitization at each call. This is not optimal.

So I decided to refactor the code that separates ASD preparation and noise generation.
See the refact_galaxy branch linked to this ticket #144
branch
https://github.com/grand-mother/grand/tree/refact_galaxy

[luckyjim]

Remark

the unit of v_amplitude_infile is : $\micro V/\sqrt{MHz}$

see, #144 (comment)

[luckyjim]

Final PSD of galactic signal through antenna provided from file with variable v_amplitude_infile