I have a problem with the right scaling of the PSD when I use an window. Let me show you an example.

I am using a box window (no window) and a hamming window for comparison.

This is the signal I produce:

f_aq = 2**18 # acquisition rate
Np = 2**22 # data points

t = np.arange(0,Np) / f_aq

f1 = 70 # frequency 1 of the signal in Hz
A1 = 6
y = A1 * np.sin(2*np.pi*f1 * t)

Now I apply the hamming window:

window = np.hamming(Np)
y_window = y * window

Y = fft(y,n=Np, norm='backward')
Y_window = fft(y_window, n=Np, norm='backward')

freq = fftfreq(Np, 1/f_aq)

and plot the one-sided FFT for the cases box and hamming window:

I scaled the hamming window result with 1/np.mean(window) and as you can see the amplitudes match.

After this, I am going to calculate the Power Spectrum:

# Power spectrum
Y_one_sided_fft_power_spectrum = 2*np.abs((Y_one_sided)/2)**2 /len(y)**2 # here dividing by the number of points which is S1 for a box filter
# S1 = np.sum(window)/len(window)
S1 = np.sum(window)
Y_one_sided_fft_power_spectrum_window = 2*np.abs((Y_one_sided_window/2))**2 /S1**2

plt.plot(freq_one_sided,Y_one_sided_fft_power_spectrum)
plt.plot(freq_one_sided+2,Y_one_sided_fft_power_spectrum_window)




plt.title("one-sided power spectrum")
plt.xlim(-1, 100)
plt.xlabel("Frequency in Hz")
plt.ylabel("Power spectrum")

obtaining the same result for both as you can see in this plot:

Here I had to scale the hamming window result with S1 being the sum of all windows. If I would use an forward FFT, this scaling factor would be the mean of the window.

Now the problem case: The PSDs of the two cases do not match. I am using the EBNW in the right way but the results are different:

# power spectral density (PSD)
Y_one_sided_fft_psd = Y_one_sided_fft_power_spectrum / (f_aq/Np)

# first is the raw fft which is an imaginary number (Y_one_sided)

S1 = np.sum(window)
S2 = np.sum(window**2)
ENBW = f_aq * S2 / S1**2

Y_one_sided_fft_psd_window_1 = 2*np.abs((Y_one_sided_window/2))**2 / (f_aq * S2)

Y_one_sided_fft_psd_window_2 = Y_one_sided_fft_power_spectrum_window / ENBW



plt.plot(freq_one_sided,Y_one_sided_fft_psd)
plt.plot(freq_one_sided+2,Y_one_sided_fft_psd_window_1)


plt.title("one-sided PSD")
plt.xlim(-1, 100)
plt.xlabel("Frequency in Hz")
plt.ylabel("PSD")
# plt.yscale('log')
energy = np.dot(y,y) / len(y)

psd = energy / (f_aq/Np)

plt.axhline(psd, color='k', alpha=0.3)

I do not understand this behaviour since I am using the right scaling factors as given in this document: https://holometer.fnal.gov/GH_FFT.pdf

I double check the result with the welch function of scipy and I obtain exactly the same result. What did I forgot? I manually calculated the psd from the data and the box window is given the right value as you can see in the horizontal line. The hamming window case exhibits a lower amplitude for the PSD.