transformada - Trazar una transformación rápida de Fourier en Python
transformada de fourier python (6)
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 el 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 asumido por 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'')
Tengo acceso a numpy y scipy y quiero crear un FFT simple de un conjunto de datos. Tengo dos listas, una con valores y y otra con marcas de tiempo para esos valores y.
¿Cuál es la forma más sencilla de alimentar estas listas en un método descuidado o nudoso y trazar la FFT resultante?
He buscado ejemplos, pero todos se basan en la creación de un conjunto de datos falsos con cierto número 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 intentado 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 fft a mi conjunto de datos y lo trazo, obtengo resultados extremadamente extraños, parece que la escala de la frecuencia puede estar desactivada. No estoy seguro
Aquí hay un pastebin de los datos que estoy tratando de FFT
http://pastebin.com/0WhjjMkb http://pastebin.com/ksM4FvZS
Cuando hago una acción en todo el asunto, solo tiene un pico enorme 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]
El pico alto que tiene se debe a la parte DC (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 la visualización, es posible que necesite trazar desde el desplazamiento 1, no desde el desplazamiento 0 de la FFT de la señal.
Modificando el ejemplo dado anteriormente 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 de registro:
Utilizando:
plt.semilogy(xf, 2.0/N * np.abs(yf[0:N/2]))
Mostrará:
Entonces ejecuto una forma funcionalmente equivalente de su código en un cuaderno de 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 una salida 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 a 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 tienes 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, un 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 a un muestreo uniforme:
https://docs.scipy.org/doc/scipy-0.14.0/reference/tutorial/interpolate.html
Cuando tenga muestras uniformes, solo tendrá que saber el tiempo delta (
t[1] - t[0]
) de sus muestras.
En este caso, puede usar 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.
Tengo que construir una función que se ocupe de trazar FFT de señales reales. En mi función existe la amplitud REAL de la señal (nuevamente, debido a la suposición de señal real, lo que significa simetría ...):
import matplotlib.pyplot as plt
import numpy as np
import warnings
def fftPlot(sig, dt=None, block=False, 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(block=block)
return sigFFTPos, freqAxisPos
if __name__ == "__main__":
dt = 1 / 1000
f0 = 1 / dt / 4
t = np.arange(0, 1 + dt, dt)
sig = np.sin(2 * np.pi * f0 * t)
fftPlot(sig, dt=dt)
fftPlot(sig)
t = np.arange(0, 1 + dt, dt)
sig = np.sin(2 * np.pi * f0 * t) + 10 * np.sin(2 * np.pi * f0 / 2 * t)
fftPlot(sig, dt=dt, block=True)
Ya hay excelentes soluciones en esta página, pero todas han asumido que el conjunto de datos se muestrea / distribuye de manera uniforme / uniforme. Trataré de proporcionar un ejemplo más general de datos muestreados al azar. También usaré este tutorial de MATLAB como ejemplo:
Agregar los módulos requeridos:
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]
Muestreo:
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 para implementar un algoritmo más canónico para obtener una transformación de Fourier de datos distribuidos de manera desigual. Puede ver el código, la descripción y el cuaderno de ejemplos de Jupyter here .