在python绘制功率谱(Plotting power spectrum in python)

我与301点的值，将其从具有301帧的影片剪辑聚集的阵列。这意味着从1帧1倍的值。影片剪辑运行以30fps，所以实际上是在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来计算与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

你永远不会看到一个功率谱的采样频率。如果你已经有偶数个样本，那么你就已经达到奈奎斯特频率，在您的案件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包含4赫兹＆7HZ正弦波，所以有以4Hz＆7HZ两个峰。

Answer 3:

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

当输入一个是时域信号，并且A = FFT（a）中，np.abs（A）是它的振幅谱和np.abs（A）** 2是它的功率谱。相位频谱由np.angle（A）获得。

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

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

[