Đề tài Thiết kế bộ lọc FIR thông dải bằng phương pháp cửa sổ

ong đó: Band [0, wp] được gọi là dải thông, và d1 là dung sai (gợn sóng) được chấp nhận trong đáp ứng dải thông lý tưởng. Band [ws, p] được gọi là dải chắn, và d2 là dung sai ở dải chắn. Band [wp, ws] được gọi là dải chuyển tiếp, và không có ràng buộc nào về đáp ứng biên độ trong dải này Các chỉ tiêu tương đối gồm có: Rp: Độ gợn sóng trong dải thông tính theo dB. As : Suy hao trong dải chắn tính theo dB. Quan hệ giữa các chỉ tiêu này như sau: (»0) (1.1) (>>1) (1.2) Các chỉ tiêu trên được đưa ra đối với bộ lọc FIR thông thấp, và tất nhiên đối với các bộ lọc khác như thông cao HPF (High Pass Filter), thông dải BPF (Band Pass Filter) đều có thể được định nghĩa tương tự. Tuy nhiên, các tham số thiết kế quan trọng nhất là các dung sai dải tần và các tần số cạnh-dải (tolerance or ripples and band-edge frequencies). Bởi vậy, trong phần 1 về cơ sở lý thuyết này chúng ta chỉ tập trung vào bộ lọc FIR thông thấp. Việc thiết kế cụ thể cho bộ lọc FIR thông dải bằng kỹ thuật cửa sổ sẽ được phát triển trên cơ sở lọc thông thấp và sẽ được mô tả chi tiết trong phần 2. Việc thiết kế và thực hiện lọc FIR có những thuận lợi sau đây: Đáp ứng pha là tuyến tính. Dễ thiết kế do không gặp các vấn đề ổn định (lọc FIR luôn ổn định). Việc thực hiện rất hiệu quả . Có thể sử dụng DFT để thực hiện Đáp ứng pha là tuyến tính (linear phase response) mang lại những thuận lợi sau: Bài toán thiết kế chỉ gồm các phép tính số học thực chứ không cần phép tính số học phức. Bộ lọc pha tuyến tính không có méo trễ nhóm và chỉ bị trễ một khoảng không đổi. Đối với bộ lọc có chiều dài M (hoặc bậc M-1) số phép toán có bậc M/2 như đã khảo sát trong thực hiện bộ lọc có pha tuyến tính. 1.2. Cấu trúc của bộ lọc FIR Một bộ lọc đáp ứng xung hữu hạn với hàm hệ thống có dạng: (1.3) Như vậy đáp ứng xung h(n) là: (1.4) Và phương trình sai phân là: (1.5) Đây chính là tích chập tuyến tính của các dãy hữu hạn. Bậc của bộ lọc là M-1, trong khi chiều dài của bộ lọc là M (bằng với số lượng các hệ số). Các cấu trúc bộ lọc FIR luôn luôn ổn định, và tương đối đơn giản hơn so với các cấu trúc bộ lọc IIR. Hơn thế nữa, các bộ lọc FIR có thể được thiết kế để có một đáp ứng pha tuyến tính và đó là điều cần thiết trong một số ứng dụng. Chúng ta sẽ xem xét lần lượt các cấu trúc của bộ lọc FIR sau đây: a. Cấu trúc dạng trực tiếp Phương trình sai phân được thực hiện bởi một dãy liên tiếp các bộ trễ do không có đường phản hồi: (1.6) Do mẫu thức bằng đơn vị nên ta chỉ có một cấu trúc dạng trực tiếp duy nhất. Cấu trúc dạng trực tiếp được cho trong hình 1.2 với M = 5: b0 z-1 b1 z-1 b2 z-1 b3 z-1 b4 y(n) x(n) Hình 1.2 Cấu trúc lọc FIR dạng trực tiếp b. Cấu trúc dạng ghép tầng: Hàm hệ thống H(z) được biến đổi thành các tích của các khâu bậc 2 với các hệ số thực. Các khâu này được thực hiện ở dạng trực tiếp và bộ lọc tổng thể có dạng ghép tầng của các khâu bậc 2. (1.7) B1,1 z-1 z-1 z-1 y(n) x(n) B2,1 B3,1 b0 B1,2 z-1 z-1 z-1 B2,2 B3,2 Hình 1.3 Cấu trúc lọc FIR dạng ghép tầng trong đó , Bk,1 và Bk,2 là các số thực đại diện cho các hệ số của các khâu bậc 2. Cấu trúc dạng ghép tầng được cho trong hình 1.3 với M = 7: c. Cấu trúc dạng pha tuyến tính: Đối với các bộ lọc chọn tần, người ta mong muốn có đáp ứng pha là hàm tuyến tính theo tần số, nghĩa là: (1.8) trong đó hoặc và là một hằng số. Đối với bộ lọc FIR nhân quả có đáp ứng xung trong khoảng [0, M-1], thì các điều kiện tuyến tính là: (1.9) (1.10) Xét phương trình sai phân được cho trong phương trình (1.5) với đáp ứng xung đối xứng trong phương trình (1.9), ta có: Sơ đồ khối thực hiện phương trình sai phân trên được mô tả trong hình 1.4 dưới đây đối với cả M lẻ và M chẵn: b0 z-1 b1 z-1 b2 x(n) z-1 z-1 z-1 y(n) b0 z-1 b1 z-1 b2 b3 y(n) x(n) z-1 z-1 z-1 z-1 M=7 M=6 Hình 1.4 Cấu trúc lọc FIR pha tuyến tính với các hệ số M chẵn và lẻ Đối với M lẻ: M = 7, còn đối với M chẵn: M = 6 Rõ ràng, với cùng một bậc của bộ lọc (cùng M) cấu trúc pha tuyến tính sẽ tiết kiệm được 50% các bộ nhân so với cấu trúc dạng trực tiếp. 1.3. Các đặc tính của bộ lọc FIR pha tuyến tính Trong phần này chúng ta sẽ thảo luận về hình dạng của đáp ứng xung, đáp ứng tần số trong hàm hệ thống của các bộ lọc FIR pha tuyến tính. Cho h(n), trong đó 0 £ n £ M – 1, là đáp ứng xung có chiều dài M thì hàm truyền hệ thống là: (1.11) có (M-1) điểm cực ở gốc (trivial poles) và M-1 điểm không nằm ở vị trí bất kỳ trên mặt phẳng z. Đáp ứng tần số là: (1.12) a. Đáp ứng xung h(n): Chúng ta có thể đưa ra ràng buộc pha tuyến tính: (1.13) trong đó: a là một hằng số trễ pha. Ta đã biết rằng h(n) phải đối xứng, nghĩa là: (1.14) Do đó h(n) là đối xứng theo a, là chỉ số đối xứng. Có hai kiểu đối xứng: Hình 1.5 Đáp ứng xung đối xứng, M lẻ M lẻ: Trong trường hợp này, là một số nguyên. Đáp ứng xung được mô tả trong hình 1.5 dưới đây: M chẵn: Trong trường hợp này, không phải là một số nguyên. Đáp ứng xung được mô tả bằng hình 1.6 dưới đây: Hình 1.6 Đáp ứng xung đối xứng, M chẵn Ta cũng có bộ lọc FIR pha tuyến tính loại hai nếu ta yêu cầu đáp ứng pha thoả mãn điều kiện: với (1.15) Đáp ứng pha là đường thẳng nhưng không đi qua gốc. Trong trường hợp này a không phải là hằng số trễ pha, nhưng: (1.16) là hằng số, chính là trễ nhóm (a là một hằng số trễ nhóm). Trong trường hợp này, các tần số được làm trễ với một tốc độ không đổi. Nhưng một số tần số có thể được làm trễ với tốc độ lớn hơn hoặc nhỏ hơn. Đối với kiểu pha tuyến tính này, có thể thấy rằng: và (1.17) Điều này có nghĩa rằng đáp ứng xung h(n) là phản đối xứng (antisymmetric). Chỉ số đối xứng vẫn là . Một lần nữa chúng ta lại có 2 kiểu, cho M lẻ và M chẵn. Hình 1.7 Đáp ứng xung phản đối xứng, M lẻ M lẻ: Trong trường hợp này, là một số nguyên. Đáp ứng xung được mô tả bằng hình 1.7 dưới đây: Lưu ý rằng mẫu h(a) tại phải bằng 0, nghĩa là, . Hình 1.8 Đáp ứng xung phản đối xứng, M chẵn M chẵn: Trong trường hợp này, không phải là một số nguyên. Đáp ứng xung được mô tả trong hình 1.8. b. Đáp ứng tần số H(ejw): Khi tổ hợp hai loại đối xứng và phản đối xứng với M chẵn và M lẻ, ta có bốn kiểu lọc FIR pha tuyến tính. Đáp ứng tần số của mỗi kiểu có biểu thức và hình dạng riêng. Để nghiên cứu các đáp ứng pha của các kiểu này, ta viết biểu thức của H(ejw) như sau: (1.18) trong đó Hr(ejw) là hàm đáp ứng độ lớn chứ không phải là hàm đáp ứng biên độ. Đáp ứng độ lớn là một hàm thực, có thể vừa dương vừa âm, không giống đáp ứng biên độ luôn luôn dương. Đáp ứng pha kết hợp với đáp ứng biên độ là một hàm không liên tục, trong khi kết hợp với đáp ứng độ lớn là một hàm tuyến tính liên tục. Bộ lọc FIR pha tuyến tính Loại-1 (Type 1): Đáp ứng xung đối xứng, M lẻ Trong trường hợp này , là một biến nguyên, và , , thì ta có thể chứng tỏ rằng: (1.19) trong đó: với mẫu ở chính giữa (1.20) với Bộ lọc FIR pha tuyến tính Loại-2 (Type 2): Đáp ứng xung đối xứng, M chẵn Trong trường hợp này , , , nhưng không phải là một biến nguyên, thì ta có thể chứng tỏ rằng: (1.21) trong đó: với (1.22) So sánh (1.21) và (1.18), ta có: (1.23) Lưu ý: Tại , ta có mà không cần quan tâm đến b(n) hoặc h(n). Do đó chúng ta không thể sử dụng loại này (h(n) đối xứng, M chẵn) đối với bộ lọc thông cao hoặc bộ lọc chắn dải. Lọc FIR pha tuyến tính Loại-3 (Type 3): Đáp ứng xung phản đối xứng, M lẻ Trong trường hợp này ta có , là một biến nguyên, , , và thì ta có thể chứng tỏ: (1.24) trong đó với (1.25) So sánh (1.24) và (1.18), ta có: (1.26) Lưu ý: Tại và , ta có mà không cần quan tâm c(n) hoặc h(n). Hơn thế nữa, , điều đó có nghĩa là là thuần ảo. Do đó, loại bộ lọc này không thích hợp đối với việc thiết kế bộ lọc thông thấp hoặc thông cao. Tuy nhiên, điều này thích hợp đối với việc xấp xỉ các bộ vi phân và bộ biến đổi Hilbert số lý tưởng. Lọc FIR pha tuyến tính Loại-4(Type 4):Đáp ứng xung phản đối xứng, M chẵn Trong trường hợp này , , , nhưng không phải là một biến nguyên, thì ta có thể chứng tỏ rằng: (1.27) trong đó: với (1.28) So sánh (1.27) và (1.18), ta có: (1.29) Lưu ý: Tại , và . Do vậy, loại này cũng thích hợp cho việc thiết kế các bộ vi phân số và bộ biến đổi Hilbert số. Bảng sau đây mô tả khả năng thích hợp trong việc thiết kế các bộ lọc và các bộ biến đổi Hilbert số, bộ vi phân số của 4 loại lọc FIR pha tuyến tính đã nêu: Type LPF HPF BPF SBF Hilbert Differentiator FIR Type 1 ü ü ü ü FIR Type 2 ü ü FIR Type 3 ü ü ü FIR Type 4 ü ü ü ü 1.4. Các kỹ thuật thiết kế cửa sổ Ý tưởng cơ bản của việc thiết kế là: chọn một bộ lọc chọn tần lý tưởng (mà đáp ứng xung luôn luôn phi nhân quả, dài vô hạn) và cắt (lấy cửa sổ - window) đáp ứng xung của nó để thu được bộ lọc FIR có pha tuyến tính và nhân quả (linear phase and causal FIR filter). Bởi vậy, điểm quan trọng trong phương pháp này là việc chọn một hàm cửa sổ thích hợp và một bộ lọc lý tưởng tương ứng. Bộ lọc thông thấp lý tưởng (ideal LPF) có tần số cắt wc < p được cho bởi: (1.30) (1.31) Đáp ứng xung hd(n) của bộ lọc lý tưởng này được cho bởi: Chú ý rằng hd(n) là đối xứng theo a, sự kiện này tiện dụng cho bộ lọc FIR có pha tuyến tính. Ngoài ra hd(n) có độ dài vô hạn và phi nhân quả. (1.32) Để thu được bộ lọc FIR có pha tuyến tính và nhân quả h(n) có độ dài M, ta cần có: Thao tác này được gọi là lấy cửa sổ (window). Tổng quát, đáp ứng xung h(n) có thể có được bằng cách lấy đáp ứng xung của bộ lọc lý tưởng hd(n) nhân với hàm cửa sổ w(n) như sau: h(n) = hd(n).w(n) (1.33) nếu 0£ n £ M-1 n khác (1.34) Tuỳ thuộc vào cách định nghĩa hàm cửa sổ w(n) (window function) chúng ta có những cửa sổ thiết kế khác nhau. Chẳng hạn nếu dùng cửa sổ chữ nhật (Rectangular), thì hàm cửa sổ được định nghĩa: Trong miền tần số, đáp ứng của lọc FIR nhân quả H(ejw) chính là tích chập vòng của đáp ứng tần số bộ lọc lý tưởng Hd(ejw) và đáp ứng tần số của hàm cửa sổ W(ejw) : H(ejw) = Hd(ejw) Ä W(ejw) (1.35) H(ejw) Hd(ejw) W(ejw) Max side-lobe height Transition bandwidth Hình 1.9 Kết quả của việc lấy cửa sổ trong miền tần số Main-lobe width Hình dạng của H(ejw) có thể được mô tả trực quan trên hình 1.9. Từ hình vẽ này chúng ta có một số nhận xét quan trọng sau đây: Do cửa sổ w(n) có chiều dài M hữu hạn, đáp ứng của nó có một búp chính (main-lobe) có độ rộng tỷ lệ với 1/M, và các búp bên (side-lobe) của nó có chiều cao thấp hơn. Tích chập tuần hoàn sinh ra một phiên bản méo của đáp ứng xung lý tưởng Hd(ejw). Búp chính (main-lobe) sinh ra một dải chuyển tiếp trong H(ejw) mà độ rộng là nguyên nhân tạo nên độ rộng dải chuyển tiếp (transition bandwidth). Độ rộng này tỷ lệ với 1/M. Độ rộng búp chính càng lớn thì độ rộng dải chuyển tiếp càng lớn. Các búp bên (side-lobes) sinh ra các gợn sóng có hình dạng như nhau trong cả dải thông và dải chắn. Với các chỉ tiêu bộ lọc đã cho, chọn chiều dài bộ lọc và một hàm cửa sổ w(n) với độ rộng main-lobe hẹp nhất và hệ số suy giảm side-lobe bé nhất có thể được. Từ nhận xét 4 nêu trên ta chú ý rằng dung sai dải thông d1 và dung sai dải chắn d2 không thể ấn định một cách độc lập. Ta lấy chung d1 = d2. Tiếp theo chúng ta xem xét các hàm cửa sổ thường được dùng, bao gồm: cửa sổ chữ nhật (Rectangular), cửa sổ tam giác (Bartlett), cửa sổ Hanning, cửa sổ Hamming, cửa sổ Blackman và cửa sổ Kaiser. a. Cửa sổ chữ nhật (Rectangular Window) w(n) = 1 0£ n £ M-1 0 Các trường hợp khác (1.36) Trong miền n cửa sổ chữ nhật được định nghĩa như sau: b. Cửa sổ tam giác (Bartlett Window) Với mục đích giảm biên độ của các đỉnh thứ cấp của cửa sổ chữ nhật, chúng ta chọn một cửa sổ khác có dạng tam giác cân, gọi là cửa sổ tam giác hay cửa sổ Bartlett Trong miền n cửa sổ tam giác được định nghĩa như sau: 0 Các trường hợp khác w(n) = 0£ n £ 2 - £ n £ M - 1 2n M-1 M-1 2 2n M-1 M-1 2 (1.37) c. Cửa sổ Hanning (Hanning Window) w(n) = 0 Các trường hợp khác 0£ n £ M (1.38) d.Cửa sổ Hamming (Hamming Window) w(n) = 0 Các trường hợp khác 0£ n £ M (1.39) e. Cửa sổ Blackman (Blackman Window) w(n) = 0 Các trường hợp khác 0£ n £ M (1.40) Hình 1.10 Hình dạng một số cửa sổ thường dùng Hình 1.11 Tóm tắt đặc tính của một số loại cửa sổ thường dùng Trên hình 1.11 cho chúng ta một sự so sánh giữa các hàm cửa sổ thường dùng về các đặc tính: độ rộng dải chuyển tiếp Dw, độ suy giảm ở dải chắn As. f. Cửa sổ Kaiser (Kaiser Window) Để đạt được độ suy giảm của dải chặn như mong muốn, các nhà thiết kế tìm một hàm cửa sổ đáp ứng được các yêu cầu của thiết kế. Nhưng các hàm cửa sổ có mức búp bên càng thấp thì độ rộng của búp chính lại càng lớn dẫn đến độ rộng dải chuyển tiếp cũng tăng. Do đó phải tăng bậc của bộ lọc để đạt được dải thông mong muốn. Cửa sổ Kaiser có thông số b có thể điều chỉnh được, do vậy có thể điều chỉnh được độ rộng búp bên so với đỉnh của búp chính. Cũng giống các hàm cửa sổ khác, độ rộng búp chính có thể thay đổi được bằng cách điều chỉnh chiều dài cửa sổ, do vậy điều chỉnh được độ rộng của dải chuyển tiếp. Với mục tiêu này, các bộ lọc số được thiết kế rất có hiệu quả khi dùng hàm cửa sổ Kaiser. Cửa sổ Kaiser được định nghĩa như sau: , với 0£ n £ M-1 (1.41) Trong đó I0[.] là hàm Bessel bậc không được hiệu chỉnh (modified zero-order Bessel function), còn b là tham số phụ thuộc vào bậc bộ lọc M và có thể chọn để có được những độ rộng dải chuyển tiếp khác nhau cũng như có được độ suy giảm ở dải chặn là gần tối ưu. Cửa sổ này có thể cung cấp những độ rộng dải chuyển tiếp khác nhau với cùng một hệ số M, đây là điều mà các hàm cửa sổ khác không làm được. Tuy nhiên, do tính phức tạp của hàm Bessel nên việc thiết kế bộ lọc với cửa sổ Kaiser không phải dễ dàng. Thật may mắn là Kaiser đã phát triển những phương trình thiết kế theo kinh nghiệm (empirical design equation) và chúng ta có thể sử dụng mà không cần phải chứng minh. Phương trình thiết kế như sau: (đối với lọc thông thấp LPF) Cho trước các chỉ tiêu của bộ lọc cần thiết kế: ws, wp, Rp và As. Độ rộng dải chuyển tiếp: 0.1102(As - 8.7), nếu As ³ 50 0.5842 (As - 21)0.4 + 0.07886(As – 21), nếu 21 < As < 50 Bậc của bộ lọc: (1.42) Tham số: b = Hình 1.12 Khảo sát đặc tính cửa sổ Kaiser Để tổng kết phần này chúng ta cùng so sánh giữa các loại cửa sổ thông dụng thường được dùng, cũng như khảo sát sự tương ứng với cửa sổ Kaiser với hệ số b khác nhau, như trên hình 1.13 dưới đây. Ta có nhận xét rằng: Ngoại trừ cửa sổ Kaiser, đối với các loại cửa sổ khác thì độ gợn sóng (ripple) ở dải thông và dải chắn là như nhau và không phụ thuộc vào bậc M của bộ lọc, và chỉ có thể thay đổi bằng cách thay đổi hình dạng của cửa sổ (nghĩa là thay đổi loại cửa sổ sử dụng). Hình 1.13 So sánh giữa các loại cửa sổ được dùng Thực hiện bằng Matlab: Matlab cung cấp sẵn các hàm (routine) để thực hiện các hàm cửa sổ chúng ta vừa khảo sát. Có thể mô tả ngắn gọn các hàm đó như sau: W=boxcar(M) : cho một hàm cửa sổ chữ nhật M điểm trong mảng W. W=triang(M) : cho một hàm cửa sổ tam giác M điểm trong mảng W. W=hanning(M) : cho một hàm cửa sổ Hanning M điểm trong mảng W. W=hamming(M): cho một hàm cửa sổ Hamming M điểm trong mảng W. W=blackman(M):cho một hàm cửa sổ Blackman M điểm trong mảng W. W=kaiser(M,beta):cho một hàm cửa sổ Kaiser với hệ số b gồm M điểm trong mảng W. Phần 2. THIẾT KẾ LỌC FIR THÔNG DẢI 2.1. Bài toán thiết kế Hãy thiết kế bộ lọc FIR thông dải pha tuyến tính theo phương pháp cửa sổ, với các chỉ tiêu bộ lọc cần thiết kế được cho như sau: Cạnh thấp của dải chắn: ws1 Cạnh thấp của dải thông: wp1 Cạnh cao của dải thông: wp2 Cạnh cao của dải chắn: ws2 Độ gợn sóng trong dải thông: Rp Suy hao trong dải chắn: As As Rp ws1 wp1 wp2 ws2 p 0 Hình 2.1 Các chỉ tiêu của bộ lọc thông dải BPF (bandpass filter) Các đại lượng này có thể được mô tả trên hình 2.1 như sau: Điều kiện: ws1 < wp1 < wp2 < ws2 Hoặc bài toán cho các chỉ tiêu d1 và d2 ta cũng có thể tính được As và Rp dựa vào quan hệ giữa chúng theo công thức (1.1) và (1.2). 2.2. Phương pháp thiết kế Bước 1. Chọn loại cửa sổ sử dụng Việc chọn loại cửa sổ sử dụng nhằm đảm bảo suy hao trong dải chắn thoả mãn chỉ tiêu thiết kế As yêu cầu của bài toán đặt ra. Có đến 6 loại cửa sổ khác nhau, trong đó ngoại trừ cửa sổ Kaiser có suy hao trong dải chắn có thể đạt được bằng cách thay đổi hệ số b, còn các cửa sổ khác có suy hao trong dải chắn là cố định. Ngoài ra, bậc M của bộ lọc phụ thuộc vào độ rộng của dải chuyển tiếp Dw và phụ thuộc vào loại cửa sổ được chọn. Với cùng 1 chỉ tiêu thiết kế, tức là cùng 1 giá trị độ rộng dải chuyển tiếp Dw, bậc của bộ lọc sẽ khác nhau nếu chọn các cửa sổ khác nhau. Bởi vậy, việc chọn loại cửa sổ nào còn phụ thuộc vào quan điểm của người thiết kế trên cơ sở dung hoà giữa việc đảm bảo được suy hao trong dải chắn theo đúng yêu cầu bài toán và đảm bảo bậc bộ lọc đủ nhỏ. Tuy nhiên, thông số được ưu tiên ở đây là suy hao As, bởi sự thay đổi của bậc bộ lọc M khi sử dụng các loại cửa sổ khác nhau là không nhiều lắm. Bước 2. Xác định bậc M của bộ lọc Nếu sử dụng cửa sổ Kaiser: M được tính theo nhóm công thức (1.42), cũng từ công thức này ta tính được hệ số b. Nếu sử dụng các cửa sổ khác: M được xác định nhờ vào quan hệ giữa M với độ rộng dải chuyển tiếp Dw và có thể tính được dựa vào bảng tóm tắt được cho trong hình 1.11. Bước 3. Tìm hàm cửa sổ w(n) Sử dụng các hàm có sẵn của Matlab với bậc bộ lọc M đã tìm được ở bước 2. Bước 4. Tìm đáp ứng xung của bộ lọc thông dải lý tưởng hd(n) Å 0 p wc1 0 p wc2 0 p wc2 wc1 - + Hình 2.2 Đáp ứng xung lọc thông dải lý tưởng từ 2 bộ lọc thông thấp lý tưởng Đáp ứng xung của bộ lọc thông dải lý tưởng có thể tìm được trên cơ sở kết hợp đáp ứng xung của 2 bộ lọc thông thấp lý tưởng theo như hình 2.2 sau đây: Trong đó wc1 và wc2 lần lượt là tần số cắt của 2 bộ lọc thông thấp lý tưởng tương ứng, có thể được tính từ các chỉ tiêu đã cho của bộ lọc thông dải thực tế: wc1 = (ws1 + wp1) / 2 wc2 = (wp2 + ws2) / 2 Đáp ứng xung của lọc thông thấp lý tưởng có thể tìm được từ công thức (1.31). (1.31) Bước 5. Tìm đáp ứng xung h(n) của bộ lọc thực tế (bộ lọc cần thiết kế) Đáp ứng xung của bộ lọc thông dải thực tế được tính bằng cách lấy đáp ứng xung lý tưởng nhân với hàm cửa sổ (công thức 1.33), đây chính là thao tác lấy cửa sổ. h(n) = hd(n).w(n) (1.33) Đến đây chúng ta đã có được bộ lọc cần thiết kế. 2.3. Thuật toán và chương trình Matlab Trong phần này sẽ thực hiện chương trình thiết kế bộ lọc thông dải bằng cách sử dụng cửa sổ Kaiser. Chương trình sẽ nhận các chỉ tiêu yêu cầu của bộ lọc cần thiết kế, sau đó thực hiện các bước thiết kế để tìm được đáp ứng xung h(n). Để khảo sát bộ lọc vừa thiết kế, chương trình cũng sẽ thực hiện tính toán và vẽ đáp ứng biên độ - tần số của bộ lọc theo dB, cũng như vẽ các đáp ứng xung lý tưởng hd(n), hàm cửa sổ w(n) và đáp ứng xung bộ lọc thực tế h(n). Chương trình được viết và chạy trên nền Matlab 6.5, với việc sử dụng một số hàm hỗ trợ có sẵn của Matlab cho xử lý tín hiệu số, và một số hàm viết thêm được tham khảo từ tài liệu [1] (các hàm dưới dạng các file .m). a. Lưu đồ thuật toán: BEGIN Nhập các chỉ tiêu w1s, w1p, w2p, w2s, As, Rp Chỉ tiêu có hợp lệ không? No Tính bậc M và hệ số b của cửa sổ Kaiser theo công thức (1.42) Tìm hàm cửa sổ w(n) (Gọi hàm của Matlab w=kaiser(M,beta) ) Tính hd(n) h(n)=hd(n).w(n) Vẽ hd(n), w(n), h(n) và đáp ứng biên độ (dB) của bộ lọc thiết kế. END Yes b. Chương trình Matlab Chương trình chính: Hàm bpf_kai(ws1,wp1,wp2,ws2,Rp,As) thực hiện thiết kế bộ lọc thông dải dùng cửa sổ Kaiser. function [h] = bpf_kai(ws1,wp1,wp2,ws2,Rp,As); %Design a bandpass filter using Kaiser window technique %------------------------------------------------------------------------ % h = bpf_kai(ws1,wp1,wp2,ws2,Rp,As) % return impulse response of design filter in row vector h % ws1: Lower Stopband Edge % wp1: Lower Passband Edge % wp2: Upper Passband Edge % ws2: Upper Passband Edge % Rp : Passband Ripple in dB % As : Stopband Attenuation in dB %------------------------------------------------------------------------ % Note: 0 < ws1 < wp1 < wp2 < ws2 < 1 % Rp < 1 dB % As > 21 dB %------------------------------------------------------------------------ % Written by: HUYNH VIET THANG % Date: November 24th, 2005 % Check number of input arguments if (nargin == 0), % if there is no input argument assigned temp=input('Lower Stopband Edge: ws1 = '); if (temp =1), error('Lower Stopband Edge ws1: 0< ws1 <1'); end ws1=temp; temp=input('Lower Passband Edge: wp1 = '); if (temp =1), error('Lower Passband Edge wp1: ws1 < wp1 < 1'); end wp1=temp; temp=input('Upper Passband Edge: wp2 = '); if (temp =1), error('Upper Passband Edge wp2: wp1 < wp2 < 1'); end wp2=temp; temp=input('Upper Stopband Edge: ws2 = '); if (temp =1), error('Upper Stopband Edge ws2: wp2 < ws2 < 1'); end ws2=temp; temp=input('Passband Ripple in dB (0 < Rp < 1): Rp = '); if (temp =1), error('Error: 0 < Rp < 1'); end Rp=temp; temp=input('Stopband Attenuation in dB (As > 21): As = '); if (temp < 0), error('Error: As must be larger than 21 dB !!'); end As=temp; elseif (nargin ~= 6), error('Number of input arguments invalid: bpf_kai(ws1,wp1,wp2,ws2,Rp,As)'); end end % Adjust the band-edge specifications ws1=ws1*pi; wp1=wp1*pi; wp2=wp2*pi; ws2=ws2*pi; % Caculate M and beta tr_width = min((wp1-ws1),(ws2-wp2)); % Transition bandwidth M = ceil((As-7.95)/(14.36*tr_width/(2*pi))+1)+1; M = 2*floor(M/2)+1 % Odd filter length n=[0:1:M-1]; if As >= 50 beta = 0.1102*(As-8.7); elseif (As 21) beta = 0.5842*(As-21)^(0.4) + 0.07886*(As-21); else error('As must be greater than 21') end beta w_kai = (kaiser(M,beta))'; % Caculate Kaiser window function wc1 = (ws1+wp1)/2; % cutoff freq of ideal LPF 1 wc2 = (wp2+ws2)/2; % cutoff freq of ideal LPF 2 hd = ideal_lp(wc2,M) - ideal_lp(wc1,M); % Caculate the ideal impulse response h = hd .* w_kai; % Actual impulse response h(n) [db,mag,pha,grd,w] = freqz_m(h,[1]); % Caculate Magnitude Response delta_w = 2*pi/1000; % disp(['Actual Passband Ripple:']) % Rp = -min(db(wp1/delta_w+1:1:wp2/delta_w)) % Actua; Passband Ripple disp(['Actual Stopband Attenuation:']) % As = -round(max(db(ws2/delta_w+1:1:501))) % Min Stopband Attenuation %================================================== % plot Ideal impulse response hd(n) subplot(1,1,1); subplot(2,2,1); stem(n,hd,'.'); title('Ideal Impulse Response hd(n)') axis([0 M-1 -0.4 0.4]); xlabel('n'); ylabel('hd(n)') %================================================== % plot Kaiser window function w_kai(n) subplot(2,2,2); stem(n,w_kai,'.');title('Kaiser Window w(n)') axis([0 M-1 -0.1 1.1]); xlabel('n'); ylabel('w(n)') %================================================== % plot design impulse response h(n) subplot(2,2,3); stem(n,h,'.');title('Actual Impulse Response h(n)') axis([0 M-1 -0.4 0.4]); xlabel('n'); ylabel('h(n)') %================================================== % plot magnitude response in dB subplot(2,2,4); plot(w/pi,db); title('Magnitude Response in dB');grid; xlabel('frequency in pi units'); ylabel('Decibels') axis([0 1 -150 10]); %================================================== % plot dash lines set(gca,'XTickMode','manual','XTick',[0,ws1/pi,wp1/pi,wp2/pi,ws2/pi,1]) set(gca,'YTickMode','manual','YTick',[-As,0]) Các hàm chính được sử dụng trong chương trình: ideal_lp(wc,M): Tính đáp ứng xung của bộ lọc thông thấp lý tưởng bậc M kaiser(M,beta): Tính cửa sổ Kaiser chiều dài M, với tham số beta. freqz_m(h,[1]): Tính đáp ứng tần số của bộ lọc FIR có đáp ứng xung h. Hai hàm ideal_lp.m và freqz_m.m được tham khảo từ tài liệu [1]. 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 %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: a=h) % a=demonitor 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); function hd = ideal_lp(wc,M); % Ideal LowPass filter computation % -------------------------------- % [hd] = ideal_lp(wc,M) % hd = ideal impulse response between 0 to M-1 % wc = cutoff frequency in radians % M = length of the ideal filter % alpha = (M-1)/2; n = [0:1:(M-1)]; m = n - alpha + eps; hd = sin(wc*m) ./ (pi*m); 2.4. Kết quả chạy chương trình thiết kế Các chỉ tiêu thiết kế được cho: ws1 = 0.3p, wp1 = 0.4p wp2 = 0.6p, ws2 = 0.7p Suy hao ở dải chắn : As = 74 dB Độ gợn sóng ở dải thông : Rp = 0.5 dB Cách gọi hàm : h=bpf_kai(0.3, 0.4, 0.6, 0.7, 0.5, 74); Kết quả chạy chương trình thiết kế: M = 95 Beta = 7.1961 Độ gợn sóng trong dải thông thực tế: Rp = 0.0029 dB Suy hao trong dải chắn thực