Problemas con la implementación de un filtro de supresión de banda en una MCU (dsPIC) utilizando aritmética de punto fijo

Estoy tratando de implementar un filtro de parada de banda (muesca) en un microcontrolador dsPIC33EP64GS506para filtrar un 100 Hzcomponente de una señal de entrada. El problema aquí es que mi microcontrolador no tiene una unidad de punto flotante, así que tengo que usar una aritmética de punto fijo, pero para ser más preciso, uso aritmética de enteros. El tiempo de muestra del sistema es$T_s=50~\mu\text{s}$. Aquí explico mi problema en detalle, y las preguntas están al final de la publicación. Encuentre adjunto un código MATLAB para ejecutar simulaciones si lo desea: enlace de descarga en mi Dropbox . Tenga en cuenta que no tiene que ser un usuario registrado de Dropbox para poder descargar el archivo.

La resolución de ADC es 12b, a la que "aumento" 15b, no para aumentar la resolución de medición (que no se puede hacer), sino para aumentar la resolución del filtro:

v_in = ADCBUF0<<3;

La función de transferencia de filtro de muesca en el dominio s se da como:

donde filtrar parámetros$\zeta_1$,$\zeta_2$, y$\omega_n$determinar completamente el filtro. Si queremos filtrar un 100 Hzcomponente, entonces establecemos$\omega_n=2\pi\cdot100~\text{s}^{-1}$. En cuanto a otros parámetros, sin explicaciones detalladas, uso$\zeta_1=0.001$y$\zeta_2=1$. Para implementar este filtro en un sistema digital, necesitamos discretizar la función de transferencia en s-domain, y para eso utilizo un método de discretización de Tustin .

La ecuación recursiva correspondiente, implementada en un sistema digital, es:

Aquí está el código de MATLAB para obtener la función de transferencia en el dominio z:

Ts = 50e-6;
zeta1 = 0.001;
zeta2 = 1;
wn = 100*(2*pi);
Gs = tf([1 2*zeta1*wn wn^2], [1 2*zeta2*wn wn^2]);
Gz = c2d(Gs, Ts, 'tustin');
bodeplot(Gz);

Consulte a continuación las características de frecuencia del filtro en el dominio z generado por MATLAB. Como se puede ver, este tipo de filtro es muy sensible en términos de respuesta de frecuencia. Incluso los cambios más pequeños en los parámetros podrían causar un comportamiento completamente diferente, por ejemplo, una ganancia inesperada, cambio de fase, etc.

Ahora, como no tengo una unidad de punto flotante a disposición, utilizo un método muy conocido llamado escalado binario. Cualquier número decimal$d \in \mathbb{R}$se puede representar como$\frac{\lfloor d\cdot2^r \rceil}{2^r}$, dónde$r\in\mathbb{N}$, y$\lfloor\cdot\rceil$es una función de redondeo al entero más cercano. El escalado binario es un método que se usa cuando queremos evitar el uso de una división, que es una operación muy costosa en términos de ciclos de CPU necesarios (normalmente 18 ciclos), ya que el desplazamiento binario puede producir los mismos resultados en muchos menos ciclos (normalmente 1 ciclo). Para este ejemplo usé$r=15$, que es la precisión más alta que podría usar considerando los bits disponibles para hacer una multiplicación (un resultado debe estar dentro de los 32 bits). La ecuación recursiva correspondiente usando aritmética de enteros se implementa de la siguiente manera:

yk0 = ( B0*uk0+B1*uk1+B2*uk2 - A1*yk1-A2*yk2 + (1<<14) ) >> 15;
uk2 = uk1; uk1 = uk0;
yk2 = yk1; yk1 = yk0;

El término (1<<14)se utiliza para garantizar que el resultado se redondee al entero más cercano después de la operación de cambio de bit. Tenga en cuenta que >>15es prácticamente una división entera por 32768, pero sabemos que no es completamente equivalente. El cambio de bit siempre redondea a menos infinito, mientras que la división de enteros siempre redondea a cero. Los valores de los parámetros de filtro son los siguientes: B0=31771, B1=-63509m B2=31769, A1=-63509, A2=30772. Esto también se puede encontrar en un código de MATLAB proporcionado.

Me sorprendió mucho cuando me di cuenta de que la variante entera del filtro no funciona en absoluto. Consulte a continuación las respuestas de tres filtros diferentes utilizados para eliminar un 100 Hzcomponente de la señal de entrada. De arriba a abajo: (1) filtro de muesca en la implementación de punto flotante, (2) filtro de muesca en la implementación aritmética de enteros usando$r=15$, (3) un filtro de promedio móvil simple.

Así es como implementé un filtro de promedio móvil en C:

uk200 = window[ind];
window[ind] = uk0;
win_sum = win_sum - uk200 + uk0;
yk0 = (((win_sum+(1<<3))>>4)*5243+(1<<15))>>16;
ind++;
if (ind==200) ind=0;

Preguntas a la comunidad.

Como no tengo mucha experiencia en filtrado digital, ¿alguien puede confirmarme si se esperan estos resultados? ¿Es realmente problemático implementar un filtro de muesca en un sistema digital usando solo aritmética de enteros ? ¿Qué filtro digital se utiliza normalmente para eliminar una determinada frecuencia ? Como veo en estos resultados, el filtro de media móvil supera a las dos implementaciones de filtro de muesca. ¿Hay alguna desventaja para el filtro de promedio móvil que deba tener en cuenta, excepto por la demanda de memoria obviamente aumentada ?