Estoy tratando de implementar un filtro de corte de banda (notch) en un microcontrolador dsPIC33EP64GS506
para filtrar un 100 Hz
componente de una señal de entrada. El problema aquí es que mi microcontrolador no tiene una unidad de punto flotante, así que tengo que utilizar una aritmética de punto fijo, pero para ser más preciso, utilizo aritmética entera. El tiempo de muestreo del sistema es \$T_s=50~\mu\text{s}\$ . Aquí explico mi problema en detalle, y las preguntas están al final del post. Adjunto un código MATLAB para realizar simulaciones si lo deseas: enlace de descarga en mi Dropbox . Ten en cuenta que no tienes que ser un usuario registrado de Dropbox para poder descargar el archivo.
La resolución del ADC es 12b
que "aumento" a 15b
no para aumentar la resolución de la medición (que no puede hacerse), sino para aumentar la resolución del filtro:
v_in = ADCBUF0<<3;
La función de transferencia del filtro notch en el dominio s viene dada por:
$$G_s(s) = \frac{s^2 + 2\zeta_1\omega_n s + \omega_n^2}{s^2 + 2\zeta_2\omega_ns + \omega_n^2},$$
donde los parámetros del filtro \$\zeta_1\$ , \$\zeta_2\$ y \$\omega_n\$ determinar completamente el filtro. Si queremos filtrar un 100 Hz
entonces fijamos \$\omega_n=2\pi\cdot100~\text{s}^{-1}\$ . En cuanto a otros parámetros, sin explicaciones detalladas, utilizo \$\zeta_1=0.001\$ y \$\zeta_2=1\$ . Para implementar este filtro en un sistema digital, necesitamos discretizar la función de transferencia en el dominio s, y para ello, utilizo un Método de discretización de Tustin .
$$G_z(z) = \frac{b_0 + b_1 z^{-1} + b_2 z^{-2}}{1 + a_1 z^{-1} + a_2 z^{-2}}.$$
La ecuación recursiva correspondiente, implementada en un sistema digital, es:
$$y(k) = b_0 u(k) + b_1 u(k-1) + b_2 u(k-2) - a_1 y(k-1) - a_2 y(k-2).$$
Este es el código 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);
A continuación se muestran las características de frecuencia del filtro en el dominio z generadas por MATLAB. Como puede verse, este tipo de filtro es muy sensible en términos de respuesta en frecuencia. Incluso los cambios más pequeños en los parámetros podrían causar un comportamiento completamente diferente, por ejemplo, una ganancia inesperada, desplazamiento de fase, etc.
Ahora bien, como no dispongo de una unidad de coma flotante, utilizo un método bien conocido llamado escalado binario. Cualquier número decimal \$ d \in \mathbb{R} \$ puede representarse como \$\frac{\lfloor d\cdot2^r \rceil}{2^r}\$ donde \$r\in\mathbb{N}\$ y \$\lfloor\cdot\rceil\$ es una función redonda al entero más cercano. El escalado binario es un método utilizado cuando queremos evitar el uso de una división, que es una operación muy costosa en términos de ciclos de CPU requeridos (normalmente 18 ciclos), ya que el desplazamiento binario puede producir los mismos resultados en muchos menos ciclos (normalmente 1 ciclo). Para este ejemplo he utilizado \$r=15\$ que es la máxima precisión que podría utilizar teniendo en cuenta los bits disponibles para realizar una multiplicación (un resultado debe estar dentro de los 32 bits). La ecuación recursiva correspondiente utilizando 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 redondea al entero más próximo después de la operación de desplazamiento de bits. Tenga en cuenta que >>15
es prácticamente una división entera por 32768
pero somos conscientes de que no es totalmente equivalente. El desplazamiento de bits siempre redondea a menos infinito, mientras que la división de enteros siempre redondea a cero. Los valores de los parámetros del filtro son los siguientes: B0=31771
, B1=-63509
m B2=31769
, A1=-63509
, A2=30772
. Esto también se puede encontrar en un código MATLAB proporcionado.
Me sorprendió mucho cuando me di cuenta de que la variante entera del filtro no funciona en absoluto. Por favor, vea a continuación las respuestas de tres filtros diferentes utilizados para eliminar un 100 Hz
componente de la señal de entrada. De arriba a abajo: (1) filtro notch en la implementación de coma flotante, (2) filtro notch en la implementación de aritmética entera utilizando \$r=15\$ (3) un filtro simple de media móvil.
He aquí cómo he implementado un filtro de media 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 estos resultados son los esperados? ¿Es realmente tan problemático implementar un filtro notch en un sistema digital utilizando sólo aritmética de enteros? ? ¿Qué filtro digital se utiliza normalmente para eliminar una frecuencia determinada? ? Como veo en estos resultados, el filtro de media móvil supera a ambas implementaciones del filtro notch. ¿Hay alguna desventaja del filtro de media móvil que deba tener en cuenta, excepto, obviamente, el aumento de la demanda de memoria? ?