Filtrado IIR de orden 4 en forma directa II

CHRISTIAN BOCERO TOLEDANO

 

  1. Introducción y Cálculos previos

Los Filtros IIR son sistemas cuya salida depende además de salidas anteriores y que, estando en reposo, al ser estimulados con una entrada impulsional su salida no vuelve al reposo, de ahí el  calificativo de filtros de respuesta impulsional infinita (IIR).

Dibuje la respuesta en frecuencias analógicas del filtro a implementar y compruebe que todos los coeficientes son menores o iguales que la unidad.

 

» fs=32000;  %frec. de muestreo del enunciado

 

» fc=7000;    %frec. de corte (donde la respuesta en frecuencia cae 3 dB)

» fd=fc/fs     %%frec. discreta normalizada

fd =   0.2188

» wn=fd*2  %Es necesario multiplicar por 2 porque para Matlab wn va entre 0 y 1 y no entre 0 y 0.5 como en TDS

wn =  0.4375

» [B,A]=butter(4,wn) %Aquí ya tenemos los coeficientes del filtro de orden 4

B =  0.0618    0.2471    0.3707    0.2471    0.0618

A =   1.0000   -0.4884    0.5615   -0.1069    0.0224

 

» [H,W]=freqz(B,A,2048);

» plot(W/(2*pi),20*log10(abs(H)),’r’),grid on

He puesto format long para calcular más decimales y haciendo

»ginput(1)

ans =   0.218811670214084  -3.017977778642544

Observamos que a -3dB la frecuencia normalizada del filtro es la que tiene que ser con lo cual los cálculos están bien hechos.

 

 

 

 

 

 

 

 

 

» plot(W/(2*pi),abs(H)),grid on

 

 

 

 

 

 

 

 

 

 

»  plot(fs*W/2/pi,(abs(H))),grid on %Otra forma de verlo

 

 

 

 

 

 

 

 

 

 

» cerosH=roots(B);

» polosH=roots(A);

» zplane(cerosH,polosH);

» max(A)  ans =   1        » max(B) ans =  0.3707

 

 

 

 

 

 

  1. Estudio teórico en Matlab, Genere el fichero con los coeficientes hexadecimales.

% Llegados a este punto tenemos que cuantificar los coeficientes con 16 bits para meterlos en el DSP ya que el DSP divide los coeficientes por una constante.

Para conseguir que el DSP guarde los coeficientes con los valores que nos han salido debemos multiplicarlos previamente por dicha constante y además redondear el resultado para un correcto funcionamiento del DSP .

» coefsB=B*(2^15)

coefsB =  1.0e+004 *

0.2025    0.8098    1.2147    0.8098    0.2025

» coefsA=A*(2^15)

coefsA = 1.0e+004 *

3.2768   -1.6004    1.8400   -0.3504    0.0733

» coefsB=round(coefsB)

coefsB =    2025        8098       12147        8098        2025

» coefsA=round(coefsA)

coefsA =   32768      -16004       18400       -3504         733

» coefA=coefsA(2:5) %porque recuerda que en los filtros no existe el coeficiente a0 pero si el b0

coefA =     -16004       18400       -3504         733

Ahora ya tenemos los coeficientes que buscábamos, aunque estos no están en Hexadecimal. El próximo paso será guardar los coeficientes en Hexadecimal. Para ello lo primero que hay que hacer es guardarlos en formato ASCII.

 

 

 

» save coefsA.txt coefA –ascii

-1.6004000e+004  1.8400000e+004 -3.5040000e+003  7.3300000e+002

» save coefsB.txt coefsB –ascii

2.0250000e+003  8.0980000e+003  1.2147000e+004  8.0980000e+003  2.0250000e+003

Ahora que los tenemos  en ASCII sólo nos falta la conversión a Hexadecimal. Utilizando el MS-DOS de Windows hacemos:

asc2hex coefsA.txt coefsA.hex

c17d

47e0

f251

02dd

asc2hex coefsB.txt coefsB.hex

07e9

1fa2

2f73

1fa2

07e9

