计算功率谱 [英] Computing a power spectrum

问题描述

我想使用Python3计算功率谱.从有关该主题的另一个主题中，我了解了基本要素.我认为应该是这样的:

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

此处t是时间，而x是光子计数.我也尝试过:

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:

应该给我580Hz左右的信号.我在这里做什么错了?

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

推荐答案

@kwinkunks '答案中我缺少一些东西:

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

1. 您提到在零处看到一个大峰值.正如我在上面的评论中所说，如果您的输入信号的均值非零，那么这是可以预期的.如果您想摆脱直流分量，那么在进行DFT之前应先消除信号趋势，例如减去均值.

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.

在进行DFT之前，应始终对信号应用窗口功能为了避免光谱泄漏.

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

尽管采用DFT的平方模可以粗略估计频谱密度，但这对信号中的任何噪声都非常敏感.一种用于噪声数据的更健壮的方法是计算信号的多个较小段的周期图，然后将其平均化.为了提高鲁棒性，这需要在频域中进行某种分辨率的折衷. 韦尔奇的方法使用了这一原理.

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.

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

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

这篇关于计算功率谱的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！