在Python中绘制功率谱

Question

我有一个包含301个值的数组，这些值来自于一个包含301帧的电影剪辑。这意味着每帧取一个值。该电影剪辑以30fps的速度运行，因此实际上长度为10秒。

现在我想获取这个“信号”的功率谱（并且需要正确的轴）。我尝试过：

 X = fft(S_[:,2]);
 pl.plot(abs(X))
 pl.show()

我也尝试了：

 X = fft(S_[:,2]);
 pl.plot(abs(X)**2)
 pl.show()

虽然我不认为这是真正的光谱。

信号： 光谱： 功率谱： 有人能帮助一下吗？我想要一个以Hz为单位的图。

- Olivier_s_j

---

为什么您认为这不是真正的频谱？ - Jakub M.

Answer 1

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

在功率谱中，你从来没有看到采样频率。如果你有偶数个样本，那么你将会达到奈奎斯特频率，在你的情况下是15赫兹（虽然numpy会计算为-15）。

- Jaime

---

在您上面的评论中，频率应该使用Hz单位而不是您使用的kHz单位吗？ - Cabbage soup

在这种情况下，x轴和y轴的标签是什么？ - FaCoffee

x轴标签将为Hz，y轴标签将为数据单位的平方。例如，如果数据的单位为m/s，则功率谱将为(m/s)^2。 - Arun

@Arun，功率谱密度的单位是SI^2 / Hz。因此，如果数据是m/s，则y轴的单位为(m/s)^2 / Hz。 - H. Vabri

Answer 2

如果rate是采样率（赫兹），那么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处有两个峰值。

- HYRY

---

一个小更正，在使用fft.rfft时：p[0] -= 6.02; p[-1] -= 6.02 (absfft2[0] /= 2; absfft2[-1] /= 2) - 请参见例如 Numerical Recipes 第 653页。 - denis

我认为最后一行应该是 pl.plot(f, p) 才能运行代码。非常感谢您的回答，很有教育意义。 - wancharle

如果您正在使用 np.fft.frrt，则相应的频率函数为 np.fft.rfftfreq。 - Itamar Mushkin

Answer 3

你也可以使用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()

[

- Noam Peled

---

可能比较fft，使用rfft进行比较会更好。 - endolith

最新文档的新链接：https://docs.scipy.org/doc/scipy/tutorial/signal.html#spectral-analysis-using-welch-s-method - undefined

Answer 4

从 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) 获得。

- glormph

---

我使用np.abs(A)**2添加了图表。但是，我该如何绘制它以便我可以看到Hz？我怀疑当我恰好有301个样本时，它是否从0到301 Hz。 :P - Olivier_s_j

FFT 只能处理等间距数据（例如在常规网格上），而不能处理物理量。 - Francesco Montesano

把结果值取log10，得到的结果以分贝（dB）为单位，这样不是更有价值吗？ - Laurent

Answer 5

由于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)

- Abhi

---

通常情况下，FFT不是对称的。 在您的情况下，np.fft.rfft应该更快地工作，并且使用更少的资源。 - kimstik