python - 计算功率谱

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 左右的信号。我在这里做错了什么？

score 5 · Accepted Answer

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

你提到在零处看到一个大峰值。正如我在上面的评论中所说，如果您的输入信号具有非零均值，这是可以预期的。如果您想摆脱直流分量，那么您应该在进行 DFT 之前去除信号的趋势，例如减去平均值。
在进行 DFT 之前，您应该始终对信号应用窗口函数，以避免频谱泄漏问题。
尽管取 DFT 的模平方可以粗略估计频谱密度，但这对信号中的任何噪声都非常敏感。对于噪声数据，一种更稳健的方法是计算信号的多个较小段的周期图，然后跨段平均这些周期图。这牺牲了频域中的一些分辨率以提高鲁棒性。韦尔奇的方法使用了这个原理。

我个人会使用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

score 1 · Accepted Answer

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

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)

我得到这个光谱：

python - 计算功率谱

2 回答 2

Related

Reference