Filtro IIR en una MCU de 16 bits (PIC24F)

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);

Simulación Respuesta de frecuencia

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 :)

Respuestas (6)

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.

No resolvió, pero ayudó a comprender mejor el código y los "trucos" del compilador. El problema estaba en los coeficientes. Debería ser a1*xn y b0*y. Un error realmente cojo de hecho eheh ¡Gracias por su ayuda!
+1: técnicamente, el estándar C no define si >> es un cambio aritmético o lógico ... o al menos eso es lo que me dijo la gente de Microchip hace unos años cuando informé que tenía problemas para usar >> para power-of -2-se divide por falta de extensiones de signo.

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 .

Aún mejor, no tiene que definir estos typedefs usted mismo si su sistema tiene un <stdint.h>: pubs.opengroup.org/onlinepubs/007904975/basedefs/stdint.h.html

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.

+1 por usar intNN_t y uintNN_t typedefs, también el enfoque de unión, los cuales uso en mis sistemas.
@JasonS - Gracias. El enfoque de unión es lo que separa a los desarrolladores integrados de los de PC. Los chicos de PC parecen no tener idea de cómo funcionan los números enteros.

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.

Debo tener que usar esta MCU (PIC24F) porque se usará en una aplicación portátil de baja potencia y dsPIC no cumple con ese requisito. El filtro es simple, es de 1er orden. Supongo que mi problema es convertir doble-> formato de punto fijo. Sin embargo, al leer el código, creo que estoy haciendo lo correcto, pero no funciona... Me pregunto por qué.
"utilizado en una aplicación portátil de baja potencia y dsPIC no cumple con ese requisito" <-- este tipo de declaración general hace sonar las alarmas en mi cabeza. El consumo de energía depende de la aplicación y el soporte de hardware a menudo significa un menor consumo de energía general. Cualquier aplicación de bajo consumo mantendrá el procesador dormido la mayor parte del tiempo, por lo que reducir el tiempo de "encendido" es el camino para reducir el consumo de energía...
eso es punto fijo o punto flotante? Flotar sería una ventaja, pero no necesita un DSP para implementar un filtro bi-quad. Necesita mucho más un software de diseño de filtros para encontrar los parámetros.

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);
}