Algoritmo de Goertzel

Técnica en procesamiento de señales digitales

El algoritmo de Goertzel es una técnica de procesamiento de señales digitales (DSP) para la evaluación eficiente de los términos individuales de la transformada discreta de Fourier (DFT). Es útil en ciertas aplicaciones prácticas, como el reconocimiento de tonos de señalización multifrecuencia de doble tono (DTMF) producidos por los botones del teclado de un teléfono analógico tradicional . El algoritmo fue descrito por primera vez por Gerald Goertzel en 1958. ^[1]

Al igual que la DFT, el algoritmo de Goertzel analiza un componente de frecuencia seleccionable de una señal discreta . ^{[2] [3] [4]} A diferencia de los cálculos DFT directos, el algoritmo de Goertzel aplica un solo coeficiente de valor real en cada iteración, utilizando aritmética de valor real para secuencias de entrada de valor real. Para cubrir un espectro completo (excepto cuando se usa para un flujo continuo de datos donde los coeficientes se reutilizan para cálculos posteriores, que tiene una complejidad computacional equivalente a la DFT deslizante ), el algoritmo de Goertzel tiene un orden de complejidad más alto que los algoritmos de transformada rápida de Fourier (FFT), pero para calcular una pequeña cantidad de componentes de frecuencia seleccionados, es numéricamente más eficiente. La estructura simple del algoritmo de Goertzel lo hace muy adecuado para procesadores pequeños y aplicaciones integradas.

El algoritmo de Goertzel también se puede utilizar "a la inversa" como una función de síntesis sinusoidal, que requiere solo 1 multiplicación y 1 resta por muestra generada. ^[5]

El algoritmo

El cálculo principal del algoritmo de Goertzel tiene la forma de un filtro digital , por lo que el algoritmo suele denominarse filtro de Goertzel . El filtro funciona sobre una secuencia de entrada en cascada de dos etapas con un parámetro , que proporciona la frecuencia que se va a analizar, normalizada a radianes por muestra. ${\displaystyle x[n]}$ ${\displaystyle \omega _{0}}$

La primera etapa calcula una secuencia intermedia : ${\displaystyle s[n]}$

La segunda etapa aplica el siguiente filtro a , produciendo la secuencia de salida : ${\displaystyle s[n]}$ ${\displaystyle y[n]}$

Se puede observar que la primera etapa del filtro es un filtro IIR de segundo orden con una estructura de forma directa . Esta estructura particular tiene la propiedad de que sus variables de estado internas son iguales a los valores de salida anteriores de esa etapa. Se supone que todos los valores de entrada para son iguales a 0. Para establecer el estado inicial del filtro de modo que la evaluación pueda comenzar en la muestra , se asignan valores iniciales a los estados del filtro . Para evitar los peligros de aliasing , la frecuencia a menudo se restringe al rango de 0 a π (consulte el teorema de muestreo de Nyquist-Shannon ); usar un valor fuera de este rango no carece de sentido, pero es equivalente a usar una frecuencia aliasing dentro de este rango, ya que la función exponencial es periódica con un período de 2π en . ${\displaystyle x[n]}$ ${\displaystyle n<0}$ ${\displaystyle x[0]}$ ${\displaystyle s[-2]=s[-1]=0}$ ${\displaystyle \omega _{0}}$ ${\displaystyle \omega _{0}}$

Se puede observar que el filtro de segunda etapa es un filtro FIR , ya que sus cálculos no utilizan ninguna de sus salidas pasadas.

Los métodos de transformada Z se pueden aplicar para estudiar las propiedades de la cascada de filtros. La transformada Z de la primera etapa de filtro dada en la ecuación (1) es

La transformada Z de la segunda etapa de filtro dada en la ecuación (2) es

La función de transferencia combinada de la cascada de las dos etapas de filtrado es entonces

Esto se puede transformar nuevamente en una secuencia equivalente en el dominio del tiempo, y los términos se pueden desenrollar nuevamente hasta el primer término de entrada en el índice : ^{[ cita requerida ]} ${\displaystyle n=0}$

Estabilidad numérica

Se puede observar que los polos de la transformada Z del filtro están ubicados en y , en un círculo de radio unitario centrado en el origen del plano complejo de la transformada Z. Esta propiedad indica que el proceso del filtro es marginalmente estable y vulnerable a la acumulación de errores numéricos cuando se calcula utilizando aritmética de baja precisión y secuencias de entrada largas. ^[6] Christian Reinsch propuso una versión numéricamente estable . ^[7] ${\displaystyle e^{+j\omega _{0}}}$ ${\displaystyle e^{-j\omega _{0}}}$

