Power spectrum, Power spectrum density and windowing functions

Question

At the moment I have some trouble understanding the different methods. I would like to calculate both the Power Spectrum and the Power Spectrum Density.

As the name suggests I thought that $PSD = PS / \Delta f$ (the power spectrum divided by the bandwidth to get Pa / Hz)

For stationary signal I use the Hanning window.

What I did so far in python:

    d = (d.T - d.mean(axis=1)).T
    # convert to torch tensor
    d = torch.from_numpy(d)
    zpd = 1
    Nt = d.shape[-1]
    zeroNt = zpd * Nt
    win = np.hanning(d.shape[-1])
    winfac = np.mean(win ** 2)**-0.5
    # torch.cuda.empty_cache()
    fft = torch.fft.rfft(d * win, n=zeroNt)
    # psd https://dsp.stackexchange.com/questions/32187/what-should-be-the-correct-scaling-for-psd-calculation-using-tt-fft/32205#32205
    psd = fft.abs() ** 2
    fs = SAMPLE_FREQUENCY
    S = np.sum(win ** 2)
    psd = 2 * psd / (fs * S)
PS = fft.abs()**2 / (S**2)

f = torch.fft.rfftfreq(zeroNt, 1 / SAMPLE_FREQUENCY).numpy()
fft[:, 1:-1] *= winfac * 2
fft = Nt / (2 * zeroNt) * np.abs(fft)**2
f = f[:fft.shape[-1]]

So what I did here was:

Split the signal (1D) into chunks (2D).
Remove the mean for each chunk
Zero padding set to 0 (total length is thus zpd=1)
Shape of the signal (Nt)
Calculate the hanning window based on shape
Calculate the window correction factor -> This is one part I do not understand
Use torch for faster computations, but results are the same as with numpy.
calculate the RFFT (real FFT) with the signal, times the widow.
calculate PSD based on What should be the correct scaling for PSD calculation using $\tt fft$
Calculate PS to be the same as scipy.signal.spectrogram (this was trial and error, but in the end it was this function which gives the same result?
Calculate Power as I was used to (where I overwrite fft)

So what I do not understand at first is how to calculate the PD. Why do one want to divided by the squared window?

For the window correction factor (for the PSD) why divide by $\sum x_i^2$? The window correction factor according to this answer is https://nl.mathworks.com/matlabcentral/answers/372516-calculate-windowing-correction-factor, which is correct for winfac, the correction factor that I normally use. S=2999.625 and winfac=1.633095.

Can one explain the connection between PSD and PS? And why another window correction factor is used?

Power spectrum, Power spectrum density and windowing functions

0 Answers0