FFT bins from exact frequencies

Question

Estoy tratando de entender algunos conceptos sobre las Transformadas de Fourier (principalmente en el contexto del procesamiento de señales).

Imaginemos una señal muestreada a 10kHz y que el tamaño de la FFT es 1000. Si 1000 muestras se procesan a través de esta FFT (solo reales, asumiendo ventana rectangular), y si tomamos la amplitud del resultado (solo la primera mitad y descartando la fase), esto da 500 bins de ancho 10Hz (excepto el primero), ya que la frecuencia máxima/Nyquist es de 5000Hz.

Así que:

• bin 0, centrado en 0Hz va desde 0Hz a 5Hz
• bin 1, centrado en 10Hz va desde 5Hz a 15Hz
• bin 2, centrado en 20Hz va desde 15Hz a 25Hz
• bin 3, centrado en 30Hz va desde 25Hz a 35Hz
• (y así sucesivamente...)

Aunque tal espectro se puede representar como una línea, tiendo a pensar en esto como un histograma, una especie de "integral" (estoy usando este término de forma flexible, siéntete libre de corregirme), representando la energía que cae dentro del rango de cada bin.

Por lo tanto, si el objetivo fuera encontrar la energía correspondiente a un bin más amplio (digamos 15Hz a 35Hz), yo "sumaría" los valores de los bins 2 y 3 (o tomaría la raíz cuadrada de la suma de los valores al cuadrado). En primer lugar, ¿tiene sentido agrupar los bins de esta manera, para obtener un único valor para un grupo contiguo de bins (o es completamente absurdo, y nada se puede decir del rango de 15Hz a 35Hz a partir de esos dos bins)?

Luego, digamos que alguien quiere saber la "energía" entre dos frecuencias que no caen en los delimitadores de los bins, digamos 5Hz a 6Hz y 151Hz a 3002Hz. He escuchado que podría ser posible usar la interpolación lineal:

• Para 5->6Hz, tomar la décima parte del valor para el bin de 5->15Hz,
• Para 151->3002Hz, tomar la cuarta parte del valor para 145->155Hz (4/10 para ir de 151 a 155), todos los bins de 155Hz a 2995Hz y siete décimas del bin de 2995->3005Hz (7/10 para ir de 2995 a 3002).

No estoy seguro de si esto es correcto. Si bien el margen de error puede ser bajo en la banda de 151->3002Hz, porque los dos bins incompletos en los bordes podrían ser despreciables en comparación con el número de bins completos en el medio, parece que cualquier valor obtenido para 5->6Hz usando este método podría estar bastante alejado del valor real que habría sido para ese bin si se hubiera usado un tamaño de FFT más grande (y por lo tanto hubiera tenido bins más estrechos que permitieran que 5->6Hz caiga en bins reales).

¿Tiene sentido este tipo de interpolación matemáticamente, especialmente para frecuencias más bajas? ¿Funcionaría para bandas aún más estrechas, por ejemplo, de 5.1Hz a 5.3Hz?

Gracias.

Answer 1

FFT generalmente requiere ventanas de tamaño de potencia de 2, así que usemos solo DFT (alternativamente use $1024$ muestras).

Una onda senoidal con una frecuencia de $6\:\mathrm{Hz}$ no es ortogonal a ninguna de las ondas de $0\:\mathrm{Hz}$, $10\:\mathrm{Hz}$,... con respecto al producto escalar $L^2$ sobre un intervalo de $100\:\mathrm{ms}$, por lo que de hecho aparecerá en todos los bins. Ese es el problema con una función de ventana rectangular: a menos que comiences con una superposición perfecta de solo los valores de frecuencia cuantizados, terminarás con una horrible propagación en todo el rango de frecuencia. Para evitar tener que introducir una función de ventana no trivial aquí, hablemos de señales de soporte compacto, como wavelets, centradas en nuestra ventana DFT. Como seguramente sabes, tales funciones siempre tienen una indeterminación de frecuencia intrínseca, lo que significa esencialmente que la transformada de Fourier (infinita) consiste no en picos definidos claramente (dirac) sino en picos en forma de campana de tipo Gaussiano. Si la confinación en el tiempo fue de $100\:\mathrm{ms}$, la indeterminación de frecuencia será más de $\frac1{100\:\mathrm{ms}}=10\:\mathrm{Hz}$. Así que como ves, la anchura de los bins de DFT no es solo un problema técnico con el algoritmo específico de la transformada de Fourier, representa la incapacidad general para definir la frecuencia de una señal "correctamente procesable" más precisamente que el ancho del bin.

