Bài tập thiết kế bộ lọc số IIR thông cao: Xử lý tín hiệu số 2

ọ ễ H và tên: Nguy n Đình Quý

ễ ướ Nguy n Ph ắ c Th ng

ầ Tr n Văn Duy Ph ướ c

ớ

L p: 11DT3 Nhóm: 40

Ộ ế ế ộ ọ ố ề ầ ỏ t k b l c s IIR thông cao th a mãn các yêu c u đ bài. : Thi N I DUNG

Ầ Ế Ế YÊU C U THI T K :

Stopband edge: 0.3 π

Stopband attenuation: As = 15 dB

Passband edge: 0.5 π

Passband ripple: Rp = 1 dB

ế ế ằ ử ụ ế ế ổ t k b ng cách s d ng b l c ChebyshevII, phép bi n đ i song tuy n

ộ ọ ể ể ừ ộ ọ ố BÀI 1. Thi tính( bilinear transform), hàm zmapping(đ chuy n t ấ b l c s thông th p sang thông cao).

ươ ế  Ph ng Pháp Thi ế ộ ọ ố t K B L c S IIR:

(cid:0) ự ấ Thi ế ế ộ ọ ươ t k b l c t ng t thông th p.

(cid:0) ụ ế ế ể ổ ượ ộ ọ ố ấ Áp d ng bi n đ i song tuy n tính đ thu đ c b l c s thông th p.

(cid:0) ụ ế ể ổ ể ừ ộ ọ ố b l c s ộ ọ ố ấ Áp d ng bi n đ i Frequencyband(dùng hàm zmapping) đ chuy n t thông th p sang b l c s thông cao.

ướ ự ệ  Các b c th c hi n:

(cid:0) ề ủ ộ ọ ầ c 1: ấ ừ ầ ế ế ộ ọ ộ ọ T yêu c u đ bài, ta có ws và wp c a b l c thông cao. Đ u tiên, ta ầ ố t k b l c thông th p có các t n s wslp và wplp đ sau đó chuy n ị ể ị ợ ằ ọ ọ ể ừ ứ ề ướ B ể ầ c n thi ổ đ i thành b l c thông cao. Ta ch n giá tr wplp b ng 1 giá tr h p lý, ch n ổ wplp = 0.2*pi. T công th c chuy n đ i Frequency Band trong mi n Z, ta có:

ượ ị ủ ị Thay wp và wplp vào công th c trên, ta thu đ c giá tr c a alpha. Ta xác đ nh ừ ứ ứ ị giá tr wslp t công th c:

ớ V i Z = exp(j*wslp) và z = exp(j*ws). Suy ra:

Wslp = angle((exp(j*ws)+alpha)/(1+alpha*exp(j*ws)))

(cid:0) ướ ế ế ộ ọ ấ Thi t k b l c analog thông th p ChebysevII B c 2:

ọ ị Ch n giá tr T = 1. a.i.1.

ị ầ ố ắ a.i.2. Prewarp các giá tr t n s c t wplp và wslp.

OmegaP = (2/T)*tan(wplp/2) OmegaS = (2/T)*tan(wslp/2) ự a.i.3. Thi ng t ứ ở ử ượ ố ớ ChebysevII v i các thông s : ẫ và m u c các đa th c t ề ở ế ế ộ ọ ươ t k b l c t OmegaP, OmegaS, As, Rp. Ta thu đ ủ c a hàm truy n ề mi n s.

(cid:0) ụ ướ c 3: ấ ố ộ ọ ể ứ ở ử c các đa th c t ộ ọ ề ủ ộ ọ ấ ở ự ấ ươ ể Áp d ng bilinear transform đ chuy n b l c thông th p t ng t B ẫ ủ ượ thành b l c thông th p s . Ta thu đ và m u c a hàm ề mi n Z. truy n c a b l c thông th p

(cid:0) ể c 4:

ử ụ ấ ộ ọ ượ ổ ừ S d ng hàm zmapping(Frequencyband transform) đ chuy n đ i t ẫ và m u ể ứ ở ử t c các đa th c

ở ướ B ộ ọ b l c thông th p thành b l c thông cao. Ta thu đ ề ủ ộ ọ ủ c a hàm truy n c a b l c thông cao ề mi n z.

