Вычисление спектра мощности в python
У меня есть массив с 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()
Хотя я не думаю, что это реальный спектр.
сигнал:
![enter image description here]()
Спектр: ![enter image description here]()
Спектр мощности:
![enter image description here]()
Может ли кто-нибудь помочь с этим? Я хотел бы иметь график в Гц.
Ответы
Ответ 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])
![enter image description here]()
Обратите внимание, что наибольшая частота, которую вы видите в вашем случае, не равна 30 Гц, но
In [7]: max(freqs)
Out[7]: 14.950166112956811
Вы никогда не видите частоту дискретизации в спектре мощности. Если бы у вас было четное количество образцов, вы бы достигли частоты Найквиста, 15 Гц в вашем случае (хотя numpy рассчитал бы его как -15).
Ответ 2
если скорость - это частота дискретизации (Гц), то 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)
![enter image description here]()
сигнал x содержит синусоидальную частоту 4 Гц и 7 Гц, поэтому есть два пика при 4 Гц и 7 Гц.
Ответ 3
На странице 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).
Ответ 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)