python – 计算功率谱

我想用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的信号.我在这做错了什么？

解决方法:

@kwinkunks的答案中有一些我觉得遗漏的东西：

>你提到看到一个大的零峰值.正如我在上面的评论中所说,如果您的输入信号具有非零均值,则可以预期这种情况.如果你想摆脱DC component那么你应该在采取DFT之前去掉你的信号,例如减去平均值.
>在进行DFT之前,应始终对信号应用window function以避免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

标签：python,numpy,fft,signal-processing,spectral-density
来源： https://codeday.me/bug/20191008/1874033.html