Ya hemos obtenido los coeficientes en Hexadecimal!

 

 

 

 

 

 

  1. Programación en DSP

3.1)  Ajuste la frecuencia de muestreo, compile el programa y ejecútelo.

Vamos a modificar el archivo fir.dsp, para adecuarlo a nuestro filtro IIR de orden 4, (obviamos las partes de código idénticas al fir y adjuntamos el archivo IIR4B.dsp donde está el código completo para probar en el DSP)

{Inicializamos los buffers y las variables}

.var/dm/ram/circ                    muestrasl[5];

.var/dm/ram/circ                    muestrasr[5];

.var/pm/ram/circ                    tapsA[4];

 

.var/pm/ram/circ                    tapsB[5];

{Cambiar “fich.hex” por el nombre del fichero con los coeficientes }

{en hexadecimal }

.init tapsA: <COEFSA.hex>;

.init tapsB: <COEFSB.hex>;

{Es fundamental darse cuenta de cambiar la frecuencia de muestreo del DSP y ponerle la frecuencia de muestreo para lo que hemos hecho los cálculos}

0xc856,     {  * LAB_TDS * Cambiar aqui la frecuencia de muestreo

modificando el £ltimo caracter hexadecimal

segun la tabla de abajo.

######## Frecuencias de muestreo ################

b0-3: 0=  8.

1=  5.5125

2= 16.

3= 11.025

4= 27.42857

5= 18.9

                       6= 32.

7= 22.05

8=   .

9= 37.8

a=   .

b= 44.1

c= 48.

d= 33.075

e=  9.6

f=  6.615

{Inicialización de los componentes y punteros en el DSP}

ax0=0;

dm(stat_flag)=ax0;

ax0=l3;

i2=^muestrasl; {Apuntamos al principio de las muestras}

i3=^muestrasr; {lo defino, pero en el primer apartado no se usa, se usa en la ampliación}

i4=^tapsA; {Apuntamos al Primer  coeficiente de A}

i5=^tapsB; {Apuntamos al Primer  coeficiente de B}

m2=1; {para aumentar en uno la dirección del puntero i2}

m4=-1; {para disminuir en uno la dirección de los punteros i4 e i5}

{Recuerda que por temas del DAG i0,i1,i2 e i3 están en el DAG1 con m0,m1,m2,m3 y i4,i5,i6 e i7 están en el DAG2 con m4,m5,m6,m7 y no puedes usar m4 con i2 por ejemplo etc.}

 

l2=%muestrasl;  {longitud del buffer de muestras del canal izq.}

l3=%muestrasr; {longitud del buffer de muestras del canal derecho}

l4=%tapsA; {longitud del buffer de los coefs de A}

l5=%tapsB; {longitud del buffer de los coefs de B}

ax0=l4;  {metemos la longitud del buffer de coefs. de A qué será de tamaño 4}

ay0=1; {cargamos un uno en ay0 para luego restárselo a ax0 y calcular las vueltas que va a dar el bucle}

ar=ax0-ay0;

dm(tapsA_1)=ar; {En tapsA_1 he calculado el número de veces que tiene que hacer el bucle para calcular la parte A del IIR}

ax1=pm(i4,m4); {Esto es Superimportante para que i4 que al principio hemos dicho que apuntaba al primer coef a1 ahora apunte a a4 hay que hacer esta línea. En ax1 se guardará el coeficiente a1 que será un número pero que no nos sirve para nada porque lo que queríamos con esta línea es colocar el puntero i4 en el coeficiente a4}

ax0=l5;

ay0=1;

ar=ax0-ay0;

dm(tapsB_1)=ar; {Hacemos lo mismo para calcular el número de vueltas que va a dar en el segundo bucle para calcular la parte B del IIR}

ax1=pm(i5,m4); {Hacemos igual que antes, colocamos el puntero i5 que al principio apunta siempre al primer sitio, en este caso a b0, pues con esta línea apuntará a b4 de tal forma que lo que se guarde en ax1, que se guardará b0 nos da igual, lo que de verdad interesa es que te quedas apuntando al coef b4}

ax0=0;

set fl1;

cntr=%muestrasl;

 

 

 

{Por último en input_samples es donde calculamos los bucles y las salidas}

input_samples:

{      i2=^muestrasl;     i2: Apunta a la memoria del canal L;

i3=^muestrasr;     i3: Apunta a la memoria del canal R;

i4=^tapsA               i4: Apunta a los coefic. de A del filtro.

i5=^tapsB                     i5: Apunta a los coefic. de B del filtro.

m2=1;                      m2: Para incrementar direcci¢n;

m4=-1;                          m4: Para decrementar dirección

i0,i1,m0,m1,l0,l1: NO se deben usar pues los usa el autobuffering de

salida  y  entrada}

{ CANAL L}

ax0 = dm (rx_buf + 1);     {Leer nueva muestra de A/D}

1à        mr=0, mx0=dm(i2,m2), my0=pm(i4,m4); {Inicializo el bucle, pongo a cero el acumulador de sumas, i2 apunta a la muestra mas antigua pero i4 apunta al coef a4. Esto es como la preparación de lo que luego va a multiplicar}

cntr=dm(tapsA_1);

do filt_l1 until ce;

2à  filt_l1:         mr=mr-mx0*my0 (ss),mx0=dm(i2,m2),my0=pm(i4,m4);

{Efectivamente, primero multiplico y luego actualizo los coeficientes y muestras en los registros mx0 y my0 y por último actualizo los punteros, de tal forma que cuando vuelva a entrar al bucle lo primero que hace es multiplicar lo anterior (mx0 y my0)}

3à       mr=mr-mx0*my0 (ss);

ay0=mr1; {guardo el resultado anterior}

ar=ax0+ay0;  {a las “sumas”(o restas en este caso) anteriores, les debemos sumar la muestra de entrada que guardé previamente en ax0}

4à       dm(i2,m2)=ar;  {Guardo en V0 el resultado anterior para empezar con los coefs. de B, i2 apunta a la muestra más antigua e i4 apuntará al coef. a4 }

 

 

 

 

 

 

 

 

 

 

 

Traza del programa hasta aquí y primer bucle:

1à                         2à                 2à                 2à                    3à                            4à

mr=0 mr=mr-a4*a4 mr=mr-V3*a3 mr=mr-V2*a2 mr=mr-V1*a1 V0=X+mr
mx0=V4 mx0=V3 mx0=V2 mx0=V1 mx0=V1
my0=a4 my0=a3 my0=a2 my0=a1 my0=a1
i2=V3 i2=V2 i2=V1 i2=V0 i2=V0 i2=V4
i4=a3 i4=a2 i4=a1 i4=a4 i4=a4 i4=a4

 

{Segundo Bucle: Coefs. B del filtro}

 

1à        mr=0, mx0=dm(i2,m2), my0=pm(i5,m4); {Inicializo el bucle, pongo a cero el acumulador de sumas, i2 apunta a la muestra más antigua (V4) pero i5 apunta al coef. b4. Esto es como la preparación de lo que luego va a multiplicar}

cntr=dm(tapsB_1); {guardo en el contador cntr el numero de vueltas del bucle}

2à        do filt_l2 until ce; {Realizamos el bucle hasta que el contador expire, hasta que cntr valga cero}

filt_l2: mr=mr+mx0*my0 (ss), mx0=dm(i2,m2),my0=pm(i5,m4);

3à        mr=mr+mx0*my0 (ss);

 

 

ax0=dm(i2,m2); {Como la longitud del vector de muestras tiene una componente más que el vector de coeficientes A, debemos aumentar el puntero del vector de muestras en 1 posición más}

dm (tx_buf + 1) = mr1; {Enviar resultado a D/A}

rti;

.endmod;

 

 

 

 

 

 

 

 

 

 

 

 

1à                          2à                2à                 2à                 2à                 3à

mr=0 mr=mr+V4*b4 mr=mr+V3*b3 mr=mr+V2*b2 mr=mr+V1*b1 mr=mr+V0*b0
mx0=V4 mx0=V3 mx0=V2 mx0=V1 mx0=V0
my0=b4 my0=b3 my0=b2 my0=b1 my0=b0
i2=V3 i2=V2 i2=V1 i2=V0 i2=V4 i2=V4
i5=b3 i5=b2 i5=b1 i5=b0 i5=b4 i5=b4

 

 

 

 

 

3.2) Proceda ahora a la medida de la respuesta en frecuencia del filtro diseñado. Para ello utilice el generador de funciones, y el osciloscopio conectados como se muestra en la figura 3. Confeccione una curva con la respuesta en frecuencia práctica y compárela con la teórica calculada con Matlab.

