在python中绘制功率谱
Posted
技术标签:
【中文标题】在python中绘制功率谱【英文标题】:Plotting power spectrum in python 【发布时间】:2013-03-01 04:43:15 【问题描述】:我有一个包含 301 个值的数组,这些值是从一个包含 301 帧的影片剪辑中收集的。这意味着来自 1 帧的 1 个值。影片剪辑以 30 fps 的速度运行,因此实际上是 10 秒长
现在我想获得这个“信号”的功率谱(使用右轴)。我试过了:
X = fft(S_[:,2]);
pl.plot(abs(X))
pl.show()
我也试过了:
X = fft(S_[:,2]);
pl.plot(abs(X)**2)
pl.show()
虽然我不认为这是真正的光谱。
信号:
频谱:
功率谱:
任何人都可以提供一些帮助吗? 我想要以赫兹为单位的绘图。
【问题讨论】:
为什么你“认为这不是真正的光谱”? 【参考方案1】:来自numpy fft页面http://docs.scipy.org/doc/numpy/reference/routines.fft.html:
当输入a为时域信号且A = fft(a)时,np.abs(A)为 它的幅度谱和 np.abs(A)**2 是它的功率谱。这 相位谱由np.angle(A)得到。
【讨论】:
我用 np.abs(A)**2 添加了情节。不过,我怎样才能绘制它以便我可以看到赫兹?我怀疑它从 0 到 301 Hz,当我正好有 301 个样本时:P 您必须自己动手:FFT 只知道等距数据(例如在规则网格上),而不知道物理量。 用 log10 的结果值来获得以 dB 为单位的结果难道不值得吗?【参考方案2】:如果rate是采样率(Hz),那么np.linspace(0, rate/2, n)
是fft中每个点的频率数组。您可以使用rfft
来计算您数据中的 fft 是真实值:
import numpy as np
import pylab as pl
rate = 30.0
t = np.arange(0, 10, 1/rate)
x = np.sin(2*np.pi*4*t) + np.sin(2*np.pi*7*t) + np.random.randn(len(t))*0.2
p = 20*np.log10(np.abs(np.fft.rfft(x)))
f = np.linspace(0, rate/2, len(p))
plot(f, p)
信号 x 包含 4Hz 和 7Hz 正弦波,因此在 4Hz 和 7Hz 有两个峰值。
【讨论】:
一个小的修正,当使用fft.rfft
: p[0] -= 6.02; p[-1] -= 6.02
(absfft2[0] /= 2; absfft2[-1] /= 2
) -- 参见例如数字食谱 p。第653章
我认为最后一行应该是pl.plot(f, p)
才能运行代码。感谢您的回答,这非常具有指导意义。
如果你使用np.fft.frrt
,频率对应的函数是np.fft.rfftfreq
【参考方案3】:
Numpy 有一个方便的函数 np.fft.fftfreq
来计算与 FFT 分量相关的频率:
from __future__ import division
import numpy as np
import matplotlib.pyplot as plt
data = np.random.rand(301) - 0.5
ps = np.abs(np.fft.fft(data))**2
time_step = 1 / 30
freqs = np.fft.fftfreq(data.size, time_step)
idx = np.argsort(freqs)
plt.plot(freqs[idx], ps[idx])
请注意,您在案例中看到的最大频率不是 30 Hz,而是
In [7]: max(freqs)
Out[7]: 14.950166112956811
您永远不会在功率谱中看到采样频率。如果您有偶数个样本,那么您将达到Nyquist frequency,在您的情况下为 15 Hz(尽管 numpy 会将其计算为 -15)。
【讨论】:
在您上面的评论中,频率是否应该使用 Hz 单位而不是您使用的 kHz 单位? 本例中的 x 轴和 y 轴标签是什么? x 轴标签是 Hz,y 轴标签是数据单位的平方。例如,如果数据的单位为 m/s,则功率谱将为 (m/s)^2。 @Arun,功率谱密度的单位是 SI^2 / Hz。所以如果数据是m/s,y的单位是(m/s)^2/Hz。【参考方案4】:由于 FFT 在其中心上是对称的,因此一半的值就足够了。
import numpy as np
import matplotlib.pyplot as plt
fs = 30.0
t = np.arange(0,10,1/fs)
x = np.cos(2*np.pi*10*t)
xF = np.fft.fft(x)
N = len(xF)
xF = xF[0:N/2]
fr = np.linspace(0,fs/2,N/2)
plt.ion()
plt.plot(fr,abs(xF)**2)
【讨论】:
【参考方案5】:您还可以使用scipy.signal.welch 使用 Welch 方法估计功率谱密度。 下面是 np.fft.fft 和 scipy.signal.welch 的比较:
from scipy import signal
import numpy as np
import matplotlib.pyplot as plt
fs = 10e3
N = 1e5
amp = 2*np.sqrt(2)
freq = 1234.0
noise_power = 0.001 * fs / 2
time = np.arange(N) / fs
x = amp*np.sin(2*np.pi*freq*time)
x += np.random.normal(scale=np.sqrt(noise_power), size=time.shape)
# np.fft.fft
freqs = np.fft.fftfreq(time.size, 1/fs)
idx = np.argsort(freqs)
ps = np.abs(np.fft.fft(x))**2
plt.figure()
plt.plot(freqs[idx], ps[idx])
plt.title('Power spectrum (np.fft.fft)')
# signal.welch
f, Pxx_spec = signal.welch(x, fs, 'flattop', 1024, scaling='spectrum')
plt.figure()
plt.semilogy(f, np.sqrt(Pxx_spec))
plt.xlabel('frequency [Hz]')
plt.ylabel('Linear spectrum [V RMS]')
plt.title('Power spectrum (scipy.signal.welch)')
plt.show()
[
【讨论】:
比较rfft
而不是fft
可能会更好以上是关于在python中绘制功率谱的主要内容,如果未能解决你的问题,请参考以下文章