Trazando una transformada rápida de Fourier en Python

95

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].

user3123955
fuente
muéstrenos lo que ha intentado, cómo falló y los ejemplos con los que está trabajando.
Paul H
Publiqué el ejemplo que probé y lo que pensé de él, creo que estoy confundido sobre cómo trazar la salida correctamente.
user3123955
ese es un gran ejemplo, pero ¿cuál es exactamente el problema? ese código funciona muy bien para mí. ¿La trama simplemente no aparece?
Paul H
es decir, qué tipo de argumentos está utilizando (necesitamos ver al menos algunos de sus datos)
Paul H
He agregado el pastebin de los ejes xey, los datos x están en segundos y los datos y son solo una lectura del sensor. Cuando coloco estas listas de datos en el ejemplo de fft, solo tiene un gran aumento en cero
user3123955

Respuestas:

99

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.

ingrese la descripción de la imagen aquí

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)

ingrese la descripción de la imagen aquí

Paul H
fuente
1
No es que el ejemplo sea incorrecto, es que no sé cómo tomarlo y aplicarlo a mis datos.
user3123955
@ user3123955, a la derecha. es por eso que necesitamos ver sus datos y cómo fallan si vamos a ayudarlo.
Paul H
He agregado el pastebin
user3123955
2
@ user3123955 entonces, ¿qué esperas que haga cualquier algoritmo FFT al respecto? necesita limpiar sus datos.
Paul H
6
@PaulH no debe la amplitud en la frecuencia 50 Hzsea 1y en la frecuencia 80 Hzsea 0.5?
Furqan Hashim
24

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.

ssm
fuente
2
He interpolado mis datos para un espaciado uniforme. ¿Puede decirme exactamente qué hace fftfreq? ¿Por qué necesita mi eje x? ¿Por qué trazas los abdominales de Y y el ángulo? ¿Es el ángulo la fase? ¿A qué se refiere la fase? Cuando hago esto con mis datos, solo tiene un pico gigante a 0Hz y se reduce muy rápidamente, pero lo estoy alimentando con datos que no tienen un desplazamiento constante (hago un paso de banda grande en los datos con bordes de 0.15 Gz a 12Hz para deshacerme del desplazamiento constante, mis datos no deberían ser mayores de 4 Hz de todos modos, por lo que la banda debería hacerme perder información).
user3123955
2
1. fftfreqle proporciona los componentes de frecuencia correspondientes a sus datos. Si grafica freq, verá que el eje x no es una función que sigue aumentando. Deberá asegurarse de tener los componentes de frecuencia correctos en el eje x. Puede consultar el manual: docs.scipy.org/doc/numpy/reference/generated/…
ssm
2
2. A la mayoría de la gente le gustará observar la magnitud y la fase del fft. Es difícil explicar en una oración lo que la información de fase te dirá, pero todo lo que puedo decir es que es significativo cuando combinas señales. Cuando combina señales de la misma frecuencia que están en fase, se amplifican, mientras que cuando están desfasadas en 180 grados, se atenúan. Esto se vuelve importante cuando diseña amplificadores o cualquier elemento que tenga retroalimentación.
ssm
2
3. Generalmente, su frecuencia más baja tendrá prácticamente fase cero, y es en referencia a esto. Cuando las señales se mueven a través de su sistema, cada frecuencia se mueve con una velocidad diferente. Esta es la velocidad de fase. El diagrama de fase le proporciona esta información. No sé con qué sistema está trabajando, así que no puedo darle una respuesta definitiva. Para tales preguntas, es mejor leer sobre control de retroalimentación, electrónica analógica, procesamiento de señales digitales, teoría de campos electromagnéticos, etc., o algo que sea más específico para su sistema.
ssm
4
4. En lugar de utilizar sus datos, ¿por qué no empieza generando sus propias señales t = linspace(0, 10, 1000); ys = [ (1.0/i)*sin(i*t) for i in arange(10)]; y = reduce(lambda m, n: m+n, ys)? Luego grafique cada uno de los ysy el total yy obtenga el fft de cada componente. Ganará confianza con su programación. Entonces puede juzgar la autenticidad de su resultado. Si la señal que está tratando de analizar es la primera de la que ha tomado el fft, siempre sentirá que está haciendo algo mal ...
ssm
12

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: Trazar la señal FFT con DC y luego al eliminarla (frecuencia de salto = 0)

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

Utilizando:

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

Mostrará: ingrese la descripción de la imagen aquí

hesham_EE
fuente
Sí, está en Hz. En el código, la definición de xfasigna los contenedores de fft a las frecuencias.
hesham_EE
1
¡Agradable! ¿Y qué hay del eje y? ¿Amplitud? Muchas gracias hesham_EE
Victor Aguiar
Sí, el eje y es el valor absoluto del complejo fft. Observe el uso denp.abs()
hesham_EE
8

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: ingrese la descripción de la imagen aquí

ewerlopes
fuente
5

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)

resultado de la gráfica FFT analítica

YoniChechik
fuente
4

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")

ingrese la descripción de la imagen aquí

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)

ingrese la descripción de la imagen aquí

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í .

Foad
fuente
No veo nada en los documentos que sugiera que resamplemaneja tiempos de muestra no uniformes. Acepta un parámetro de tiempo (que no se usa en el ejemplo), pero eso parece asumir también tiempos de muestra uniformes.
user2699
@ user2699 este ejemplo podría ayudar
Foad
@ user2699 Editó el código. Le agradecería si pudiera echar un vistazo y decirme si está bien ahora.
Foad
1
'scipy.signal.resample` utiliza el método FFT para volver a muestrear los datos. No tiene sentido usarlo para volver a muestrear datos no uniformes para obtener una FFT uniforme.
user2699
1
Hay ventajas y desventajas en todos los métodos que ha proporcionado (aunque tenga en cuenta que sklearn.utils.resampleno realiza interpolación). Si desea discutir las opciones disponibles para encontrar frecuencias de una señal muestreada irregularmente, o los méritos de diferentes tipos de interpolación, inicie otra pregunta. Ambos son temas interesantes, pero más allá del alcance de las respuestas sobre cómo trazar una FFT.
user2699
4

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.

bousof
fuente