Cálculos DFT

Para el caso importante de calcular un término DFT, se aplican las siguientes restricciones especiales.

• El filtrado termina en el índice , donde es el número de términos en la secuencia de entrada de la DFT. ${\displaystyle n=N}$ ${\displaystyle N}$
• Las frecuencias elegidas para el análisis de Goertzel están restringidas a la forma especial

• El número de índice que indica el "bin de frecuencia" de la DFT se selecciona del conjunto de números de índice. ${\displaystyle k}$

Haciendo estas sustituciones en la ecuación (6) y observando que el término , la ecuación (6) toma entonces la siguiente forma: ${\displaystyle e^{+j2\pi k}=1}$

Podemos observar que el lado derecho de la ecuación (9) es extremadamente similar a la fórmula de definición para el término DFT , el término DFT para el número índice , pero no exactamente igual. La suma que se muestra en la ecuación (9) requiere términos de entrada, pero solo los términos de entrada están disponibles cuando se evalúa una DFT. Un expediente simple pero poco elegante es extender la secuencia de entrada con un valor artificial más . ^[8] Podemos ver a partir de la ecuación (9) que el efecto matemático en el resultado final es el mismo que eliminar el término de la suma, entregando así el valor DFT deseado. ${\displaystyle X[k]}$ ${\displaystyle k}$ ${\displaystyle N+1}$ ${\displaystyle N}$ ${\displaystyle x[n]}$ ${\displaystyle x[N]=0}$ ${\displaystyle x[N]}$

Sin embargo, existe un enfoque más elegante que evita el paso de filtro adicional. De la ecuación (1), podemos observar que cuando se utiliza el término de entrada extendido en el paso final, ${\displaystyle x[N]=0}$

De esta forma el algoritmo se puede completar de la siguiente manera:

• terminar el filtro IIR después de procesar el término de entrada , ${\displaystyle x[N-1]}$
• aplicar la ecuación (10) para construir a partir de los resultados anteriores y , ${\displaystyle s[N]}$ ${\displaystyle s[N-1]}$ ${\displaystyle s[N-2]}$
• aplicar la ecuación (2) con el valor calculado y con el producido por el cálculo directo final del filtro. ${\displaystyle s[N]}$ ${\displaystyle s[N-1]}$

Las dos últimas operaciones matemáticas se simplifican combinándolas algebraicamente:

Tenga en cuenta que al detener las actualizaciones del filtro en el término y aplicar inmediatamente la ecuación (2) en lugar de la ecuación (11) se pierden las actualizaciones finales del estado del filtro, lo que produce un resultado con una fase incorrecta. ^[9] ${\displaystyle N-1}$

La estructura de filtrado particular elegida para el algoritmo de Goertzel es la clave para sus cálculos DFT eficientes. Podemos observar que solo se utiliza un valor de salida para calcular la DFT, por lo que se omiten los cálculos para todos los demás términos de salida. Dado que no se calcula el filtro FIR, los cálculos de la etapa IIR , etc. se pueden descartar inmediatamente después de actualizar el estado interno de la primera etapa. ${\displaystyle y[N]}$ ${\displaystyle s[0],s[1]}$

Esto parece generar una paradoja: para completar el algoritmo, la etapa de filtrado FIR debe evaluarse una vez utilizando las dos salidas finales de la etapa de filtrado IIR, mientras que para lograr eficiencia computacional, la iteración del filtro IIR descarta sus valores de salida. Aquí es donde se aplican las propiedades de la estructura de filtrado de forma directa. Las dos variables de estado internas del filtro IIR proporcionan los dos últimos valores de la salida del filtro IIR, que son los términos necesarios para evaluar la etapa de filtrado FIR.

Aplicaciones

Términos de espectro de potencia

Examinando la ecuación (6), un pase final del filtro IIR para calcular el término usando un valor de entrada suplementario aplica un multiplicador complejo de magnitud 1 al término anterior . En consecuencia, y representan la potencia de señal equivalente. Es igualmente válido aplicar la ecuación (11) y calcular la potencia de señal a partir del término o aplicar la ecuación (2) y calcular la potencia de señal a partir del término . Ambos casos conducen a la siguiente expresión para la potencia de señal representada por el término DFT : ${\displaystyle y[N]}$ ${\displaystyle x[N]=0}$ ${\displaystyle y[N-1]}$ ${\displaystyle y[N]}$ ${\displaystyle y[N-1]}$ ${\displaystyle y[N]}$ ${\displaystyle y[N-1]}$ ${\displaystyle X[k]}$