(cid:0) ứ ở ử ể ổ ẫ ủ ề ừ ướ Chuy n đ i các đa th c t và m u c a hàm truy n t Direct form c 5:

B sang Cascade form.

(cid:0) ướ ủ ộ ọ ố ứ Tìm hàm đáp ng xung c a b l c s thông cao. B c 6:

(cid:0) ộ ươ ứ ộ ố ở Tính các hàm đáp ng biên đ , đáp ng biên đ t ng đ i dB, đáp c 7: B ứ ứ ủ ộ ọ ộ ễ ướ ng pha, đ tr nhóm c a b l c.

(cid:0) ướ ủ ộ ọ ự ế Tính As và Rp c a b l c th c t . B c 8:

(cid:0) ướ ẽ ồ ị ộ ọ ự ế ớ ế ề ố V đ th và so sánh, đ i chi u b l c th c t ầ v i yêu c u đ bài. B c 9:

 Thi ế ế t k :

ử ụ S d ng các Function:

(cid:0) Hàm u_chb1ap:

function [b,a] = u_chb2ap(N,Rp,Omegac); % Unnormalized Chebyshev2 Analog Lowpass Filter Prototype % % [b,a] = u_chb1ap(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] = cheb2ap(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;

(cid:0) Hàm afd_chb2:

function [b,a] = afd_chb1(Wp,Ws,Rp,As); % Analog Lowpass Filter Design: Chebyshev2 % % [b,a] = afd_chb2(Wp,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 edge') end if (Rp <= 0) | (As < 0) error('PB ripple and/or SB attenuation ust be larger than 0') end

ep = sqrt(10^(Rp/10)1); A = 10^(As/20); OmegaC = Wp; OmegaR = Ws/Wp; g = sqrt(A*A1)/ep; N = ceil(log10(g+sqrt(g*g1))/log10(OmegaR+sqrt(OmegaR*OmegaR1))); fprintf('\n*** Chebyshev2 Filter Order = %2.0f \n',N) [b,a]=u_chb2ap(N,Rp,OmegaC);

(cid:0) Hàm zmapping:

function [bz,az] = zmapping(bZ,aZ,Nz,Dz) % Frequency band Transformation from Zdomain to zdomain % % [bz,az] = zmapping(bZ,aZ,Nz,Dz) % performs: % b(z) b(Z)| % ---- = ----| N(z) % a(z) a(Z)|@Z = ---- % D(z) % bzord = (length(bZ)1)*(length(Nz)1); azord = (length(aZ)1)*(length(Dz)1); bz = zeros(1,bzord+1); for k = 0:bzord pln = [1]; for l = 0:k1 pln = conv(pln,Nz); end pld = [1]; for l = 0:bzordk1 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:k1 pln = conv(pln,Nz); end pld = [1]; for l = 0:azordk1 pld = conv(pld,Dz);

end az = az+aZ(k+1)*conv(pln,pld); end az1 = az(1); az = az/az1; bz = bz/az1;

(cid:0) Hàm freqz_m:

function [db,mag,pha,grd,w] = freqz_m(b,a);

radians

% Modified version of freqz subroutine % % [db,mag,pha,grd,w] = freqz_m(b,a); % db = Relative magnitude in dB computed over 0 to pi % mag = absolute magnitude computed over 0 to pi radians % pha = Phase response in radians over 0 to pi radians % grd = Group delay over 0 to pi radians % w = 501 frequency samples between 0 to pi radians % b = numerator polynomial of H(z) (for FIR: b=h) % a = denominator polynomial of H(z) (for FIR: a=[1]) % [H,w] = freqz(b,a,1000,'whole'); H = (H(1:1:501))'; w = (w(1:1:501))'; mag = abs(H); db = 20*log10((mag+eps)/max(mag)); pha = angle(H); grd = grpdelay(b,a,w);

(cid:0) Hàm dir2cas:

function [b0,B,A] = dir2cas(b,a); % DIRECTform to CASCADEform conversion (cplxpair version) % % [b0,B,A] = dir2cas(b,a) % b0 = gain coefficient % B = K by 3 matrix of real coefficients containing bk's % A = K by 3 matrix of real coefficients containing ak's % b = numerator polynomial coefficients of DIRECT form % a = denominator polynomial coefficients 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,NM)]; elseif M > N a = [a zeros(1,MN)]; 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

(cid:0) ươ Code ch ng trình: