Estoy tratando de implementar un filtro IIR de primer orden simple en una MCU (PIC24FJ32GA002), sin éxito hasta ahora. El filtro es un filtro de seguimiento de CC (filtro de paso bajo) cuyo propósito es rastrear el componente de CC de una señal de 1,5 Hz. La ecuación de diferencia se tomó de una nota de aplicación de TI:
y(n)=K x(n)+y(n-1) (1-K)
con K = 1/2^8
Hice un script de MATLAB para probarlo y funciona bien en la simulación. Código utilizado:
K=1/2^8
b = K
a = [1 -(1-K)]
Fs=200; // sampling frequency
Ts=1/Fs;
Nx=5000; // number of samples
nT=Ts*(0:Nx-1);
fin=1.5; // signal frequency
randn('state',sum(100*clock));
noise=randn(1,Nx);
noise=noise-mean(noise);
xin=200+9*(cos(2*pi*fin*nT));
xin=xin+noise;
out = filter(b,a,xin);
Sin embargo, no puedo implementarlo en un microcontrolador PIC24F. Estoy representando los coeficientes en formato Q15 (1.15), almacenándolos en variables cortas y usando una larga para multiplicaciones. Aquí está el código:
short xn;
short y;
short b0 = 128, a1 = 32640; // Q15
long aux1, aux2;
// (...)
while(1){
xn = readADC(adc_ch);
aux1 = ((long)b0*xn) << 1;
aux2 = ((long)a1*y) << 1;
y = ((aux1 + aux2) >> 16);
delay_ms(5);
}
Long cast se usa para extender la señal para que la operación de multiplicación se realice correctamente. Después de cada multiplicación, cambio un bit a la izquierda para eliminar el bit de señal extendida. Al sumar, cambio 16 bits a la derecha para obtener y en formato Q15.
Estoy depurando la MCU con Pickit2 y la ventana "Ver-> Ver" (MPLAB IDE 8.53) y probando el filtro con una señal de CC (cambio la señal de CC con un potenciómetro para probar diferentes valores). El ADC tiene una resolución de 10 bits y la MCU se suministra con 3,3 V. Algunos resultados:
1V --> xn = 312 (correcto), yn = 226 (incorrecto)
1.5V --> xn = 470 (correcto), yn = 228 (completamente incorrecto)
¿Qué estoy haciendo mal? ¿Alguna sugerencia sobre cómo implementar este filtro IIR en una MCU de 16 bits?
Muchas gracias de antemano :)
No me sumergí demasiado en el diseño de su filtro, pero con solo mirar el código fuente surgen un par de cosas. Por ejemplo, estas líneas:
aux1 = ((long)b0*xn) << 1;
aux2 = ((long)a1*y) << 1;
y = ((aux1 + aux2) >> 16);
El primer problema que veo es el ((long)b0*xn). Me encontré con compiladores que compilarían esto incorrectamente como ((long)(b0*xn)), lo cual es completamente incorrecto. Solo para estar seguro, escribiría esto como (((long)b0)*((long)xn)). Sin duda, esto es programación paranoica, pero...
Luego, cuando haces "<<1", NO es lo mismo que "*2". Para la mayoría de las cosas, está cerca, pero no para DSP. Tiene que ver con cómo la MCU/DSP maneja las condiciones de desbordamiento y firma las extensiones, etc. Pero incluso si funcionó como un *2, está eliminando un bit de resolución que no necesita eliminar. Si realmente tiene que hacer un * 2, entonces haga un * 2 y deje que el compilador descubra si podría sustituir un << 1 en su lugar.
El >>16 también es problemático. En mi cabeza, no sé si va a hacer un cambio lógico o aritmético. Quieres un cambio aritmético. Los cambios aritméticos manejarán el bit de signo correctamente, mientras que un cambio lógico insertará ceros para los nuevos bits. Además, puede ahorrar bits de resolución eliminando el <<1 y cambiando el >>16 por >>15. Bueno, y cambiar todo esto a multiplicaciones y divisiones normales.
Así que aquí está el código que usaría:
aux1 = ((long)b0) * ((long)xn);
aux2 = ((long)a1) * ((long)y);
y = (short)((aux1+aux2) / 32768);
Ahora, no pretendo que esto resuelva su problema. Puede o no, pero mejora su código.
Le sugiero que reemplace "corto" con "int16_t" y "largo" con "int32_t". Guarde "corto", "int", "largo", etc. para lugares donde el tamaño realmente no importa, como los índices de bucle.
Luego, el algoritmo funcionará igual si compila en su computadora de escritorio. Será mucho más fácil depurar en el escritorio. Cuando esté funcionando satisfactoriamente en el escritorio, muévalo al microprocesador. Para hacer esto más fácil, incluya la complicada rutina en su propio archivo. Probablemente querrá escribir una función main() simple para el escritorio que ejerza la rutina y salga de cero (para éxito) o distinto de cero (para falla).
Luego integre la construcción y la ejecución de la prueba en el sistema de construcción de su proyecto. Cualquier sistema de desarrollo que valga la pena considerar le permitirá hacer esto, pero es posible que deba leer el manual. Me gusta GNU Make .
David Kessner mencionó que el >> 16 no es confiable, lo cual es cierto. He sido mordido por eso antes, así que ahora siempre tengo una afirmación estática en mi código que comprueba que >> funciona como lo espero.
Sugirió una división real en su lugar. Pero he aquí cómo hacer una aritmética adecuada >> 16 si resulta que la tuya no lo es.
union
{
int16_t i16s[2];
uint16_t i16u[2];
int32_t i32s;
uint32_t i32u;
}union32;
union32 x;
x.i32s = -809041920;
// now to divide by 65536
x.i16u[0] = x.i16u[1]; // The shift happens in one go.
// Now we do the sign extension
if (x.i16u[0] & 0x8000) // Was this a negative number?
x.i16u[1] = 0xFFFF; // Then make sure the final result is still negative
else
x.i16u[1] = 0x0000; // Otherwise make sure the final result is still positive
// Now x.i32s = -12345
Verifique que su compilador esté generando un buen código para (x.i16u[0] y 0x8000). Los compiladores de Microchip no siempre aprovechan las instrucciones BTFSC de PIC en casos como este. Especialmente los compiladores de 8 bits. A veces me veo obligado a escribir esto en asm cuando realmente necesito el rendimiento.
Hacer el cambio de esta manera es mucho más rápido en los PIC que solo pueden hacer cambios de un bit a la vez.
El uso de un dsPIC en lugar del PIC24 podría facilitar mucho las cosas (son compatibles con pines). El compilador dsPIC viene con una biblioteca DSP que incluye filtros IIR y el hardware admite directamente la aritmética de punto fijo.
Parece que nunca le estás dando valor antes de usarlo. Utiliza y en el cálculo de aux2, lo cual es correcto porque es la muestra anterior filtrada, pero no tiene idea de cómo comenzará en esta situación.
Le sugiero que tome una muestra de ADC antes de ingresar a su ciclo while (1). Asigne ese valor a y, retrase 5 ms, luego ingrese su ciclo. Esto permitirá que su filtro actúe normalmente.
También te recomendaría que sigas adelante y cambies tus pantalones cortos para que sean variables más largos. A menos que se esté quedando sin espacio, tendrá un código mucho más simple si no tiene que escribir cast en cada ciclo.
Muchas de las primeras MCU ni siquiera tienen una instrucción de multiplicación, y mucho menos una división. El filtro IIR con el valor mágico con K = 1/2 ^ 8 fue diseñado específicamente para tales MCU.
// filter implements
// y(n)=K*x(n)+y(n-1)*(1-K)
// with K = 1/2^8.
// based on ideas from
// http://techref.massmind.org/techref/microchip/math/filter.htm
// WARNING: untested code.
// Uses no multiplies or divides.
// Eventually y converges to a kind of "average" of the x values.
short average; // global variable
// (...)
while(1){
short xn = readADC(adc_ch);
// hope that we don't overflow the subtract --
// -- perhaps use saturating arithmetic here?
short aux1 = xn - average;
short aux2 = aux1 / 256;
// compiler should use an arithmetic shift,
// or perhaps an unaligned byte read, not an actual divide
average += aux2;
delay_ms(5);
}
rasgo
jason s