El algoritmo de Goertzel permite calcular un coeficiente aislado de la DFT sobre un conjunto de N muestras con una complejidad temporal de O(n) y una complejidad espacial de O(1), además N no tiene por qué ser potencia de 2. Este algoritmo se basa en la aplicación de una ecuación en diferencias finitas (un filtro IIR).
A mayor N, mayor resolución en frecuencia y también mayor latencia en la detección de los tonos. Es necesario, por tanto encontrar un compromiso. La resolución en hercios del algoritmo de Goertzel viene dada por:
$$R = {f_s \over N}$$
Siendo $f_s$ la frecuencia de muestreo y $N$ el número de muestras que se evalúan en cada pasada.
1. Se calcula el índice del coeficiente correspondiente
Lo primero que hay que hacer es calcular el índice del coeficiente asociado a la DFT a partir de la frecuencia que queremos detectar:
$$k = N{f \over f_s}$$
Siendo $f$ la frecuencia que queremos detectar, $f_s$ la frecuencia de muestreo y $N$ el número de muestras que procesamos cada vez. La máxima frecuencia que podemos detectar será la mitad de la frecuencia de muestreo.
2. Para cada conjunto de N muestras
2.1. Se aplica la ecuación en diferencias
Para cada una de las muestras que van llegando, vamos aplicando la siguiente ecuación en diferencias finitas:
$$s[n] = x[n] + 2\cos\left({2 \pi k \over N}\right)s[n-1] - s[n-2]$$
$$n = 0..N$$
Con codiciones iniciales $s[-1] = s[-2] = 0$. Se trata, como se puede ver, de un sencillo filtro IIR de segundo orden.
2.2. Se obtiene del coeficiente k-esimo de la DFT
Se puede demostrar que el coeficiente k-ésimo de la DFT de tamaño N es:
$$X(k) = s[N] - W_N^ks[N-1]$$
Siendo:
$$W_N = e^{-j\left({2 \pi \over N}\right)}$$
Por tanto:
$$X(k) = s[N] - e^{-j\left({2 \pi \over N}\right)k}s[N-1]$$
Si desarrollamos la exponencial compleja mediante la fórmula de Euler tenemos que:
$$e^{-j\left({2 \pi \over N}\right)k} = \cos\left({2 \pi \over N}\right) - j \sin\left({2 \pi k \over N}\right)$$
Y, por tanto:
$$X(k)_{real} = s[N] - \cos\left({2 \pi k \over N}\right)s[N-1]$$
$$X(k)_{imag} = \sin\left({2 \pi k \over N}\right)s[N-1]$$
Para calcular la magnitud de la banda de frecuencia correspondiente calculamos el módulo de $X(k)$:
$$M^2 = \left|X(k)\right|^2 = X(k)_{real}^2 + X(k)_{imag}^2$$
De esta forma podemos medir la magnitud de la banda de frecuencia correspondiente al coeficiente k-ésimo de la DFT o, lo que es lo mismo, la magnitud de la banda de frecuencia correspondiente a la frecuencia $f$.
$$k = N{f \over f_s} \Rightarrow f = {k f_s \over N}$$
Implementación
A continuación se puede ver una implementación del algoritmo de Goertzel en Arduino:
const int ANALOG_INPUT = A0; const int SAMPLE_RATE_HZ = 3000; const int SAMPLE_PERIOD_US = 1000000 / SAMPLE_RATE_HZ; const int F_HZ = 440; const int N = 3000; const int K = 440; // N * F_HZ / SAMPLE_RATE_HZ; unsigned long tPrev; struct goertzelFilter { float s1, s2; int k, n; float sinv, cosv, cosv2; int nextIteration; }; struct goertzelFilter filter; void goertzelFilterReset(struct goertzelFilter &f) { f.s1 = 0; f.s2 = 0; f.nextIteration = 0; } void goertzelFilterInit(struct goertzelFilter &f, int k, int n) { f.k = k; f.n = n; f.cosv = cos(2 * PI * k / n); f.sinv = sin(2 * PI * k / n); f.cosv2 = 2 * f.cosv; goertzelFilterReset(f); } boolean goertzelFilterFinished(struct goertzelFilter &f) { return (f.nextIteration == (f.n + 1)); } void goertzelFilterRun(struct goertzelFilter &f, float input) { float s = input + (f.cosv2 * f.s1) - f.s2; f.s2 = f.s1; f.s1 = s; f.nextIteration++; } float goertzelFilterGetMagnitude(struct goertzelFilter &f) { float real = f.s1 - (f.cosv * f.s2); float imag = f.sinv * f.s2; return ((real * real) + (imag * imag)); } void setup() { Serial.begin(9600); tPrev = micros(); goertzelFilterInit(filter, K, N); goertzelFilterReset(filter); } void loop() { unsigned long t = micros(); if ((t - tPrev) >= SAMPLE_PERIOD_US) { int v = analogRead(ANALOG_INPUT); if (goertzelFilterFinished(filter)) { float magnitude = goertzelFilterGetMagnitude(filter); Serial.println(magnitude); goertzelFilterReset(filter); } else goertzelFilterRun(filter, ((float) v - 512) / 512); tPrev = t; } }
Se ha elegido una frecuencia de muestreo baja (3000Hz) para poder trabajar cómodamente con tipos float. Utilizando aritmética de punto fijo podríamos incremenentar la frecuencia de muestreo y que la detección sea más precisa.
La entrada de audio se toma de la entrada analógica A0 a la que se conecta un sencillo circuito amplificador:
Elegir el valor de N
A la hora de elegir la N lo ideal es escogerla lo más grande posible, que nos permita una latencia razonable y que se cumpla que:
$$k = N{f \over f_s} \in \mathbb{N}$$
Por ejemplo, para detectar un tono de 1Khz sobre una señal muestreada a 6KHz lo ideal sería que la N valiese: 96, 102, 114… Ya que para todos estos valores se cumple que $k$ es número natural.
Magnitudes medidas para diferentes tonos
En ausencia de señal de entrada: 0.30, 0.37, 0.30, 0.35, 0.44, 0.32, 0.37, 0.28...
Con una señal de entrada de 200Hz: 2.56, 1.73, 3.56, 0.81, 0.58, 1.67, 4.71, 6.81...
Con una señal de entrada de 440Hz (la frecuencia del detector): 138.41, 87.29, 441.14, 185.20, 233.03, 762.27, 80.62, 330.98...
Con una señal de entrada de 500Hz: 25.70, 9.60, 22.50, 2.76, 16.62, 18.75, 23.56, 35.58...
[ añadir comentario ] ( 2532 visualizaciones ) | [ 0 trackbacks ] | enlace permanente
Tweet | |
( 3.1 / 3185 )Hace un tiempo encontré en el siguiente enlace un interesantísimo artículo explicando las diferentes técnicas que existen para la detección de BPM (tempo) sobre pistas de audio.
http://www.flipcode.com/misc/BeatDetect ... rithms.pdf
En particular me llamó la atención la técnica basada en el cálculo de la correlación de la señal con un tren de pulsos. En la sección soft puede descargarse una implementación sencilla en C de este algoritmo. Da muy buenos resultados con música pop, rock y electrónica. No he probado el código con otros estilos musicales.
Código fuente aquí.
[ añadir comentario ] ( 1160 visualizaciones ) | [ 0 trackbacks ] | enlace permanente
Tweet | |
( 3 / 1848 )El algoritmo HPS (Harmonic Product Spectrum) permite detectar la frecuencia fundamental (la altura) de una nota tocada por cualquier tipo de instrumento armónico.
Asume que el espectro generado por el instrumento está formado por frecuencias múltiplas enteras de la frecuencia fundamental, lo cual es cierto en el 99% de los instrumentos musicales convencionales, incluyendo la voz.
En la sección soft he colgado una sencilla implementación en C de este algoritmo y la he acompañado de dos ejemplos de tonos generados por un saxo. Uno de los ejemplos es una nota LA que el algoritmo detecta de 440.0848 Hz (440 Hz es la frecuencia de la nota LA en la 4ª octava) y el otro es de una nota RE que el algoritmo detecta de 293.3899 Hz (293.665 Hz es la frecuencia de la nota RE en la 4ª octava).
Más info sobre el HPS y otros algoritmos de detección de pitch en http://cnx.org/content/m11714/latest/.
EDITADO: He corregido un bug en el código que provocaba que la detección estuviera desplazada una octava.
Sección soft.
[ añadir comentario ] ( 1176 visualizaciones ) | [ 0 trackbacks ] | enlace permanente
Tweet | |
( 3 / 3463 )Sonority es un motor de síntesis de sonido mediante modelado analógico implementado en ANSI C que utiliza exclusivamente aritmética de punto fijo (formato Q16.16, 16 bits de parte entera y 16 bits de parte fraccionaria, 32 bits en total).
En esta nueva versión se han añadido dos modos de portamento (ALWAYS y LEGATTO), dos modos de disparo de envolvente (ALWAYS y STACATTO) y tablas precalculadas (tanto la tabla de frecuencias como las tablas de ondas).
El fichero test.c contiene un ejemplo de utilización del motor.
+info y descargas aquí.
[ añadir comentario ] ( 1369 visualizaciones ) | [ 0 trackbacks ] | enlace permanente
Tweet | |
( 3 / 2372 )He encontrado un artículo muy interesante que habla sobre la FHT (Fast Hartley Transform), para los neófitos como yo: algo así como la FFT pero con números reales :-).
http://www.embedded.com/2000/0009/0009feat3.htm
El artículo incluye una implementación en C para los que tienen prisa ;-).
Saludos y feliz Navidad a todos.
[ añadir comentario ] ( 1145 visualizaciones ) | [ 0 trackbacks ] | enlace permanente
Tweet | |
( 3 / 2341 )
Calendario