En el pseudocódigo a continuación, los datos de entrada de valor real se almacenan en la matriz x y las variables sprevy sprev2almacenan temporalmente el historial de salida del filtro IIR. Ntermses el número de muestras en la matriz y Ktermcorresponde a la frecuencia de interés, multiplicada por el período de muestreo.

Nterms definidos aquíKterm seleccionado aquíω = 2 × π × Ktérmino / Ntérminos;coeficiente := 2 × cos(ω)difundir := 0sprev2 := 0para cada índice n en el rango 0 a Nterms-1 hacer s := x[n] + coef × sprev - sprev2 sprev2 := sprev difundir := sfinpotencia := sprev 2 + sprev2 2 - (coeficiente × sprev × sprev2)

Es posible ^[10] organizar los cálculos de modo que las muestras entrantes se entreguen individualmente a un objeto de software que mantiene el estado del filtro entre actualizaciones, y se acceda al resultado de potencia final después de que se realice el resto del procesamiento.

Término DFT único con aritmética de valor real

El caso de datos de entrada con valores reales surge con frecuencia, especialmente en sistemas embebidos donde los flujos de entrada resultan de mediciones directas de procesos físicos. Cuando los datos de entrada tienen valores reales, las variables de estado internas del filtro sprevtambién sprev2pueden observarse como valores reales, por lo que no se requiere ninguna aritmética compleja en la primera etapa de IIR. La optimización para la aritmética con valores reales suele ser tan simple como aplicar tipos de datos con valores reales apropiados para las variables.

Una vez finalizados los cálculos con el término de entrada y las iteraciones de filtro, se debe aplicar la ecuación (11) para evaluar el término DFT. El cálculo final utiliza aritmética de valores complejos, pero se puede convertir en aritmética de valores reales separando los términos reales e imaginarios: ${\displaystyle x[N-1]}$

En comparación con la aplicación de espectro de potencia, la única diferencia es el cálculo utilizado para finalizar:

(Los mismos cálculos de filtro IIR que en la implementación de potencia de señal)XKreal = sprev * cr - sprev2;XKimag = sprev * ci;

Detección de fase

Esta aplicación requiere la misma evaluación del término DFT , como se explicó en la sección anterior, utilizando un flujo de entrada de valor real o de valor complejo. Luego, la fase de la señal se puede evaluar como ${\displaystyle X[k]}$

tomando las precauciones adecuadas para singularidades, cuadrantes, etc. al calcular la función tangente inversa.

Señales complejas en aritmética real

Dado que las señales complejas se descomponen linealmente en partes reales e imaginarias, el algoritmo de Goertzel se puede calcular en aritmética real por separado sobre la secuencia de partes reales, obteniendo , y sobre la secuencia de partes imaginarias, obteniendo . Después de eso, los dos resultados parciales de valor complejo se pueden recombinar: ${\displaystyle y_{\text{r}}[n]}$ ${\displaystyle y_{\text{i}}[n]}$

Complejidad computacional

• Según la teoría de la complejidad computacional , calcular un conjunto de términos DFT utilizando aplicaciones del algoritmo de Goertzel en un conjunto de datos con valores con un "costo por operación" de tiene complejidad . ${\displaystyle M}$ ${\displaystyle M}$ ${\displaystyle N}$ ${\displaystyle K}$ ${\displaystyle O(KNM)}$

Para calcular un único intervalo DFT para una secuencia de entrada compleja de longitud , el algoritmo de Goertzel requiere multiplicaciones y sumas/restas dentro del bucle, así como 4 multiplicaciones y 4 sumas/restas finales, para un total de multiplicaciones y sumas/restas. Esto se repite para cada una de las frecuencias. ${\displaystyle X(f)}$ ${\displaystyle N}$ ${\displaystyle 2N}$ ${\displaystyle 4\ N}$ ${\displaystyle 2N+4}$ ${\displaystyle 4N+4}$ ${\displaystyle M}$

• Por el contrario, utilizar una FFT en un conjunto de datos con valores tiene complejidad . ${\displaystyle N}$ ${\displaystyle O(KN\log _{2}(N))}$

Esto es más difícil de aplicar directamente porque depende del algoritmo FFT utilizado, pero un ejemplo típico es una FFT de base 2, que requiere multiplicaciones y adiciones/sustracciones por contenedor DFT , para cada uno de los contenedores. ${\displaystyle 2\log _{2}(N)}$ ${\displaystyle 3\log _{2}(N)}$ ${\displaystyle N}$