Para crear la curva de respuesta en frecuencia de forma práctica. En primer lugar hemos conectado la salida del generador de funciones a la entrada del canal izquierdo del DSP y la salida del canal izquierdo del DSP a la entrada del canal uno del osciloscopio y hemos ido variando la frecuencia con el generador de funciones y comprobando con el osciloscopio como se va reduciendo la amplitud de la señal conforme aumentaba la frecuencia. Una vez anotados los datos en una hoja, hemos representado la gráfica en Matlab con los siguientes comandos:

 

» ap = [0.29 0.29 0.29 0.29 0.28 0.26 0.2 0.18 0.17 0.15 0.14 0.12 0.1 0.09 0.08 0.07 0.06 0.05 0.04 0.04 0.035 0.03 0.006 0.003 0 0 0 0];

» f = [1 2 3 4 5 6 7 7.2 7.4 7.6 7.8 8 8.2 8.4 8.6 8.8 9 9.2 9.4 9.6 9.8 10 11 12 13 14 15 16];

>> f=f*1e3;

» H1=ap./ap(1);

» plot(f,H1),grid on

» hold on

» plot(fs*W/2/pi,(abs(H)),’m’),grid on

»  hold off

 

 

 

 

 

 

 

 

 

 

 

Así pues, obtenemos finalmente una gráfica donde representamos la curva obtenida de forma práctica (azul) y la curva obtenida de forma teórica (magenta).

