Aplicación de FFT (Fast Fast Fourier Transform)

Importación de dependencias

import numpy as np import matplotlib.pyplot as plt plt.rcParams['figure.figsize'] = [16, 8] plt.rcParams.update({'font.size' : 14})

Ejemplo 1: Filtrado de señales con FFT

Creación de la señal

Ts = 0.001 # Se define un tiempo de muestreo Fs=1/Ts # Se define una frecuencia de muestreo print(f"Frecuencia de muestreo: {Fs} [Hz]")

w1 = 2*np.pi*50 # Se define una frecuencia de 50 Hz para la señal 1 w2 = 2*np.pi*120 # Se define una frecuencia de 120 Hz para la señal 2 n = Ts*np.arange(0, 1000) N = len(n) # N es la cantidad de muestras x = np.sin(w1*n) + np.sin(w2*n) # Se construye la señal compuesta

ruido = 2.5 * np.random.randn(N) # Se adiciona ruido a la señal compuesta x_ruido = x + ruido

Gráfica: señal vs señal con ruido

plt.plot(n,x_ruido,color='r',lw=1.5,label='Señal con ruido') plt.plot(n,x,color='b',lw=1.5,label='Señal limpia') plt.xlim(n[0],n[-1]) plt.xlabel('Tiempo [s]') plt.ylabel('Amplitud') plt.legend() plt.show()

Cómputo de la FFT

X_FFT = np.fft.fft(x_ruido, N) # Se calcula la FFT. El segundo parámetro corresponde al tamaño de la salida mag_X = abs(X_FFT) # Se calcula la Magnitud de la FFT

freq = Fs*np.arange(0, N)/N # Se define el vector de frecuencias

Espectro de la señal

plt.plot(freq, mag_X) plt.xlabel('Frecuencia [Hz]', fontsize='14') plt.ylabel('Magnitud FFT', fontsize='14') plt.xlim(freq[0],freq[-1]) plt.show()

L = np.arange(0, np.floor(N/2), dtype='int') # Para plotear la primera mitad del espectro

plt.plot(freq[L], mag_X[L]) plt.xlabel('Frecuencia [Hz]', fontsize='14') plt.ylabel('Magnitud FFT', fontsize='14') plt.xlim(freq[L[0]],freq[L[-1]]) plt.show()

PSD (Densidad espectral de potencia)

PSD = X_FFT * np.conj(X_FFT) / N # Densidad espectral de potencia (Power spectrum, power per frecuency)

plt.plot(freq[L],PSD[L],color='r',lw=1.5,label='PSD') plt.xlim(freq[L[0]],freq[L[-1]]) plt.xlabel('Frecuencia (Hz)') plt.ylabel('$|C_{k}|^{2}$') plt.legend() plt.show()

Filtrado a través de la PSD

umbral = 100 indices = PSD > umbral # Arroja las potencias superiores al umbral arbitrario PSD_filtrado = PSD * indices # Hace cero a las menores que el umbral en el PSD X_FFT = indices * X_FFT # Hace cero a los correspondientes en los coeficientes de Fourier x_filtrada = np.fft.ifft(X_FFT) # FFT inversa para filtrar la señal de tiempo

fig,axs = plt.subplots(2,1) plt.sca(axs[0]) plt.plot(freq[L],PSD[L],color='r',lw=1.5,label='PSD') plt.xlim(freq[L[0]],freq[L[-1]]) plt.legend() plt.sca(axs[1]) plt.plot(freq[L],PSD_filtrado[L],color='g',lw=1.5,label='PSD Filtrado') plt.xlim(freq[L[0]],freq[L[-1]]) plt.legend() plt.show()

fig,axs = plt.subplots(3,1) plt.sca(axs[0]) plt.plot(n,x,color='b',lw=1.5,label='Señal compuesta') plt.xlim(n[0],n[-1]) plt.legend() plt.sca(axs[1]) plt.plot(n,x_ruido,color='r',lw=1.5,label='Señal con ruido') plt.xlim(n[0],n[-1]) plt.legend() plt.sca(axs[2]) plt.plot(n,x_filtrada,color='g',lw=1.5,label='Señal Filtrada') plt.xlim(n[0],n[-1]) plt.legend() plt.show()

Ejemplo 2 (con audio)

Importación de dependencias

from scipy.io import wavfile # Permite leer y grabar audio from IPython.display import Audio # Permite reproducir el audio

Carga del archivo de audio

AudioName = "hornero.wav" # Archivo de Audio fs, Audiodata = wavfile.read(AudioName) duracion = Audiodata.shape[0]/fs print(f'Duracion = {duracion} , Frecuencia de Muestreo = {fs} [Muestras/Seg]' \ f', Wav format = {Audiodata.dtype}')

Audio(AudioName)

dt = 1/fs # Tiempo entre muestras t = np.arange(0, duracion, dt) # Se genera el vector tiempo x = Audiodata[:,0] # Se crea la señal de sonido

Gráfica del audio en el tiempo

plt.plot(t,x,color='r',lw=1,label='Audio') plt.xlim(t[0],t[-1]) plt.xlabel('Tiempo [s]') plt.xlabel('Amplitud') plt.legend() plt.show()

N = len(t) # Cantidad de muestras totales AUDIO_FFT = np.fft.fft(x,N) # Computo de la FFT mag_AUDIO_FFT = abs(AUDIO_FFT) # Magnitud de los coefiecientes PSD = AUDIO_FFT * np.conj(AUDIO_FFT) / N # Densidad espectral de potencia (Power spectrum, power per frecuency) freq = (1/(dt*N)) * np.arange(N) # Eje x de frecuencias L = np.arange(0, np.floor(N/2), dtype='int') # Para plotear la primera mitad del espectro

Espectro del audio

plt.plot(freq[L], mag_AUDIO_FFT[L]) plt.xlabel('Frecuencia [Hz]') plt.ylabel('Magnitud FFT') plt.xlim(freq[L[0]],freq[L[-1]]) plt.show()

PSD del audio

plt.plot(freq[L],PSD[L],color='r',lw=2,label='PSD') plt.xlim(freq[L[0]],freq[L[-1]]) plt.xlabel('Frecuencia [Hz]') plt.ylabel('Potencia') plt.legend() plt.show()

L = np.arange(0, np.floor(N/10), dtype='int')

plt.plot(freq[L],PSD[L],color='r',lw=2,label='PSD') plt.xlim(freq[L[0]],freq[L[-1]]) plt.xlabel('Frecuencia [Hz]') plt.ylabel('Potencia') plt.legend() plt.show()