Probablemente ya sepas esto. De todos modos, echemos un vistazo ahora a un wavelet con una frecuencia centrada alrededor de $60\:\mathrm{Hz}$, como obtendrías al aplicar una función de ventana al zumbido principal. Suponiendo que la indeterminación de la frecuencia no sea mayor de lo necesario, esto te dará un pico bastante pronunciado en el bin de $55$ a $65\:\mathrm{Hz}$, con solo pequeños valores en los bins vecinos, por lo que podemos aproximar la energía total, es decir, la norma $L^2$ de nuestra señal, solo por el cuadrado del valor en este bin (esto se debe a la igualdad de Bessel). Del mismo modo, si estuvieras interesado en la energía entre $45$ y $105\:\mathrm{Hz}$, solo sumarías los valores al cuadrado de esos bins y obtendrías, correctamente, la energía total del wavelet. Donde se vuelve interesante es cuando deseas saber sobre la energía en el rango de $61$ a $63\:\mathrm{Hz}$. Según tu propuesta, esto debería calcularse como $\tfrac15$ del valor al cuadrado en el bin de $55$ a $65\:\mathrm{Hz}$, es decir, $\tfrac15$ de la energía total de nuestro wavelet. ¡Y eso es bastante bueno, de hecho, porque como dijimos la energía de este wavelet está realmente dispersa sobre un intervalo de $10\:\mathrm{Hz}$ alrededor de $60\:\mathrm{Hz}$, así que es una aproximación bastante razonable decir que $\tfrac15$ de ella está en el rango de $61$ a $63\:\mathrm{Hz}!

¿Qué pasa con un wavelet centrado alrededor de $65\:\mathrm{Hz}$? Si DFT esto, aparecerá en los bins de $55$ a $65\:\mathrm{Hz}$ y de $65$ a $75\:\mathrm{Hz}$, con valores de cada uno de $\sqrt{\tfrac12}$ de la amplitud total. Si ahora calculas la energía entre $45$ y $105\:\mathrm{Hz}$, obtendrás $$ 0 + \sqrt{\tfrac12}^2E + \sqrt{\tfrac12}^2E + 0 + 0 = E $$ así que nuevamente esa es la energía total, correctamente. Si deseas la energía entre $55$ y $65\:\mathrm{Hz}$, obtienes $$ \sqrt{\tfrac12}^2E = \frac{E}2 $$ lo cual es bastante razonable, porque de hecho solo alrededor de la mitad de la energía de la señal se encuentra en esta banda.

Pero donde comienzas a obtener resultados extraños es cuando calculas la energía entre $55$ y $56\:\mathrm{Hz}$, lo cual resulta en $$ \frac1{10}\sqrt{\tfrac12}^2E = \frac{E}{20} $$ y lo comparas con la energía entre $65$ y $66\:\mathrm{Hz}$, para la cual obviamente obtendrías el mismo resultado. ¡Pero entonces, el wavelet real no tiene realmente ningún componente de frecuencia notable en $55$ a $56\:\mathrm{Hz} en absoluto, mientras que $65$ a $66\:\mathrm{Hz}$ es donde es más fuerte!

En conclusión: tiene sentido hacer esta interpolación, pero debe manejarse con cuidado.

---

Como acabo de darme cuenta, lo que haces de hecho no es una interpolación lineal, sino simplemente una extensión de dominio de orden $0$. Una interpolación lineal o de orden superior sufriría menos del problema que acabo de explicar.

Answer 2

Cuando tomas la DFT (transformada de Fourier discreta, con un número finito de muestras) obtienes una transformada de Fourier discreta-finita, en algunos valores discretos de la frecuencia, digamos, $F_d(\omega_k)$. Estás preguntando sobre la relación entre estos coeficientes de transformada y la "verdadera" transformada de Fourier de la señal discreta $F(\omega)$ (la DTFT que se define para todas las frecuencias, desde 0 hasta la frecuencia de Nyquist): por ejemplo, si $F_d(0)$ (el "primer bin") es igual a $F(0)$, o alguna integral de $F(w)$ alrededor de cero, o qué - y si hay una forma de interpolar $F_d(\omega_k)$ para obtener $F(w)$ para cualquier valor (continuo) de frecuencia.

Para responder a esto, uno debería estudiar la relación precisa de $F_d(\omega_k)$ y $F(\omega)$ en las frecuencias correspondientes. $F_d(\omega_k)$ es equivalente a la transformada "exacta" (DTFT) de la señal original ventaneada, es decir, multiplicada por una señal rectangular. Entonces, $F_d(\omega_k)$ se puede ver como obtenido por el siguiente procedimiento: toma la verdadera transformada de Fourier $F(\omega)$ de la señal original, haz una convolución con una función sinc (con un ancho igual al ancho del "bin"), y muestreala en puntos discretos equiespaciados. Observa que la convolución con un sinc es conceptualmente similar a un desenfoque, o una integración que 'promedia' los valores de las frecuencias vecinas. Pero esta integración no se limita a las frecuencias vecinas (la del 'bin'), de hecho incluye un poco de todas las frecuencias. Este efecto (indeseable) se conoce como Fuga Espectral. (Para aliviar esto, solo parcialmente, se aplica frecuentemente una ventana suave a la señal muestreada)

Entonces, tu intuición es solo parcialmente correcta. $F_d(0)$ corresponde, sí, a una integración de los valores reales de la transformada alrededor de la frecuencia 0, principalmente para tu cero 'bin', pero, tristemente, no se limita a esas frecuencias.

Dado que no hay una relación simple entre $F_d(\omega_k)$ y la verdadera transformada en los puntos $\omega_k$, aún menos se puede decir para valores de frecuencia intermedios. Creo que se podría hacer una interpolación (un sinc sería el ideal, pero típicamente se usa alguna aproximación lineal o cuadrática) para otras frecuencias, pero de todos modos eso nos daría los valores interpolados de $F_d(\cdot)$ (la transformada de la señal ventaneada, no la señal real).