Học liệu
Thiết kế bộ lọc số thông dải sử dụng bộ lọc Chebychev-I
1. Bài toán thiết kế
Thiết kế bộ lọc thông dải sử dụng bộ lọc Chebychev-I và phương pháp biến đổi bất biến xung, có các đặc tính sau:
- Biên dải chắn dưới: ωs1=0.3*pi
- Biên dải chắn trên: ωs2=0.4*pi
- Biên dải thông dưới: ωp1=0.6*pi
- Biên dải thông trên: ωp2=0.7*pi
Điều kiện bộ lọc thông dải: ωs1 < ωp1 < ωp2< ωs2
- Độ gợn dải thông: Rp=1
- Suy hao dải chắn: As=16
2. Chương trình thiết kế
function [b,a] = u_chb1ap(N,Rp,Omegac);
%unnormalized chebyshev-1 analog Lowpass filter prototype
%--------------------------------------------------
%[b,a] = u_chblap(N,Rp,Omegac);
% b = numerator polynomial coefficients
% a = denominator polynomial coefficients
% N = Order of the Elliptic Filter
% Rp = Passband ripple in dB,Rp > 0
%Omegac = Cutoff frequency in radians/sec
%
[z,p,k] = cheblap(N,Rp);
a = real(poly(p));
aNn = a(N+1);
p=p*Omegac;
a=real(poly(p));
aNu = a(N+1);
k = k*aNu/aNn;
b0 = k;
B = real(poly(z));
b = k*B;
function [b,a] =afd_chb1(Wp,Ws,Rp,As);
% Analog Lowpass Filter Design : Butterworth
% ------------------------------------------
% [b,a] = afd_butt(Wa,Ws,Rp,As);
% b = Numerator coefficients of Ha(s)
% a = Denominator coefficients of Ha(s)
% Wp = Passband edge frequency in rad/sec ; Wp>0
% Ws = Stopband edge frequency in rad/sec ; Ws>Wp>0
% Rp = Passband ripple in +dB (Rp>0)
% As = Stopband attenuation in +dB ;( As >0)
%
if Wp <= 0
error ('Passband edge must be larger than 0 ')
end
if Ws <= Wp
error (' Stopband edge must be larger than Passband egde ')
end
if (Rp <= 0) | (As <0 )
error (' Pb ripple and/or SB attenuation ust ba laeger than 0 ')
end
ep=sqrt(10^(Rp/10)-1);
A=10^(As/20);
OmegaC=Wp;
OmegaR=Ws/Wp;
g=sqrt(A*A-1)/ep;
N = ceil(log10(g+sqrt(g*g-1))/log10(OmegaR+sqrt(OmegaR*OmegaR-1)));
fprintf('\n*** Chebyshev Filter Order = %2.Of \n',N);
[b,a]=u_chb1ap(N,Rp,OmegaC);
function [b,a]=imp_invr(c,d,T)
% Impulse Invariance Transformation from Analog to Digital Filter
%---------------------------------------------------------
% [b,a]= imp_invr(c,d,T)
% b= Numerator polynomial in z ^(-1)of the digital filter
% a= Denominator polynomial in z ^(-1)of the digital filter
% c= Numerator polynomial in s of the analog filter
% d= Denominator polynomial in s of the analog filter
% T=Sampling (transformation)parameter
%
[R,p,k]=residue(c,d);
p=exp(p*T);
[b,a]=residuez(R,p,k);
b=real(b'); a=real(a')
function [bz,az] = zmapping(bZ, aZ, Nz, Dz)
% frequency band transformation from Z-domain to z-domain
%-------------------------------------------------------------
% [bz, az] = zmapping(bZ, aZ, Nz, Dz)
%performs:
% b(z) b(Z)
% ----- -----
% a(z) a(Z) Z= N(z)/D(z)
bzord = (length(bZ) – 1)*(lenth(Nz) - 1);
azord = (length(aZ) – 1)*(length(Dz) – 1);
bz = zeros(1,bzord + 1);
for k = 0:bzord
pln = [1];
for l = 0:k – 1
pln = conv(pln,Nz);
end
pld = [1];
for l = 0:bzord – k – 1
pld = conv(pld,Dz);
end
bz = bz+bZ(k+1)*conv(pln,pld);
end
az = zeros(1,azord + 1);
for k = 0:azord
pln = [1];
for l = 0:k – 1
pln = conv(pln,Nz);
end
pld = [1];
for l = 0:azord – k – 1
pld = conv(pld,Dz);
end
az = az+aZ(k+1)*conv(pln,pld);
end
az1 = az(1); az = az/az1; bz=bz/az1;
function [b0,B,A] = dir2cas(b,a);
%DIRECT-form to CASCADE-form conversion (cplxpair version)
%.........................................................
%[b0,B,A]=dir2cas(b,a)
%b0=gain coefficient
%B= K by 3 matrix of real coefficient
%A= K by 3 matrix of real coefficient
%b = numerator polynomial coefficient of DIRECT form
%a = numerator polynomial coefficient of DIRECT form
% compute gain coefficient b0
b0 = b(1); b = b/b0;
a0 = a(1); a = a/a0;
b0 = b0/a0;
%
M = length(b);N = length(a);
if N > M
b=[b zeros(1,N-M)];
elseif M > N
a=[a zeros(1,M-N)];N=M;
else
NM=0;
end
%
K=floor(N/2);B = zeros(K,3);A = zeros(K,3)
if K*2 == N;
b=[b 0];
a=[a 0];
end
%
broots = cplxpair(roots(b));
aroots = cplxpair(roots(a));
for i=1:2:2*K
Brow = broots(i:1:i+1,:);
Brow = real(poly(Brow));
B(fix((i+1)/2),:) = Brow;
Arow = aroots(i:1:i+1,:);
Arow = real(poly(Arow));
A(fix((i+1)/2),:) = Arow;
end
% Chuyen bo loc thogn thap sang thong dai
% su dung ham ZMAPPING
ws1 = 0.3*pi ;
wp1 = 0.4*pi;
wp2 = 0.6*pi;
ws2 = 0.7*pi;
T=1;
Rp = 1; % Do gon dai thong(dB)
As = 16; % Suy hao dai chan(dB)
wplp = 0.2*pi; % tần số cắt dải thông của DLP
beta = (cos((wp2 + wp1)/2))/(cos((wp2 - wp1)/2));
k = cot((wp2 - wp1)/2)*tan(wplp/2);
alpha1 = -(2*beta*k)/(k+1);
alpha2 = (k-1)/(k+1);
% TÍNH wslp
wslp = angle(-(exp(-j*ws1)*exp(-j*ws1)-alpha1*exp(-j*ws1) +
alpha2)/(alpha2*exp(-j*ws1)*exp(-j*ws1)-alpha1*exp(-j*ws1)+1));
% tần số cắt chắn của DLP
OmegaP = wplp/T; % tần số cắt dải thông ALP
OmegaS = wslp/T; % tần số cắt dải chắn ALP
ep = sqrt(10^(Rp/10)-1);
Ripple = sqrt(1/(1+ep*ep)); % passband Ripple
Attn = 1/(10^(As/20)); % stopband Attenuation
%Tinh toan bo loc tuong tu Chebyshev
[cs,ds] = afd_chb1(OmegaP,OmegaS,Rp,As);
% Bien doi bat bien xung:
[blp,alp] = imp_invr(cs,ds,T);
% Chuyen bo loc thong thap sang thong dai
Nz = -[alpha2,-alpha1,1]
Dz = [1,-alpha1,alpha2]
[bbp,abp] = zmapping(blp,alp,Nz,Dz) ;
[C,B,A] = dir2cas(bbp,abp) % Cascade realization of the digital bandpass filter :
[db1,mag1,pha1,grd1,w1] = freqz_m(bbp,abp);
subplot(2,2,1); plot(w/pi,magl); title('MagnitudeResponse of Chebyshev-I Digital Lowpass Filter')
xlabel('frequency in pi units'); ylabel('|H|'); axis([0,1,0,1]);
set(gca,'XTickMode','manual','XTick',[0;ws1/pi;wp1/pi;wp2/pi;ws2/pi;1],); set(gca,'YTickMode','manual','YTick',[0,Ripple,1]); grid
subplot(2,2,2); plot(w1/pi,pha1/pi); title('Phase Response')
xlabel('frequency in pi units'); ylabel('pi units');
axis([0,1,-1,1.1]);
set(gca,'XTickMode','manual','XTick',[0,wplp/pi,wslp/pi,1]);
set(gca,'YTickmode','manual','YTick',[-1,0,1]); grid
subplot(2,2,3);plot(w1/pi,db1); title('Magnitude in dB')
xlabel('frequency in pi units'); ylabel('Decibels'); axis([0,1,-(As+10),5]);
set(gca,'XTickMode','manual','XTick',[0,wplp/pi,wslp/pi,1]);
set(gca,'YTickmode','manual','YTick',[-As,-Rp,0]); grid
subplot(2,2,4);plot(w1/pi,grd1); title('Group Delay')
xlabel('frequency in pi units'); ylabel('Samples');
axis([0,1,0,15]);
set(gca,'XTickMode','manual','XTick',[0,wplp/pi,wslp/pi,1]);
set(gca,'YTickmode','manual','YTick',[0,max(grd1)]); grid