En la foto siguiente se pretende enseñar el montaje y el funcionamiento del DSP con una señal sinusoidal de 1Khz,  el filtro la deja pasar ocupando en el osciloscopio 4 divisiones verticales, como se puede observar haciendo zoom a esta foto.

 

 

 

 

 

 

 

 

 

 

 

Como es lógico cuando meta la señal de 7Khz que es la frec. de corte del filtro, en el osciloscopio la señal ocupará 4divisiones*0.7= 2.8 divisiones, ya que en la frecuencia de corte del filtro la amplitud de la señal se atenúa 3dB (0.7 en lineal porque es      10^(-3/20)) como se puede observar en la siguiente foto.

 

 

 

 

 

 

 

 

 

En la siguiente foto, vemos como al ir aumentando la frecuencia de la señal de entrada, el filtro va atenuando la amplitud de la señal de esta.

 

 

 

 

 

 

 

 

 

 

Para los más puristas, se puede demostrar que sale correctamente haciendo un rápido cálculo en Matlab.

» plot(fs*W/(2*pi),20*log10(abs(H)),’r’),grid on

>> ginput(1)

ans =  1.0e+003 *

9.0000   -0.0139

>> 10^(-13.9/20) %paso a lineal

ans =    0.2018

Y efectivamente 4 divisiones del original multiplicado por 0.2018 igual a 0.8073 divisiones, que es lo que salía en la foto anterior.

 

 

 

 

 

  1. Ampliación.

4.1) Determine la amplitud máxima de entrada Amax para que no se produzca saturación. Aplique un factor de ganancia a las muestras de entrada de forma que se guarden en el buffer circular ya multiplicadas por ese factor. Tenga en cuenta que el MAC divide por 2^15 al realizar la multiplicación.

Tenemos que determinar el valor máximo de señal de entrada para no se produzca saturación

>> hxa=impz(1,A); %Calculo la respuesta al impulso de 1/A

>> Sxa=sum(abs(hxa))

