Vẽ phổ tín hiệu bằng MATLAB

Successfully reported this slideshow.

Your SlideShare is downloading. ×

Working at Trường Đại Học Bách Khoa - ĐH Đà Nẵng

More Related Content

1. 1. BỘ GIÁO DỤC VÀ ĐÀO TẠO ĐẠI HỌC CÔNG NGHỆ TP.HCM THỰC HÀNH XỬ LÝ TÍN HIỆU SỐ Biên soạn: ThS. Phạm Hùng Kim Khánh www.hutech.edu.vn
2. 2. THỰC HÀNH XỬ LÝ TÍN HIỆU SỐ Ấn bản 2014
3. 3. MỤC LỤC I MỤC LỤC MỤC LỤC ................................................................................................................... 1 HƯỚNG DẪN .............................................................................................................. 3 BÀI 1: PHẦN MỀM MATLAB.......................................................................................... 5 1.1 KHỞI ĐỘNG MATLAB.............................................................................................. 5 1.2 CÁC VẤN ĐỀ CƠ BẢN .............................................................................................. 8 1.2.1 Các phép toán và toán tử ................................................................................. 8 1.2.2 Khai báo biến ................................................................................................. 9 1.2.3 Các lệnh thường dùng...................................................................................... 9 1.3 LẬP TRÌNH TRONG MATLAB................................................................................... 10 1.3.1 Các phát biểu điều kiện if, else, elseif .............................................................. 10 1.3.2 Switch ......................................................................................................... 10 1.3.3 While........................................................................................................... 10 1.3.4 For .............................................................................................................. 11 1.3.5 Break: ......................................................................................................... 11 1.4 MA TRẬN ............................................................................................................ 11 1.4.1 Các thao tác trên ma trận .............................................................................. 11 1.4.2 Vector.......................................................................................................... 15 1.4.3 Đa thức........................................................................................................ 15 1.5 ĐỒ HOẠ ............................................................................................................. 16 1.5.1 Các lệnh vẽ .................................................................................................. 16 1.5.2 Tạo hình vẽ .................................................................................................. 16 1.5.3 Kiểu đường vẽ .............................................................................................. 16 1.5.4 Vẽ với hai trục y............................................................................................ 17 1.5.5 Vẽ đường cong 3-D........................................................................................ 18 1.5.6 Vẽ nhiều trục toạ độ ...................................................................................... 18 1.5.7 Đặt các thông số cho trục............................................................................... 19 1.5.8 Đồ hoạ đặc biệt............................................................................................. 23 1.5.9 Đồ hoạ 3D.................................................................................................... 25 1.5.10 Thực hành vẽ đồ thị..................................................................................... 26 1.6 CÁC FILE VÀ HÀM ................................................................................................ 28 1.6.1 Script file (file kịch bản) ................................................................................. 28 1.6.2 File hàm....................................................................................................... 28 1.6.3 Các hàm toán học cơ bản ............................................................................... 29 1.6.4 Các phép toán trên hàm toán học.................................................................... 30
4. 4. TRANG 2 > THỰC HÀNH XỬ LÝ TÍN HIỆU SỐ 1.6.5 Thực hành trên script và function .....................................................................31 BÀI 2: TÍN HIỆU RỜI RẠC THEO THỜI GIAN ............................................................. 38 2.1 CÁC TÍN HIỆU SƠ CẤP ...........................................................................................38 2.2 CÁC PHÉP TOÁN ...................................................................................................39 2.3 KIỂM TRA TÍNH CHẤT TUYẾN TÍNH VÀ BẤT BIẾN .......................................................39 2.4 HỆ LTI.................................................................................................................41 BÀI 3: BIẾN ĐỔI Z .................................................................................................... 44 3.1 CÁC ĐIỂM CỰC VÀ ĐIỂM KHÔNG .............................................................................44 3.2 PHÂN TÍCH DÙNG PHƯƠNG PHÁP THẶNG DƯ ............................................................45 3.3 BIẾN ĐỔI Z VÀ Z NGƯỢC .......................................................................................46 BÀI 4: BIẾN ĐỔI FOURIER RỜI RẠC.......................................................................... 49 4.1 TÍNH DTFT...........................................................................................................49 4.2 FFT VÀ CÁC TÍNH CHẤT .........................................................................................50 BÀI 5: BỘ LỌC SỐ FIR ............................................................................................... 53 5.1 CÁC LOẠI BỘ LỌC .................................................................................................53 5.2 PHƯƠNG PHÁP CỬA SỔ..........................................................................................55 5.3 PHƯƠNG PHÁP LẤY MẪU TẦN SỐ.............................................................................58 5.4 PHƯƠNG PHÁP LẶP ...............................................................................................60 BÀI 6: BỘ LỌC SỐ IIR ............................................................................................... 65 6.1 THIẾT KẾ BỘ LỌC TƯƠNG TỰ..................................................................................65 6.2 THIẾT KẾ BỘ LỌC SỐ.............................................................................................69 TÀI LIỆU THAM KHẢO .................................................................................................73
5. 5. HƯỚNG DẪN II I HƯỚNG DẪN MÔ TẢ MÔN HỌC Thực hành xử lý tín hiệu số là môn học hỗ trợ cho môn Xử lý tín hiệu số của chuyên ngành Kỹ thuật Điện tử Truyền thông. Môn học này dựa trên MATLAB để kiểm chứng các lý thuyết đã học trong môn Xử lý tín hiệu số. NỘI DUNG MÔN HỌC  Bài 1. Phần mềm MATLAB: cơ bản về MATLAB, cách lập trình cũng như cách xử lý ma trận, vẽ đồ thị trong MATLAB.  Bài 2: Tín hiệu rời rạc theo thời gian: cách biểu diễn tín hiệu và hệ thống rời rạc theo thời gian, các tính chất và đáp ứng xung của hệ LTI.  Bài 3: Biến đổi z và z ngược: chuyển tín hiệu từ miền thời gian sang miền z, các tính chất của biến đổi z và chuyển tín hiệu hữu tỷ trên miền z sang miền thời gian.  Bài 4: Biến đổi Fourier rời rạc: chuyển tín hiệu rời rạc trên miền thời gian sang miền tần số, dùng các thuật toán FFT để xác định biến đổi Fourier rời rạc.  Bài 5: Bộ lọc số FIR: thiết kế bộ lọc FIR theo yêu cầu cho trước.  Bài 6: Bộ lọc số IIR: thiết kế bộ lọc IIR theo yêu cầu cho trước. KIẾN THỨC TIỀN ĐỀ Môn học Thực hành Xử lý tín hiệu số đòi hỏi sinh viên có nền tảng về Xử lý tín hiệu số. YÊU CẦU MÔN HỌC Người học phải dự học đầy đủ các buổi lên lớp và làm bài tập đầy đủ. CÁCH TIẾP NHẬN NỘI DUNG MÔN HỌC
6. 6. TRANG 4 > THỰC HÀNH XỬ LÝ TÍN HIỆU SỐ Để học tốt môn này, người học cần thực hành theo hướng dẫn, làm các bài tập; đọc trước bài mới và tìm thêm các thông tin liên quan đến bài học. PHƯƠNG PHÁP ĐÁNH GIÁ MÔN HỌC Môn học được đánh giá gồm:  Điểm quá trình: 30%. Hình thức và nội dung do giảng viên quyết định, phù hợp với quy chế đào tạo và tình hình thực tế tại nơi tổ chức học tập.  Điểm thi: 70%. Hình thức bài thi trên máy tính trong 60 phút.
7. 7. BÀI 1: PHẦN MỀM MATLAB 5 BÀI 1: PHẦN MỀM MATLAB Sau khi học xong bài này, người học có thể:  Sử dụng phần mềm MATLAB.  Thực hiện tạo các script file hay function và lưu trữ trên MATLAB.  Biết được các công cụ cơ bản trên MATLAB. 1.1 KHỞI ĐỘNG MATLAB MATLAB (Matrix laboratory) là phần mềm dùng để giải các bài toán kỹ thuật, đặc biệt là các bài toán liên quan đến ma trận. MATLAB cung cấp các toolboxes, tức các hàm mở rộng môi trường MATLAB để giải quyết các vấn đề đặc biệt như xử lý tín hiệu số, hệ thống điều khiển, mạng neuron, fuzzy logic, mô phỏng v.v. Cửa sổ biểu tượng của chương trình MATLAB: Hình 1.1 - Cửa sổ khởi động của MATLAB
8. 8. 6 BÀI 1: PHẦN MỀM MATLAB Cửa sổ làm việc của MATLAB: Hình 1.2 – Cửa sổ làm việc của MATLAB  Cửa sổ lệnh (Command window): Là cửa sổ giao tiếp chính của MATLAB bởi đây là nơi nhập giá trị các biến, hiển thị giá trị, tính toán giá trị của biểu thức, thực thi các hàm có sẵn trong thư viện hoặc các hàm do người dùng lập trình ra trong M-files. Các lệnh được nhập sau dấu nhắc ‘>>‘ và thực thi lệnh bằng phím Enter. Để mở chương trình soạn thảo trong MATLAB, gõ lệnh: >>edit Hình 1.3 – Cửa sổ edit để soạn script file hay function Thư mục hiện hành Cửa sổ lệnh Danh sách các file hiện có trong thư mục hiện hành Workspace: các biến trong bộ nhớ Các lệnh đã thực hiện Nút Start: chứa các toolbox
9. 9. BÀI 1: PHẦN MỀM MATLAB 7 Sau đó nhập vào đoạn chương trình sau: % Chuong trinh trong M-file x= 0:pi/6:2*pi; y=sin(x); plot(x,y); Lưu chương trình với tên file plot_sin.m bằng cách nhấn Ctrl+S hay nhấn vào biểu tượng Save . Giải thích đoạn chương trình trên: Dòng 1 là dòng chú thích, các chuỗi ở phía sau dấu % sẽ không được dịch. Dòng 2 định nghĩa vector x trong khoảng từ 0 đến 2 và mỗi giá trị cách nhau một khoảng /6. Dòng 3 gán biến y = sin(x) và dòng 4 vẽ đồ thị trong đó x là trục hoành và y là trục tung. Hình 1.4 – Lưu file trong cửa số Edit Thực thi chương trình trên trong Command window bằng dòng lệnh sau: >>plot_sin  Cửa sổ Command History: Các dòng đã nhập trong Command window (các dòng này có thể là dòng nhập biến, có thể là dòng lệnh) được giữ lại trong cửa sổ Command History và cửa sổ này cho phép ta sử dụng lại những lệnh đó bằng cách nhấp đúp chuột lên các lệnh hay biến đó  Cửa sổ Workspace:
10. 10. 8 BÀI 1: PHẦN MỀM MATLAB Là cửa sổ thể hiện tên các biến bạn sử dụng cùng với kích thước vùng nhớ (số bytes), kiểu dữ liệu (lớp), các biến được giải phóng sau mỗi lần tắt chương trình. Cửa sổ Workspace cho phép thay đổi giá trị của biến bằng cách nhấn phím chuột phải lên các biến và chọn Edit.  Nút Start:  Cửa sổ Edit: Là cửa sổ dùng để soạn thảo chương trình ứng dụng, được khởi động bằng lệnh edit trong Command Window. Chương trình soạn thảo trong cửa sổ edit có 2 dạng: + Dạng Script file : tập hợp các câu lệnh viết dưới dạng liệt kê, không có biến dữ liệu vào và biến ra. + Dạng function: có biến dữ liệu vào và biến ra. 1.2 CÁC VẤN ĐỀ CƠ BẢN 1.2.1 Các phép toán và toán tử Các phép toán: + , - , * , / , (chia trái) , ^ (mũ) , ‘ (chuyển vị hay số phức liên hiệp). Các toán tử quan hệ : < , <= , > , >= , == , ~= Các toán tử logic : & , | (or) , ~ (not) Các hằng: pi 3.14159265 Hình 1.5 – Nút Start
11. 11. BÀI 1: PHẦN MỀM MATLAB 9 i, j số ảo eps sai số 2-52 realmin số thực nhỏ nhất 2-1022 realmax số thực lớn nhất 21023 inf vô cùng lớn NaN Not a number Chú ý : + Các lệnh kết thúc bằng dấu chấm phẩy, MATLAB sẽ không thể hiện kết quả trên màn hình. + Các chú thích được đặt phía sau dấu %. + Trong quá trình nhập nếu các phần tử trên một hàng dài quá ta có thể xuống dòng bằng toán tử ba chấm (. . .) 1.2.2 Khai báo biến - Phân biệt chữ hoa và chữ thường. - Không cần phải khai báo kiểu biến. - Tên biến phải bắt đầu bằng ký tự và không được có khoảng trắng. - Không đặt tên trùng với các tên đặc biệt của MATLAB. - Để khai báo biến toàn cục (sử dụng được trong tất cả chương trình con), phải dùng thêm từ khoá global phía trước. 1.2.3 Các lệnh thường dùng >>help tên_hàm % tham khảo help của hàm >>lookfor ‘chuỗi’ %Tìm kiếm chuỗi Ví dụ: >>lookfor ‘filter’ % Tìm các hàm có liên quan đến mạch lọc >>clc % Xoá màn hình >>clear tên_biến % Xoá biến >>clear all %Xoá tất cả các biến >> clf %Xoá figure >>save % Lưu các biến hiện có trong bộ nhớ >>load % Lấy nội dung các biến đã lưu >>who % liệt kê các biến trong bộ nhớ >>whos % liệt kê chi tiết các biến trong bộ nhớ >>which % Xác định vị trí của hàm hay file
12. 12. 10 BÀI 1: PHẦN MỀM MATLAB Ví dụ: >>which plot %Xác định vị trí của hàm plot >>what % Liệt kê các file có trong một thư mục 1.3 LẬP TRÌNH TRONG MATLAB 1.3.1 Các phát biểu điều kiện if, else, elseif Cú pháp của if: if end Nếu cho kết quả đúng thì phần lệnh trong thân của if được thực hiện. Các phát biểu else và elseif cũng tương tự. 1.3.2 Switch Cú pháp của switch như sau: switch case n1 case n2 . . . . . . . . . . . . . . . case nn Otherwise end 1.3.3 While Vòng lặp while dùng khi không biết trước số lần lặp. Cú pháp của nó như sau: while end
13. 13. BÀI 1: PHẦN MỀM MATLAB 11 1.3.4 For Vòng lặp for dùng khi biết trước số lần lặp. Cú pháp như sau: for =:: 1.3.5 Break: Phát biểu break để kết thúc vòng lặp for hay while mà không quan tâm đến điều kiện kết thúc vòng lặp đã thoả mãn hay chưa. 1.4 MA TRẬN 1.4.1 Các thao tác trên ma trận 1.4.1.1 Nhập ma trận Ma trận là một mảng có m hàng và n cột. Trường hợp ma trận chỉ có một phần tử (ma trận 1x1) ta có một số. Ma trận chỉ có một cột hay một hàng được gọi là một vector. Ta có thể nhập ma trận vào MATLAB bằng nhiều cách: • Nhập một danh sách các phần tử từ bàn phím. • Nạp ma trận từ file. • Tạo ma trận nhờ các hàm có sẵn trong MATLAB. • Tạo ma trận nhờ hàm tự tạo. Khi nhập ma trận từ bàn phím ta phải tuân theo các quy định sau: • Ngăn cách các phần tử của ma trận bằng dấu “,” hay khoảng trắng. • Dùng dấu “;” để kết thúc một hàng. • Bao các phần tử của ma trận bằng cặp dấu ngoặc vuông [ ]. 1.4.1.2 Chỉ số Phần tử ở hàng i cột j của ma trận có ký hiệu là A(i,j). Tuy nhiên, ta cũng có thể tham chiếu tới phần tử của mảng nhờ một chỉ số, ví dụ A(k). Trong trường hợp này, ma trận được xem là một cột dài tạo từ các cột của ma trận ban đầu. Như vậy viết A(8) có nghĩa là tham chiếu phần tử A(4, 2) (nếu ma trận có 4 hàng). Lưu ý rằng các chỉ số của ma trận thường bắt đầu từ 1. 1.4.1.3 Toán tử “:”
14. 14. 12 BÀI 1: PHẦN MỀM MATLAB Toán tử “:” là một toán tử quan trọng của MATLAB. Nó xuất hiện ở nhiều dạng khác nhau. Biểu thức 1:10 là một vector hàng chứa 10 số nguyên từ 1 đến 10. >>1:10 >>100:‐7:50 %tạo dãy số từ 100 đến 51, cách đều nhau 7 >>0: pi/4: pi %tạo một dãy số từ 0 đến , cách đều nhau /4 Các biểu thức chỉ số có thể tham chiếu tới một phần của ma trận. A(1:k,j) xác định k phần tử đầu tiên của cột j. Ngoài ra toán tử “:” tham chiếu tới tất cả các phần tử của một hàng hay một cột. >>A(:,3) >>A(3, :) >>B = A(:, [1 3 2 4]) %tạo ma trận B từ ma trận A bằng cách đổi thứ tự các cột từ [1 2 3 4] thành [1 3 2 4] 1.4.1.4 Tạo ma trận bằng hàm có sẵn MATLAB cung cấp một số hàm để tạo các ma trận cơ bản: - zeros tạo ra ma trận mà các phần tử đều là 0. >>z = zeros(2, 4) - ones tạo ra ma trận mà các phần tử đều là 1. >>x = ones(2, 3) >>y = 5*ones(2, 2) - rand tạo ra ma trận mà các phần tử ngẫu nhiên phân bố đều. >>d = rand(4, 4) - randn tạo ra ma trận mà các phần tử ngẫu nhiên phân bố chuẩn. >>e = randn(3, 3) - magic(n) tạo ra ma trận cấp n gồm các số nguyên từ 1 đến n2 với tổng các hàng bằng tổng các cột và bằng tổng các đường chéo (n  3). - pascal(n) tạo ra tam giác Pascal. >>pascal(4) - eye(n) tạo ma trận đơn vị >>eye(3) - eye(m,n) tạo ma trận đơn vị mở rộng
15. 15. BÀI 1: PHẦN MỀM MATLAB 13 >>eye(3,4) 1.4.1.5 Nối ma trận Ta có thể nối các ma trận có sẵn thành một ma trận mới. Ví dụ: >>a = ones(3, 3); >>b = 5*ones(3, 3); >>c = [a+2; b] 1.4.1.6 Xoá hàng và cột Ta có thể xoá hàng và cột từ ma trận bằng dùng dấu []. >>b(:, 2) = [] ; %xoá cột thứ 2 >>b(1:2:5) = []; % xoá các phần tử bắt đầu từ 1 đến 5 và cách 2 (1,3,5) rồi sắp xếp lại ma trận. 1.4.1.7 Các lệnh xử lý ma trận Cộng : X = A + B Trừ : X = A - B Nhân : X = A*B X = A.*B nhân các phần tử tương ứng với nhau, yêu cầu 2 ma trận A và B phải có cùng kích thước. Chia : X = A/B lúc đó X = A * inv(B) X = AB lúc đó X = inv(A) * B X=A./B chia các phần tử tương ứng với nhau, 2 ma trận A và B có cùng kích thước. Luỹ thừa : X = A^2 X = A.^2 Nghịch đảo: X = inv(A) Định thức: d = det(A) Để tạo ma trận trong MATLAB ta chỉ cần liệt các phần tử của ma trận trong cặp dấu ngoặc vuông ([…]). Các phần tử trên cùng hàng được phân biệt bởi dấu phẩy (,) hoặc khoảng trắng (space). Các hàng của ma trận, phân cách nhau bởi dấu chấm phẩy (;). Ví dụ, nhập ma trận A có 4 hàng, 4 cột như sau: >> A=[16 3 2 13; 5 10 11 8; 9 6 7 12; 4 15 14 1]
16. 16. 14 BÀI 1: PHẦN MỀM MATLAB >> size(A) Để truy xuất đến từng phần tử của ma trận ta dùng chỉ số phần tử tương ứng. Ví dụ, phần tử ở hàng thứ 2, cột thứ 3 của A là A(2,3). >> A(2,3) Bài 1.1. Cho ma trận A = [2 4 1; 6 7 2; 3 5 9], sinh viên dùng các lệnh cần thiết để: a. Lấy dòng đầu tiên của ma trận A. b. Tạo ma trận B bằng 2 dòng cuối cùng của A. c. Tính tổng các phần tử trên các cột của A. (gợi ý: tính tổng các phần tử trên cột 1: sum(A(:,1))). d. Tính tổng các phần tử trên các dòng của A. Bài 1.2. Cho ma trận A = [2 7 9 7; 3 1 5 6; 8 1 2 5], sinh viên giải thích kết quả của các lệnh sau: a. A' b. A(:,[1 4]) c. A([2 3],[3 1]) d. reshape(A,2,6) e. A(:) f. [A;A(end,:)] g. A(1:3,:) h. [A ; A(1:2,:)] i. sum(A) j. sum(A') k. [ [ A ; sum(A) ] [ sum(A,2) ; sum(A(:)) ] ] Bài 1.3. Giải hệ phương trình sau: 2x1 + 4x2 + 6x3 – 2x4 = 0 x1 + 2x2 + x3 + 2x4 = 1 2x2 + 4x3 + 2x4 = 2 3x1 - x2 + 10x4 = 10. Gợi ý: A.X = B  X = inv(A).B
17. 17. BÀI 1: PHẦN MỀM MATLAB 15 Bài 1.4. Chứng tỏ rằng (A+B)C=AC+BC, với: 1.4.2 Vector Vector thực chất là ma trận có kích thước n x 1 hay 1 x n, nên ta có thể tạo ra vector như cách tạo ra ma trận. Ngoài ra, có thể dùng một số cách sau: >>x=0:0.1:1 >>y=linspace(1, 10, 20) % vector 20 phan tu cach deu nhau tu 1 den 10 >>z=rand(10,1) ; tạo 10 số ngẫu nhiên phân bố đều Bài 1.5. Cho vector x = [3 1 5 7 9 2 6], giải thích kết quả của các lệnh sau: a. x(3) b. x(1:7) c. x(1:end) d. x(1:end-1) e. x(6:-2:1) f. x([1 6 2 1 1]) g. sum(x) Bài 1.6. Tạo một vector x có 100 phần tử, sao cho: x(n) = (-1)n+1 /(2n+1) với n = 0 – 99. 1.4.3 Đa thức Các đa thức trong MATLAB được mô tả bằng các vector hàng với các phần tử của vector chính là các hệ số của đa thức, xếp theo thứ tự số mũ giảm dần. Ví dụ, đa thức m = s4 -s3 +4s2 -5s-1 được biểu diễn là: >>m=[1 -1 4 -5 -1] Để xác định giá trị của đa thức, ta dùng hàm polyval. Ví dụ, xác định giá trị của đa thức tại điểm s=2:
18. 18. 16 BÀI 1: PHẦN MỀM MATLAB >>polyval(m,2) Để xác định nghiệm của đa thức, ta dùng hàm roots. Ví dụ: >>roots(m) Bài 1.7. Cho phương trình ax2 +bx+c=0, giải phương trình dùng hàm roots. >>m=[a b c]; >>x=roots(m) Thay đổi các giá trị khác nhau của a, b và c tương ứng trong 2 cách giải trên. Bài 1.8. Giải phương trình x3 - 2x2 +4x+5=0. Kiểm chứng kết quả thu được bằng hàm polyval. Sinh viên có nhận xét gì về kết quả kiểm chứng. Bài 1.9. Lặp lại bài 1.8 cho phương trình x7 -2=0. 1.5 ĐỒ HOẠ 1.5.1 Các lệnh vẽ MATLAB cung cấp một loạt hàm để vẽ biểu diễn các vector cũng như giải thích và in các đường cong này. plot: đồ họa 2-D với số liệu 2 trục vô hướng và tuyến tính plot3: đồ họa 3-D với số liệu 2 trục vô hướng và tuyến tính loglog: đồ hoạ với các trục x, y ở dạng logarit semilogx: đồ hoạ với trục x logarit và trục y tuyến tính semilogy: đồ hoạ với trục y logarit và trục x tuyến tính 1.5.2 Tạo hình vẽ Hàm plot có các dạng khác nhau phụ thuộc vào các đối số đưa vào. Ví dụ nếu y là một vector thì plot(y) tạo ra một đường quan hệ giữa các giá trị của y và chỉ số của nó. Nếu ta có 2 vector x và y thì plot(x,y) tạo ra đồ thị quan hệ giữa x và y. >>t = [0:pi/100:2*pi]; >>y = sin(t); >>plot(t,y); % Vẽ hình sin từ 0->2 >>grid on 1.5.3 Kiểu đường vẽ Ta có thể dùng các kiểu đường vẽ khác nhau khi vẽ hình. >>t = [0:pi/100:2*pi];
19. 19. BÀI 1: PHẦN MỀM MATLAB 17 >>y = sin(t); >>plot(t,y,’.‘) % vẽ bằng đường chấm chấm Để xác định màu và kích thước đường vẽ, ta dùng các tham số sau: LineWidth : độ rộng đường thẳng, tính bằng số điểm MarkerEdgeColor : màu của các cạnh của khối đánh dấu MarkerFaceColor : màu của khối đánh dấu MarkerSize : kích thước của khối đánh dấu Màu được xác định bằng các tham số: r: red m magenta g: green y: yellow b: blue k: black c: cyan w: white Các dạng đường thẳng xác định bằng: - đường liền -- đường đứt nét : đường chấm chấm -. đường chấm gạch Các dạng điểm đánh dấu xác định bằng: + dấu cộng . điểm o vòng tròn x chữ thập * dấu sao s hình vuông d hạt kim cương v tam giác hướng xuống ^ tam giác hướng lên < tam giác sang trái > tam giác sang phải h lục giác p ngũ giác >>x = -pi : pi/10 : pi; >>y = tan(sin(x)) - sin(tan(x)); >>plot(x,y,ʹ‐‐rs’,ʹLineWidthʹ,2,ʹMarkerEdgeColorʹ,ʹkʹ,... ʹMarkerFaceColorʹ,ʹgʹ,ʹMarkerSizeʹ,10) Để vẽ hai hàm trên cùng một đồ thị, ta dùng lệnh: >>hold on 1.5.4 Vẽ với hai trục y Hàm plotyy cho phép tạo một đồ thị có hai trục y. Ta cũng có thể dùng plotyy để cho giá trị trên hai trục y có kiểu khác nhau nhằm tiện so sánh. >>t = 0:900; >>A = 1000;
20. 20. 18 BÀI 1: PHẦN MỀM MATLAB >>b = 0.005; >>a = 0.005; >>z2 = sin(b*t); >>z1 = A*exp(-a*t); >>[haxes, hline1, hline2] = plotyy(t,z1,t,z2,ʹsemilogyʹ,ʹplotʹ); 1.5.5 Vẽ đường cong 3-D Nếu x, y, z là 3 vector có cùng độ dài thì plot3 sẽ vẽ đường cong 3D. >>t = 0:pi/50:10*pi; >>plot3(sin(t),cos(t),t) >>axis square; >>grid on 1.5.6 Vẽ nhiều trục toạ độ Dùng hàm subplot để vẽ nhiều trục toạ độ. >>subplot(2,3,5) %2,3: xác định có 2 hàng, 3 cột % 5: chọn trục thứ 5 (đếm từ trái sang phải, trên xuống dưới) >>x=linspace(0,2*pi); >>y1=sin(x); >>y2=cos(x) >>y3=2*exp(-x).*sin(x); >>x1=linspace(-2*pi,2*pi); >>y4=sinc(x1); >>subplot(221); >>plot(x,y1); >>title(‘Ham y = sinx’); >>subplot(222); >>plot(x,y2); >>title(‘Ham y = cosx’); >>subplot(223); >>plot(x,y3);
21. 21. BÀI 1: PHẦN MỀM MATLAB 19 >>title(‘Ham y = 2e^{-x}sinx’); >>subplot(224); >>plot(x1,y4); >>title(‘Ham y = $${sin pi x over pi x}$$','interpreter','latex'); 1.5.7 Đặt các thông số cho trục Khi ta tạo một hình vẽ, MATLAB tự động chọn các giới hạn trên trục toạ độ và khoảng cách đánh dấu dựa trên dữ liệu dùng để vẽ. Tuy nhiên ta có thể mô tả lại phạm vi giá trị trên trục và khoảng cách đánh dấu theo ý riêng. Ta có thể dùng các hàm sau: axis đặt lại các giá trị trên trục toạ độ. axes tạo một trục toạ độ mới với các đặc tính được mô tả. get và set cho phép xác định và đặt các thuộc tính của trục toạ độ đang có. gca trở về trục toạ độ cũ. 1.5.7.1 Giới hạn của trục và chia vạch trên trục MATLAB chọn các giới hạn trên trục toạ độ và khoảng cách đánh dấu dựa trên số liệu dùng để vẽ. Dùng lệnh axis có thể đặt lại giới hạn này. Cú pháp của lệnh: axis[xmin , xmax , ymin , ymax] >>x = 0:0.025:pi/2; >>plot(x,tan(x),ʹ‐roʹ) >>axis([0 pi/2 0 5]) MATLAB chia vạch trên trục dựa trên phạm vi dữ liệu và chia đều. Ta có thể mô tả cách chia nhờ thông số xtick và ytick bằng một vector tăng dần. >>x = -pi:.1:pi; >>y = sin(x); >>plot(x,y) >>set(gca,ʹxtickʹ,‐pi:pi/2:p); >>set(gca,ʹxticklabelʹ,{ʹ‐piʹ,ʹ‐pi/2ʹ,ʹ0ʹ,ʹpi/2ʹ,ʹpiʹ}) 1.5.7.2 Ghi nhãn lên các trục toạ độ MATLAB cung cấp các lệnh ghi nhãn lên đồ hoạ gồm:
22. 22. 20 BÀI 1: PHẦN MỀM MATLAB title thêm nhãn vào đồ hoạ. xlabel thêm nhãn vào trục x. ylabel thêm nhãn vào trục y. zlabel thêm nhãn vào trục z. legend thêm chú giải vào đồ thị. text hiển thị chuỗi văn bản ở vị trí nhất định. gtext đặt văn bản lên đồ hoạ nhờ chuột. >>x = -pi:.1:pi; >>y = sin(x); >>plot(x,y) >>xlabel(ʹt = 0 to 2piʹ,ʹFontsizeʹ,16) >>ylabel(ʹsin(t)ʹ,ʹFontsizeʹ,16) >>title(ʹit{Gia tri cua sin tu zero đến 2 pi}ʹ,ʹFontsizeʹ,16) Ta có thể thêm văn bản vào bất kỳ chỗ nào trên hình vẽ nhờ hàm text. >>text(3*pi/4,sin(3*pi/4),ʹleftarrowsin(t)=0.707ʹ,ʹFontSizeʹ,1 2) Ta có thể sử dụng đối tượng văn bản để ghi chú các trục ở vị trí bất kỳ. MATLAB định vị văn bản theo đơn vị dữ liệu trên trục. Ví dụ để vẽ hàm y = Aet với A = 0.25, t = 0 đến 900 và = 0.005: >>t = 0:900; >>plot(t,0.25*exp(‐0.005*t)) Để thêm ghi chú tại điểm t = 300, ta viết: >>text(300,.25*exp(‐.005*300),... ’bulletleftarrowfontname{times}0.25{ite}^(‐0.005{itt}} at,... {itt}=300’,ʹFontSize’,14) Tham số HorizontalAlignment và VerticalAlignment định vị văn bản so với các toạ độ x, y, z đã cho.
23. 23. BÀI 1: PHẦN MỀM MATLAB 21 Để thêm các ký tự đặc biệt, ta dùng format dạng Tex:
24. 24. 22 BÀI 1: PHẦN MỀM MATLAB bf — Bold font it — Italic font sl — Oblique font (rarely available) rm — Normal font fontname{fontname} fontsize{fontsize} color(colorSpec) Các chỉ số trên và dưới thực hiện bằng ^ và _. >>title('ite^{iomega_0tau} = cos(omega_0tau) + i sin(omega_0tau)') Kết quả: 𝑒 𝑖𝜔0 𝜏 = cos( 𝜔0 𝜏) + 𝑖𝑠𝑖𝑛(𝜔0 𝜏) Để thêm các công thức toán học, ta dùng dạng LaTeX. Một số ví dụ: >>text('units','inch', 'position',[.2 5], ... 'fontsize',14, 'interpreter', 'latex', 'string' ,... ['$$hbox {magic(3) is } left( {matrix{ 8 & 1 & 6 cr'... '3 & 5 & 7 cr 4 & 9 & 2 } } right)$$']); >>text('units','inch', 'position',[.2 4], ... 'fontsize',14, 'interpreter','latex', 'string',... ['$$left[ {matrix{cos(phi) & -sin(phi) cr'... 'sin(phi) & cos(phi) cr}} right]'... 'left[ matrix{x cr y} right]$$']); >>text('units','inch', 'position',[.2 3], ... 'fontsize',14, 'interpreter','latex', 'string',... ['$$L{f(t)} equiv F(s) = int_0^infty!!{e^{-st}'... 'f(t)dt}$$']); >>text('units','inch', 'position',[.2 2], ... 'fontsize',14, 'interpreter','latex', 'string',... '$$e = sum_{k=0}^infty {1 over {k!} } $$'); >>text('units','inch', 'position',[.2 1], ...
25. 25. BÀI 1: PHẦN MỀM MATLAB 23 'fontsize',14, 'interpreter','latex', 'string',... ['$$m ddot y = -m g + C_D cdot {1 over 2}'... 'rho {dot y}^2 cdot A$$']); >>text('units','inch', 'position',[.2 0], ... 'fontsize',14, 'interpreter','latex', 'string',... '$$int_{0}^{infty} x^2 e^{-x^2} dx = frac{sqrt{pi}}{4}$$'); Các khai báo cụ thể tham khảo tại: http://www.latex-project.org/ 1.5.8 Đồ hoạ đặc biệt 1.5.8.1 Khối và vùng Đồ hoạ khối và vùng biểu diễn số liệu là vector hay ma trận. MATLAB cung cấp các hàm đồ hoạ khối và vùng: bar hiển thị các cột của ma trận m*n như là m nhóm, mỗi nhóm có n bar. barh hiển thị các cột của ma trận m*n như là m nhóm, mỗi nhóm có n bar nằm ngang. bar3 hiển thị các cột của ma trận m*n như là m nhóm, mỗi nhóm có n bar dạng 3D. bar3h hiển thị các cột của ma trận m*n như là m nhóm, mỗi nhóm có n bar dạng 3D nằm ngang. Mặc định, mỗi phần tử của ma trận được biểu diễn bằng một bar. >>y = [5 2 1 6 7 3 8 6 3 5 5 5 1 5 8]; >>bar(y) Sau đó nhập lệnh bar3(y) ta có đồ thị 3D. 1.5.8.2 Xếp chồng đồ thị Ta có thể xếp chồng số liệu trên đồ thị thanh bằng cách tạo ra một trục khác trên cùng một vị trí và như vậy ta có một trục y độc lập với bộ số liệu khác.
26. 26. 24 BÀI 1: PHẦN MỀM MATLAB >>TCE = [515 420 370 250 135 120 60 20]; >>nhdo = [29 23 27 25 20 23 23 27]; >>ngay = 0:5:35; >>bar(ngay,nhdo) >>xlabel(ʹNgayʹ) >>ylabel(ʹNhiet do (^{o}C)ʹ) Để xếp chồng một số liệu lên một đồ thị thanh ở trên, có trục thứ 2 ở cùng vị trí như trục thứ nhất ta viết : >>h2 = gca; và tạo trục thứ 2 ở vị trí trục thứ nhất trước nhất vẽ bộ số liệu thứ 2 >>h2 = axes(ʹPositionʹ,get(h2,ʹPositionʹ)); >>plot(days,TCE,ʹLineWidthʹ,3) Để trục thứ 2 không gây trở ngại cho trục thứ nhất ta viết : >>set(h2,ʹYAxisLocationʹ,ʹrightʹ,ʹColorʹ,ʹnoneʹ,ʹXTickLabelʹ,[]) >>set(h2,ʹXLimʹ,get(h2,ʹXLimʹ),ʹLayerʹ,ʹtopʹ) Để ghi chú lên đồ thị ta viết: >>text(11,380,ʹMat doʹ,ʹRotationʹ,‐‐55,ʹFontSizeʹ,16) >>ylabel(ʹTCE Mat do (PPM)ʹ) >>title(ʹXep chong do thiʹ,ʹFontSizeʹ,16) 1.5.8.3 Đồ hoạ vùng Hàm area hiển thị đường cong tạo từ một vector hay từ một cột của ma trận. Nó vẽ các giá trị của một cột của ma trận thành một đường cong riêng và tô đầy vùng không gian giữa các đường cong và trục x. >>Y = [5 1 2 8 3 7 9 6 8 5 5 5 4 2 3]; >>area(Y)
27. 27. BÀI 1: PHẦN MỀM MATLAB 25 1.5.8.4 Đồ thị pie Đồ thị pie hiển thị theo tỉ lệ phần trăm của một phần tử của một vector hay một ma trận so với tổng các phần tử. pie và pie3 tạo ra đồ thị 2D và 3D. >>X = [19.3 22.1 51.6; 34.2 70.3 82.4; 61.4 82.9 90.8; 50.5 54.9 59.1; 29.4 36.3 47.0]; >>x = sum(X); >>explode = zeros(size(x)); >>[c,offset] = max(x); >>explode(offset) = 1; >>h = pie(x,explode); >>h = pie3(x,explode); 1.5.9 Đồ hoạ 3D 1.5.9.1 Các hàm cơ bản Hàm mesh và surf tạo ra mặt 3D từ ma trận dữ liệu. Gọi ma trận dữ liệu là z mà mỗi phần tử của nó z(i,j) xác định tung độ của mặt thì mesh(z) tạo ra một lưới có màu thể hiện mặt z còn surf(z) tạo ra một mặt có màu z. 1.5.9.2 Đồ thị các hàm hai biến Bước thứ nhất để thể hiện hàm 2 biến z=f(x,y) là tạo ma trận x và y chứa các toạ độ trong miền xác định của hàm. Hàm meshgrid sẽ biến đổi vùng xác định bởi 2 vector x và y thành ma trận x và y. Sau đó ta dùng ma trận này để đánh giá hàm. >>[x,y] = meshgrid(-8:.5:8); >>r = sqrt(x.^2 + y.^2); Ma trận r chứa khoảng cách từ tâm của ma trận. Tiếp theo ta dùng hàm mesh để vẽ hàm. >>z = sin(r)./r; >>mesh(z)
28. 28. 26 BÀI 1: PHẦN MỀM MATLAB 1.5.10 Thực hành vẽ đồ thị Bài 1.10. Vẽ đồ thị hàm số y1=sinx.cos2x và hàm số y2=sinx2 trong [0-2], trên cùng hệ trục tọa độ, ta lần lượt thực hiện như sau: >>x=0:0.01:2*pi; >>y1=sin(x).*cos(2*x); %nhan tuong ung tung phan tu >>plot(x,y1) >>grid on %hien thi luoi Sau khi thu được đồ thị hàm y1, để vẽ y2 trên cùng đồ thị, ta thực hiện: >>hold on %giu hinh, mac nhien la hold off >>y2=sin(x.^2); %luy thua tung phan tu >>plot(x,y2,’k’) %duong ve co mau den >>axis([0 4*pi –1.25 1.25]) %dinh lai toa do hien thi Ta có thể đặt nhãn cho các trục cũng như tiêu đề cho đồ thị: >>xlabel(‘Time’) >>ylabel(‘Amplitude’) >>title(‘y1=sinx.cos2x and y2=sin(x^2)’) >>legend(‘sinx.cos2x’,’sinx^2’) Bài 1.11. Thực hiện như trên nhưng dùng các hàm semilogx, semilogy, loglog thay thế cho plot. Bài 1.12. Thực hiện như trên cho hàm số y = 𝑒−𝑥 2𝑒−𝑥+2 Bài 1.13. Vẽ hàm số r = sin (5) trong toạ độ cực: >>theta=0:0.05:2*pi; >>r=sin(5*theta); >>polar(theta,r) Bài 1.14. Vẽ hàm số r = 2sin() + 3cos() Bài 1.15. Vẽ hàm số 2x2 + y2 = 10 ở dạng toạ độ cực. Gợi ý: x = rcos, y = rsin Bài 1.16. Dùng hàm text xuất các công thức sau ra trục toạ độ:
29. 29. BÀI 1: PHẦN MỀM MATLAB 27 ∫ 𝑒−𝑥 + 𝑒 𝑥 √2𝑥 + 1 sin(𝜔( 𝑥 − 𝜏)) 𝑑𝑥 ∞ −∞ ∑ 𝑧−𝑛 𝑥(𝑛) ∞ 𝑛=−1 𝑦̈ = 𝑥 [ 𝑥 𝑥2 𝑥3 𝑥2 𝑥3 𝑥 𝑥3 𝑥 𝑥2 ] - MATLAB cung cấp nhiều hàm vẽ đồ thị 3D, chẳng hạn: plot3 - dùng để vẽ các đường trong không gian 3 chiều; mesh và surf - dùng để vẽ vật thể 3D. Bài 1.17. Vẽ đồ thị 3D bằng hàm plot3: >>t=0:pi/50:10*pi; >>x=sin(t); >>y=cos(t); >>z=t; >>subplot(121), plot3(x,y,z) >>grid on >>subplot(122), plot3(x,y,t.^2) >> grid on Bài 1.18. Vẽ mặt paraboloid z=x2 +y2 trong không gian 3 chiều: >>close all >>t=-5:0.1:5; >> [x,y]=meshgrid(t); %dinh luoi ve >>z=x.^2+y.^2; >> subplot(2,2,1), mesh(z) >> title('mesh(z)') >> subplot(2,2,2), meshc(z) >> title('meshc(z)') >> subplot(2,2,3), meshz(z) >> title('meshz(z)') >> subplot(2,2,4), waterfall(z)
30. 30. 28 BÀI 1: PHẦN MỀM MATLAB >> title('waterfall(z)') Bài 1.19. Nhận xét về các hàm vẽ trên. Bài 1.20. Vẽ mặt 𝑧 = 𝑠𝑖𝑛(√𝑥2+𝑦2) 2(√𝑥2+𝑦2) dùng hàm surf và mesh. 1.6 CÁC FILE VÀ HÀM 1.6.1 Script file (file kịch bản) Kịch bản là M‐file đơn giản nhất, không có đối số. Nó dùng khi thi hành một loạt lệnh MATLAB theo một trình tự nhất định. Ta xét ví dụ tạo ra các số Fibonacci nhỏ hơn 1000. f = [1 1]; i = 1; while(f(i)+f(i+1))<1000 f(i + 2)= f(i) +f(i+1); i = i + 1; end plot(f) Ta lưu đoạn mã lệnh này vào một file tên là fibo.m. Đây chính là một script file. Để thực hiện các mã chứa trong file fibo.m từ cửa sổ lệnh ta nhập: >>fibo Và nhấn enter. Kết quả MATLAB sẽ vẽ ra đồ thị của chuỗi Fibonacci. 1.6.2 File hàm Hàm là M-file có chứa các đối số. Ta có một ví dụ về hàm: function y = tb(x) %Tinh tri trung binh cua cac phan tu [m,n] = size(x); if m == 1 m = n; end y = sum(x)/m; Từ ví dụ trên ta thấy một hàm M-file gồm các phần cơ bản sau :
31. 31. BÀI 1: PHẦN MỀM MATLAB 29 • Một dòng định nghĩa hàm: function y = tb(x) gồm từ khoá function, đối số trả về y, tên hàm tb và đối số vào x. • Dòng kế tiếp là dòng trợ giúp đầu tiên. Vì đây là dòng văn bản nên nó phải đặt sau %. Nó xuất hiện khi ta nhập lệnh help . Phần văn bản này giúp người dùng hiểu tác dụng của hàm. • Thân hàm chứa mã MATLAB. • Các lời giải thích dùng để cho chương trình rõ ràng. Nó được đặt sau dấu %. Cần chú ý là tên hàm phải bắt đầu bằng ký tự và cùng tên với file chứa hàm. Tên hàm là tb thì tên file cũng là tb.m. Từ cửa sổ MATLAB ta gõ lệnh: >>z = 1:99; >>tb(z) Các biến khai báo trong một hàm của MATLAB là biến địa phương. Các hàm khác không nhìn thấy và sử dụng được biến này. Muốn các hàm khác dùng được biến nào đó của hàm ta cần khai báo nó là global. Nếu hàm có nhiều thông số ngõ vào và ngõ ra thì khai báo như sau: function [y1,y2,y3] = tb(x1,x2,x3) Lưu ý rằng trong một file .m có thể có nhiều hàm nhưng hàm đầu tiên phải có tên trùng với tên file. 1.6.3 Các hàm toán học cơ bản exp(x) hàm mũ cơ số e sqrt(x) căn bậc hai của x log(x) logarit cơ số e log10(x) logarit cơ số 10 abs(x) module của số phức x (giá trị tuyệt đối của số thực) angle(x) argument của số phức a conj(x) số phức liên hợp của x imag(x) phần ảo của x real(x) phần thực của x
32. 32. 30 BÀI 1: PHẦN MỀM MATLAB sign(x) dấu của x cos(x) tính cos sin(x) tính sin tan(x) tính tang acos(x) tính cos-1 asin(x) tính sin-1 atan(x) tính tang-1 cosh(x) tính 𝑒 𝑥+𝑒−𝑥 2 coth(x) tính cosh(x)/sinh(x) sinh(x) tính 𝑒 𝑥−𝑒−𝑥 2 tanh(x) tính sinh(x)/cosh(x) acosh(x) tính cosh-1 acoth(x) tính coth-1 asinh(x) tính sinh-1 atanh(x) tính tanh-1 1.6.4 Các phép toán trên hàm toán học a. Biểu diễn hàm toán học MATLAB biểu diễn các hàm toán học bằng cách dùng các biểu thức đặt trong M-file. Ví dụ để khảo sát hàm: 𝑓( 𝑥) = 2 𝑥 + 1 + 1 𝑥2 + 2 − 1 Ta tạo ra một file, đặt tên là ab.m có nội dung: function y = ab(x) y = 2./(x + 1) + 1./(x.^2 + 2) ‐ 1 ; Cách thứ hai để biểu diễn một hàm toán học trên dòng lệnh là tạo ra một đối tượng inline từ một biểu thức chuỗi. Ví dụ ta có thể nhập từ dòng lệnh như sau: >>f = inline(‘2./(x + 1) + 1./(x.^2 + 2) ‐ 1’); Ta có thể tính trị của hàm tại x = 2 như sau: >>f(2)
33. 33. BÀI 1: PHẦN MỀM MATLAB 31 b. Vẽ đồ thị của hàm Hàm fplot vẽ đồ thị hàm toán học giữa các giá trị đã cho. >>fplot(@(x)[tan(x),sin(x),cos(x)], 2*pi*[-1 1 -1 1]) c. Tìm cực tiểu của hàm Cho một hàm toán học một biến, ta có thể dùng hàm fminbnd của MATLAB để tìm cực tiểu địa phương của hàm trong khoảng đã cho. >>f=inline(ʹ1./((x‐0.3).^2+0.01)+1./(x.^2+0.04)‐6ʹ); >>x = fminbnd(f,0.3,1) Hàm fminsearch tương tự hàm fminbnd dùng để tìm cực tiểu địa phương của hàm nhiều biến. Ta có hàm three_var.m: function b = three_var(v) x = v(1); y = v(2); z = v(3); b = x.^2 + 2.5*sin(y) ‐ z^2*x^2*y^2; Và bây giờ tìm cực tiểu đối với hàm này bắt đầu từ x = -0.6, y = -1.2; z = 0.135 >>v = [-0.6 -1.2 0.135]; >>a = fminsearch(ʹthree_varʹ,v) d. Tìm điểm zero Hàm fzero dùng để tìm điểm không của hàm một biến. Ví dụ để tìm giá trị không của hàm lân cận giá trị -0.2, ta viết: >>f=inline(ʹ1./((x‐0.3).^2+0.01)+1./(x.^2+0.04)‐6ʹ); >>a = fzero(f,-0.2) 1.6.5 Thực hành trên script và function 1.6.5.1 Script Tập hợp các dòng lệnh của MATLAB được sắp xếp theo một cấu trúc nào đó và lưu thành file có phần mở rộng *.m được gọi là script file. Ta có thể chạy file này từ cửa sổ lệnh giống hệt như các lệnh của MATLAB. Cấu trúc của một script file như sau: % Phần viết sau dấu ‘%’ ở đây dùng cho lệnh help
34. 34. 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;
35. 35. BÀI 1: PHẦN MỀM MATLAB 33 pause(2); % dung trong 2 giay % ---- Hien thi logo cua MATLAB --------------- logo clear all close % ket thuc script file Sau khi lưu file này, từ cửa sổ lệnh của MATLAB, nhập: >>help vd Để thi hành script file vừa soạn, nhập: >>vd 1.6.5.2 Sử dụng các tool xây dựng sẵn MATLAB hỗ trợ một thư viện hàm rất phong phú, xây dựng trên các giải thuật nhanh và có độ chính xác cao. Ngoài các hàm cơ bản của MATLAB, tập hợp các hàm dùng để giải quyết một ứng dụng chuyên biệt nào đó gọi là Toolbox, ví dụ: Xử lý số tín hiệu (Digital Signal Processing), Điều khiển tự động (Control), mạng neural (Neural networks), … help % chuc nang toolbox >>help control % liet ke ham cua control toolbox Ta có thể tìm kiếm các hàm liên quan bằng cách cung cấp cho hàm lookfor của MATLAB một từ khóa: lookfor >>lookfor filter % tìm các hàm liên quan đến mạch lọc 1.6.5.3 Xây dựng hàm Xây dựng hàm cũng được thực hiện tương tự như script file. Tuy nhiên, đối với hàm ta cần quan tâm đến các tham số truyền cho hàm và các kết quả trả về sau khi thực hiện. Có 3 điểm cần lưu ý: - Tên hàm phải được đặt trùng với tên file lưu trữ. - Phải có từ khóa function ở dòng đầu tiên. - Trong một hàm có thể xây dựng nhiều hàm con (điều này không có trong script file). function [out1,out2,…]=tenham(in1,in2,…)
36. 36. 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ả:
37. 37. BÀI 1: PHẦN MỀM MATLAB 35 >>help gptb2 >>[x1,x2]=gptb2(1,6,-7) >>[x1,x2]=gptb2(2,7,14) >>[x1,x2]=gptb2(0,4,3) >>[x1,x2]=gptb2(1,6) Bài 1.22. Cho biết ý nghĩa của từ khóa nargin? Bài 1.23. Viết lại hàm này để kết quả chỉ trả về nghiệm số thực. Bài 1.24. Xây dựng hàm vdcongdb(a,m,method) để vẽ một số đường cong trong hệ tọa độ cực, với a là bán kính và m là số đường cong vẽ trên cùng trục tọa độ. Trường hợp này hàm không trả về giá trị. Tuỳ theo giá trị của tham số ‘method’ mà ta vẽ đồ thị tương ứng: Nếu method = ’Becnulli’: Vẽ đường Becnulli: 𝑟 = 𝑎√|2𝑐𝑜𝑠2𝜃| Nếu method = ’Astroit’: Vẽ đường Astroit: 𝑟 = 𝑎√|1 − 𝑠𝑖𝑛3𝜃 4 | Nếu method = ‘Xoanoc’: Vẽ đường xoắn ốc: r = cosθ +1 Nội dung hàm như sau: function vdcongdb(a,m,method) % Ve duong cong trong toa do cuc: vdcongdb(a,m,method) % method = 'Becnulli' - Ve duong Lemniscat Becnulli: % r=a*sqrt(abs(2*cos(2*theta))) % 'Astroit' - Ve duong Astroit: % r=a*sqrt(abs(1-sin(3*theta)/4)) % 'Xoanoc' - Ve duong xoan oc: % r=a*cos(theta)+1 % Voi: a-ban kinh; m-so duong cong ve tren cung he truc if nargin<3 error('Vui long nhap du 3 thong so cua ham') else theta=0:0.01:2*pi; method=upper(method); switch method
38. 38. 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’)
39. 39. BÀI 1: PHẦN MỀM MATLAB 37 Bài 1.25. Xây dựng hàm dudoan() để dự đoán kết quả sau mỗi lần tung một xúc xắc đồng nhất, 6 mặt. Nội dung hàm như sau: function dudoan() % Du doan ket qua sau moi lan tung ngau nhien mot xuc xac 6 mat % Chuong trinh lap lai cho den khi nguoi su dung khong doan tiep % tiep = 'y'; sai=0; dung=0; disp('Chao mung ban den voi chuong trinh du doan xuc xac!') while(lower(tiep)=='y') doan=input('Moi ban du doan ket qua (1-6):'); kqua=tungxx; if (doan ~= kqua) disp('Xin loi, ban da doan sai!') sai=sai+1; else disp('Xin chuc mung!') dung=dung+1; end tiep=input('Ban muon choi tiep(''y''/''n''):'); end disp(['Dung ' num2str(dung) ' trong tong so ' num2str(sai+dung) ' lan doan']) % subfunction -------------- function mat = tungxx() mat=floor(6*rand(1,1))+1; % end Kết luận về sự khác nhau giữa script file và hàm không có tham số vào. Bài 1.26. Viết chương trình in tam giác Pascal n dòng trong màn hình đồ họa với n được nhập từ bàn phím.
40. 40. 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)n u(n), 3n u(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)n u(n) trong khoảng (-10,10); (0,1000); (0,1e6).
41. 41. BÀI 2: TÍN HIỆU RỜI RẠC THEO THỜI GIAN 39 Bài 2.4. Tính năng lượng và công suất tín hiệu x(n) = (1 - 2n )u(n) trong khoảng (- 10,10); (0,1000); (0,1e6). Năng lượng của tín hiệu trong khoảng [-N,N]: 𝐸 𝑁 = ∑ | 𝑥( 𝑛)|2 𝑁 𝑛=−𝑁 Năng lượng tín hiệu: 𝐸 = lim 𝑁→∞ 𝐸 𝑁 Công suất trung bình của tín hiệu: 𝑃 = lim 𝑁→∞ 𝐸 𝑁 2𝑁 + 1 2.2 CÁC PHÉP TOÁN - Dịch: x(n)  x(n – k) - Ảnh gương: x(n)  x(-n) - Co trên miền thời gian: x(n)  x(n) - Cộng: y(n) = x1(n) + x2(n) - Nhân: y(n) = x1(n)x2(n) - Co biên độ: y(n) = Ax(n) Dịch phải xung đơn vị 5 mẫu: >>delta=[zeros(1,15) 1 zeros(1,5)]; >>stem(n,delta); Bài 2.5. Viết chương trình và vẽ dạng tín hiệu hàm u(n – 3), u(n + 2). Bài 2.6. Viết chương trình và vẽ dạng tín hiệu hàm x(n) = 2u(n – 3) + (n – 2) trong khoảng (-10,10). Từ đó vẽ các tín hiệu x(-n), 2x(n), x(2n). Bài 2.7. Viết function thực hiện cộng và nhân 2 tín hiệu. Ví dụ thực hiện cho 2 tín hiệu x1(n) = {1,-1,2,3,-2} và x2(n) = {-2,-2,1,1,-4}.   2.3 KIỂM TRA TÍNH CHẤT TUYẾN TÍNH VÀ BẤT BIẾN Hệ thống H bất biến theo thời gian nếu và chỉ nếu: y(n) = H[x(n)]  y(n – k) = H[x(n – k)]
42. 42. 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: Tính kết hợp: Tính phân phối: 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. h2(n) h2(n)x(n) y(n) h2(n) h2(n)x(n) y(n) h2(n) h2(n) h2(n) * h2(n)x(n) y(n) y(n)x(n) h2(n) h2(n) h2(n) + h2(n)x(n) y(n) y(n)x(n)
43. 43. BÀI 2: TÍN HIỆU RỜI RẠC THEO THỜI GIAN 41 >>clear all >>clf >>n = -10:10; >>x1 = randn(size(n)); %Tín hiệu x1 ngẫu nhiên >>x2 = randn(size(n)); %Tín hiệu x2 ngẫu nhiên >>a1 = 3; a2 = -2; %a1, a2 tùy ý >>y1 = n.*x1; >>y2 = n.*x2; >>y = n.*(a1*x1 + a2*x2); >>subplot(211), stem(n,a1*y1+a2*y2); >>title(‘a_1y_1(n)+a_2y_2(n)’); >>subplot(212), stem(n,y); >>title(‘H[a_1x_1(n)+a_2x_2(n)]’); Kết luận về tính chất tuyến tính của y(n) = nx(n). Bài 2.11. Xác định tính chất tuyến tính của hệ thống có phương trình y(n) = x2 (n) và y(n) = x(n2 ) 2.4 Hệ LTI Phương trình sai phân: 𝑦(𝑛) = − ∑ 𝑎 𝑘 𝑦(𝑛 − 𝑘) 𝑁 𝑘=1 + ∑ 𝑏 𝑘 𝑥(𝑛 − 𝑘) 𝑀 𝑘=0 Bài 2.12. Xét hệ thống có phương trình sai phân: y(n) = 0.3x(n) + 0.2x(n – 1) – 0.3x(n – 2) -0.9y(n – 1) + 0.9y(n – 2). Xác định đáp ứng xung đơn vị của hệ thống. >>N = 40; >>num = [0.3 0.2 -0.3]; >>den = [1 0.9 -0.9]; >>h = impz(num,den,N);% N: số lượng mẫu tính toán >>stem(h); Xác định ngõ ra khi biết đáp ứng xung và ngõ vào: >>x = randn(1,10);
44. 44. 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: >>h2 = [1 2 -2 -3]; Hệ thống 1 >>h2 = [-2 0 3 1]; Hệ thống 2 >>h = conv(h2,h2); >>N = 30; >>x = randn(1,N); >>y11 = conv(x,h2); >>y1 = conv(y11,h2); >>y21 = conv(x,h2); >>y2 = conv(y21,h2); >>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:
45. 45. BÀI 2: TÍN HIỆU RỜI RẠC THEO THỜI GIAN 43 Bài 2.16. Xác định đáp ứng xung tương đương của hệ thống sau: h2 = {1,0,-1,3}; h2 = {2,-2,1}; h3 = {3,4,-1,1}; h4 = {-3,5,6,-1,1}     Viết chương trình xác định ngõ ra của hệ thống khi ngõ vào là x(n) = (-2)n u(n) (tính toán cho giá trị n từ -20 đến 20): a. Dùng theo sơ đồ trên. b. Dùng đáp ứng xung tương đương. h2(n) h2(n) h3(n) x(n) y(n) h4(n)
46. 46. 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: 𝑋( 𝑧) = 𝑁(𝑧) 𝐷(𝑧) = 𝑏0 𝑧−𝑀 𝑎0 𝑧−𝑁 𝑧 𝑀 + 𝑏1 𝑏0 𝑧 𝑀−1 + ⋯ + 𝑏 𝑀 𝑏0 𝑧 𝑁 + 𝑎1 𝑎0 𝑧 𝑀−1 + ⋯ + 𝑎 𝑀 𝑎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);
47. 47. BÀI 3: BIẾN ĐỔI Z 45 Ta có thể vẽ các điểm cực và điểm không nếu đã biết điểm cực và điểm không bằng cách đưa thông số vào hàm zplane ở dạng vector cột: >>zero = [-1 1+j*1]; >>pole = [j*2 -1+j]; >>zplane(zero’,pole’); Để xác định điểm cực và không, ta dùng hàm tf2zp: [z,p,k] = tf2zp(num,den) trong đó z, p là các điểm cực và không lưu dạng vector hàng, k là hệ số khuếch đại: >>num = [1 2 3]; % Tử số >>den = [2 4 7]; % Mẫu số >>[z,p,k] = tf2zp(num,den) Nếu đã cho điểm cực và điểm không, ta có thể xác định lại biểu thức của biến đổi z bằng hàm zp2tf: [num,den] = zp2tf(z,p,k) (z, p ở dạng vector cột) >>zero = [-1 1+j*1]; >>pole = [j*2 -1+j]; >>k = 2; >>[num,den] = zp2tf(zero’,pole’,k) Bài 3.2. Xác định và vẽ điểm cực, điểm không của các hàm hệ thống sau: 𝐻( 𝑧) = 2 + 5𝑧−1 + 4𝑧−2 + 5𝑧−3 + 3𝑧−4 5 + 45𝑧−1 + 2𝑧−2 + 𝑧−3 + 𝑧−4 Từ đó xác định các miền hội tụ có thể có. So sánh với lý thuyết. Bài 3.3. Xác định biểu thức của biến đổi z có các điểm cực 0.5; 0.75; 1+j0.5; 1-j0.5 và các điểm không 0.3; 0.1; 2-j2; 2+j2 với hệ số khuếch đại k = 0.7 3.2 PHÂN TÍCH DÙNG PHƯƠNG PHÁP THẶNG DƯ Phân tích thành các thừa số theo phương pháp thặng dư: 𝑁1(𝑧) 𝐷(𝑧) = 𝐴1 1 − 𝑝1 𝑧−1 + ⋯ + 𝐴 𝑁 1 − 𝑝 𝑁 𝑧−1 Bài 3.4. Xác định các hệ số của biểu thức biến đổi z bằng hàm residuez: >>num = [1 2 3]; % Tử số >>den = [2 4 7]; % Mẫu số
48. 48. 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: 𝑋( 𝑧) = 𝑁(𝑧) 𝐷(𝑧) = ∑ 𝑐 𝑘 𝑧−𝑘 𝑀−𝑁 𝑘=0 + 𝑁1(𝑧) 𝐷(𝑧) (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ư:
49. 49. BÀI 3: BIẾN ĐỔI Z 47 𝑁1(𝑧) 𝐷(𝑧) = 𝐴1 1 − 𝑝1 𝑧−1 + ⋯ + 𝐴 𝑁 1 − 𝑝 𝑁 𝑧−1 Nếu các giá trị pj = … = pm thì chuyển các số hạng từ Aj đến Am thành: 𝐴𝑗 1 − 𝑝 𝑁 𝑧−1 + 𝐴𝑗+1 (1 − 𝑝 𝑁 𝑧−1)2 + ⋯ + 𝐴𝑗+𝑚−1 (1 − 𝑝 𝑁 𝑧−1) 𝑚 Áp dụng kết quả: 𝐴 1 − 𝑎𝑧−1 𝑧 ↔ 𝐴𝑎 𝑛 𝑢( 𝑛), 𝑅𝑂𝐶: | 𝑧| > | 𝑎| Bài 3.7. Dùng hàm ztrans để biến đổi z ở dạng công thức: >>syms n x >>x = 2^n; >>ztrans(x) >>x = (-1/2)^n; >>ztrans(x) Bài 3.8. Xác định biến đổi z của các hàm sau: a. x(n) = (-2)n-1 u(n) b. x(n) = n3n u(n) c. x(n) = n2 4n u(n) Bài 3.9. Biến đổi z ngược theo giá trị bằng hàm impz. >>num = [1 1 2]; >>den = [1 -1 2]; >>L = 50; %Số lượng mẫu cần tính >>x = impz(num,den,L) % x là biến đổi z ngược >>impz(num,den,L); % Vẽ trên đồ thị Bài 3.10. Xác định và vẽ 100 mẫu đầu tiên của biến đổi z ngược của hàm: 𝑋( 𝑧) = 0.9 + 0.7𝑧−1 + 0.1𝑧−2 − 𝑧−3 + 0.5𝑧−4 1 + 0.5𝑧−1 − 0.2𝑧−3 + 2𝑧−4 + 𝑧−5 Bài 3.11. Thực hiện lại bài 3.10 với kết quả ở dạng công thức. Để xác định biến đổi z ngược ở dạng công thức, ta dùng hàm residuez để phân tích thành dạng biểu thức đơn giản và dùng các kết quả biến đổi ngược để xác định. Bài 3.12. Ta cũng có thể xác định biến đổi z ngược bằng cách dùng hàm iztrans.
50. 50. 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 + 1 2 𝑧−1) 2 (1 + 2𝑧−1) 𝑌(𝑧) = 2𝑧−1 (3 + 4𝑧−1 + 𝑧−2)(1 + 𝑧−1)
51. 51. BÀI 4: BIẾN ĐỔI FOURIER RỜI RẠC 49 BÀI 4: BIẾN ĐỔI FOURIER RỜI RẠC Sau khi học xong bài này, người học có thể:  Chuyển tín hiệu rời rạc trên miền thời gian sang miền tần số bằng cách dùng biến đổi Fourier rời rạc.  Tính toán biến đổi Fourier bằng một số thuật toán. 4.1 TÍNH DTFT Biến đổi Fourier rời rạc (DTFT) mô tả như sau: 𝑋( 𝑒 𝑗𝜔) = ∑ 𝑥( 𝑛) 𝑒−𝑗𝑛𝜔 ∞ 𝑛=−∞ = |𝑋( 𝑒 𝑗𝜔)|𝑒 𝑗Φ(𝜔) Bài 4.1. Tính và vẽ DTFT có dạng: 𝑋( 𝑒−𝑗𝜔) = 1 + 𝑒−𝑗𝜔 − 2𝑒−𝑗2𝜔 1 + 0.5𝑒−𝑗𝜔 >>w = linspace(-4*pi,4*pi,512); % Tạo 512 giá trị từ -4 đến 4 >>num = [1 1 -2]; >>den = [1 0.5]; >>h = freqz(num,den,w); >>subplot(211),plot(w/pi,abs(h)); >>xlabel(‘omega/pi’); >>ylabel(‘Bien do’); >>title(‘Pho bien do |X(e^{jomega})|’); >>subplot(212),plot(w/pi,angle(h)); >>xlabel(‘omega/pi’); >>ylabel(‘Pha [rad]’);
52. 52. 50 BÀI 4: BIẾN ĐỔI FOURIER RỜI RẠC >>title(‘Pho pha arg(X(e^{jomega}))’); Để 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);
53. 53. BÀI 4: BIẾN ĐỔI FOURIER RỜI RẠC 51 >>subplot(321),stem(abs(x)); >>title(‘x(n)’);xlabel(‘n’); >>subplot(323),stem(abs(y)); >>title([‘FFT(n,’ N ‘)’]);xlabel(‘k’);ylabel(‘Bien do’); >>subplot(324),stem(angle(y)); >>title([‘FFT(n,’ N ‘)’]); xlabel(‘k’);ylabel(‘Pha’); >>subplot(325),stem(abs(x1)); >>title([‘IFFT(FFT(n,’ N ‘))’]); >>xlabel(‘k’);ylabel(‘Bien do’); >>subplot(326),stem(angle(x1)); >>title([‘IFFT(FFT(n,’ N ‘))’]); >>xlabel(‘k’);ylabel(‘Pha’); Bài 4.9. Xác định và vẽ FFT 16 điểm của x(n) cho tùy ý. Bài 4.10. Tạo function cshift để dịch vòng một chuỗi m giá trị: function out = cshift(x,m) m0 = m; if abs(m0) > length(x) m0 = rem(m0,length(x)); end while (m0<0) m0 = m0 + length(x); end out= [x(length(x)-m0+1:length(x)) x(1:length(x)-m0)]; Khảo sát tính chất dịch vòng của DFT-N điểm: >>N = 20; >>m = 3; >>x = randn(1,N); >>y = fft(x,N); >>x1 = cshift(x,m); %Có thể thay bằng hàm circshift như sau: x1 = circshift(x’,m)’;
54. 54. 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.
55. 55. BÀI 5: BỘ LỌC SỐ FIR 53 BÀI 5: BỘ LỌC SỐ FIR Sau khi học xong bài này, người học có thể:  Thiết kế một bộ lọc số FIR dựa theo các thông số cho trước. 5.1 CÁC LOẠI BỘ LỌC Điều kiện đối xứng và phản đối xứng của mạch lọc: h(n) =  h(M – 1 – n) Dựa trên tính chất đối xứng hay phản đối xứng của chuỗi đáp ứng xung và chiều dài N của chuỗi đáp ứng xung, người ta phân loại bộ lọc FIR làm 4 loại sau: • Bộ lọc FIR loại 1: h(n) đối xứng, M lẻ, β = 0,  = (M – 1)/2 • Bộ lọc FIR loại 2: h(n) đối xứng, M chẵn, β = 0,  = (M – 1)/2 • Bộ lọc FIR loại 3: h(n) phản đối xứng, M lẻ, β = /2,  = (M – 1)/2 • Bộ lọc FIR loại 4: h(n) phản đối xứng, M chẵn, β = /2,  = (M – 1)/2 Đáp ứng tần số của FIR cho từng loại: • Loại 1: H(ω) = [∑ a( 𝑛)cosωn 𝑀−1 2 𝑛=0 ] e−j (M−1) 2 ω với { 𝑎(0) = ℎ ( 𝑀−1 2 ) 𝑎( 𝑛) = 2ℎ ( 𝑀−1 2 − 𝑛) 𝑣ớ𝑖 1 ≤ 𝑛 ≤ 𝑀−1 2 • Loại 2: H(ω) = [∑ b( 𝑛)cosω (n − 1 2 ) 𝑀 2 𝑛=1 ] e−j (M−1) 2 ω với 𝑏( 𝑛) = 2ℎ ( 𝑀 2 − 𝑛) 𝑣ớ𝑖 1 ≤ 𝑛 ≤ 𝑀 2 • Loại 3: H(ω) = [∑ c( 𝑛)sinωn 𝑀−1 2 𝑛=1 ] e−j( π 2 − (M−1) 2 ω) với 𝑐( 𝑛) = 2ℎ ( 𝑀−1 2 − 𝑛) 𝑣ớ𝑖 1 ≤ 𝑛 ≤ 𝑀−1 2 • Loại 4: H(ω) = [∑ d( 𝑛)sinω (n − 1 2 ) 𝑀 2 𝑛=1 ] e−j( π 2 − (M−1) 2 ω) với 𝑑( 𝑛) = 2ℎ ( 𝑀 2 − 𝑛) 𝑣ớ𝑖 1 ≤ 𝑛 ≤ 𝑀 2
56. 56. 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 h2 = [1.5 -2.5 3 -2.5 1.5]: >>h2 = [1.5 -2.5 3 -2.5 1.5]; >> [a,w,L,Hr]=FIR_t1(h2); 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)’);
57. 57. BÀI 5: BỘ LỌC SỐ FIR 55 b. Các hệ số của bộ lọc: >>[a,w,L,Hr] = FIR_t1(h); >>subplot(222);stem(0:L,a); >>title(‘Cac he so a(n)’);xlabel(‘n’);ylabel(‘a(n)’); c. Đáp ứng tần số: >>w = linspace(0,2*pi,100)'; >>subplot(223);plot(w,Hr); >>title(‘Dap ung tan so’);xlabel(‘omega’);ylabel(‘H(omega)’); d. Phân bố cực - không: >>subplot(224); zplane(h,1); >>title(‘Bieu do cuc – khong’);xlabel(‘Thuc’);ylabel(‘Ao’); Bài 5.6. Thực hiện lại bài 5.5 với đáp ứng xung h = [-1 2 1.3 -2.2 0.6 3 3 0.6 -2.2 1.3 2 -1]. Bài 5.7. Thực hiện lại bài 5.5 với đáp ứng xung h = [-1 2 1.3 -2.2 0.6 3 -0.6 2.2 -1.3 -2 1]. Bài 5.8. Thực hiện lại bài 5.5 với đáp ứng xung h = [-1 2 1.3 -2.2 0.6 3 -3 -0.6 2.2 -1.3 -2 1]. 5.2 PHƯƠNG PHÁP CỬA SỔ Bài 5.9. Tạo hàm ideal_lp xác định đáp ứng xung của bộ lọc thông thấp lý tưởng theo tần số cắt c và chiều dài chuỗi đáp ứng xung. function hd = ideal_lp(wc,M) alpha = (M-1)/2; n = [0:1:(M-1)]; m = n - alpha + eps; hd = sin(wc*m)./(pi*m); Bài 5.10. Tạo hàm freqz_m tính toán độ lớn và pha của đáp ứng tần số, hàm trễ nhóm. function [db,mag,pha,grd,w] = freqz_m(b,a)
58. 58. 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: Rp = −20𝑙𝑜𝑔 1−𝛿 𝑝 1+𝛿 𝑝
59. 59. BÀI 5: BỘ LỌC SỐ FIR 57 • Độ suy giảm dải chắn theo dB: A 𝑆 = −20𝑙𝑜𝑔 𝛿 𝑠 1+𝛿 𝑝 Các hàm cửa sổ trong MATLAB: boxcar hay rectwin, bartlett, hanning, hamming, blackman, Kaiser. >>wp = 0.2*pi; ws =0.3*pi; >>tr_width = ws - wp; >>M = ceil(6.6*pi/tr_width) + 1; >>n = [0:1:M-1]; >>wc = (ws+wp)/2; >>hd = ideal_lp(wc,M); >>w_ham = (hamming(M))'; >>h = hd.*w_ham; >>[db,mag,pha,grd,w] = freqz_m(h,[1]); >>delta_w = 2*pi/1000; >>Rp = -(min(db(1:1:wp/delta_w+1))) >>As = -round(max(db(ws/delta_w+1:1:501))) a. Đáp ứng xung: >>subplot(221); stem(n,hd); >>axis([0,M-1,-0.1,0.3]); >>title('Dap ung xung ly tuong');xlabel('n'); ylabel('h_d(n)'); b. Cửa sổ Hamming: >>subplot(222); stem(n,w_ham); >>axis([0,M-1,0,1.1]); >>title('Cua so Hamming');xlabel('n'); ylabel('w(n)'); c. Đáp ứng tần số: >>subplot(223); stem(n,h); >>axis([0,M-1,-0.1,0.3]);
60. 60. 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à: H() = ∑ ℎ(𝑛)𝑒−𝑗𝜔𝑛𝑀−1 𝑛=0 Ta chỉ định một tập tần số như sau: k = 2𝜋 𝑀 ( 𝑘 + 𝛼), { 𝑘 = 0, 1, … , 𝑀−1 2 𝑛ế𝑢 𝑀 𝑙ẻ 𝑘 = 0, 1, … , 𝑀 2 − 1 𝑛ế𝑢 𝑀 𝑐ℎẵ𝑛 𝛼 = 0 ℎ𝑎𝑦 1 2 Khi đó: H(k + α) = H( 2𝜋 𝑀 (𝑘 + 𝛼)) = ∑ ℎ(𝑛)𝑒−𝑗 2𝜋 𝑀 (𝑘+𝛼)𝑛𝑀−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).
61. 61. BÀI 5: BỘ LỌC SỐ FIR 59 Bài 5.13. Thiết kế bộ lọc thông thấp theo phương pháp lấy mẫu tần số với các tham số như sau: p = 0.2; s = 0.3; Rp = 0.25 dB; As = 50 dB Chọn đáp ứng xung có chiều dài 60 ứng với 60 mẫu tần số trong khoảng [0,2). Dải thông có độ rộng là 0.2 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.5925 và T2 = 0.1099. Mẫu các tần số được cho như sau: H(k) = [1,1,1,1,1,1,1,T1,T2,0,0,…,0,0,T2,T1,1,1,1,1,1,1] (43 số 0) >>M = 60; alpha = (M-1)/2; L = 0:M-1; wl = (2*pi/M)*L; >>Hk = [ones(1,7),0.5925,0.1099,zeros(1,43), 0.1099, 0.5925, ones(1,6)]; % Dap ung tan so mau ly tuong >>Hdr = [1,1,0,0]; wdl = [0,0.2,0.3,1]; % Dap ung tan so ly tuong de bieu dien do thi >>k1 = 0:floor((M-1)/2); k2 = floor((M-1)/2)+1:M-1; >>angH = [-alpha*(2*pi)/M*k1, alpha*(2*pi)/M*(M-k2)]; >>H = Hk.*exp(j*angH); >>h = real(ifft(H,M)); >>[db,mag,pha,grd,w] = freqz_m(h,1); >>[a,ww,L,Hr] = FIR_t2(h); a. Chuỗi mẫu tần số: >>subplot(221); plot(wl(1:31)/pi,Hk(1:31),'o',wdl,Hdr); >>axis([0,1,-0.1,1.1]); >>title('Cac mau tan so: M=60, T2 = 0.5925, T1 = 0.1099'); >>xlabel(f[*pi]'); ylabel('Hr(k)'); b. Đáp ứng xung của bộ lọc thực tế: >>subplot(222); stem(L,h);
62. 62. 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:
63. 63. BÀI 5: BỘ LỌC SỐ FIR 61 𝐸(𝜔) = 𝑊(𝜔)[𝐻 𝑑(𝜔) − 𝐻(𝜔)] = 𝑊(𝜔)𝑄(𝜔) [ 𝐻 𝑑(𝜔) 𝑄(𝜔) − 𝑃(𝜔)] Ta định nghĩa: 𝑊̂ (𝜔) = 𝑊(𝜔)𝑄(𝜔) 𝐻̂ 𝑑(𝜔) = 𝐻 𝑑(𝜔) 𝑄(𝜔) Khi đó: 𝐸(𝜔) = 𝑊̂ (𝜔)[𝐻̂ 𝑑(𝜔) − 𝑃(𝜔)] Parks và McClellan đã đưa ra giải pháp sử dụng thuật toán Remez để tìm ra đáp ứng xung của bộ lọc tối ưu, tức là gần đúng theo nghĩa Chebyshev đối với một bộ lọc lý tưởng, cho giá trị M là chiều dài của chuỗi đáp ứng xung với các điều kiện ràng buộc về độ gợn sóng ở dải thông và dải chắn như sau: 1. Xác định loại bộ lọc, tính giá trị R và xây dựng các hàm W(ω), Q(ω). Loại bộ lọc R P() Q() FIR loại 1 𝑀 − 1 2 ∑ 𝑎̅(𝑛)𝑐𝑜𝑠𝜔𝑛 𝑅 𝑛=0 1 FIR loại 2 𝑀 2 − 1 ∑ 𝑏̅(𝑛)𝑐𝑜𝑠𝜔𝑛 𝑅 𝑛=0 cos(/2) FIR loại 3 𝑀 − 1 2 − 1 ∑ 𝑐̅( 𝑛)𝑐𝑜𝑠𝜔𝑛 𝑅 𝑛=0 sin FIR loại 4 𝑀 2 − 1 ∑ 𝑑̅( 𝑛)𝑐𝑜𝑠𝜔𝑛 𝑅 𝑛=0 sin(/2) 2. Xây dựng bài toán gần đúng bằng cách xác định các hàm 𝑊̂ (𝜔), 𝐻 𝑑(𝜔). 3. Sử dụng thuật toán trao đổi Remez để tìm ra hàm tối ưu P(ω).  Chọn lấy R+2 điểm rời rạc, giả sử đó là các cực trị của hàm sai số.  Trên cơ sở tại R+2 điểm rời rạc nói trên, hàm E(ω) luân phiên đổi dấu và có trị tuyệt đối bằng một giá trị δ nào đó, tính nội suy lại giá trị δ và hàm P(ω), từ đó tính ra hàm sai số E(ω), tính được cực trị thực của hàm sai số.
64. 64. 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à: 𝑀0 = −20𝑙𝑜𝑔√𝛿1 𝛿2 − 13 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];
65. 65. BÀI 5: BỘ LỌC SỐ FIR 63 h = firpm(M-1,f,m,weights); [db,mag,pha,grd,w] = freqz_m(h,[1]); Asd = -max(db(wsi:1:501)) while Asd
66. 66. 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
67. 67. BÀI 6: BỘ LỌC SỐ IIR 65 BÀI 6: BỘ LỌC SỐ IIR Sau khi học xong bài này, người học có thể:  Thiết kế một bộ lọc IIR dựa theo các thông số cho trước. 6.1 THIẾT KẾ BỘ LỌC TƯƠNG TỰ Các yêu cầu thiết kế cho bộ lọc tương tự dựa trên điều kiện giới hạn các chỉ tiêu kỹ thuật cho bình phương biên độ của hàm đáp ứng tần số. Hàm truyền Ha(s) của một hệ thống tuyến tính bất biến tương tự là biến đổi Laplace hàm đáp ứng xung ha(t) của hệ thống và Ha(s) chính bằng tỷ số giữa biến đổi Laplace của tín hiệu đầu ra với biến đổi Laplace của tín hiệu đầu vào. Hàm đáp ứng tần số Ha() của một hệ thống tuyến tính bất biến là hàm thể hiện đáp ứng của hệ thống với đầu vào là các giá trị tần số khác nhau. Do đó, hàm đáp ứng tần số là biến đổi Fourier hàm đáp ứng xung ha(t) của hệ thống và hàm Ha() cũng chính là hàm hệ thống Ha(s) đánh giá trên trục ảo. Đối với các hệ thống thực hiện được về mặt vật lý, đáp ứng xung bao giờ cũng là một hàm thực. Do đó, hàm truyền luôn đối xứng qua trục thực. Mặt khác, một hệ thống tương tự là nhân quả và ổn định nếu và chỉ nếu tất cả các điểm cực của hàm truyền đạt nằm ở nửa bên trái của mặt phẳng s hoặc nằm ở gốc toạ độ. Khi ta xét đến bình phương biên độ của hàm hệ truyền đạt, các điểm cực của hàm số này sẽ phân bố trên tất cả các góc phần tư của mặt phẳng S. Lúc này, việc xét đáp ứng tần số của hệ thống cũng thuận tiện và tất cả các điểm cực nằm bên trái mặt phẳng S của hàm |Ha(s)|2 . Yêu cầu về chỉ tiêu kỹ thuật của bộ lọc thông thấp tương tự thường được cho dưới dạng như sau: 1 1 + 𝜖2 ≤ |𝐻 𝑎(𝜔)| ≤ 1, |𝜔| ≤ 𝜔 𝑝 0 ≤ |𝐻 𝑎(𝜔)| ≤ 1 𝐴2 , |𝜔| ≥ 𝜔𝑠 với p và s lần lượt là các tần số cắt dải thông và tần số cắt dải chắn,  là tham số gợn sóng và A là tham số suy giảm. Quan hệ giữa  và A với Rp và As hay δp và δs: 𝑅 𝑝 = −10𝑙𝑜𝑔 1 1 + 𝜖2 → 𝜖 = √10 𝑅 𝑝 10 − 1
68. 68. 66 BÀI 6: BỘ LỌC SỐ IIR 𝐴 𝑠 = −10𝑙𝑜𝑔 1 𝐴2 → 𝐴 = 10 𝐴 𝑠 20 1 − 𝛿 𝑝 1 + 𝛿 𝑝 = √ 1 1 + 𝜖2 → 𝜖 = 2√ 𝛿 𝑝 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
69. 69. BÀI 6: BỘ LỌC SỐ IIR 67 % w = Cac mau tan so tren doan tu 0 den wmax % b = Cac he so da thuc tu so cua Ha(s) % a = Cac he so da thuc mau so cua Ha(s) % wmax = Tan so cuc dai theo don vi rad/sec tren doan % tan so mong muon tim dap ung tan so w = [0:1:500]*wmax/500; H = freqs(b,a,w); mag = abs(H); db = 20*log10((mag+eps)/max(mag)); pha = angle(H); Hàm cheb1ap trả về danh sách các điểm không, điểm cực và độ lợi của hàm truyền đạt cho thiết kế bộ lọc dạng Chebyshev I, tần số cắt đã được chuẩn hoá. Đáp ứng tần số của bộ lọc thông thấp Chebyshev I bậc N: |𝐻 𝑎(𝜔)| = 1 1 + 𝜖2 𝑇𝑁 2 ( 𝜔 𝜔𝑐 ) Trong đó  là tham số gợn sóng dải thông, c là tần số cắt và TN(x) là đa thức Chebyshev bậc N: 𝑇 𝑁(𝑥) = { cos[𝑁𝑐𝑜𝑠−1(𝑥)] −1 ≤ 𝑥 ≤ 1 cosh[𝑁𝑐𝑜𝑠ℎ−1(𝑥)] |𝑥| > 1 Giá trị thích hợp của N cho như sau: 𝑁 = ⌈ log(𝑔 + √𝑔2 − 1) log(𝜔𝑟 + √𝜔𝑟 2 − 1 ⌉ Với: 𝑔 = √ 𝐴2 − 1 𝜖2 𝜔𝑟 = 𝜔𝑠 𝜔 𝑝 Bài 6.2. Tạo hàm u_chb1ap nhằm trả về hệ số của các đa thức tử số và đa thức mẫu số của hàm truyền đạt cho thiết kế bộ lọc dạng Chebyshev I có tần số cắt tuỳ ý: function [b,a] = u_chb1ap(N,Rp,Omegac) % Bo loc thong thap dang Chebyshev-1 % tan so cat khong duoc chuan hoa
70. 70. 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
71. 71. BÀI 6: BỘ LỌC SỐ IIR 69 error('Tan so cat dai thong phai lon hon 0') end if Ws <= Wp error('Tan so cat dai thong phai lon hon tan so cat dai chan') end if (Rp <= 0) | (As<0) error('Do gon dai thong va/hoac Do suy giam dai chan phai lon hon 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*** Bac cua bo loc Chebyshev-1 = %2.0fn',N); [b,a] = u_chb1ap(N,Rp,OmegaC); Bài 6.4. Thiết kế bộ lọc thông thấp tương tự, định dạng Chebyshev-I với các tham số đầu vào như sau: p = 0.2; s = 0.3; Rp = 1 dB; As = 16 dB Viết chương trình tính và biểu diễn trên đồ thị: a. Độ lớn của đáp ứng tần số b. Đáp ứng pha của bộ lọc c. Độ lớn tương đối tính theo dB của đáp ứng tần số d. Đáp ứng xung của bộ lọc tương tự 6.2 THIẾT KẾ BỘ LỌC SỐ Biến đổi song tuyến tính là phép ánh xạ từ mặt phẳng s sang mặt phẳng z định nghĩa như sau:
72. 72. 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');
73. 73. BÀI 6: BỘ LỌC SỐ IIR 71 c. Hàm đáp ứng pha: subplot(223); plot(w/pi,pha/pi); axis([0,1,-1,1]); grid title('Dap ung pha'); xlabel('f[*pi]'); ylabel('Angle(H_r(omega))'); d. Trễ nhóm theo tần số: subplot(224); plot(w/pi,grd); axis([0,1,0,15]); grid title('Group Delay'); xlabel(' f[*pi]'); ylabel('Samples'); Gõ lệnh: >>bai605 tại cửa sổ lệnh để kiểm tra kết quả. Phương pháp bất biến xung thực hiện lấy mẫu bộ lọc tương tự để tạo ra bộ lọc số: h(n) = ha(nT) Đáp ứng tần số của bộ lọc số có quan hệ với đáp ứng tần số của bộ lọc tương tự như sau: 𝐻(𝜔) = 1 𝑇 ∑ 𝐻 𝑎 [ ( 𝜔 − 2𝜋𝑘) 𝑇 ] ∞ 𝑘=−∞ Bài 6.6. Thực hiện yêu cầu của bài 6.5 theo phương pháp bất biến xung, dùng hàm impinvar của MATLAB. So sánh kết quả thu được với câu trên. Bài 6.7. Tạo hàm zmapping thực hiện việc chuyển đổi băng tần số, trả về hàm truyền của bộ lọc mới với tham số đầu vào là hàm truyền đạt của bộ lọc thông thấp, hàm đa thức thể hiện phép đổi biến số độc lập: % function [bz,az] = zmapping(bZ,aZ,Nz,Dz) 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];
74. 74. 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ố.
75. 75. TÀI LIỆU THAM KHẢO 1. Phạm Hùng Kim Khánh. Xử lý tín hiệu số. ĐH Công nghệ TPHCM. 2. Sanjit K. Mitra. Digital Signal Processing Laboratory using MatLab. 3. Đinh Đức Anh Vũ, Vũ Tuấn Thanh, Lê Trọng Nhân, Tôn Thất Đại Hải. Giáo trình Thực hành Xử lý tín hiệu số. Bộ môn Kỹ thuật Máy tính, ĐH Bách Khoa TPHCM. 4. Hồ Văn Sung. Thực hành Xử lý tín hiệu số với MatLab. NXB Khoa học và Kỹ thuật.