Computing a power spectrum

问题:

I would like to compute a power spectrum using Python3. From another thread about this topic I got the basic ingredients. I think it should be something like:

ps = np.abs(np.fft.fft(x))**2
timeres = t[1]-t[0]
freqs = np.fft.fftfreq(x.size, timeres)
idx = np.argsort(freqs)
plt.plot(freqs[idx], ps[idx])
plt.show()

Here t are the times and x is the photon count. I have also tried:

W = fftfreq(x.size, timeres=t[1]-t[0])
f_x = rfft(x)
plt.plot(W,f_x)
plt.show()

But both mostly just give me a peak around zero (though they are not the same). I am trying to compute the power spectrum from this:

Which should give me a signal around 580Hz. What am I doing wrong here?

回答1:

There are a few things I feel are missing from @kwinkunks' answer:

1. You mentioned seeing a large peak at zero. As I said in the comments above, this is expected if your input signal has non-zero mean. If you want to get rid of the DC component then you should detrend your signal before taking the DFT, for example by subtracting the mean.

2. You should always apply a window function to your signal before taking the DFT in order to avoid the problem of spectral leakage.

3. Although taking the modulus squared of the DFT will give you a rough estimate of the spectral density, this will be quite sensitive to any noise in your signal. A more robust method for noisy data is to compute the periodograms for multiple smaller segments of your signal, then average these across segments. This trades some resolution in the frequency domain for improved robustness. Welch's method uses this principle.

Personally I would use scipy.signal.welch, which addresses all of the points I mentioned above:

from scipy.signal import welch

f, psd = welch(x,
               fs=1./(t[1]-t[0]),  # sample rate
               window='hanning',   # apply a Hanning window before taking the DFT
               nperseg=256,        # compute periodograms of 256-long segments of x
               detrend='constant') # detrend x by subtracting the mean

回答2:

For what it's worth, here's how I do it:

from scipy.fftpack import fft, fftfreq
import matplotlib.pyplot as plt

dt = 0.001
X = fft(x)
freq = fftfreq(x.size, d=dt)

# Only keep positive frequencies.
keep = freq>=0
X = X[keep]
freq = freq[keep]

ax1 = plt.subplot(111)
ax1.plot(freq, np.absolute(X)/3000.)
ax1.set_xlim(0,60)
plt.show()

For example, using this signal...

T = 10
f = 2
f1 = 5

t = np.linspace(0, T, T/dt)
x = 0.4 * np.cos(2*np.pi*f*t) + np.cos(2*np.pi*f1*t)

I get this for the spectrum: