Trazando una transformada rápida de Fourier en Python

Tengo acceso a NumPy y SciPy y quiero crear una FFT simple de un conjunto de datos. Tengo dos listas, una de yvalores y la otra de marcas de tiempo para esos yvalores.

¿Cuál es la forma más sencilla de introducir estas listas en un método SciPy o NumPy y trazar la FFT resultante?

He buscado ejemplos, pero todos se basan en la creación de un conjunto de datos falsos con una cierta cantidad de puntos de datos, frecuencia, etc. y realmente no muestran cómo hacerlo con solo un conjunto de datos y las marcas de tiempo correspondientes. .

He probado el siguiente ejemplo:

from scipy.fftpack import fft

# Number of samplepoints
N = 600

# Sample spacing
T = 1.0 / 800.0
x = np.linspace(0.0, N*T, N)
y = np.sin(50.0 * 2.0*np.pi*x) + 0.5*np.sin(80.0 * 2.0*np.pi*x)
yf = fft(y)
xf = np.linspace(0.0, 1.0/(2.0*T), N/2)
import matplotlib.pyplot as plt
plt.plot(xf, 2.0/N * np.abs(yf[0:N/2]))
plt.grid()
plt.show()

Pero cuando cambio el argumento de ffta mi conjunto de datos y lo trazo, obtengo resultados extremadamente extraños y parece que la escala de la frecuencia puede estar desactivada. No estoy seguro

Aquí hay un pastebin de los datos que estoy intentando FFT

http://pastebin.com/0WhjjMkb http://pastebin.com/ksM4FvZS

Cuando lo uso fft()en todo, solo tiene un gran pico en cero y nada más.

Aquí está mi código:

## Perform FFT with SciPy
signalFFT = fft(yInterp)

## Get power spectral density
signalPSD = np.abs(signalFFT) ** 2

## Get frequencies corresponding to signal PSD
fftFreq = fftfreq(len(signalPSD), spacing)

## Get positive half of frequencies
i = fftfreq>0

##
plt.figurefigsize = (8, 4));
plt.plot(fftFreq[i], 10*np.log10(signalPSD[i]));
#plt.xlim(0, 100);
plt.xlabel('Frequency [Hz]');
plt.ylabel('PSD [dB]')

El espaciado es igual a xInterp[1]-xInterp[0].

Entonces ejecuto una forma funcionalmente equivalente de su código en un cuaderno IPython:

%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
import scipy.fftpack

# Number of samplepoints
N = 600
# sample spacing
T = 1.0 / 800.0
x = np.linspace(0.0, N*T, N)
y = np.sin(50.0 * 2.0*np.pi*x) + 0.5*np.sin(80.0 * 2.0*np.pi*x)
yf = scipy.fftpack.fft(y)
xf = np.linspace(0.0, 1.0/(2.0*T), N/2)

