0
$\begingroup$

In order to compute the power spectral density of a signal, we can use np.fft.ftt in python in the following formula:

def direct_fft(f,f_s):
    '''
    - f: signal
    - f_s: sampling frequency
    '''
    return  np.fft.fft(f) / f_s

def PSD_fct(f, f_s):
    Fourier = direct_fft(f,f_s)
    T_obs   = np.real(len(f))/f_s
    return np.real(Fourier  * np.conjugate(Fourier) ) / T_obs

Then if I want to plot this PSD, i need to calculate the corresponding frequencies with np.fft.freq.

But now my problem is that I have multiple datasets with different lengths. And I would like to calculate the average PSD of those datasets. But I can't do a simple sum/N because the frequencies of those PSD are different, right ? Because the array have different length.. How can I calculate the average, I'm confused...

Improve this question
$\endgroup$
10
  • $\begingroup$ Try docs.scipy.org/doc/scipy/reference/generated/…. In most cases taking the FFT of the whole thing will not give you a good estimate of the underlying power spectrum. I suggest reading up on "windowing" and "spectral leakage" $\endgroup$
    – Hilmar
    Commented Mar 28 at 13:07
  • $\begingroup$ Parseval's theorem says that what you're doing is pointless: you can just as well calculate the average square magnitude in the time domain. Alos, why np.real(len(…))? a length of a vector is always already a real number (an integer, to be specific, both in the mathematical as well as in the Python sense). $\endgroup$
    – Marcus Müller
    Commented Mar 28 at 13:17
  • 1
    $\begingroup$ You have to interpolate each dataset to a uniform frequency grid, typically via zero-padding, and then normalize by FFT length, which is the periodogram normalization. Because of the periodogram variance, you may not get the most accurate power estimates from the interpolation. $\endgroup$
    – Baddioes
    Commented Mar 28 at 20:05
  • $\begingroup$ @MarcusMüller and why would that be pointless to calculate the average square magnitude in the time domain ? Indeed, np.real in no needed here. $\endgroup$
    – Apinorr
    Commented Apr 12 at 9:07
  • $\begingroup$ @Himar why on Earth In most cases taking the FFT of the whole thing will not give you a good estimate of the underlying power spectrum ? Okay I'll look into it thank you $\endgroup$
    – Apinorr
    Commented Apr 12 at 9:10

2 Answers 2

Reset to default
0
$\begingroup$

To compute the power spectral density with Python use the Welch method as given by scipy.welch. The function provided in all of these tools properly compensates for all the parameters (window used, fft length) to provide an accurate power spectral density.

With that I recommend that the OP compute the PSD for each dataset using Welch directly and then average those results. The number of samples returned by the Welch function is a parameter and can be set to be the same for each result. Another option is to concatenate the different sets and let Welch do the averaging, but this will then be affected by the discontinuities at each boundary (which can be countered with windowing, which would then modify the PSD result, etc so easier in my opinion to do my first suggestion).

The Welch method in simplest explanation provides a noise reduced estimate of the power spectral density by averaging the result of several overlapped shorter FFTs for a longer dataset. This results in significantly less noise in the estimate than one long FFT. The trade however is longer resolution bandwidth in the measurement since the resolution bandwidth is a reciprocal of the time duration of the signal in each FFT block. (This is what causes the noise floor in an FFT to go up and down; the resulting power from each bin represents the total power within the resolution bandwidth of the measurement). The Welch function provided by Python will accurately compensate for the resolution bandwidth and provide a result as a density in power/Hz units.

Because of that scaling compensation, the Welch method is not appropriate for computing the power of individual tones, or any signal that has an occupied bandwidth that is less than the resolution bandwidth of the measurement (as set by the window used and number of points in the FFT block).

I detail further how the Welch Method compares to a direct FFT spectrum in DSP.SE #79255 and the preferred accuracy of the FFT spectrum for the processing of single tones in DSP.SE #87723.

Improve this answer
$\endgroup$
0
$\begingroup$

So after a lot of thinking and searching, I found a way of computing the average power spectrum of datasets that have different length and hence, different frequencies.

  1. I defined a common frequency grid that covers all datasets

  2. I interpolated Each Power Spectrum onto the common frequency grid with scipy.interpolate.interp1d

  3. I averaged the Interpolated Power Spectra

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.