Filtros FIR Ejemplos Ejemplos de cálculo 1. Método de las ventanas. 1. Diseñar un filtro FIR pasa-bajos con las características enunciadas posteriormente: •
Frecuencia de corte de la banda pasante: 1000 Hz
•
Frecuencia de corte de la banda suprimida: 1600 Hz
•
Frecuencia de muestreo y reconstrucción: 8 kHz
•
Atenuación en la banda suprimida: 50 dB
Solución: Las ventanas que permiten cumplir con las especificaciones de rizado en la banda suprimida son la de Hamming y Blackman porque permiten una atenuación mayor a 50 dB. Las frecuencias y pulsaciones digitales de corte para la banda suprimida y la pasante son
fs
=
fp
=
Fs 1.6kHz = = 0.2 → ωs = 2πfs = 0.4 π Fm 8kHz Fp 1kHz = = 0.125 → ωp = 2πfp = 0.25 π Fm 8kHz El ancho de transición puede calcularse como
∆ω = ωs − ωp = 0.15π La frecuencia de corte es
ωc =
ωs + ωp 2
= 0.325π
El orden de la ventana se determina con la fórmula de la tabla de ventanas en función del ancho de transición ∆ω. En caso de elegir la ventana de Hamming 8π M= = 53.55 0.15π Para que cumpla la especificación de banda pasante se requiere una respuesta al impulso simétrica (el caso antisimétrico corresponde a los transformadores de Hilbert y
diferenciadores). Si el orden fuera impar (respuesta al impulso Tipo I) se podría implementar cualquier banda pasante (LP, HP, BP, BS). Incrementando en dos unidades para que N sea un entero impar M=55 La respuesta espectral deseada del filtro ideal tiene la forma
( )
Hd e jω = e
− jω
M −1 2
− jω M −1 2 ω e Π = 2ωc 0
si si
ω ≤ ωc 0< ω
En el ejemplo analizado la fase es α=(M-1)/2=27, y la frecuencia de corte ωc=0.325π, por lo cual la respuesta es e − j 27ω Hd e jω = 0
( )
si si
ω ≤ 0.325π 0< ω
El filtro ideal tendría una respuesta impulsiva con una cantidad infinita de muestras
hd [n ] =
M − 1 ωc ωc senc n − π 2 π
Con los valores del ejemplo hd n = 0.325senc [(n − 27 )0.325] =
sen[0.325π (n − 27 )] π (n − 27 )
La ventana de Hamming es una función dada por 2πn 0.54 − 0.64 cos si 0 ≤ n ≤ M − 1 w [n ] = M − 1 0 cc Con los valores del ejemplo se obtiene una ventana de 55 muestras 2πn 0.54 − 0.64 cos si 0 ≤ n ≤ 54 w [n ] = 54 0 cc
Así, la respuesta al impulso del filtro enventanado resulta de multiplicar los 55 valores de la ventana de Hamming por los valores de la respuesta impulsiva ideal resultando 2πn 0.325senc ((n − 27 )0.325 ) 0.54 − 0.64 cos si h[n ] = hd [n ]w [n ] = 54 0
0 ≤ n ≤ 54 cc
Los coeficientes resultantes son hd[n] w[n] h[n] … 0 0 h [ 0 ] = h [ 54 ] = 0.00765650867949 * 0.08000000000000 = 0.00061145191874 h [ 1 ] = h [ 53 ] = 0.01209196012615 * 0.08311035543871 = 0.00100321355619 h [ 2 ] = h [ 52 ] = 0.00487247679202 * 0.09239935953328 = 0.00044942816560 h [ 3 ] = h [ 51 ] = -0.00779574403157 * 0.10774139443848 = -0.00083845876477 h [ 4 ] = h [ 50 ] = -0.01379689752327 * 0.12892898545123 = -0.00177571617112 h [ 5 ] = h [ 49 ] = -0.00656862110457 * 0.15567560675005 = -0.00102078980605 h [ 6 ] = h [ 48 ] = 0.00791983136504 * 0.18761955616527 = 0.00148332250057 h [ 7 ] = h [ 47 ] = 0.01591549430919 * 0.22432884658038 = 0.00356407472526 h [ 8 ] = h [ 46 ] = 0.00875349782452 * 0.26530704781672 = 0.00231831241629 h [ 9 ] = h [ 45 ] = -0.00802831468337 * 0.31000000000000 = -0.00248443493158 h [ 10 ] = h [ 44 ] = -0.01866639076678 * 0.35780330762199 = -0.00666724247944 h [ 11 ] = h [ 43 ] = -0.01169361604736 * 0.40807051295290 = -0.00476349364161 h [ 12 ] = h [ 42 ] = 0.00812079465337 * 0.46012183827321 = 0.00373003512099 h [ 13 ] = h [ 41 ] = 0.02245649737714 * 0.51325337870118 = 0.01150576187740 h [ 14 ] = h [ 40 ] = 0.01590197956510 * 0.56674662129882 = 0.00899666763830 h [ 15 ] = h [ 39 ] = -0.00819693035903 * 0.61987816172679 = -0.00507223221200 h [ 16 ] = h [ 38 ] = -0.02813772351730 * 0.67192948704710 = -0.01887357642469 h [ 17 ] = h [ 37 ] = -0.02250790790393 * 0.72219669237801 = -0.01622677336412 h [ 18 ] = h [ 36 ] = 0.00825644079990 * 0.77000000000000 = 0.00634636640681 h [ 19 ] = h [ 35 ] = 0.03784133643203 * 0.81469295218328 = 0.03077527703870 h [ 20 ] = h [ 34 ] = 0.03457781948052 * 0.85567115341962 = 0.02953561646482 h [ 21 ] = h [ 33 ] = -0.00829910612703 * 0.89238044383473 = -0.00739303749136 h [ 22 ] = h [ 32 ] = -0.05881599776824 * 0.92432439324995 = -0.05427020089651 h [ 23 ] = h [ 31 ] = -0.06437952685006 * 0.95107101454877 = -0.06112266373281 h [ 24 ] = h [ 30 ] = 0.00832476861040 * 0.97225860556152 = 0.00807970515553 h [ 25 ] = h [ 29 ] = 0.14180809265171 * 0.98760064046672 = 0.13980539297808 h [ 26 ] = h [ 28 ] = 0.27140379367128 * 0.99688964456129 = 0.27008753708905 h [ 27 ] = 0.32500000000000 * 100.000.000.000.000 = 0.32443291372751 Un procedimiento similar se debe seguir en caso de elegir la ventana de Blackman.
86
2. Método del muestreo en frecuencia. 2. Diseñar el filtro del ejercicio anterior haciendo el muestreo en frecuencia de la respuesta espectral ideal con 32 muestras. Solución: Primero debe establecerse la relación entre los valores de ω y k. El límite inferior de la escala de frecuencias angulares digitales ω=0 corresponde a un valor de k=0 y el límite superior ω=2π corresponde al valor de la última muestra, k=M. Por regla de 3 simple se pueden calcular los valores de k que corresponden a una frecuencia angular digital cualquiera
2π __________ ____ M
ω __________ _____ k
⇒k =
Mω 2π
Así la frecuencia de corte de la banda pasante corresponde a kp =
ωpM 0.25π 32 = =4 2π 2π Igualmente, la frecuencia de corte de la banda suprimida corresponde a
ks =
ωsM 0.4π 32 64 = = = 6.4 ≅ 6 2π 2π 10
El valor de k debe ser entero. En este último caso, para ser conservativos en el cálculo, conviene considerar la zona de transición más angosta de lo exigido en las especificaciones, por lo que ks se aproxima con el entero inferior. De lo anterior, se desprende que la muestra k=5 se encuentra en la zona de transición. En general no interesa qué valor exacto toma la atenuación del filtro en esa zona, sino simplemente importa el ancho en frecuencia de esa región. Es decir que H[5] podría fijarse arbitrariamente. El valor que optimiza el rizo del lóbulo lateral máximo en la banda de rechazo puede buscarse en las Tablas de Coeficientes de Transición de Rabiner et al. (reimpresas en el libro Tratamiento Digital de Señales; J. Proakis, D. Manolakis; Apéndice C). En este caso hay sola muestra en la zona de transición (un solo Ti) y el ancho de banda pasante es BW=5 muestras (k=0, 1, ...,4). Con α=0 se obtiene un valor de atenuación T1=0.37897949. Entonces, para la muestra k=5, el valor de H[k]=T1 y la atenuación en esa muestra es de H[k]dB=-40.85183477 dB. De esta forma, los valores de atenuación H[k] correspondientes a cada intervalo de k son
1 si 2π Hr k = T 1 si 32 0 si
k = 0;1;...;4 k =5 k = 6;7;...;15
con H [k ] = H [32 − k ]
La función auxiliar G[k] es (− 1)k 2πk G[k ] = (− 1) Hr = −T1 M 0 k
si
k = 0;1;...;4
si
k =5
si
k = 6;7;...;15
con G[k ] = −G[M − k ]
Con esta función auxiliar se obtiene la respuesta al impulso h[n], es decir, los coeficientes del filtro 1 h[n ] = M
U 2πk G[0] + 2∑ G[k ]cos k =1 M
15 2πk 1 1 n + = 1 + 2∑ G[k ]cos 2 32 k =1 32
1 n + 2
Calculando con una computadora los 32 valores anteriores se obtienen los coeficientes h[ 0 ]=h[31]= 0.00739435497045 h[ 1 ]=h[30]= 0.00535455868240 h[ 2 ]=h[29]=-0.00088103104021 h[ 3 ]=h[28]=-0.01048388541473 h[ 4 ]=h[27]=-0.01633937096262 h[ 5 ]=h[26]=-0.00968048242689 h[ 6 ]=h[25]= 0.00979906225005 h[ 7 ]=h[24]= 0.02854450588269 h[ 8 ]=h[23]= 0.02684227563082 h[ 9 ]=h[22]=-0.00339733474223 h[10]=h[21]=-0.04510122574620 h[11]=h[20]=-0.06103186085330 h[12]=h[19]=-0.01811476589717 h[13]=h[18]= 0.08499188112202 h[14]=h[17]= 0.20891457190105 h[15]=h[16]= 0.29318874664386 Con estos coeficientes puede implementarse el filtro inmediatamente en la forma directa. Sin embargo hay una estructura más eficiente y para el cálculo de la misma, hace uso de las definiciones de las transformadas y antitransformadas. La transferencia de un sistema causal cuya respuesta impulsiva tiene M muestras es H [z ] =
1 M
h[n ] =
M −1
∑ H [k ]e
j
2πnk M
k =0
Usando estas definiciones, se puede escribir la transferencia de la siguiente forma H [z ] =
M −1
M −1
1
M −1
n =0
n =0
k =0
∑ h[n ]z −n = ∑ M ∑ H [k ]e
j
2πkn M
−n z
Operando con el orden de las sumatorias y sumando sobre el índice n se obtiene
1 H [z ] = ∑ H [k ] M k =0 M −1
M −1
∑e
j
2πkn M
n =0
z −1
−n
=
1 − z −M M
M −1
∑
k =0
H [k ] 1− e
j
2πkn M
z −1
Esto puede verse como el producto de la transferencia de un sistema H1[1] “todos ceros” y otro H2[1] “todo polos”
H [z ] =
1 − z −M M
M −1
∑
k =0
H [k ] 1− e
j
2πkn M
= H 1[z ]H 2[z ] z
−1
Desarrollando la sumatoria para el ejemplo de 32 muestras, donde H[k]=1 para k=0;..;4 y por simetría k=27;...;31, y también H[k]=T1 para k=5 y por simetría para k=26 H [z ] =
1 − z −32 32
31
∑
k =0
H [k ] 1− e
j
2πkn 32
z −1
La sumatoria de la transferencia del sistema “todo polos” muestra que éste es un sistema en paralelo conformado por tantas ramas como muestras no nulas tenga H[k]. Es decir que la estructura será eficiente cuando la mayoría de las muestras sean cero. El único inconveniente que presenta esta transferencia para la implementación es que posee coeficientes complejos. Pero haciendo uso de la simetría de H[k] pueden agruparse los pares de polos complejos conjugados para conformar sistemas de segundo orden con coeficientes reales. Utilizando dos funciones auxiliares A[k] y B[k]
M −1
∑ h[n ]z n
k =0
Y la respuesta impulsiva puede evaluarse como la antitransformada discreta de Fourier
A[k ] = H [k ] + H [M − k ] y
B[k ] = H [k ]e
−j
2πk M
+ H [M − k ]e
j
2πk M
Con esa función y por la simetría, H2[k] puede expresarse como
87
M −1 2 H [0] A[k ] + B[k ]z −1 + ∑ −1 2πk −1 1− z −2 k =1 1 − 2 cos z + z M H 2[z ] = M H M2 −1 A[k ] + B[k ]z −1 0 [ ] H 2+ + ∑ 1 1 − − 1 − z 2πk −1 1+ z −2 k =1 1 − 2 cos z + z M
p/M
impar
p/M
par
En el ejemplo que se analiza, pueden evaluarse A[k] y B[k] ó k = 27;...;31 2 si k = 0;...;4 A[k ] = 2T 1 si k =5 ó k = 26 0 si k = 6;...;25 2πk ó k = 27;...;31 2 cos si k = 0;...;4 M 2πk B[k ] = 2 cos k =5 ó k = 26 T 1 si M 0 si k = 6;...;25 Y la transferencia resultante es la siguiente
1 − z − 32 H [z ] = 32
... +
A[1] + B[1]z −1 A[2] + B[2]z −1 1 + + ... 1 − z −1 + π −1 2π −1 −2 −2 − + − + 1 2 cos z z 1 2 cos z z 16 16
A[3] + B[3]z A[4] + B[4]z + 3π −1 4π −1 −2 −2 1 − 2 cos 1 − 2 cos z + z z + z 16 16 −1
−1
A[5] + B[5]z + 5π −1 −2 1 − 2 cos z + z 16 −1
donde se observa que, para todos los bloques de segundo orden, a0=1 y a2=1. A su vez, para el k-ésimo bloque, a1=2cos(2πk/M).
88
Ejemplos de simulación 1. Método de las ventanas. 1.1. Analizar la variación del rizado en la respuesta espectral para el caso de un filtro FIR pasabajos con las características enunciadas posteriormente, al utilizar diferentes ventanas manteniendo fijo el orden •
Frecuencia de muestreo y reconstrucción: 20 kHz
•
Frecuencia de corte: 5000 Hz
•
Orden del filtro: 20
•
Atenuación en la banda suprimida para el caso de Kaiser: 90 dB
Solución: %Ejemplo 7.1.1. %Comparacion de las ventanas en el diseno de filtro FIR clc, clear, close all %Especificaciones Fm=20000; %frecuencia de muestreo Fc=5000; %frecuencia de corte M=20; %orden del filtro A=90; %rizado en dB (Kaiser) puntos=512; %cantidad de puntos para el grafico %Calculo de parametros wc=2*pi*Fc/Fm; %pulsacion digital de corte vc=2*Fc/Fm; %frecuencia relativa a Fm/2 usada por Matlab de corte %parametro beta de Kaiser if(A<=21) beta=0; elseif (A>=50) beta=0.1102*(A-8.7); else beta=0.5842*(A-21)^0.4+0.07886*(A-21); end %Coeficientes de los filtros enventanados bboxcar=fir1(M-1,vc,boxcar(N)); bhamming=fir1(N-1,vc,hamming(N)); bhann=fir1(N-1,vc,hanning(N));
bbartlett=fir1(N-1,vc,bartlett(N)); bblackman=fir1(N-1,vc,blackman(N)); bkaiser=fir1(N-1,vc,kaiser(N,beta)); %Respuesta espectral de los filtros enventanados [Hboxcar,w]=freqz(bboxcar,1,puntos); [Hhamming,w]=freqz(bhamming,1,puntos); [Hhann,w]=freqz(bhann,1,puntos); [Hbartlett,w]=freqz(bbartlett,1,puntos); [Hblackman,w]=freqz(bblackman,1,puntos); [Hkaiser,w]=freqz(bkaiser,1,puntos); %Respuesta espectral del filtro ideal Hideal=rectpuls(w,2*wc); %Graficos de la respuesta espectral figure(1) plot(w,Hideal,'y-'); hold on plot(w,abs(Hboxcar),'r--') plot(w,abs(Hhamming),'m:') plot(w,abs(Hhann),'g-') plot(w,abs(Hbartlett),'b--') plot(w,abs(Hblackman),'c:') plot(w,abs(Hkaiser),'k-') legend('Ideal','Boxcar','Hamming','Hann','Bartlett','Blackman','Kaiser') xlabel('\omega,rad/muestra') ylabel('Abs(H(e^j^\omega))') disp('Pulsar una tecla para ver detalle'),pause; %Zoom para detalle del mismo grafico axis([0 1.1*wc 0.9 1.1]) disp('Pulsar una tecla para ver respuesta impulsiva'),pause; %Graficos de la respuesta impulsiva figure(2) stem(bboxcar,'r.-') hold on stem(bhamming,'m.:') stem(bhann,'g.-') stem(bbartlett,'b.-') stem(bblackman,'c.:') stem(bkaiser,'ko:') legend('Boxcar','Hamming','Hann','Bartlett','Blackman','Kaiser') xlabel('n') ylabel('h[n]')
89
1.2. Analizar la variación del ancho de transición en la respuesta espectral para el caso de un filtro FIR pasabajos con ventana de Kaiser, bajo las especificaciones enunciadas posteriormente, al modificar el rizado en la banda suprimida manteniendo fijo el orden •
Frecuencia de muestreo y reconstrucción: 20 kHz
•
Frecuencia de corte: 5000 Hz
•
Orden del filtro: 20
•
Atenuación en la banda suprimida: 10; 30; 50; 70 y 90 dB
Solución: %Ejemplo 7.1.2. %Filtro FIR con ventana de Kaiser %Analisis del efecto de la variacion del rizado en la banda suprimida clc, clear, close all %Especificaciones Fm=20000; %frecuencia de muestreo Fc=5000; %frecuencia de corte M=20; %orden del filtro puntos=512; %cantidad de puntos para el grafico A=10:20:90; %atenuacion de la banda suprimida en dB %Calculo de parametros wc=2*pi*Fc/Fm; %pulsacion digital de corte vc=2*Fc/Fm; %frecuencia relativa a Fm/2 usada por Matlab de corte deltaw=(A-8)/(2.285*N); %ancho de banda de transicion %Respuesta impulsiva y espectral del filtro enventanado %Calculo para los 5 valores de A for(i=1:length(A)) %parametro beta para la ventana de Kaiser if(A(i)<=21) beta=0; elseif (A(i)>=50) beta=0.1102*(A(i)-8.7); else beta=0.5842*(A(i)-21)^0.4+0.07886*(A(i)-21); end %Coeficientes del filtro enventanado bkaiser=fir1(M-1,vc,kaiser(M,beta)); %Respuesta espectral del filtro enventanado [Hkaiser(:,i),w]=freqz(bkaiser,1,puntos);
90
%Ancho de banda de transicion disp('Para una atenuacion de rizado de'),disp(A(i)),... disp('el ancho de la banda de transicion es'),disp(deltaw(i)) end %respuesta espectral en dB del filtro HdBkaiser=20*log10(abs(Hkaiser)); %graficos plot(w,HdBkaiser(:,1),'b-',w,HdBkaiser(:,2),'c:',... w,HdBkaiser(:,3),'g--',w,HdBkaiser(:,4),'r:',... w,HdBkaiser(:,5),'k-') hold on title('Filtro FIR con ventana de Kaiser') legend('A=10','A=30','A=50','A=70','A=90') grid xlabel('\omega,rad/muestra') ylabel('Abs(H(e^j\omega)),dB')
1.3. Analizar la variación del rizado en la banda suprimida en la respuesta espectral para el caso de un filtro FIR pasabajos con ventana de Kaiser, bajo las especificaciones enunciadas posteriormente, al modificar el ancho de transición manteniendo fijo el orden •
Frecuencia de muestreo y reconstrucción: 20 kHz
•
Frecuencia de corte: 5000 Hz
•
Atenuación en la banda suprimida: 90 dB
•
Ancho de la transición: 500; 1000; 1500; 2000 y 2500 rad/s
Solución: %Ejemplo 7.1.3. %Filtro FIR con ventana de Kaiser %Analisis del efecto de la variacion %del ancho de la transicion clear, close all, clc %Especificaciones Fm=20000; %frecuencia de muestreo Fc=5000; %frecuencia de corte puntos=512; %cantidad de puntos para el grafico A=90; %atenuacion de la banda suprimida en dB deltaW=500:500:2500; %ancho de banda analogico de transicion
%Calculo de parametros wc=2*pi*Fc/Fm; %pulsacion digital de corte vc=2*Fc/Fm; %frecuencia relativa a Fm/2 usada por Matlab de corte deltaw=deltaW/Fm; %ancho de banda digital de transicion %parametro beta para la ventana de Kaiser if(A<=21) beta=0; elseif (A>=50) beta=0.0012*(A-8.7); else beta=0.5842*(A-21)^0.4+0.07886*(A-21); end %orden del filtro M=(A-7.95)./(14.36*deltaw); M=ceil(M); %redondeo hacia arriba N=M+1; %orden del filtro
2. Método de diseño óptimo con rizado constante.
%Respuesta impulsiva y espectral del filtro enventanado %Calculo para los 5 valores de N for(i=1:length(N)) %Coeficientes del filtro bkaiser=fir1(N(i)-1,vc,kaiser(N(i),beta)); %Respuesta espectral del filtro [Hkaiser(:,i),w]=freqz(bkaiser,1,puntos); %Orden del filtro disp('Para un ancho de banda de transicion de'),disp(deltaw(i)),... disp('el orden del filtro es'),disp(N(i)) end
%Caso Fs variable %Especificaciones Fs=210:10:300; %frecuencia de corte de la banda suprimida A1=3; %rizado en dB de la banda pasante A2=-60; %rizado en dB en la banda suprimida
%respuesta espectral en dB HdBkaiser=20*log10(abs(Hkaiser)); %graficos plot(w,HdBkaiser(:,1),'b-',w,HdBkaiser(:,2),'c:',... w,HdBkaiser(:,3),'g--',w,HdBkaiser(:,4),'r:',w,HdBkaiser(:,5),'k-') hold on title('Filtro FIR con ventana de Kaiser') legend('\delta\omega=500','\delta\omega=1000',... '\delta\omega=1500','\delta\omega=2000','\delta\omega=2500') grid xlabel('\omega,rad/muestra') ylabel('Abs(H(e^j\omega)),dB')
2.1. Verificar el efecto de la variación del orden en función de la variación del rizado en las bandas y del ancho de banda de transición Solución: %Filtro pasa-bajos por el Metodo de Optimizacion de Errores %Efecto de la variacion de los parametros clc, clear, close all %Especificaciones F0=200; %frecuencia de corte Fm=1000; %frecuencia de muestreo puntos=512; %puntos para el grafico del espectro nivel=[1 0]; %vector de nivel en cada banda
%Estimacion de parametros delta1=(10^(A1/20)-1)/(10^(A1/20)+1); %rizado en la banda pasante delta2=10^(A2/20); %rizado en la banda suprimida Fp=2*F0-Fs; %frecuencia de corte de la banda pasante DeltaF=Fs-Fp; %ancho de banda de transicion Fc=[Fp' Fs']; %Vector de frecuencias de corte delta=[delta1 delta2]; %Vector de atenuaciones para la funcion remez for i=1:length(Fs) [M(i),vc,niveles,W]=remezord(Fc(i,:),nivel,delta,Fm); %parametros remez end %Grafico orden vs anchura de la banda de transicion figure(1) plot(DeltaF,M,'k*',DeltaF,M,'k') grid xlabel('Anchura de la banda de transicion, Hz') ylabel('Orden del filtro'); title('Orden del filtro en funcion de la anchura de la banda de transicion') %A1 variable %Especificaciones clear M A1=1:1:5; %rizado en dB de la banda pasante
91
A2=-60; %rizado en dB de la banda suprimida DeltaF=20; %ancho de banda de transicion %Estimacion de parametros Fp=(2*F0-DeltaF)/2; %frecuencia de corte de la banda pasante Fs=(2*F0+DeltaF)/2; %frecuencia de corte de la banda suprimida delta2=10^(A2/20); %rizado en la banda suprimida delta1=(10.^(A1/20)-1)./(10.^(A1/20)+1); %rizado en la banda pasante Fc=[Fp Fs]; %Vector de frecuencias de corte delta=[delta1;ones(1,length(delta1))*delta2]; %Vector de rizado for i=1:length(A1) [M(i),vc,niveles,W]=remezord(Fc,nivel,delta(:,i),Fm); %parametros remez end %Grafico del orden en funcion del rizado en dB en la banda pasante figure(2) plot(A1,M,'k*',A1,M,'k') grid xlabel('Rizado en la banda pasante, dB') ylabel('Orden del filtro'); title('Orden del filtro en funcion del Rizado en dB en la banda pasante') %A2 variable %Especificaciones clear M A1=3; %rizado en dB de la banda pasante A2=-20:-10:-100; %rizado en dB de la banda suprimida DeltaF=20; %ancho de banda de transicion Fp=(2*F0-DeltaF)/2; %frecuencia de corte de la banda pasante Fs=(2*F0+DeltaF)/2; %frecuencia de corte de la banda suprimida delta2=10.^(A2/20); %rizado en la banda suprimida delta1=(10^(A1/20)-1)/(10^(A1/20)+1); %rizado en la banda pasante Fc=[Fp Fs]; %Vector de frecuencias de corte delta=[delta1*ones(1,length(delta2));delta2]; %Vector de rizados for i=1:length(A2) [M(i),vc,niveles,W]=remezord(Fc,nivel,delta(:,i),Fm); %parametros remez end %Grafico del orden en funcion del rizado en dB en la banda suprimida figure(3) plot(A2,M,'k*',A2,M,'k') grid xlabel('Rizado en la banda suprimida, dB') ylabel('Orden del filtro'); title('Orden del filtro en funcion del Rizado en dB en la banda suprimida')
92
%Diseno del filtro por el Metodo de Optimizacion del Error bremez=remez(M(5),vc,niveles,W); %coeficientes del filtro [Bremez,w]=freqz(bremez,1,puntos); %respuesta espectral del filtro BdBremez=20*log10(abs(Bremez)); %amplitud en dB de la respuesta espectral %Grafico figure(4) plot(w,BdBremez,'k-') grid xlabel('\omega,rad/muestra') ylabel('Abs(H(e^(j*\omega))),dB'); title('Amplitud de la respuesta espectral del filtro en dB caso intermedio')
Ejemplos de laboratorio 1. Método de las ventanas.
M
Escribir un programa que implemente el algoritmo de filtrado en la forma directa de un FIR declarando los coeficientes en una tabla. Utilizar como base el programa de muestreo y reconstrucción.
Solución: ;-------------------------------------------------------------------------; ffir.asm ; Programa con algoritmo generico para implementar filtros FIR ; mediante una tabla de coeficientes. ; ; Como ejemplo se programo un filtro con las siguientes caracteristicas: ; Orden: 61 ; Rizado en la banda pasante: 3db ; Atenuacion minima en la banda rechazada: -60db ; Frecuencia de corte de la banda pasante: 1500Hz ; Frecuencia de corte de la banda rechazada: 1650Hz ; Frecuencia de muestreo y reconstruccion:9796.27kHz ; ; El codigo del filtro fue tomado de la Guia de Usuario del TMS320C3x. ; Los parámetros del AIC y control fueron diseñados para trabajar ; con el TMS320C31 DSK ;-------------------------------------------------------------------------; Definicion constantes usadas en el programa ;-------------------------------------------------------------------------; Valores de los registros del AIC que definen las ; frecuencias de conversion del DAC (TB) y del ADC (RB) ; frec. de conversion = frec. del SCF / contador B = 9796.24Hz ; frec. del SCF = Master Clock / (2 * Contador A) = 3.55KHz ; Master Clock = Timer cero del DSP = 6.25MHz TA .set 11 TB .set 29 RA .set 11 RB .set 29 ; Constante utilizada para poner en 1 el bit GIE ; del registro ST que habilita las interrupciones GIE .set 0x2000 .include "C3XMMRS.ASM" .brstart "AICTEST",64 ; inicio de ensamblado .sect "AICTEST" ; tabla en un contorno de 2^N ; hasta 64 lugares de largo
;-------------------------------------------------------------------------; Buffer de datos de 32 bits en punto flotante de datos entrantes del AIC ; La ubicación en memoria deberia ser en un contorno de 2^N ; y el tamaño de la tabla de coeficientes ; deberia ser el mismo que de la de datos ;-------------------------------------------------------------------------ADC_recv .float 0.0 ; 0 .float 0.0 ; 1 .float 0.0 ; 2 .float 0.0 ; 3 .float 0.0 ; 4 .float 0.0 ; 5 .float 0.0 ; 6 .float 0.0 ; 7 .float 0.0 ; 8 .float 0.0 ; 9 .float 0.0 ; 10 .float 0.0 ; 11 .float 0.0 ; 12 .float 0.0 ; 13 .float 0.0 ; 14 .float 0.0 ; 15 .float 0.0 ; 16 .float 0.0 ; 17 .float 0.0 ; 18 .float 0.0 ; 19 .float 0.0 ; 20 .float 0.0 ; 21 .float 0.0 ; 22 .float 0.0 ; 23 .float 0.0 ; 24 .float 0.0 ; 25 .float 0.0 ; 26 .float 0.0 ; 27 .float 0.0 ; 28 .float 0.0 ; 29 .float 0.0 ; 30 .float 0.0 ; 31 .float 0.0 ; 32 .float 0.0 ; 33 .float 0.0 ; 34 .float 0.0 ; 35 .float 0.0 ; 36 .float 0.0 ; 37 .float 0.0 ; 38 .float 0.0 ; 39
93
.float 0.0 ; 40 .float 0.0 ; 41 .float 0.0 ; 42 .float 0.0 ; 43 .float 0.0 ; 44 .float 0.0 ; 45 .float 0.0 ; 46 .float 0.0 ; 47 .float 0.0 ; 48 .float 0.0 ; 49 .float 0.0 ; 50 .float 0.0 ; 51 .float 0.0 ; 52 .float 0.0 ; 53 .float 0.0 ; 54 .float 0.0 ; 55 .float 0.0 ; 56 .float 0.0 ; 57 .float 0.0 ; 58 .float 0.0 ; 59 .float 0.0 ; 60 ;-------------------------------------------------------------------------; Tabla de coeficientes del Filtro FIR ;-------------------------------------------------------------------------FIR_coef .float -9.50540294E-03 ; 0 .float -9.35886553E-03 ; 1 .float -1.23368528E-04 ; 2 .float 9.91606486E-03 ; 3 .float 1.10832484E-02 ; 4 .float 1.54810237E-03 ; 5 .float -1.02919941E-02 ; 6 .float -1.31031575E-02 ; 7 .float -3.33973261E-03 ; 8 .float 1.06304609E-02 ; 9 .float 1.55434458E-02 ; 10 .float 5.64951200E-03 ; 11 .float -1.09289976E-02 ; 12 .float -1.86197869E-02 ; 13 .float -8.74282943E-03 ; 14 .float 1.11854199E-02 ; 15 .float 2.27360857E-02 ; 16 .float 1.31313007E-02 ; 17 .float -1.13978462E-02 ; 18 .float -2.87412365E-02 ; 19 .float -1.99426028E-02 ; 20 .float 1.15647137E-02 ; 21
94
.float .float .float .float .float .float .float .float .float .float .float .float .float .float .float .float .float .float .float .float .float .float .float .float .float .float .float .float .float .float .float .float .float .float .float .float .float .float .float
3.87635101E-02 3.22412721E-02 -1.16847926E-02 -6.00486194E-02 -6.22930183E-02 1.17571967E-02 1.43340077E-01 2.69586301E-01 3.21551942E-01 2.69586301E-01 1.43340077E-01 1.17571967E-02 -6.22930183E-02 -6.00486194E-02 -1.16847926E-02 3.22412721E-02 3.87635101E-02 1.15647137E-02 -1.99426028E-02 -2.87412365E-02 -1.13978462E-02 1.31313007E-02 2.27360857E-02 1.11854199E-02 -8.74282943E-03 -1.86197869E-02 -1.09289976E-02 5.64951200E-03 1.55434458E-02 1.06304609E-02 -3.33973261E-03 -1.31031575E-02 -1.02919941E-02 1.54810237E-03 1.10832484E-02 9.91606486E-03 -1.23368528E-04 -9.35886553E-03 -9.50540294E-03
; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ;
22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
END_coef ;-------------------------------------------------------------------------BufSz .set END_coef - FIR_coef ;-------------------------------------------------------------------------SIZE .word BufSz ; Tamano del filtro ADC_first .word ADC_recv ADC_end .word FIR_coef
ADC_last .word ADC_recv FIR_coefx .word FIR_coef ;-------------------------------------------------------------------------; Definicion de algunas constantes de almacenamiento de datos ;-------------------------------------------------------------------------; Si la primer palabra envia al AIC tiene los dos bits menos significativos ; en uno indica que la proxima es una palabra secundaria. ; Luego, de acuerdo a los dos bits menos significativos ; de la palabra secundaria se realiza una tarea diferente ; (cf manual USER'S GUIDE STARTER KIT pag 4-17). ; Se acomodan los valores de los registros TA, RA, TB y RB ; de acuerdo al formato de la palabra secundaria (cf. manual citado) A_REG .word (TA<<9)+(RA<<2)+0; registro A B_REG .word (TB<<9)+(RB<<2)+2; registro B C_REG .word 00000111b ; registro de control ; Valores que se cargaran en los registros ; de control de la puerta serie 0 S0_gctrl_val .word 0x0E970300 S0_xctrl_val .word 0x00000111 S0_rctrl_val .word 0x00000111 ;************************************************************************** ; Inicio del codigo principal ;************************************************************************** main or GIE,ST ; activa interrupciones ldi 0xF4,IE ; habilita XINT/RINT/INT2 b main ; repeticion ;-------------------------------------------------------------------------DAC2 push ST ; rutina servicio interrupcion DAC push R0 pushf R0 push R2 pushf R2 push AR0 push AR1 ldi @ADC_last,AR1 ldi @FIR_coefx,AR0 ldi @SIZE,BK FIR mpyf3 *AR0++,*AR1++(1)%,R0 ldf 0.0,R2 ldi @SIZE,RC subi 2,RC rptb FIR2 mpyf3 *AR0++,*AR1++(1)%,R0 FIR2 || addf3 R0,R2,R2 addf R2,R0 fix R0,R0
mpyi 4,R0 andn 3,R0 sti R0,@S0_xdata ; Saca el nuevo valor del DAC pop AR1 pop AR0 popf R2 pop R2 popf R0 pop R0 pop ST reti ;-------------------------------------------------------------------------ADC2 push ST push R3 pushf R3 push AR0 ldi @S0_rdata,R3 lsh 16,R3 ash -16,R3 ldi float stf
@ADC_last,AR0 R3,R3 R3,*AR0++
cmpi ldige sti
@ADC_end,AR0 @ADC_first,AR0 AR0,@ADC_last
pop AR0 popf R3 pop R3 pop ST reti ;************************************************************************** ; Esta parte del programa se usa solamente durante la inicializacion ; y se puede sobreescribir sin problemas por el stack o datos ;************************************************************************** .entry ST_STUB ; inicio de la depuracion ; Carga del puntero pagina de datos DP con los 8 bits mas significativos ; del T0_ctrl (=808020h) o sea 80, definiendo asi que se use el DP ; y el stack del Kernel ST_STUB ldp T0_ctrl ; Carga del SP con la direccion del lugar asignado con la directiva ; que se usa mas adelante, que indica que el stack comienza ahí. ; Se programa el TIMER0 que se usa para generar el Master clock. ; El TIMER0 se ha programado en modo pulso con el periodo en 1.
95
; La frecuencia de salida es fsp=(50mhz/4)/2=12.5MHz/2=6.25MHz ldi @stack,SP ldi 0,R0 ; Halt TIM0 & TIM1 sti R0,@T0_ctrl sti R0,@T0_count ; Ajusta la cuenta a 0 ldi 1,R0 ; Ajusta los periodos a 1 sti R0,@T0_prd ; GO=1 hace el reset y comienzo de la cuenta y HLDneg=1 quita el hold ; (en la palabra de control del TIMER0=T0_ctrl) ldi 0x3C1,R0 ; Reinicializa ambos temp sti R0,@T0_ctrl ;-------------------------------------------------------------------------; Programacion del registro del control de la transmision ldi @S0_xctrl_val,R0 ; Control transmision: Con FSX/DX/CLKX en modo de transmision ; de puerta serie (podrian usarse como simples pines de E/S). sti R0,@S0_xctrl ; Programacion del registro de control de la recepcion ; definiendo los pines FSR/DR/CLKR para comunicacion serie (no E/S) ldi @S0_rctrl_val,R0 sti R0,@S0_rctrl ; control de recepcion ldi 0,R0 sti R0,@S0_xdata ; DXR data value ; Programa el Global Control Register para la puerta serie. ; El 'C31 le envia al conversor (AIC) los datos en serie. ; Se le carga una palabra definida al inicio del programa (0x0E970300) ; y se define fuente externa de reloj para la Tx y Rx ; (en este caso la genera el AIC). ldi @S0_gctrl_val,R0 ; Setup serial port sti R0,@S0_gctrl ; Global control ;========================================================================== ; Esta seccion inicializa el AIC ;========================================================================== AIC_INIT LDI 0x10,IE ; Habilita solo XINT int andn 0x34,IF ; Deshabilita int, int2, ; xint0, rint0 ldi 0,R0 ; transmite un cero sti R0,@S0_xdata ; Por el pin XF0 (controlado por el registro IOF) ; saca un cero y luego un uno, ; escalon necesario para resetear el AIC ; y luego ponerlo en RUN rpts 0x040 ldi 2,IOF ; XF0=0 resets AIC rpts 0x40 ldi 6,IOF ; XF0=1 runs AIC
96
;-------------------------------------------------------------------------ldi @C_REG,R0 ; Setup control register call prog_AIC ldi 0xfffc,R0 ; Program AIC ; to be real slow call prog_AIC ldi 0xfffc|2,R0 call prog_AIC ldi @B_REG,R0 ; Bump up the Fs ; to final rate call prog_AIC ; (smallest divisor ; should be last) ldi @A_REG,R0 call prog_AIC b main ; teh DDR before going ; to the main loop ;-------------------------------------------------------------------------prog_AIC ldi @S0_xdata,R1 ; Use original DXR data ; during 2 ndy sti R1,@S0_xdata idle ldi @S0_xdata,R1 ; Use original DXR data ; during 2 ndy or 3,R1 ; Request 2 ndy XMIT sti R1,@S0_xdata idle sti R0,@S0_xdata ; Send register value idle andn 3,R1 sti R1,@S0_xdata ; Leave with original safe ; value in DXR ;-------------------------------------------------------------------------ldi @S0_rdata,R0 ; Fix the receiver ; underrun by reading rets ; the DRR before ; going to the main loop stack .word $ ; Put stack here ;************************************************************************** ; Install the XINT/RINT ISR handler directly * ; into the vector RAM location it will be used for * ;************************************************************************** .start "SP0VECTS",0x809FC5 .sect "SP0VECTS" B DAC2 ; XINT0 B ADC2 ; RINT
