Bài tập môn Xử lý tín hiệu số 2 - Nội dung: Thiết kế bộ lọc số IIR thông cao thỏa mãn các yêu cầu đề bài - Nguyễn Đình Quý

BÀI TẬP 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 NỘI DUNG: Thiết kế bộ lọc số IIR thông cao thỏa mãn các yêu cầu đề bài. 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 BÀI 1. Thiết kế bằng cách sử dụng bộ lọc Chebyshev-II, phép biến đổi song tuyến 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: Thiết kế bộ lọc tương tự thông thấp. Áp dụng biến đổi song tuyến tính để thu được bộ lọc số thông thấp. Áp dụng biến đổi Frequency-band(dùng hàm zmapping) để chuyển từ bộ lọc số thông thấp sang bộ lọc số thông cao. Các bước thực hiện: Bướ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 cần thiết kế bộ lọc thông thấp có các tần số wslp và wplp để sau đó chuyển đổ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))) Bước 2: Thiết kế bộ lọc analog thông thấp Chebysev-II Chọn giá trị T = 1. 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) Thiết kế bộ lọc tương tự Chebysev-II với các thông số: OmegaP, OmegaS, As, Rp. Ta thu được các đa thức ở tử và mẫu của hàm truyền ở miền s. Bước 3: Áp dụng bilinear transform để chuyển bộ lọc thông thấp tương tự thành bộ lọc thông thấp số. Ta thu được các đa thức ở tử và mẫu của hàm truyền của bộ lọc thông thấp ở miền Z. Bước 4: Sử dụng hàm zmapping(Frequency-band transform) để chuyển đổi từ bộ lọc thông thấp thành bộ lọc thông cao. Ta thu được các đa thức ở tử và mẫu của hàm truyền của bộ lọc thông cao ở miền z. Bước 5: Chuyển đổi các đa thức ở tử và mẫu của hàm truyền từ Direct form sang Cascade form. Bước 6: Tìm hàm đáp ứng xung của bộ lọc số thông cao. Bước 7: Tính các hàm đáp ứng biên độ, đáp ứng biên độ tương đối ở dB, đáp ứng pha, độ trễ nhóm của bộ lọc. Bước 8: Tính As và Rp của bộ lọc thực tế. Bước 9: Vẽ đồ thị và so sánh, đối chiếu bộ lọc thực tế với yêu cầu đề bài. Thiết kế: Sử dụng các Function: Sử dụng các Function: Hàm u_chb1ap: function [b,a] = u_chb2ap(N,Rp,Omegac); % Unnormalized Chebyshev-2 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; Hàm afd_chb2: function [b,a] = afd_chb1(Wp,Ws,Rp,As); % Analog Lowpass Filter Design: Chebyshev-2 % ----------------------------------------- % [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*A-1)/ep; N = ceil(log10(g+sqrt(g*g-1))/log10(OmegaR+sqrt(OmegaR*OmegaR-1))); fprintf('\n*** Chebyshev-2 Filter Order = %2.0f \n',N) [b,a]=u_chb2ap(N,Rp,OmegaC); Hàm zmapping: 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)| % ---- = ----| 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: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; Hàm freqz_m: function [db,mag,pha,grd,w] = freqz_m(b,a); % Modified version of freqz subroutine % ------------------------------------ % [db,mag,pha,grd,w] = freqz_m(b,a); % db = Relative magnitude in dB computed over 0 to pi radians % 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); Hàm dir2cas: 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 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,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 Code chương trình: clc; clear; close all; wp = 0.5*pi; ws = 0.3*pi; As = 15; Rp = 1; % Determine digital lowpass cutoff frequencies(wp & ws) wplp = 0.2*pi; alpha = -(cos((wplp+wp)/2))/(cos((wplp-wp)/2)); wslp = angle(-(exp(-j*ws)+alpha)/(1+alpha*exp(-j*ws))); T = 1; Fs = 1/T; OmegaP = (2/T)*tan(wplp/2); OmegaS = (2/T)*tan(wslp/2); %Chebysev Prototype Lowpass Filter % ve do thi thong thap analog [b1,a1]=afd_chb2(wplp,wslp,Rp,As); [db1,mag1,pha1,w1]=freqs_m(b1,a1,0.5*pi) [ha1,x1,t1]=impulse(b1,a1) figure(1) subplot(2,2,1); plot(w1/pi,mag1);title('Magnitude Response');grid; xlabel('Frequency in pi units');ylabel('|H|'); subplot(2,2,2); plot(w1/pi,db1);title('Magnitude Response in db');grid; xlabel('Frequency in pi units');ylabel('decibels'); subplot(2,2,3); plot(w1/pi,pha1);title('Phase Response');grid; xlabel('Frequency in pi units');ylabel('Radians'); subplot(2,2,4); plot(t1,ha1);title('Impulse Response');grid; xlabel('n');ylabel('h(n)'); %Bilinear transformation [cs,ds] = afd_chb2(OmegaP,OmegaS,Rp,As); [blp,alp] = bilinear(cs,ds,Fs); %Transform digital lowpass into highpass filter Nz = -[alpha,1]; Dz = [1,alpha]; [b,a] = zmapping(blp,alp,Nz,Dz); %From direct to cascade form [C,B,A] = dir2cas(b,a); %Impulse response [h,t] = impz(b,a); %Get Magnitude/Phase Response in db, Group Delay [db,mag,pha,grd,w] = freqz_m(b,a); %Calculate actual Rp,As delta_w = 2*pi/1000; Rp = -(min(db(wp/delta_w+1:1:501))) As = -round(max(db(1:1:ws/delta_w+1))) %Plot figure(3) subplot(3,2,1); plot(w/pi,mag);axis tight;title('Magnitude Response');grid; xlabel('Frequency in pi units');ylabel('|H|'); subplot(3,2,2); plot(w/pi,db);axis([0,1,-30,0]);title('Magnitude Response in db');grid; xlabel('Frequency in pi units');ylabel('decibels'); subplot(3,2,3); plot(w/pi,pha);title('Phase Response');grid; xlabel('Frequency in pi units');ylabel('Radians'); subplot(3,2,4); plot(t,h);title('Impulse Response');grid; xlabel('n');ylabel('h(n)'); subplot(3,2,5); plot(w/pi,grd);axis([0 1 0 10]);title('Group Delay');grid; xlabel('Frequency in pi units');ylabel('Samples'); Kết quả: Đồ thị: