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

Chia sẻ: Hoàng Văn Thuận | Ngày: | Loại File: DOCX | Số trang:9

Thêm vào BST

Báo xấu

276
lượt xem 22
download

Download Vui lòng tải xuống để xem tài liệu đầy đủ

Nhằm giúp các bạn có thêm tài liệu phục vụ nhu cầu học tập và nghiên cứu, mời các bạn cùng tham khảo nội dung bài tập Xử lý tín hiệu số 2 "Thiết kế bộ lọc số IIR thông cao" dưới đây. Hy vọng nội dung tài liệu phục vụ hữu ích nhu cầu học tập và nghiên cứu.

Chủ đề:

Bình luận(0) Đăng nhập để gửi bình luận!

Lưu

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

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 ChebyshevII, 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 Frequencyband(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 ChebysevII a.i.1. Chọn giá trị T = 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ết kế bộ lọc tương tự ChebysevII 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(Frequencyband 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 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; 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
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); 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; 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); % 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 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ị: