Cálculo de un espectro de potencia.

Me gustaría calcular un espectro de potencia utilizando Python3. De otro hilo sobre este tema obtuve los ingredientes básicos. Creo que debería ser algo como:

ps = np.abs(np.fft.fft(x))**2 timeres = t[1]-t[0] freqs = np.fft.fftfreq(x.size, timeres) idx = np.argsort(freqs) plt.plot(freqs[idx], ps[idx]) plt.show() 

Aquí están los tiempos y x es el recuento de fotones. También he intentado:

 W = fftfreq(x.size, timeres=t[1]-t[0]) f_x = rfft(x) plt.plot(W,f_x) plt.show() 

Pero ambos, en su mayoría, solo me dan un pico alrededor de cero (aunque no son lo mismo). Estoy tratando de calcular el espectro de potencia a partir de esto:

datos

Lo que debería darme una señal de alrededor de 580Hz. ¿Qué estoy haciendo mal aquí?

Hay algunas cosas que siento que faltan en la respuesta de @kwinkunks :

  1. Usted mencionó ver un gran pico en cero. Como dije en los comentarios anteriores, esto se espera si su señal de entrada tiene una media distinta de cero. Si desea deshacerse del componente de CC, debe desviar su señal antes de tomar la DFT, por ejemplo, restando la media.

  2. Siempre debe aplicar una función de ventana a su señal antes de tomar la DFT para evitar el problema de la fuga espectral .

  3. Aunque tomar el módulo al cuadrado de la DFT le dará una estimación aproximada de la densidad espectral, esto será bastante sensible a cualquier ruido en su señal. Un método más robusto para datos ruidosos es calcular los periodogtwigs para múltiples segmentos más pequeños de su señal, luego promediarlos entre segmentos. Esto intercambia alguna resolución en el dominio de la frecuencia para mejorar la robustez. El método de Welch utiliza este principio.

Personalmente usaría scipy.signal.welch , que aborda todos los puntos que mencioné anteriormente:

 from scipy.signal import welch f, psd = welch(x, fs=1./(t[1]-t[0]), # sample rate window='hanning', # apply a Hanning window before taking the DFT nperseg=256, # compute periodograms of 256-long segments of x detrend='constant') # detrend x by subtracting the mean 

Para lo que vale, así es como lo hago:

 from scipy.fftpack import fft, fftfreq import matplotlib.pyplot as plt dt = 0.001 X = fft(x) freq = fftfreq(x.size, d=dt) # Only keep positive frequencies. keep = freq>=0 X = X[keep] freq = freq[keep] ax1 = plt.subplot(111) ax1.plot(freq, np.absolute(X)/3000.) ax1.set_xlim(0,60) plt.show() 

Por ejemplo, usando esta señal …

 T = 10 f = 2 f1 = 5 t = np.linspace(0, T, T/dt) x = 0.4 * np.cos(2*np.pi*f*t) + np.cos(2*np.pi*f1*t) 

Consigo esto para el espectro:

Ejemplo de espectro simple