Вычисление спектра мощности в 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()

Хотя я не думаю, что это реальный спектр.

сигнал:

Спектр:

Спектр мощности:

Может ли кто-нибудь помочь с этим? Я хотел бы иметь график в Гц.

Ответ 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).

Ответ 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)

сигнал 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)