[python] Plotting power spectrum in python

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

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.

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 rfft to calculate the fft in your data is real values:

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.

Since FFT is symmetric over it's centre, half the values are just enough.

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)

From the numpy fft page 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).

You can also use scipy.signal.welch to estimate the power spectral density using Welch’s method. Here is an comparison between np.fft.fft and 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()

[