用python绘制功率谱

Question

我有一个值为301的数组，它是从一个有301帧的电影剪辑中收集的。这意味着1帧中有1个值。电影剪辑的速度是每秒30帧，实际上是10秒长

现在我想得到这个“信号”的功率谱（右轴）。我试过：

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

我也试过：

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

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

信号：

光谱：

功率谱：

有人能帮忙吗？我想要一个赫兹的曲线图。

Answer 1

如果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处有两个峰值。

Answer 2

从numpy fft页http://docs.scipy.org/doc/numpy/reference/routines.fft.html：

When the input a is a time-domain signal and A = fft(a), np.abs(A) is its amplitude spectrum and np.abs(A)**2 is its power spectrum. The phase spectrum is obtained by np.angle(A).

Answer 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赫兹，而是

In [7]: max(freqs)
Out[7]: 14.950166112956811

你永远看不到功率谱中的采样频率。如果你有偶数个样本，那么在你的例子中，你将达到Nyquist frequency，15hz（尽管numpy将其计算为-15）。