Sxa =    2.5388

>> Amax_a=1/Sxa

Amax_a =   0.3939

>> hxy=impz(B,A);  %Calculo la respuesta al impulso de B/A

>> Sxy=sum(abs(hxy))

Sxy =    1.4839

>> Amax_y=1/Sxy

Amax_y =    0.6739

>> Amax=min(Amax_a,Amax_y); %El Amax más restrictivo es el mínimo de los 2.

Amax =

0.3939

 

Así pues, lo único que nos queda es aplicar el factor de ganancia a las muestras de entrada para que se guarden en el buffer circular ya multiplicadas. Para hacer esto escribimos en el código las siguientes líneas:

 

.const amax= 12907; {declaramos como constante el valor de amplitud máximo ya multiplicado por 2^15}

mx0 = dm (rx_buf + 1); {Leer nueva muestra de A/D}

my0 = dm (amax);

mr= mx0*my0 (ss);

ax0=mr1; {Guardo la muestra de entrada multiplicada por la constante}

mr=0, mx0=dm(i2,m2), my0=pm(i4,m4); cntr=dm(tapsA_1);

do filt_l1 until ce;

filt_l1: mr=mr-mx0*my0 (ss),mx0=dm(i2,m2),my0=pm(i4,m4);

mr=mr-mx0*my0 (ss);

ay0=mr1; {Guardo el resultado anterior}

ar=ax0+ay0;

dm(i2,m2)=ar; {Igual que antes, acabado el bucle guardo en V0 el resultado}

 

{Segundo Bucle: Coefs. B del filtro}

 

mr=0, mx0=dm(i2,m2), my0=pm(i5,m4);

cntr=dm(tapsB_1);

do filt_l2 until ce;

filt_l2: mr=mr+mx0*my0 (ss),mx0=dm(i2,m2),my0=pm(i5,m4);

mr=mr+mx0*my0 (ss);

ax0=dm(i2,m2);

dm (tx_buf + 1) = mr1;

rti;

.endmod;

 

 

 

 

4.2) Ampliación Albiol. Para que en vez de que salga sólo por el canal izquierdo salga también por el derecho (de mono a stereo)

input_samples:

{      i2=^muestrasl;     i2: Apunta a la memoria del canal L;

i3=^muestrasr;     i3: Apunta a la memoria del canal R;

i4=^tapsA          i4: Apunta a los coefic.del numerador del filtro.

i5=^tapsB           i5: Apunta a los coefic. del denominador del filtro.

m2=1;              m2: Para incrementar direcci¢n;

m4=-1;                                   m4: Para decrementar dirección

i0,i1,m0,m1,l0,l1: NO se deben usar pues los usa el autobuffering de

salida  y  entrada}

{Copiamos sólo el canal R que iría a continuación del canal L }

{ CANAL R}

ax0 = dm (rx_buf + 2);     {Leer nueva muestra de A/D Derecho y guardarla en registro ax0}

mx0=dm(i3,m2), my0=pm(i4,m4);

cntr=dm(tapsA_1);

do filt_r1 until ce;

filt_r1:         mr=mr-mx0*my0 (ss), mx0=dm(i3,m2),my0=pm(i4,m4);

mr=mr-mx0*my0 (ss);

ay0=mr1;

ar=ax0+ay0;

dm(i3,m2)=ar; {Guardar nueva muestra en memoria}

 

 

 

 

 

 

 

 

{Segundo Bucle: Coefs. B del filtro}

 

mr=0, mx0=dm(i3,m2), my0=pm(i5,m4);

cntr=dm(tapsB_1);

do filt_r2 until ce;

filt_r2: mr=mr+mx0*my0 (ss), mx0=dm(i3,m2),my0=pm(i5,m4);

mr=mr+mx0*my0 (ss);

 

ax0=dm(i3,m2);

 

dm (tx_buf + 2) = mr1; {Enviar resultado a D/A DERECHO}

rti;

.endmod;

 

Deja un comentario

Tu dirección de correo electrónico no será publicada. Los campos obligatorios están marcados con *