在python中绘制功率谱
声明：本页面是StackOverFlow热门问题的中英对照翻译，遵循CC BY-SA 4.0协议，如果您需要使用它，必须同样遵循CC BY-SA许可，注明原文地址和作者信息，同时你必须将它归于原作者（不是我）：StackOverFlow
原文地址: http://stackoverflow.com/questions/15382076/
Warning: these are provided under cc-by-sa 4.0 license. You are free to use/share it, But you must attribute it to the original authors (not me):
StackOverFlow
Plotting power spectrum in python
提问by Olivier_s_j
I have an array with 301 values, which were gathered from a movie clip with 301 frames. This means 1 value from 1 frame. The movie clip is running at 30 fps, so is in fact 10 sec long
我有一个包含 301 个值的数组，这些值是从具有 301 帧的影片剪辑中收集的。这意味着来自 1 帧的 1 个值。影片剪辑以 30 fps 运行，因此实际上是 10 秒长
Now I would like to get the power spectrum of this "signal" ( with the right Axis). I tried:
现在我想获得这个“信号”的功率谱（使用正确的轴）。我试过：
X = fft(S_[:,2]);
pl.plot(abs(X))
pl.show()
I also tried:
我也试过：
X = fft(S_[:,2]);
pl.plot(abs(X)**2)
pl.show()
Though I don't think this is the real spectrum.
虽然我不认为这是真正的频谱。
the signal:
信号：
The spectrum:
频谱：
The power spectrum :
功率谱：
Can anyone provide some help with this ? I would like to have a plot in Hz.
任何人都可以提供一些帮助吗？我想在 Hz 中有一个情节。
采纳答案by Jaime
Numpy has a convenience function, np.fft.fftfreqto compute the frequencies associated with FFT components:
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])
Note that the largest frequency you see in your case is not 30 Hz, but
请注意，您在案例中看到的最大频率不是 30 Hz，而是
In [7]: max(freqs)
Out[7]: 14.950166112956811
You never see the sampling frequency in a power spectrum. If you had had an even number of samples, then you would have reached the Nyquist frequency, 15 Hz in your case (although numpy would have calculated it as -15).
您永远不会在功率谱中看到采样频率。如果您有偶数个样本，那么您将达到奈奎斯特频率，在您的情况下为 15 Hz（尽管 numpy 会将其计算为 -15）。
回答by glormph
From the numpy fft page http://docs.scipy.org/doc/numpy/reference/routines.fft.html:
从 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).
当输入 a 为时域信号且 A = fft(a) 时，np.abs(A) 为其幅度谱，np.abs(A)**2 为其功率谱。相位谱由 np.angle(A) 获得。
回答by HYRY
if rate is the sampling rate(Hz), then np.linspace(0, rate/2, n)is the frequency array of every point in fft. You can use rfftto calculate the fft in your data is real values:
如果 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)
signal x contains 4Hz & 7Hz sin wave, so there are two peaks at 4Hz & 7Hz.
信号 x 包含 4Hz 和 7Hz 正弦波，因此在 4Hz 和 7Hz 处有两个峰值。
回答by Abhi
Since FFT is symmetric over it's centre, half the values are just enough.
由于 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)
回答by Noam Peled
You can also use scipy.signal.welchto estimate the power spectral density using Welch's method. Here is an comparison between np.fft.fft and scipy.signal.welch:
您还可以使用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()
[
[