En las expresiones de orden de complejidad, cuando el número de términos calculados es menor que , la ventaja del algoritmo de Goertzel es clara. Pero debido a que el código FFT es comparativamente complejo, el factor de "costo por unidad de trabajo" es a menudo mayor para una FFT, y la ventaja práctica favorece al algoritmo de Goertzel incluso para valores varias veces mayores que . ${\displaystyle M}$ ${\displaystyle \log N}$ ${\displaystyle K}$ ${\displaystyle M}$ ${\displaystyle \log _{2}(N)}$

Como regla general para determinar si una FFT de base 2 o un algoritmo de Goertzel es más eficiente, ajuste la cantidad de términos en el conjunto de datos hacia arriba hasta la potencia exacta de 2 más cercana, llamando a esto , y es probable que el algoritmo de Goertzel sea más rápido si ${\displaystyle N}$ ${\displaystyle N_{2}}$

${\displaystyle M\leq {\frac {5N_{2}}{6N}}\log _{2}(N_{2})}$

Las implementaciones de FFT y las plataformas de procesamiento tienen un impacto significativo en el rendimiento relativo. Algunas implementaciones de FFT ^[11] realizan cálculos internos de números complejos para generar coeficientes sobre la marcha, lo que aumenta significativamente su "costo K por unidad de trabajo". Los algoritmos FFT y DFT pueden utilizar tablas de valores de coeficientes precalculados para lograr una mejor eficiencia numérica, pero esto requiere más accesos a los valores de coeficientes almacenados en la memoria externa, lo que puede generar una mayor contención de caché que contrarreste parte de la ventaja numérica.

Ambos algoritmos ganan aproximadamente un factor de 2 en eficiencia cuando utilizan datos de entrada de valor real en lugar de datos de valor complejo. Sin embargo, estas ganancias son naturales para el algoritmo de Goertzel pero no se lograrán para la FFT sin utilizar ciertas variantes del algoritmo ^{[ ¿cuáles? ]} especializadas para transformar datos de valor real .

Véase también

• Algoritmo FFT de Bluestein (chirp-Z)
• Modulación por desplazamiento de frecuencia (FSK)
• Modulación por desplazamiento de fase (PSK)

Referencias

1. Goertzel, G. (enero de 1958), "Un algoritmo para la evaluación de series trigonométricas finitas", American Mathematical Monthly , 65 (1): 34–35, doi :10.2307/2310304, JSTOR  2310304
2. Mock, P. (21 de marzo de 1985), "Añadir generación y decodificación DTMF a diseños DSP-μP" (PDF) , EDN , ISSN  0012-7515; también se encuentra en Aplicaciones DSP con la familia TMS320, Vol. 1, Texas Instruments, 1989.
3. Chen, Chiouguey J. (junio de 1996), Algoritmo de Goertzel modificado en la detección de DTMF utilizando el DSP TMS320C80 (PDF) , Informe de aplicación, Texas Instruments, SPRA066
4. Schmer, Gunter (mayo de 2000), Generación y detección de tonos DTMF: una implementación utilizando el TMS320C54x (PDF) , Informe de aplicación, Texas Instruments, SPRA096a
5. Cheng, Eric; Hudak, Paul (enero de 2009), Procesamiento de audio y síntesis de sonido en Haskell (PDF) , archivado desde el original (PDF) el 28 de marzo de 2017
6. Gentleman, WM (1 de febrero de 1969). "Análisis de errores del método de Goertzel (Watt) para calcular coeficientes de Fourier". The Computer Journal . 12 (2): 160–164. doi : 10.1093/comjnl/12.2.160 .
7. Stoer, J.; Bulirsch, R. (2002), Introducción al análisis numérico , Springer, ISBN 9780387954523
8. "Algoritmo de Goertzel". Cnx.org. 12 de septiembre de 2006. Consultado el 3 de febrero de 2014 .
9. "Electronic Engineering Times | Conectando a la comunidad electrónica global". EE Times . Consultado el 3 de febrero de 2014 .
10. Elmenreich, Wilfried (25 de agosto de 2011). «Detección eficiente de una frecuencia mediante un filtro de Goertzel» . Consultado el 16 de septiembre de 2014 .
11. Prensa; Flannery; Teukolsky; Vetterling (2007), "Capítulo 12", Recetas numéricas, el arte de la computación científica , Cambridge University Press

Lectura adicional

• Proakis, JG; Manolakis, DG (1996), Procesamiento de señales digitales: principios, algoritmos y aplicaciones , Upper Saddle River, NJ: Prentice Hall, págs. 480–481, Bibcode :1996dspp.book.....P

Enlaces externos

• Algoritmo de Goertzel en Wayback Machine (archivado el 28 de junio de 2018)
• Un algoritmo DSP para el análisis de frecuencia
• El algoritmo de Goertzel de Kevin Banks