transformada - señales de audio en python

Trazado del espectro de potencia en python

Question

Tengo una matriz con 301 valores, que se recopilaron de un clip de película con 301 fotogramas. Esto significa 1 valor de 1 marco. El clip de película se está ejecutando a 30 fps, por lo que en realidad dura 10 segundos

Ahora me gustaría obtener el espectro de potencia de esta "señal" (con el eje derecho). Lo intenté:

 X = fft(S_[:,2]);
 pl.plot(abs(X))
 pl.show()

También intenté:

 X = fft(S_[:,2]);
 pl.plot(abs(X)**2)
 pl.show()

Aunque no creo que este sea el espectro real.

la señal:

El espectro:

El espectro de potencia:

¿Alguien puede proporcionar alguna ayuda con esto? Me gustaría tener una parcela en Hz .

Answer 1

Como la FFT es simétrica en su centro, la mitad de los valores son suficientes.

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 2

Desde la página de números de fpy http://docs.scipy.org/doc/numpy/reference/routines.fft.html :

Cuando la entrada a es una señal de dominio de tiempo y A = fft (a), np.abs (A) es su espectro de amplitud y np.abs (A) ** 2 es su espectro de potencia. El espectro de fase se obtiene mediante np.angle (A).

Answer 3

Numpy tiene una función de conveniencia, np.fft.fftfreq para calcular las frecuencias asociadas con los componentes 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])

Tenga en cuenta que la frecuencia más grande que ve en su caso no es de 30 Hz, pero

In [7]: max(freqs)
Out[7]: 14.950166112956811

Nunca se ve la frecuencia de muestreo en un espectro de potencia. Si hubiera tenido un número par de muestras, entonces habría alcanzado la frecuencia de Nyquist , 15 Hz en su caso (aunque Numpy lo habría calculado como -15).

Answer 4

Si la frecuencia es la frecuencia de muestreo (Hz), entonces np.linspace(0, rate/2, n) es la matriz de frecuencia de cada punto en fft. Puede usar rfft para calcular la FFT en sus datos como valores reales:

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)

la señal x contiene una onda sin de 4Hz y 7Hz, así que hay dos picos a 4Hz y 7Hz.

Answer 5

También puede usar scipy.signal.welch para estimar la densidad espectral de potencia utilizando el método de Welch. Aquí hay una comparación entre np.fft.fft y 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()

El