Thực hành xử lý tín hiệu số - Phạm Hùng Kim Khánh

Nội dung text: Thực hành xử lý tín hiệu số - Phạm Hùng Kim Khánh

1. 32 BÀI 1: PHẦN MỀM MATLAB % Thông thường phần này mô tả chức năng, cách sử dụng, % ví dụ minh họa hay những lưu ý đặc biệt mà tác giả mong muốn trợ % giúp cho người sử dụng. [global tênbiến1, tênbiến2, ] % Khai báo biến toàn cục % (nếu có) % phần trình bày câu lệnh Khởi động MATLAB Editor: >>edit Tạo một script file có tên vd.m, với nội dung như sau: % Doan script file nay hien thi loi chao trong 2s. Sau do hien thi logo cua MATLAB roi thoat close all % Tao mot cua so do hoa figure('Color',[0 0 0], 'Name','Welcome to MATLAB Experiments', 'NumberTitle','off', 'MenuBar','none'); % Hien thi loi chao text( 'String','Welcome to MATLAB', 'Color',[.25 .25 .25], 'Position',[0.01 .501], 'Fontsize',32, 'FontAngle','italic'); text( 'String','Welcome to MATLAB', 'Color','w', 'Position',[0 .5], 'Fontsize',32, 'FontAngle','italic'); axis off;
2. 34 BÀI 1: PHẦN MỀM MATLAB % % Hiển thị khi người sử dụng dùng lệnh help tenham % [global ] %khai báo biến toàn cục (nếu có) out1=kết quả 1 %kết quả trả về của hàm out2=kết quả 2 % Các hàm con (nếu có) function [subout1,subout2, ]=tenhamcon(subin1,subin2, ) Bài 1.21. Xây dựng hàm gptb2 để giải phương trình bậc hai ax2+bx+c=0. Nội dung hàm như sau: function [x1,x2]=gptb2(a,b,c) % Giai phuong trinh bac hai ax^2+bx+c=0 % [x1,x2]=gptb2(a,b,c) % Trong do: x1,x2 la nghiem % a, b, c la 3 he so cua phuong trinh if nargin<3 error('Error! Nhap 3 he so cua phuong trinh') elseif a==0 x1=-c/b; x2=[]; else delta = b^2 - 4*a*c; x1 = (-b+sqrt(delta))/(2*a); x2 = (-b-sqrt(delta))/(2*a); end Lưu tên file là gptb2.m và kiểm tra kết quả:
3. 36 BÀI 1: PHẦN MỀM MATLAB case 'BECNULLI' r=a*sqrt(abs(2*cos(2*theta))); case 'ASTROIT' r=a*sqrt(abs(1-sin(3*theta)/4)); case 'XOANOC' r=a*cos(theta)+1; otherwise error('Chon: ''Becnuli'', ''Astroit'' hoac ''Xoanoc''') end % end of switch % ve do thi close all; figure('Color','w'); for k=1:m hold on r1=r*k; mau=[rand(1,1) rand(1,1) rand(1,1)]; h=polar(theta,r1); set(h,'color',mau,'LineWidth',2); axis equal; end % end of for hold off; axis off end % end of if Kiểm tra lại hoạt động của hàm, ví dụ: >>help vdcongdb >>vdcongdb(1,5,’Becnulli’) >>vdcongdb(1,5,’ Astroit’) >>vdcongdb(1,5,’Xoanoc’) >> vdcongdb(1,5,’saikieu’) >> vdcongdb(5,’becnulli’)
4. 38 BÀI 2: TÍN HIỆU RỜI RẠC THEO THỜI GIAN BÀI 2: TÍN HIỆU RỜI RẠC THEO THỜI GIAN Sau khi học xong bài này, người học có thể: Biết được các tín hiệu sơ cấp. Thực hiện các phép toán đơn giản. Tính năng lượng của tín hiệu. Xác định các tính chất của hệ rời rạc. 2.1 CÁC TÍN HIỆU SƠ CẤP Hàm xung đơn vị: 1 푛 = 0 훿(푛) = { 0 푛 ≠ 0 >>n = -10:10; >>delta=[zeros(1,10) 1 zeros(1,10)]; >>stem(n,delta); Bài 2.1. Viết chương trình và vẽ dạng tín hiệu hàm bước đơn vị u(n). Hàm bước đơn vị: 1 푛 ≥ 0 (푛) = { 0 푛 < 0 Bài 2.2. Viết chương trình và vẽ dạng tín hiệu hàm mũ (1/2)nu(n), 3nu(n). Hàm mũ: 푛 푛 ≥ 0 (푛) = { = 푛 (푛) 0 푛 < 0 Bài 2.3. Tính năng lượng tín hiệu x(n) = (1/2)nu(n) trong khoảng (-10,10); (0,1000); (0,1e6).
5. 40 BÀI 2: TÍN HIỆU RỜI RẠC THEO THỜI GIAN Hệ thống là tuyến tính nếu và chỉ nếu: H[a1x1(n) + a2x2(n)] = a1H[x1(n)] + a2H[x2(n)] Tính giao hoán: x(n) y(n) x(n) y(n) h1(n) h2(n) h2(n) h1(n) Tính kết hợp: x(n) y(n) x(n) y(n) h1(n) h2(n) h1(n) * h2(n) Tính phân phối: h1(n) y(n) x(n) y(n) x(n) h1(n) + h2(n) h2(n) Bài 2.8. Xét hệ thống y(n) = nx(n). >>n = -10:10; >>x = randn(size(n)); %Tín hiệu x ngẫu nhiên >>y = n.*x; %y(n) = nx(n) >>ynk = [0 0 0 0 y]; %Dịch phải y(n) 4 mẫu -> y(n – 4) >>x1 = [0 0 0 0 x]; %Dịch phải x(n) 4 mẫu >>n1 = [n 11:14]; % Bổ sung them giá trị cho n >>yn = n1.*x1; % yn = H[x(n – 4)] >>subplot(211), stem(n1,ynk), title(‘y(n – k)’); >>subplot(212), stem(n1,yn), title(‘H[x(n – k)]’); Kết luận về tính chất bất biến theo thời gian của y(n) = nx(n). Bài 2.9. Xác định tính chất bất biến theo thời gian của hệ thống có phương trình y(n) = x(-n) và y(n) = x(n)cos(0.5n). Bài 2.10.
6. 42 BÀI 2: TÍN HIỆU RỜI RẠC THEO THỜI GIAN >>y = conv(x,h); >>subplot(311),stem(x); >>subplot(312),stem(h); >>subplot(313),stem(y); Bài 2.13. Kiểm tra tính giao hoán và kết hợp: >>h1 = [1 2 -2 -3]; Hệ thống 1 >>h2 = [-2 0 3 1]; Hệ thống 2 >>h = conv(h1,h2); >>N = 30; >>x = randn(1,N); >>y11 = conv(x,h1); >>y1 = conv(y11,h2); >>y21 = conv(x,h2); >>y2 = conv(y21,h1); >>y = conv(x,h); >>subplot(311),stem(y1); >>title(‘y(n) = (x*h_1(n))*h_2(n)’); >>subplot(312),stem(y2); >>title(‘y(n) = (x*h_2(n))*h_1(n)’); >>subplot(313),stem(y); >>title(‘y(n) = x*(h_1(n)*h_2(n)’); Bài 2.14. Kiểm tra tính giao hoán và kết hợp của hai hệ thống ghép liên tầng sau: Hệ thống 1: y(n) = 2x(n) – 0.5x(n – 1) + 0.5x(n – 3) + 0.1y(n – 1) Hệ thống 2: y(n) = 0.3x(n) + 0.2x(n – 2) - 0.1y(n – 2) Bài 2.15. Xác định ngõ ra của hệ thống sau:
7. 44 BÀI 3: BIẾN ĐỔI Z BÀI 3: BIẾN ĐỔI Z Sau khi học xong bài này, người học có thể: Biết cách chuyển đổi tín hiệu trên miền thời gian rời rạc sang miền z. Biết các tính chất của biến đổi z. Biểu diễn hàm hệ thống của LTI có quan hệ vào – ra là phương trình sai phân hệ số hằng bằng biến đổi z hữu tỉ. 3.1 CÁC ĐIỂM CỰC VÀ ĐIỂM KHÔNG Biến đổi z của tín hiệu rời rạc x(n): ∞ ( ) = ∑ (푛) −푛 푛=−∞ X(z) là hàm hữu tỉ: ( ) ∑ − ( ) =0 = = − ( ) ∑ =0 Giả sử a0 0 và b0 0: 1 −1 − + + ⋯ + ( ) ( ) 0 0 0 = = − ( ) 0 + 1 −1 + ⋯ + 0 0 Do N(z) và D(z) là các đa thức theo z nên có thể biểu diễn như sau: ∏ ( ) − =1 − ( ) = ∏ =1( − ) Để biểu diễn trên đồ thị, điểm cực được đánh dấu bằng x và điểm không được đánh dấu bằng o. Bài 3.1. Xác định điểm cực và không dựa vào hàm zplane: >>num = [1 2 3]; % Tử số >>den = [2 4 7]; % Mẫu số >>zplane(num,den);
8. 46 BÀI 3: BIẾN ĐỔI Z >>[A,p,k] = residuez(num,den); Ta cũng có thể dùng hàm residuez để xác định lại tử số và mẫu số: >>[num,den] = residuez(A,p,k); Ghi lại công thức biến đổi và so sánh với kết quả ban đầu. Bài 3.5. Phân tích biểu thức sau dùng phương pháp thặng dư: 1 − 4.2 −1 + 0.8 −2 ( ) = 1 − 2.5 −1 + 3 −2 − −3 1 + −1 ( ) = (1 + 2 −1)2(1 − −1) Tính toán lại kết quả theo lý thuyết. Bài 3.6. Cho hệ thống có phương trình vào / ra là phương trình sai phân hệ số hằng: y(n) = x(n) – 2x(n – 2) + 0.81y(n – 1). Xác định H(z), từ đó viết chương trình: a. Xác định và vẽ các điểm cực, không. b. Phân tích dùng phương pháp thặng dư. 3.3 BIẾN ĐỔI Z VÀ Z NGƯỢC Xét hệ LTI biểu diễn bằng phương trình sai phân hệ số hằng: (푛) = − ∑ (푛 − ) + ∑ (푛 − ) =1 =0 Hàm hệ thống của hệ LTI biểu diễn bằng phương trình sai phân hệ số hằng: 푌( ) ∑ − ( ) =0 = = − ( ) 1 + ∑ =1 Biến đổi z ngược: − ( ) ( ) ( ) = = ∑ − + 1 ( ) ( ) =0 (Nếu bậc của tử số nhỏ hơn bậc của mẫu số) − − − ∑ ↔ ∑ 훿(푛 − ) =0 =0 Để tính thành phần còn lại, ta phân tích thành các thừa số theo phương pháp thặng dư:
9. 48 BÀI 3: BIẾN ĐỔI Z >>syms F z >>F = 2*z^(-1)/(1-3*z^(-1)); >>iztrans(F) Bài 3.13. Xác định biến đổi z ngược của các hàm sau: −1 + −2 ( ) = 1 2 −1 ( −1) (1 + 2 ) 1 + 2 2 −1 푌( ) = (3 + 4 −1 + −2)(1 + −1)
10. 50 BÀI 4: BIẾN ĐỔI FOURIER RỜI RẠC >>title(‘Pho pha arg(X(e^{j\omega}))’); Để khảo sát chuỗi x(n) hữu hạn, ta cho thông số thứ 2 của hàm freqz là 1. >>x = [1 2 3 4 5 6 7]; >>h = freqz(x,1,w); Bài 4.2. Tính và vẽ DTFT trong khoảng [- , ]: 0.8 + 3푒−푗휔 − 0.2푒−푗2휔 + 0.1푒−푗3휔 (푒−푗휔) = 1 + 0.1푒−푗휔 − 0.4푒−푗2휔 + 0.7푒−푗3휔 Bài 4.3. Khảo sát DTFT của x(n) = [1 -2 2 -3 3 4 0 -1] trong khoảng [- , ] Tính chất dịch thời gian: 퐹 (푛 − ) ↔ 푒−푗휔 (푒푗휔) Bài 4.4. Khảo sát tính chất dịch thời gian: Tính và vẽ DTFT trong khoảng [- , ] của x(n - 3) với x(n) cho như bài 4.3. Tính chất dịch tần số: 퐹 푒−푗휔0푛 (푛) ↔ (푒푗(휔−휔0)) Bài 4.5. Khảo sát tính chất dịch tần số: Tính và vẽ DTFT trong khoảng [- , ] của x(n)e-j3n với x(n) cho như bài 4.3. Tính chất đảo thời gian: 퐹 (−푛) ↔ (푒−푗휔) Bài 4.6. Khảo sát tính chất đảo thời gian: Tính và vẽ DTFT trong khoảng [- , ] của x(-n) với x(n) cho như bài 4.3. Dùng hàm fliplr để chuyển x(n) thành x(-n) và nhân thêm hệ số ej(L-1) khi biến đổi. Bài 4.7. Thực hiện lại từ 4.3 đến 4.6 với x(n) cho như bài 4.2. 4.2 FFT VÀ CÁC TÍNH CHẤT Bài 4.8. Dùng hàm fft và ifft để tính DFT và IDFT của x(n): >>N = 32; >>x = randn(1,N); >>y = fft(x,N); >>x1 = ifft(y,N);
11. 52 BÀI 4: BIẾN ĐỔI FOURIER RỜI RẠC >>y1 = fft(x1,N); >>k = 0:N-1; >>y2 = exp(-j*2*pi*k*m/N).*y; >>subplot(221),stem(abs(y));title(‘y1’); >>subplot(222),stem(angle(y)); title(‘y1’); >>subplot(223),stem(abs(y2));title(‘y2’); >>subplot(224),stem(angle(y2)); title(‘y2’); Bài 4.11. Viết chương trình khảo sát tính chất dịch vòng trên miền tần số. Bài 4.12. Viết chương trình khảo sát tính chất chập vòng. Dùng hàm cconv để tính tích chập vòng. Bài 4.13. Viết chương trình khảo sát tính chất đảo trên miền thời gian.
12. 54 BÀI 5: BỘ LỌC SỐ FIR Bài 5.1. Xác định đáp ứng tần số của bộ lọc FIR loại 1 từ chuỗi đáp ứng xung. Tạo function FIR_t1 như sau: function [a,w,L,Hr] = FIR_t1(h) M = length(h); L = (M-1)/2; a = [h(L+1) 2*h(L:-1:1)]; n = [0:1:L]; w = linspace(0,2*pi,100)’; Hr = cos(w*n)*a'; stem(Hr); Lưu file với tên FIR_t1.m. Thực hiện tính toán với đáp ứng xung h1 = [1.5 -2.5 3 -2.5 1.5]: >>h1 = [1.5 -2.5 3 -2.5 1.5]; >> [a,w,L,Hr]=FIR_t1(h1); Bài 5.2. Xác định đáp ứng tần số cho bộ lọc FIR loại 2: Viết function FIR_t2. Thực hiện tính toán với đáp ứng xung h2 = [1.5 -2.5 3 3 -2.5 1.5]. Bài 5.3. Xác định đáp ứng tần số cho bộ lọc FIR loại 3: Viết function FIR_t3. Thực hiện tính toán với đáp ứng xung h3 = [1.5 -2.5 3 2.5 -1.5]. Bài 5.4. Xác định đáp ứng tần số cho bộ lọc FIR loại 4: Viết function FIR_t4. Thực hiện tính toán với đáp ứng xung h4 = [1.5 -2.5 3 -3 2.5 -1.5]. Bài 5.5. Cho đáp ứng xung của bộ lọc FIR như sau: h = [-1 2 1.3 -2.2 0.6 3 0.6 -2.2 1.3 2 -1]. Đáp ứng xung h(n) đối xứng với M lẻ nên đây là bộ lọc FIR loại 1. a. Biểu diễn đáp ứng xung: >>h = [-1 2 1.3 -2.2 0.6 3 0.6 -2.2 1.3 2 -1]; >>M = length(h); >>n = 0:M-1; >>subplot(221); stem(n,h); >>title(‘Dap ung xung’);xlabel(‘n’);ylabel(‘h(n)’);
13. 56 BÀI 5: BỘ LỌC SỐ FIR % db = Do lon tuong doi theo dB tren doan tu 0 den pi % mag = Do lon tuyet doi tren doan tu 0 den pi % pha = Dap ung pha tren doan tu 0 den pi % grd = Tre nhom tren doan tu 0 den pi % w = Cac mau tan so doan tu 0 den pi % b = Cac he so da thuc tu so cua H(z) (voi FIR: b=h) % a = Cac he so da thuc mau so cua H(z)(voi 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); Bài 5.11. Thiết kế bộ lọc thông thấp theo phương pháp cửa sổ Hamming với các tham số như sau: p = 0.2 ; s = 0.3 ; Rp = 0.25 dB; As = 50 dB Việc thiết kế bộ lọc là quá trình tìm ra các tham số, hay chuỗi đáp ứng xung của bộ lọc, thoả mãn các yêu cầu chỉ tiêu kỹ thuật cho trước, cụ thể là một số hoặc tất cả các tham số tuyệt đối (absolute specification) sau: • Tần số cắt dải thông ω p • Tần số cắt dải thông ω s • Bề rộng dải quá độ Δω • Độ gợn sóng dải thông δ p • Độ gợn sóng dải chắn δ s Các tham số thường được cho dưới dạng đơn vị dB: • Độ gợn sóng dải thông theo dB: 1−훿 Rp = −20푙표 1+훿
14. 58 BÀI 5: BỘ LỌC SỐ FIR >>title('Dap ung tan so');xlabel('n'); ylabel('h(n)'); d. Đáp ứng tần số theo dB: >>subplot(224); plot(w,db); >>axis([0,1,-100,10]); >>title('Dap ung tan so theo dB');xlabel('\omega'); ylabel('dB'); Bài 5.12. Thiết kế bộ lọc thông dải theo phương pháp cửa sổ Blackman với các tham số như sau: s1 = 0.2 ; p1 = 0.35 ; p2 = 0.65 ; s2 = 0.8 ; Rp = 1 dB; As = 60 dB Quá trình thực hiện như bài 5.11. Về lý thuyết, đáp ứng xung lý tưởng của bộ lọc thông dải lý tưởng là hiệu đáp ứng xung của hai bộ lọc thông thấp lý tưởng. Dùng hàm idea_lp như bài 5.9 để xác định đáp ứng xung của bộ lọc thông dải lý tưởng trong đó các tần số cắt có thể chọn là trung bình của các tần số cắt dải thông và dải chắn. Thông số độ rộng dải chuyển tiếp để tính chiều dài chuỗi đáp ứng xung có thể chọn là giá trị nhỏ nhất của độ rộng hai dải chuyển tiếp, từ dải chắn lên dải thông [s1,p1] và từ dải thông xuống dải chắn [p2,s2]. 5.3 PHƯƠNG PHÁP LẤY MẪU TẦN SỐ Xét đáp ứng xung h(n) của hệ thống LTI có đáp ứng tần số là: −1 −푗휔푛 H() = ∑푛=0 ℎ(푛)푒 Ta chỉ định một tập tần số như sau: −1 = 0, 1, , 푛ế 푙ẻ 2 2 k = ( + 훼), = 0, 1, , − 1 푛ế ℎẵ푛 2 1 훼 = 0 ℎ { 2 Khi đó: 2 2 −푗 ( +훼)푛 H(k + α) = H( ( + 훼)) = ∑ −1 ℎ(푛)푒 푛=0 Tập hợp các giá trị {H(k + α)} gọi là các mẫu tần số của H(). Trong trường hợp α = 0, {H(k)} tương ứng với DFT M điểm của h(n).
15. 60 BÀI 5: BỘ LỌC SỐ FIR >>axis([-1,M,-0.1,0.3]); >>title('Dap ung xung'); >>xlabel('n'); ylabel('h(n)'); c. Biên độ của đáp ứng tần số: >>subplot(223); plot(ww/pi,Hr,wl(1:31)/pi,Hk(1:31),'o'); >>axis([0,1,-0.2,1.2]); >>title('Bien do cua dap ung tan so'); >>xlabel('f[*\pi]'); ylabel('Hr(w)'); d. Biên độ của đáp ứng tần số theo dB: >>subplot(224); plot(w/pi,db); >>axis([0,1,-100,10]); grid >>title('Bien do cua dap ung tan so '); >>xlabel('f[*\pi]'); ylabel('dB'); Bài 5.14. Thiết kế bộ lọc thông dải theo phương pháp lấy mẫu tần số với các tham số như sau: s1 = 0.2 ; p1 = 0.35 ; p2 = 0.65 ; s2 = 0.8 ; Rp = 1 dB; As = 60 dB Chọn đáp ứng xung có chiều dài 40 ứng với 40 mẫu tần số trong khoảng [0,2 ). Dải thông có độ rộng là 0.3 tương đương với 7 mẫu nhận giá trị 1. Giả sử quá trình tối ưu hoá chỉ ra nên chọn dải chuyển tiếp 2 mẫu nhận các giá trị T1 = 0. 109021 và T2 = 0.59417456. Mẫu các tần số được cho như sau: H(k) = [0,0,0,0,0,T1,T2,1,1,1,1,1,1,1,T2,T1,0, ,0,T1,T2,1,1,1,1,1,1,1,T2,T1,0,0,0,0] (9 số 0) Quá trình thực hiện tương tự bài 5.13. 5.4 PHƯƠNG PHÁP LẶP Đáp ứng tần số của 4 loại bộ lọc FIR: H() = P()Q() Hàm sai số giữa bộ lọc thực tế và bộ lọc lý tưởng:
16. 62 BÀI 5: BỘ LỌC SỐ FIR  Xem xét xem các giá trị rời rạc được chọn ban đầu có thực sự là các điểm mà hàm sai số E(ω) đạt cực trị và có trị tuyệt đối bằng nhau hay không. Nếu không, tìm các điểm tại đó E(ω) đạt cực trị.  Trong các điểm cực trị của E(ω) lấy ra R+2 điểm và quay về lặp lại từ bước 2.  Lặp lại các bước 2, 3, và 4 cho đến khi tập hợp các điểm rời rạc hội tụ.  Từ tập các điểm rời rạc cuối cùng, tính ra hàm P(ω), từ đó tính ra các hệ số của P(ω). 4. Tính các giá trị của chuỗi đáp ứng xung h(n). Khi chọn giá trị M càng chuẩn thì kết quả là thu được bộ lọc có hàm đáp ứng tần số càng gần với yêu cầu bài toán. Nếu như với giá trị M nào đó mà chưa thoả mãn được yêu cầu thì phải tăng giá trị M đến khi nào thoả mãn các điều kiện ràng buộc cho δp và δs (hay As và Rp). Một công thức lựa chọn ban đầu cho chiều dài M của đáp ứng xung là: −20푙표 √훿 훿 − 13 휔 − 휔 = 1 2 푣ớ𝑖 Δ = 푠 0 14.6Δ 2 Trong MATLAB, tìm đáp ứng xung của bộ lọc tối ưu với giá trị M và hàm đáp ứng tần số lý tưởng cho trước được thực hiện bởi hàm firpm. Bài 5.15. Tạo script file sau để biểu diễn bộ lọc trên đồ thị: wp = 0.2*pi; ws =0.3*pi; Rp = 0.25; As = 50; delta_w = 2*pi/1000; wsi = ws/delta_w+1; delta1 = (10^(Rp/20)-1)/(10^(Rp/20)+1); delta2 = (1+delta1)*(10^(-As/20)); deltaH = max(delta1,delta2); deltaL = min(delta1,delta2); weights = [delta2/delta1 1]; deltaf = (ws-wp)/(2*pi); M = ceil((-20*log10(sqrt(delta1*delta2))-13)/(14.6*deltaf)+1) f = [0 wp/pi ws/pi 1]; m = [1 1 0 0];
17. 64 BÀI 5: BỘ LỌC SỐ FIR title('Dap ung tan so cua loi'); xlabel('f[*\pi]'); ylabel('Er(\omega)'); Bài 5.16. Thực hiện như bài 5.15 cho bộ lọc thông dải theo phương pháp lấy mẫu tần số với các tham số như sau: s1 = 0.2 ; p1 = 0.35 ; p2 = 0.65 ; s2 = 0.8 ; Rp = 1 dB; As = 60 dB
18. 66 BÀI 6: BỘ LỌC SỐ IIR 1 푠 푠 = −10푙표 → = 1020 2 1 − 훿 1 2√훿 √ = 2 → 휖 = 1 + 훿 1 + 휖 1 − 훿 훿 1 1 + 훿 푠 = → = 1 + 훿 훿푠 Có 4 định dạng cơ bản thường được vận dụng trong quá trình thiết kế bộ lọc tương tự là: bộ lọc Butterworth, bộ lọc Chebyshev-1, bộ lọc Chebyshev-2 và bộ lọc Elliptic. Các hàm MATLAB có thể sử dụng cho bài thực hành này là: freqs: trả về đáp ứng tần số của một hệ thống tương tự khi biết hàm truyền dưới dạng phân thức hữu tỷ. impulse: trả về đáp ứng xung của một hệ thống tương tự khi biết hàm truyền dưới dạng phân thức hữu tỷ buttap, cheb1ap, cheb2ap, ellipap: trả về các điểm không, điểm cực, và độ lợi trong thiết kế của một hàm truyền bộ lọc thông thấp bậc N, tần số cắt đã được chuẩn hoá bằng 1 với các định dạng lần lượt là Butterworth, Chebyshev-I, Chebyshev-II, và Elliptic. impinvar, bilinear: trả về các hệ số của đa thức tử số và đa thức mẫu số hàm truyền đạt của hệ thống số xuất phát từ hệ thống tương tự qua các phương pháp chuyển đổi bất biến xung và song tuyến tính. butter, cheby1, cheby2, ellip: trả về các hệ số của đa thức tử số và đa thức mẫu số hàm truyền của bộ lọc số dựa trên tham số đầu vào là các tần số cắt, phương pháp chuyển đổi được sử dụng trong các hàm này là phương pháp biến đổi song tuyến tính. Hàm freqs của MATLAB trả về đáp ứng tần số của một hệ thống tương tự khi biết trước hệ số của đa thức tử số và đa thức mẫu số của hàm truyền đạt Ha(s). Trong nhiều trường hợp, để thuận tiện ta cần tìm thêm các thông số: hàm độ lớn của đáp ứng tần số, hàm pha của đáp ứng tần số, hàm trễ nhóm, thể hiện độ lớn theo thang decibels. Bài 6.1. Tạo hàm tính đáp ứng tần số freqs_m nhằm tính các thông số trên: function [db,mag,pha,w] = freqs_m(b,a,wmax); % db = Do lon tuong doi theo dB tren doan tu 0 den wmax % mag = Do lon tuyet doi tren doan tu 0 den wmax % pha = Dap ung pha tren doan tu 0 den wmax
19. 68 BÀI 6: BỘ LỌC SỐ IIR % [b,a] = u_chb1ap(N,Rp,Omegac) % b = cac he so da thuc tu so cua Ha(s) % a = cac he so da thuc mau so cua Ha(s) % N = Bac cua bo loc Chebyshev-I % Rp = Do gon dai thong theo don vi dB; Rp > 0 % Omegac = tan so cat theo don vi radians/sec [z,p,k] = cheb1ap(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; B = real(poly(z)); b0 = k; b = k*B; Hàm số mô tả ở trên trả về hàm truyền với bậc N cho trước. Bậc của bộ lọc có thể lựa chọn cho phù hợp tối ưu với các chỉ tiêu kỹ thuật yêu cầu đầu vào. Bài 6.3. Tạo hàm afd_chb1 trả về thiết kế bộ lọc thông thấp tương tự, định dạng Chebyshev có bậc tối ưu: function [b,a] = afd_chb1(Wp,Ws,Rp,As) % Analog Lowpass Filter Design: Chebyshev-1 % [b,a] = afd_chb1(Wp,Ws,Rp,As) % b = cac he so da thuc tu so cua Ha(s) % a = cac he so da thuc mau so cua Ha(s) % Wp = tan so cat dai thong theo don vi rad/sec; Wp >0 % Ws = tan so cat dai chan theo don vi rad/sec; Ws>Wp >0 % Rp = Do gon dai thong theo don vi dB; (Rp > 0) % As = Do suy giam dai chan theo don vi +dB; (Ap > 0) if Wp <= 0
20. 70 BÀI 6: BỘ LỌC SỐ IIR 2 1 − −1 푠 = 1 + −1 Với bộ lọc tương tự cho trước có hàm hệ thống Ha(s), bộ lọc số được thiết kế như sau: 2 1 − −1 ( ) = ( ) 1 + −1 Bài 6.5. Chuyển đổi bộ lọc với các tham số đã cho ở bài 6.4 sang bộ lọc số bằng phương pháp biến đổi song tuyến. Hàm bilinear cho phép thực hiện việc chuyển đổi này. Tạo script file bai605.m như sau: wp =0.2*pi; % digital Passband freq in Hz ws =0.3*pi; % digital Stopband freq in Hz Rp = 1; % Passband ripple in dB As = 15; % Stopband attenuation in dB T = 1; Fs =1/T; % Dat T=1 OmegaP = (2/T)*tan(wp/2); OmegaS = (2/T)*tan(ws/2); % Tinh toan bo loc tuong tu: [cs, ds] = afd_chb1(OmegaP,OmegaS,Rp,As); % Bien doi song tuyen tinh: [b,a] = bilinear(cs,ds,Fs); [db,mag,pha,grd,w] = freqz_m(b,a); a. Đáp ứng tần số: subplot(221); plot(w/pi,mag); axis([0,1,0,1.2]); grid title('Dap ung tan so'); xlabel('f[*\pi]'); ylabel('|H_r(\omega)|'); b. Hàm độ lớn tương đối tính theo dB: subplot(222); plot(w/pi,db); axis([0,1,-30,10]); grid title('Dap ung tan so theo dB'); xlabel('f[*\pi]'); ylabel('dB');
21. 72 BÀI 6: BỘ LỌC SỐ IIR 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; Bài 6.8. Viết chương trình chuyển đổi từ bộ lọc thông thấp theo thiết kế của bài 6.5 sang bộ lọc thông cao có tần số cắt ωc=0.6 . Tính và biểu diễn trên đồ thị: a. Độ lớn của đáp ứng tần số b. Hàm đáp ứng pha của bộ lọc c. Hàm độ lớn tương đối tính theo dB của đáp ứng tần số d. Trễ nhóm theo tần số.