计算功率谱

Question

我想使用 Python3 计算功率谱。从关于这个主题的另一个线程中，我得到了基本的成分。我认为它应该是这样的：

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()

这里 t 是时间，x 是光子计数。我也试过：

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

但两者大多只是给我一个零附近的峰值（尽管它们不一样）。我正在尝试从中计算功率谱：

这应该给我一个大约 580Hz 的信号。我在这里做错了什么？

Answer 1

对于它的价值，这是我的做法：

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()

例如，使用这个信号...

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)

我得到这个频谱：

Answer 2

我觉得 @kwinkunks 的回答中缺少一些东西：

您提到在零处看到一个大峰值。正如我在上面的评论中所说，如果您的输入信号具有非零均值，则这是预期的。如果你想摆脱 DC component 那么你应该在进行 DFT 之前去除信号趋势，例如减去平均值。
您应该始终应用 window function to your signal before taking the DFT in order to avoid the problem of spectral leakage。
尽管采用 DFT 的模平方可以粗略估计频谱密度，但这对信号中的任何噪声都非常敏感。对于噪声数据，一种更可靠的方法是计算信号的多个较小段的周期图，然后对这些段进行平均。这会牺牲频域中的一些分辨率来提高鲁棒性。 Welch's method就是利用了这个原理。

我个人会使用 scipy.signal.welch，它解决了我上面提到的所有要点：

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

计算功率谱

Computing a power spectrum

python

signal-processing

numpy

fft

spectral-density