Học liệu

Thiết kế bộ lọc số thông dải sử dụng bộ lọc Chebychev-I

  • 17/04/2019
  • Học liệu

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: ωs=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

Các tin khác