Discussion:
numpy y fft para obtener frecuencias
Marcos Vanetta
2010-02-23 16:30:29 UTC
Permalink
Hola amigxs,
Estoy renegando con esto hace rato.
Necesito encontrar la frecuencia fundamental de una ecuación:
Es decir...

#!/usr/bin/env python
#-*- coding:utf-8 -*-

import numpy as np
import matplotlib.pyplot as plt

sample_rate = 100.0 #Hz
t = np.arange(0, 10, 0.001)

#función
y = np.sin(2 * np.pi * 50 * t) + 2 * np.sin(2 * np.pi *70 *t + np.pi/4)
freqs = np.fft.fft(y)

Una función cualquiera, que debería presentar picos en 50 y 70Hz
Pero me estoy complicando mucho con las escalas... alguno me podría dar una
mano?
Es decir, una vez que tengo el array freqs (lleno de valores complejos)
puedo identifical el
máximo (calculo que primero tendría que buscar los módulos de cada valor
complejo) y a
partir de ahí buscar la dirección del máximo?

Si alguien me puede iluminar...
muchas gracias
malev
Claudio Freire
2010-02-23 16:47:46 UTC
Permalink
Post by Marcos Vanetta
Una función cualquiera, que debería presentar picos en 50 y 70Hz
Pero me estoy complicando mucho con las escalas... alguno me podría dar una
mano?
Es decir, una vez que tengo el array freqs (lleno de valores complejos)
puedo identifical el
máximo (calculo que primero tendría que buscar los módulos de cada valor
complejo) y a
partir de ahí buscar la dirección del máximo?
Si el sample rate es R, la última frecuencia en la transformada FFT, si mal
no recuerdo, es R/2.
La primera es 0hz. (dc le llaman)

Así que sólo tenés que hacer regla de tres y encontrar el rango de
frecuencias al que responde el slot del máximo (ojo, no es una frecuencia
exacta, sino un rango, además hay ripple y otras yerbas).

Me parece que para encontrar la fundamental hay otras técnicas mejores,
aunque no se qué tan sencillo es implementarlas con scipy.
Manuel Naranjo
2010-02-23 17:13:41 UTC
Permalink
Marcos,
Post by Marcos Vanetta
#!/usr/bin/env python
#-*- coding:utf-8 -*-
import numpy as np
import matplotlib.pyplot as plt
sample_rate = 100.0 #Hz
t = np.arange(0, 10, 0.001)
#función
y = np.sin(2 * np.pi * 50 * t) + 2 * np.sin(2 * np.pi *70 *t + np.pi/4)
freqs = np.fft.fft(y)
Según el teorema de muestreo de nyquist la frecuencia de sampleo tiene
que ser mayor a 2 veces la máxima frecuencia que queres samplear, en tu
caso tendría que ser sample_rate >= 140 Hz

luego te conviene definir a t como: t=np.arange(0,10,1/sample_rate)
Post by Marcos Vanetta
Una función cualquiera, que debería presentar picos en 50 y 70Hz
Pero me estoy complicando mucho con las escalas... alguno me podría
dar una mano?
Es decir, una vez que tengo el array freqs (lleno de valores
complejos) puedo identifical el
máximo (calculo que primero tendría que buscar los módulos de cada
valor complejo) y a
partir de ahí buscar la dirección del máximo?
el tema es que freqs lo tenes normalizado en función de la frecuencia de
sampleo, f=F/Fs y ya no me acuerdo más jeje :S, hace como 3 años que di
esa matería ya.

Pero sí, podes calcular el módulo y plotearlo, y vas a encontrar la fft
de tu salida normalizada en frecuencia. Vas a notar un gran pico en f=0
que es igual a la energía de la señal y la cual se debe al proceso de
aproximación que representa la fft, sino tendrías que ver sólo 4
impulsos: -70, -50, 50, 70 (la fft es compleja y par)

Te recomiendo busques info en http://www.fceia.unr.edu.ar/tesys/

Manuel
Claudio Freire
2010-02-23 17:41:41 UTC
Permalink
Pero sí, podes calcular el módulo y plotearlo, y vas a encontrar la fft de
tu salida normalizada en frecuencia. Vas a notar un gran pico en f=0 que es
igual a la energía de la señal y la cual se debe al proceso de aproximación
que representa la fft, sino tendrías que ver sólo 4 impulsos: -70, -50, 50,
70 (la fft es compleja y par)
En realidad depende de la API usada para la transformación. Cuando
transformás una señal real, las librerías suelen usar una representación más
compacta reconociendo la paridad, y esa representación cambia todo.

Vas a encontrar más que suficientes detalles acá:
http://docs.scipy.org/doc/numpy/reference/routines.fft.html
Yo creo que vos necesitarías usar rfft, a menos que tu señal sea compleja.
Loading...