绘制蟒蛇的功率谱

Question

我有一个301值的数组，它是从带有301帧的电影剪辑中收集的。这意味着从一个帧中得到一个值。电影剪辑的速度是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

Numpy有一个方便的函数，np.fft.fftfreq用于计算与快速傅立叶变换相关的频率：

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

你永远看不到功率谱中的采样频率。如果你有一个偶数的样本，那么你会达到奈奎斯特频率，在你的情况下15赫兹(虽然numpy会计算它为-15)。

Answer 2

如果采样率是采样率(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 3

您也可以使用scipy.signal.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()

[

​