fig, ax = plt.subplots()
ax.plot(xf, 2.0/N * np.abs(yf[:N//2]))
plt.show()

Obtengo lo que creo que es un resultado muy razonable.

Ha pasado más tiempo de lo que me gustaría admitir desde que estaba en la escuela de ingeniería pensando en el procesamiento de señales, pero los picos de 50 y 80 son exactamente lo que esperaría. Entonces, ¿cuál es el problema?

En respuesta a los datos sin procesar y los comentarios que se publican

El problema aquí es que no tiene datos periódicos. Siempre debe inspeccionar los datos que ingresa en cualquier algoritmo para asegurarse de que sean apropiados.

import pandas
import matplotlib.pyplot as plt
#import seaborn
%matplotlib inline

# the OP's data
x = pandas.read_csv('http://pastebin.com/raw.php?i=ksM4FvZS', skiprows=2, header=None).values
y = pandas.read_csv('http://pastebin.com/raw.php?i=0WhjjMkb', skiprows=2, header=None).values
fig, ax = plt.subplots()
ax.plot(x, y)

Lo importante de fft es que solo se puede aplicar a datos en los que la marca de tiempo es uniforme ( es decir , muestreo uniforme en el tiempo, como lo que ha mostrado anteriormente).

En caso de muestreo no uniforme, utilice una función para ajustar los datos. Hay varios tutoriales y funciones para elegir:

https://github.com/tiagopereira/python_tips/wiki/Scipy%3A-curve-fitting http://docs.scipy.org/doc/numpy/reference/generated/numpy.polyfit.html

Si el ajuste no es una opción, puede usar directamente alguna forma de interpolación para interpolar datos en un muestreo uniforme:

https://docs.scipy.org/doc/scipy-0.14.0/reference/tutorial/interpolate.html

Cuando tenga muestras uniformes, solo tendrá que preocuparse por el tiempo delta ( t[1] - t[0]) de sus muestras. En este caso, puede utilizar directamente las funciones fft

Y    = numpy.fft.fft(y)
freq = numpy.fft.fftfreq(len(y), t[1] - t[0])

pylab.figure()
pylab.plot( freq, numpy.abs(Y) )
pylab.figure()
pylab.plot(freq, numpy.angle(Y) )
pylab.show()

Esto debería solucionar tu problema.

El pico alto que tiene se debe a la porción de CC (no variable, es decir, frecuencia = 0) de su señal. Es una cuestión de escala. Si desea ver contenido de frecuencia que no sea de CC, para visualización, es posible que deba trazar desde el desplazamiento 1, no desde el desplazamiento 0 de la FFT de la señal.

Modificando el ejemplo dado arriba por @PaulH

import numpy as np
import matplotlib.pyplot as plt
import scipy.fftpack

# Number of samplepoints
N = 600
# sample spacing
T = 1.0 / 800.0
x = np.linspace(0.0, N*T, N)
y = 10 + np.sin(50.0 * 2.0*np.pi*x) + 0.5*np.sin(80.0 * 2.0*np.pi*x)
yf = scipy.fftpack.fft(y)
xf = np.linspace(0.0, 1.0/(2.0*T), N/2)

plt.subplot(2, 1, 1)
plt.plot(xf, 2.0/N * np.abs(yf[0:N/2]))
plt.subplot(2, 1, 2)
plt.plot(xf[1:], 2.0/N * np.abs(yf[0:N/2])[1:])

Las gráficas de salida:

Otra forma es visualizar los datos en escala logarítmica:

Utilizando:

plt.semilogy(xf, 2.0/N * np.abs(yf[0:N/2]))

Mostrará:

Como complemento a las respuestas ya dadas, me gustaría señalar que a menudo es importante jugar con el tamaño de los contenedores para la FFT. Tendría sentido probar un montón de valores y elegir el que tenga más sentido para su aplicación. A menudo, tiene la misma magnitud que el número de muestras. Esto fue lo que supusieron la mayoría de las respuestas dadas y produce resultados excelentes y razonables. En caso de que uno quiera explorar eso, aquí está mi versión de código:

%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
import scipy.fftpack

fig = plt.figure(figsize=[14,4])
N = 600           # Number of samplepoints
Fs = 800.0
T = 1.0 / Fs      # N_samps*T (#samples x sample period) is the sample spacing.
N_fft = 80        # Number of bins (chooses granularity)
x = np.linspace(0, N*T, N)     # the interval
y = np.sin(50.0 * 2.0*np.pi*x) + 0.5*np.sin(80.0 * 2.0*np.pi*x)   # the signal

# removing the mean of the signal
mean_removed = np.ones_like(y)*np.mean(y)
y = y - mean_removed

# Compute the fft.
yf = scipy.fftpack.fft(y,n=N_fft)
xf = np.arange(0,Fs,Fs/N_fft)

##### Plot the fft #####
ax = plt.subplot(121)
pt, = ax.plot(xf,np.abs(yf), lw=2.0, c='b')
p = plt.Rectangle((Fs/2, 0), Fs/2, ax.get_ylim()[1], facecolor="grey", fill=True, alpha=0.75, hatch="/", zorder=3)
ax.add_patch(p)
ax.set_xlim((ax.get_xlim()[0],Fs))
ax.set_title('FFT', fontsize= 16, fontweight="bold")
ax.set_ylabel('FFT magnitude (power)')
ax.set_xlabel('Frequency (Hz)')
plt.legend((p,), ('mirrowed',))
ax.grid()

##### Close up on the graph of fft#######
# This is the same histogram above, but truncated at the max frequence + an offset. 
offset = 1    # just to help the visualization. Nothing important.
ax2 = fig.add_subplot(122)
ax2.plot(xf,np.abs(yf), lw=2.0, c='b')
ax2.set_xticks(xf)
ax2.set_xlim(-1,int(Fs/6)+offset)
ax2.set_title('FFT close-up', fontsize= 16, fontweight="bold")
ax2.set_ylabel('FFT magnitude (power) - log')
ax2.set_xlabel('Frequency (Hz)')
ax2.hold(True)
ax2.grid()

plt.yscale('log')

las gráficas de salida:

He creado una función que se ocupa de trazar FFT de señales reales. La ventaja adicional de mi función en relación con los mensajes anteriores es que obtienes la amplitud REAL de la señal. Además, debido a la suposición de una señal real, la FFT es simétrica, por lo que solo podemos trazar el lado positivo del eje x:

import matplotlib.pyplot as plt
import numpy as np
import warnings


def fftPlot(sig, dt=None, plot=True):
    # here it's assumes analytic signal (real signal...)- so only half of the axis is required

    if dt is None:
        dt = 1
        t = np.arange(0, sig.shape[-1])
        xLabel = 'samples'
    else:
        t = np.arange(0, sig.shape[-1]) * dt
        xLabel = 'freq [Hz]'

    if sig.shape[0] % 2 != 0:
        warnings.warn("signal prefered to be even in size, autoFixing it...")
        t = t[0:-1]
        sig = sig[0:-1]

    sigFFT = np.fft.fft(sig) / t.shape[0]  # divided by size t for coherent magnitude

    freq = np.fft.fftfreq(t.shape[0], d=dt)

    # plot analytic signal - right half of freq axis needed only...
    firstNegInd = np.argmax(freq < 0)
    freqAxisPos = freq[0:firstNegInd]
    sigFFTPos = 2 * sigFFT[0:firstNegInd]  # *2 because of magnitude of analytic signal

    if plot:
        plt.figure()
        plt.plot(freqAxisPos, np.abs(sigFFTPos))
        plt.xlabel(xLabel)
        plt.ylabel('mag')
        plt.title('Analytic FFT plot')
        plt.show()

    return sigFFTPos, freqAxisPos


if __name__ == "__main__":
    dt = 1 / 1000

    # build a signal within nyquist - the result will be the positive FFT with actual magnitude
    f0 = 200  # [Hz]
    t = np.arange(0, 1 + dt, dt)
    sig = 1 * np.sin(2 * np.pi * f0 * t) + \
        10 * np.sin(2 * np.pi * f0 / 2 * t) + \
        3 * np.sin(2 * np.pi * f0 / 4 * t) +\
        7.5 * np.sin(2 * np.pi * f0 / 5 * t)

    # res in freqs
    fftPlot(sig, dt=dt)
    # res in samples (if freqs axis is unknown)
    fftPlot(sig)

Ya existen excelentes soluciones en esta página, pero todos han asumido que el conjunto de datos se muestrea / distribuye de manera uniforme / uniforme. Intentaré proporcionar un ejemplo más general de datos muestreados aleatoriamente. También usaré este tutorial de MATLAB como ejemplo:

Añadiendo los módulos necesarios:

import numpy as np
import matplotlib.pyplot as plt
import scipy.fftpack
import scipy.signal

Generando datos de muestra:

N = 600 # number of samples
t = np.random.uniform(0.0, 1.0, N) # assuming the time start is 0.0 and time end is 1.0
S = 1.0 * np.sin(50.0 * 2 * np.pi * t) + 0.5 * np.sin(80.0 * 2 * np.pi * t) 
X = S + 0.01 * np.random.randn(N) # adding noise

Ordenar el conjunto de datos:

order = np.argsort(t)
ts = np.array(t)[order]
Xs = np.array(X)[order]

Remuestreo:

T = (t.max() - t.min()) / N # average period 
Fs = 1 / T # average sample rate frequency
f = Fs * np.arange(0, N // 2 + 1) / N; # resampled frequency vector
X_new, t_new = scipy.signal.resample(Xs, N, ts)

trazar los datos y los datos remuestreados:

plt.xlim(0, 0.1)
plt.plot(t_new, X_new, label="resampled")
plt.plot(ts, Xs, label="org")
plt.legend()
plt.ylabel("X")
plt.xlabel("t")

ahora calculando el fft:

Y = scipy.fftpack.fft(X_new)
P2 = np.abs(Y / N)
P1 = P2[0 : N // 2 + 1]
P1[1 : -2] = 2 * P1[1 : -2]

plt.ylabel("Y")
plt.xlabel("f")
plt.plot(f, P1)

PD : finalmente tuve tiempo de implementar un algoritmo más canónico para obtener una transformada de Fourier de datos distribuidos de manera desigual. Puede ver el código, la descripción y el ejemplo del cuaderno de Jupyter aquí .

Escribo esta respuesta adicional para explicar los orígenes de la difusión de los picos al usar fft y especialmente discutir el tutorial de scipy.fftpack con el que no estoy de acuerdo en algún momento.

En este ejemplo, el tiempo de grabación tmax=N*T=0.75. La señal es sin(50*2*pi*x)+0.5*sin(80*2*pi*x). La señal de frecuencia debe contener 2 picos en las frecuencias 50y 80con amplitudes 1y 0.5. Sin embargo, si la señal analizada no tiene un número entero de periodos, puede aparecer difusión debido al truncamiento de la señal:

• Pike 1: 50*tmax=37.5=> la frecuencia 50no es un múltiplo de 1/tmax=> Presencia de difusión debido al truncamiento de la señal en esta frecuencia.
• Pike 2: 80*tmax=60=> frecuencia 80es un múltiplo de 1/tmax=> No hay difusión debido al truncamiento de la señal en esta frecuencia.

Aquí hay un código que analiza la misma señal que en el tutorial ( sin(50*2*pi*x)+0.5*sin(80*2*pi*x)) pero con las ligeras diferencias:

1. El ejemplo original de scipy.fftpack.
2. El ejemplo original de scipy.fftpack con un número entero de períodos de señal (en tmax=1.0lugar de 0.75evitar la difusión del truncamiento).
3. El ejemplo original de scipy.fftpack con un número entero de periodos de señal y donde las fechas y frecuencias se toman de la teoría FFT.

El código:

import numpy as np
import matplotlib.pyplot as plt
import scipy.fftpack

# 1. Linspace
N = 600
# sample spacing
tmax = 3/4
T = tmax / N # =1.0 / 800.0
x1 = np.linspace(0.0, N*T, N)
y1 = np.sin(50.0 * 2.0*np.pi*x1) + 0.5*np.sin(80.0 * 2.0*np.pi*x1)
yf1 = scipy.fftpack.fft(y1)
xf1 = np.linspace(0.0, 1.0/(2.0*T), N//2)

# 2. Integer number of periods
tmax = 1
T = tmax / N # sample spacing
x2 = np.linspace(0.0, N*T, N)
y2 = np.sin(50.0 * 2.0*np.pi*x2) + 0.5*np.sin(80.0 * 2.0*np.pi*x2)
yf2 = scipy.fftpack.fft(y2)
xf2 = np.linspace(0.0, 1.0/(2.0*T), N//2)

# 3. Correct positionning of dates relatively to FFT theory (arange instead of linspace)
tmax = 1
T = tmax / N # sample spacing
x3 = T * np.arange(N)
y3 = np.sin(50.0 * 2.0*np.pi*x3) + 0.5*np.sin(80.0 * 2.0*np.pi*x3)
yf3 = scipy.fftpack.fft(y3)
xf3 = 1/(N*T) * np.arange(N)[:N//2]

fig, ax = plt.subplots()
# Plotting only the left part of the spectrum to not show aliasing
ax.plot(xf1, 2.0/N * np.abs(yf1[:N//2]), label='fftpack tutorial')
ax.plot(xf2, 2.0/N * np.abs(yf2[:N//2]), label='Integer number of periods')
ax.plot(xf3, 2.0/N * np.abs(yf3[:N//2]), label='Correct positionning of dates')
plt.legend()
plt.grid()
plt.show()

Salida:

Como puede ser aquí, incluso con el uso de un número entero de períodos, todavía queda cierta difusión. Este comportamiento se debe a un mal posicionamiento de fechas y frecuencias en el tutorial scipy.fftpack. Por lo tanto, en la teoría de las transformadas discretas de Fourier:

• la señal debe evaluarse en fechas t=0,T,...,(N-1)*Tdonde T es el período de muestreo y la duración total de la señal es tmax=N*T. Tenga en cuenta que nos detenemos en tmax-T.
• las frecuencias asociadas son f=0,df,...,(N-1)*dfdonde df=1/tmax=1/(N*T)está la frecuencia de muestreo. Todos los armónicos de la señal deben ser múltiplos de la frecuencia de muestreo para evitar la difusión.

En el ejemplo anterior, puede ver que el uso de en arangelugar de linspacepermite evitar una difusión adicional en el espectro de frecuencias. Además, el uso de la linspaceversión también conduce a un desplazamiento de los picos que se encuentran en frecuencias ligeramente más altas de lo que deberían ser, como se puede ver en la primera imagen donde los picos están un poco a la derecha de las frecuencias 50y 80.

Concluiré que el ejemplo de uso debe reemplazarse por el siguiente código (que es menos engañoso en mi opinión):

import numpy as np
from scipy.fftpack import fft
# Number of sample points
N = 600
T = 1.0 / 800.0
x = T*np.arange(N)
y = np.sin(50.0 * 2.0*np.pi*x) + 0.5*np.sin(80.0 * 2.0*np.pi*x)
yf = fft(y)
xf = 1/(N*T)*np.arange(N//2)
import matplotlib.pyplot as plt
plt.plot(xf, 2.0/N * np.abs(yf[0:N//2]))
plt.grid()
plt.show()

Salida (el segundo pico ya no se difunde):

Creo que esta respuesta todavía trae algunas explicaciones adicionales sobre cómo aplicar correctamente la transformada discreta de Fourier. Obviamente, mi respuesta es demasiado larga y siempre hay cosas adicionales que decir (@ewerlopes habló brevemente sobre aliasing, por ejemplo, y se puede decir mucho sobre las ventanas ) así que me detendré. Creo que es muy importante comprender profundamente los principios de la transformada discreta de Fourier al aplicarla porque todos conocemos a mucha gente que agrega factores aquí y allá al aplicarla para obtener lo que quieren.