Bài thí nghiệm 1 Matlab

BÀI THÍ NGHIỆM 1 MATLAB I. TÓM TẮT LÝ THUYẾT 1. Các toán tử Có 2 loại toán tử cơ bản là toán tử unary(toán tử một ngôi) và binary (toán tử hai ngôi). Ví dụ kí hiệu “-” được hiểu toán tử một ngôi là giá trị âm và hiểu toán tử hai ngôi là toán tử trừ . - Dưới đây là một vài toán tử phổ biến được sử dụng với các biểu thức số học : + toán tử cộng - giá trị âm , toán tử trừ * toán tử nhân / toán tử chia lấy tử (9/3 = 3) \ toán tử chia lấy mẫu (3/12 = 4) ^ toán tử lũy thừa. 2. Vectơ và ma trận - Vector và ma trận được sử dụng để lưu tập các giá trị có cùng kiểu dữ liệu . Một vector có thể là vector hàng hoặc vector cột . Một ma trận có thể được hình thành từ một bảng giá trị. Ma trận cỡ mxn là bảng số chữ nhật gồm mxn số được sắp xếp thành m hàng và n cột : + Nếu m = 1 thì ma trận chuyển thành vector hàng + Nếu n = 1 thì ma trận chuyển thành vector cột + Nếu m = n = 1 thì ma trận trở thành đại lượng vô hướng. - Một vector trong MATLAB được coi là mảng 1 chiều trong các ngôn ngữ khác. Một ma trận được coi là mảng 2 chiều . Vì vậy trong MATLAB , thực hiện toán tử với vector và ma trận thì được coi là các toán tử mảng . 2.1 Tạo vectơ hàng Tạo 1 vectơ v = [1,2,3,4] hoặc v = [1 2 3 4] Dùng colon operator: v = 1:2:9 với 1: giá trị đầu, 9:giá trị cuối, 2:bước nhảy ta dùng v= 1:9 trường hợp bước nhảy bằng 1. Dùng linspace: v = linspace(1,9,5) với 1: giá trị đầu, 9: giá trị cuối, 5: số phần tử. 2.2 Tạo vectơ cột Ta dùng dấu ; để ngăn cách giữa các phần tử c = [1;2;3;4] hoặc dùng dấu ’ sau vectơ hàng để được vectơ cột c = v' 2.3 Tạo ma trận Kết hợp cách tạo vectơ hàng và vectơ cột mat = [1 2 3; 4 5 6] mat = [1:3;4:6] mat = [linspace(3,15,5);1:2:9] Tạo ma trận ngẫu nhiên mat = randint(2,3,[10,20]) với 2 : số hàng, 3 : số cột, [10,20] : phạm vi giá trị. Ma trận không mat = zeros(2,4) với 2: số hàng, 4 : số cột trường hợp ma trận không vuông mat = zeros(3) với 3 : số hàng, cột 3. Plot và subplot Lệnh plot trong MATLAB cho phép vẽ đồ thị từ dữ liệu là vector. Ví dụ nếu có biến t là thời gian và y là tín hiệu , gõ lệnh plot(t,y); lập tức xuất hiện giao diện đồ thị của t và y. Các hàm phổ biến cho plot: title ('') : chú thích tên đồ thị. xlabel ('') : chú thích tên trục x của đồ thị. ylabel('') : chú thích tên trục y của đồ thị. grid on : hiển thị lưới các ô vuông. - Để hiển thị nhiều đồ thị trên một màn hình , sử dụng hàm subplot. Hàm này gồm 3 đối số (m,n,p) , m và n chia màn hình ra thành m hàng và n cột . Đối số p xác định đồ thị nào xuất hiện trước . subplot (,,): hiện thị vùng muốn vẽ. stem (signal): vẽ chuỗi dữ liệu rời rạc. angle(x): argument của số phức a. real(x): phần thực của x ; imag(x): phần ảo của x. conj(x): số phức liên hợp của x. abs(x) :modun của số phức x. 4. Điều kiện if và if – else Cấu trúc if if điều khiện câu lệnh end Cấu trúc if - else if điều khiện câu lệnh 1 else câu lệnh 2 end 5. Vòng lặp for Cấu trúc: for điều khiện câu lệnh end II. BÀI TẬP Bài tập 1 : ftemp=randint(2,6,[80,100]) ; F=ftemp C=(F-32).*5/9 ; ctemp=C Bài tập 2 : 3:6 1:0.5:3 5:-1:2 Bài tập 3 : linspace(4,8,3) linspace(-3,-15,5) linspace(9,5,3) Bài tập 4 : [7:-1:5;3:2:7] Bài tập 5 : x=-2:0.1:2; plot(x, exp(x)) title('DO THI HAM e^x') xlabel('Truc x') ylabel('e^x') Bài tập 6 : x1=linspace(0,10, 10); subplot(1,2,1) plot(x1,sin(x1)) title('Do thi ham sin(x)') xlabel('x1') ylabel('sin(x1)') x2=linspace(0,10, 100); subplot(1,2,2) plot(x2,sin(x2)) title('Do thi ham sin(x)') xlabel('x2') ylabel('sin(x2)') Bài tập 7 : Scrip: if wind<=38 fprintf('A storm is a tropical depression') else if 38<wind&&wind<=73 fprintf('A storm is a tropical storm') else fprintf('A storm is a hurricane') end end Bài tập 8: clc;close all; x=randint(3,5,[0,20]) %max cua ma tran for j=1:5 for i=1:3 if max<x(i,j) max=x(i,j); end end end max %max cua moi cot trong ma tran for j=1:5 maxc=x(1,j); for i=1:3 if maxc<x(i,j) maxc=x(i,j); end end maxc end %max cua moi hang trong ma tran for i=1:3 maxh=x(i,1); for j=1:5 if maxh<x(i,j) maxh=x(i,j); end end maxh end BÀI THÍ NGHIỆM 2 FOURIER SERIES & FOURIER TRANSFORMS (CHUỖI FOURIER VÀ BIẾN ĐỔI FOURIER) I. TÓM TẮT LÝ THUYẾT 1. Chuỗi Fourier Bất kỳ tín hiệu định kỳ x(t) với chu kỳ T0 đều có thể biểu diễn dưới dạng 𝑥(𝑡) =∑𝑥𝑛𝑒 −𝑗2𝜋𝑛𝑡 ∞ −∞ Các hệ số của chuỗi Fourier được tính bởi công thức: 𝑥𝑛 = 1 𝑇0 ∫ 𝑥(𝑡)𝑒−𝑗2𝜋𝑛𝑡/𝑇0𝑑𝑡 𝑇0 0 2. Một số lệnh conj(x): tính lượng liên hợp phức abs(x): tính biên độ angle(x): pha length(x): tính chiều dài của vectơ stem : vẽ chuỗi dữ liệu rời rạc plot : vẽ chuỗi dữ liệu liên tục II. BÀI TẬP Bài tập 1 : clc; t = 0:pi/32:pi; x = [exp(-t/2)]; subplot(4,1,1); plot(t,x) title('Tin hieu goc') Xn = fft(x,32)/32; subplot(4,1,2); Xn = [conj(Xn(32:-1:2)), Xn]; stem(Xn) title('Lay doi xung qua truc tung') Xnmag = abs(Xn); Xnangle = angle(Xn); k=-32/2+1:32/2-1 subplot(4,1,3); stem(k, Xnmag(32/2+1:length(Xn)-32/2)); title('Bien do trong mien tan so') subplot(4,1,4); stem(k,Xnangle(32/2+1:length(Xn)-32/2)); title('Pha') Bài tập 2.1 : clc;clear all; t=-5:0.1:5; for k=1:length(t); if ((-1<t(k))&&(t(k)<=0)) x(k)=t(k)+1; else if ((t(k)>0)&&(t(k)<=1)) x(k)=1; else x(k)=0; end end end fs=100; subplot(3,2,1) plot(t,x); title('Tin hieu goc x1 trong mien thoi gian') X = fft(x); X = fftshift(X); Xmag = abs(X); Xmag = Xmag/max(Xmag); Xangle = angle(X); Xangle = Xangle/max(Xangle); F = [-length(X)/2:(length(X)/2)-1]*fs/length(X); subplot(3,2,3) stem(F, Xmag); title('Bien do x1 trong mien tan so') subplot(3,2,5) stem(F, Xangle); title('Pha x1 trong mien tan so') for k=1:length(t); if ((0<t(k))&&(t(k)<=1)) x2(k)=t(k); else if ((t(k)>1)&&(t(k)<=2)) x2(k)=1; else x2(k)=0; end end end subplot(3,2,2) plot(t,x2); title('Tin hieu goc x2 trong mien thoi gian') X2=fft(x2) X2 = fftshift(X2); Xmag2 = abs(X2); Xmag2 = Xmag2/max(Xmag2); Xangle2 = angle(X2); Xangle2 = Xangle2/max(Xangle2); F2 = [-length(X2)/2:(length(X2)/2)-1]*fs/length(X2); subplot(3,2,4) stem(F2, Xmag2); title('Bien do x2 trong mien tan so') subplot(3,2,6) stem(F2, Xangle2); title('Pha x2 trong mien tan so') Bài tập 2.2 clc;clear all; t=-5:0.1:5; for k=1:length(t); if ((-3<+t(k))&&(t(k)<=3)) x(k)=1; else x(k)=0; end end fs=100; subplot(3,2,1) plot(t,x); title('Tin hieu goc x1 trong mien thoi gian') X = fft(x); X = fftshift(X); Xmag = abs(X); Xmag = Xmag/max(Xmag); Xangle = angle(X); Xangle = Xangle/max(Xangle); F = [-length(X)/2:(length(X)/2)-1]*fs/length(X); subplot(3,2,3) stem(F, Xmag); title('Bien do x1 trong mien tan so') subplot(3,2,5) stem(F, Xangle); title('Pha x1 trong mien tan so') for k=1:length(t); if ((-1<=t(k))&&(t(k)<=1)) x2(k)=1; else x2(k)=0; end end subplot(3,2,2) plot(t,x2); title('Tin hieu goc x2 trong mien thoi gian') X2=fft(x2); X2 = fftshift(X2); Xmag2 = abs(X2); Xmag2 = Xmag2/max(Xmag2); Xangle2 = angle(X2); Xangle2 = Xangle2/max(Xangle2); F2 = [-length(X2)/2:(length(X2)/2)-1]*fs/length(X2); subplot(3,2,4) stem(F2, Xmag2); title('Bien do x2 trong mien tan so') subplot(3,2,6) stem(F2, Xangle2); title('Pha x2 trong mien tan so') BÀI THÍ NGHIỆM 3 FILTERING PERIODIC SIGNALS (LỌC TÍN HIỆU) I. TÓM TẮT LÝ THUYẾT 1. Một số lệnh fft : biến đổi Fourier thuận (tức là miền thời gian sang miền tần số) ifft : biến đổi Fourier ngược (tức là biến đổi ngược lại, miền tần số qua miền thời gian). fftshift: dịch chuyển tần số zero vào trung tâm phổ. sound : âm thanh nghe chưa định tỉ lệ (clip đầu vào [-1,1] soundsc: âm thanh thu nhỏ (quy mô/ bình thường hóa đầu vào [-1,1] wavread: đọc, định dạng file WAV, tỉ lệ lấy mẫu của file WAV cũng có thể được lấy ra bằng câu lệnh như sau:[x,Fs] = wavread('filename.wav') (Với x là vector âm thanh và Fs là tần số lấy mẫu). load: tải file 2. Chuyển đổi tín hiệu trong miền tần số và ngược lại - Trong Matlab thì biến đổi Fourier tín hiệu liên tục miền thời gian không thể thực hiện một cách chính xác. Những sự khác biệt quan trọng giữa biến đổi Fourier rời rạc trong máy tính và biến đổi Fourier liên tục là : dải tần số hữu hạn và các tần số mẫu rời rạc. Độ mịn dải tần số của tín hiệu sau khi biến đổi sẽ phụ thuộc vào tần số lấy mẫu của tín hiệu. Tần số lấy mẫu càng lớn thì phổ càng mịn. - Ta xét ví dụ trong ở trong Lap 3: Tín hiệu có tần số lấy mẫu là Fs = 8000. Có nghĩa là dải tần số từ [-4000,4000]Hz. Độ phân giải tần số phụ thuộc vào độ dài của các tín hiệu (mà cũng là chiều dài của biểu diễn tần số). - Vd: >>load fall %Tải file >>x = fall; >>X = fft(x) - Lệnh fft cho ra tần số không nằm ở trung tâm.Để có thể hiển thị kết quả trong miền tần số : một là hiển thị một nửa kết quả ( phần tần số dương), hai là hiển thị cả 2 vùng tần số bằng lệnh fftshift. Lệnh fftshift để dịch phổ vào trung tâm, đối xứng tại 0. Ta sử dụng abs để hiển thị biên độ của phổ tín hiệu. Lưu ý rằng vì biến đổi Fourier của tín hiệu là giá trị phức, abs được sử dụng để xem độ lớn (biên độ) của phổ. - Đoạn code sau thực hiện cả 2 cách: >> N = length(x); %Xác định kích thước của x >>pfreq = [0:N/2]* Fs/N; %Khoảng tần số >>Xpos = X(1:N/2+1) ; %Khoảng giá trị x chạy >>plot(pfreq,abs(Xpos)) ; %Vẽ đồ thị >>figure; >>freq = [-(N/2+1) :N/2]*Fs/N; %Tần số đồ thị >>plot(freq,abs(fftshift(X))); - Sau khi thay đổi tín hiệu trong miền tần số,lệnh ifft được dùng để chuyển về lại miền thời gian nếu mong muốn. >>xnew = real(ifft(X)); + Chú ý:Chúng ta cần dùng lệnh “real” vì biến đổi Fourier ngược cho ra giá trị hỗn hợp. 3. Bộ lọc thông thấp - Bộ lọc(filter):là một quá trình lọc những dữ liệu không cần thiết, loại bỏ đi những dữ liệu này và chỉ giữ lại những thứ mà mình muốn. - Bộ lọc thông thấp là bộ lọc chỉ cho phép tín hiệu hoạt động từ 0 Hz -> tần số cắt (cutoff) của nó, các tín hiệu tần số cao sẽ được lọc đi. - Một bộ lọc thông thấp lý tưởng loại bỏ các tần số cao hoàn toàn : - Một bộ lọc thông thấp thực tế thường có giá trị thấp nhưng không bằng 0 với H(w) ở tần số cao và dần dần giảm độ lớn khi tần số tăng. Bộ lọc thông thấp đơn giản nhất(ít hiệu quả nhất) được cho bởi: - Bộ lọc thông thấp có thể được dùng trong Matlab bằng cách sử dụng phép biến đổi Fourier. Phép nhân trong miền tần số là các hoạt động tương tựnhư chập trong miền thời gian. Nếu ta có cả hai tín hiệu và bộ lọc ở miền tần số thì ta có thể nhân chúng để tìm được đầu ra của hệ thống: - Ví dụ về cách sử dụng Matlab để thực hiện một bộ lọc thông thấp với tín hiệu vào là x,các tần số cắt của bộ lọc được biểu diễn bằng hằng số a. Ta cần xác định H trong miền tần số. >> load fall % tải tín hiệu >>x = fall; >>X = fft(x); % biến đổi tín hiệu x sang miền tần số >>N = length(X); >>a = 100*2*pi; >>w = (-N/2+1:(N/2)); % Vector tần số trung tâm (rad / s) >>H = a ./ (a + i*w); % Tạo lấy mẫu trung tâm của H >>plot(w*Fs/N,abs(H)) % w chuyển đổi trở lại Hz. - Đồ thị sẽ chỉ ra các dạng đáp ứng tần số của hệ thống mà ta đã sử dụng. >>Hshift = fftshift(H); %H không nằm ở trung tâm >> Y = X .* Hshift'; %Bộ lọc tín hiệu - Chú ý:Khi nhân hai vector cần chú ý kích thước của hai vector có bằng nhau hay chưa. - Ta sẽ có đầu ra của hệ thống trong miền tần số, dùng Fourier ngược để chuyển về miền thời gian. Nghe âm thanh gốc và sửa đổi và nghe sự khác biệt. Sử dụng tần số lấy mẫu Fs. >> y = real(ifft(Y)); >>sound(x, Fs) ) % âm thanh gốc >> sound(y, Fs) % âm thanh được lọc qua thông thấp - Bộ lọc làm giảm âm thanh của tín hiệu mà khi nghe chỉ có thể sử dụng lệnh sound mà không dùng được lệnh soundsc tự động chia tỉ lệ. - Chú ý: Đôi khi, bạn có thể khuếch đại tín hiệu để nó có chiều cao giống như tín hiệu ban đầu. >> y = y * (max(abs(x))/max(abs(y))); 4. Bộ lọc thông cao - Bộ lọc thông cao là bộ lọc chỉ cho phép tín hiệu hoạt động từ tần số cắt của nó đến vô cùng các tín hiệu tần số thấp nhỏ hơn tần số cắt sẽ được lọc đi. - Một bộ lọc thông cao lý tưởng loại bỏ hoàn toàn các tần số thấp. - Bộ lọc thông cao thực tế thường có giá trị thấp nhưng không bằng 0 với H(w) ở tần số thấp và dần dần tăng độ lớn khi tần số tăng. Mạch lọc thông cao được cho bởi: - Mạch lọc này có thể thực hiện giống như mạch lọc thông thấp. Với a là tần số cắt. II. BÀI TẬP Bài tập 1 : close all; clc; clear all; [x,Fs]=wavread('castanets44m.wav'); sound(x,Fs); X=fft(x); N=length(X); a=500*2*pi; w=(-N/2+1:(N/2))*Fs/N; %vector tần số trung tâm H=a./(a+i*w); %H nam o trung tam Hshift=fftshift(H); %H khong nam o trung tam Y=X.*Hshift'; % loc tin hieu y=real(ifft(Y)); subplot(2,1,1) plot(w,abs(fftshift(X))); title('Pho tin hieu') subplot(2,1,2) plot(w,abs(fftshift(Y))); title('Loc thong thap') Bài tập 2 : clc; clear all; Fs=8000; x = repmat([zeros(1, 99) 1], 1, 5) X=fft(x); N = length(X); a=20*2*pi; %tần số lọc w = (-N/2+1:(N/2))*Fs/N; H = a ./ (a + i*w); Hshift=fftshift(H); Y=X.*Hshift; y=real(ifft(Y)); y = y * (max(abs(x))/max(abs(y))) subplot(2,1,1) plot(x) title('Do thi x(t)') xlabel('t') ylabel('x(t)') subplot(2,1,2) plot(y) title('Do thi y(t)') xlabel('t') ylabel('y(t)') Bài tập 3 : clc; clear all; [x,Fs]=wavread('castanets44m.wav'); sound(x,Fs); X=fft(x); N=length(X); a=2000*2*pi; w=(-N/2+1:(N/2))*Fs/N; H=1-a./(a+i*w); Hshift=fftshift(H); Y=X.*Hshift'; y=real(ifft(Y)); subplot(2,1,1) plot(w,abs(fftshift(X))); title('Pho tin hieu') subplot(2,1,2) plot(w,abs(fftshift(Y))); title('Loc thong cao') sound(y,Fs) Bài tập 4 : close all;clc;clear all; [x,Fs]=wavread('mixed.wav'); %sound(x,Fs); X=fft(x); N=length(X); a=100*2*pi; b=2000*2*pi; w=(-N/2+1:(N/2))*Fs/N; L=a./(a+i*w); Lshift=fftshift(L); Y=X.*Lshift'; y=real(ifft(Y)); H=1-b./(b+i*w); Hshift=fftshift(H); Z=X.*Hshift'; z=real(ifft(Z)); subplot(3,1,1) plot(w,abs(fftshift(X))); title('Tin hieu goc') subplot(3,1,2) plot(w,abs(fftshift(Y))); title('Loc thong thap') subplot(3,1,3) plot(w,abs(fftshift(Z))); title('Loc thong cao') Bài tập 5 : close all;clc;clear all; load shake;load rainstick; x=shake+ 10*rainstick; Fs=8000; X=fft(x); N=length(X); a=100*2*pi; b=2000*2*pi; w=(-N/2+1:(N/2))*Fs/N; L=a./(a+i*w); Lshift=fftshift(L); Y=X.*Lshift'; y=real(ifft(Y)); H=1-b./(b+i*w); Hshift=fftshift(H); Z=X.*Hshift'; z=real(ifft(Z)); subplot(3,1,1) plot(w,abs(fftshift(X))); title('Tin hieu vao') subplot(3,1,2) plot(w,abs(fftshift(Y))); title('Loc thong thap') subplot(3,1,3) plot(w,abs(fftshift(Z))); title('Loc thong cao') BÀI THÍ NGHIỆM 4 INVESTIGATION OF ALIASING EFFECTS (LẤY MẪU) Bài tập 1 : clc;clear all;close all; f=1000; f1=1000; t=0:1/(10000):20/f; x=2*cos(2*pi*f1*t); subplot(2,1,1) plot(t,x) title('f=1000') X = fft(x); F = [-length(X)/2:(length(X)/2)-1]*1000/length(X); subplot(2,1,2) plot(F, abs(fftshift(X))); figure(2) f2=8000; t=0:1/(10000):20/f; x=cos(2*pi*f2*t); subplot(2,1,1) plot(t,x) title('f=8000') X = fft(x); F = [-length(X)/2:(length(X)/2)-1]*1000/length(X); subplot(2,1,2) plot(F, abs(fftshift(X))); figure(3) f3=500; t=0:1/(10000):20/f; x=2*cos(2*pi*f3*t); subplot(2,1,1) plot(t,x) title('f=500') X = fft(x); F = [-length(X)/2:(length(X)/2)-1]*1000/length(X); subplot(2,1,2) plot(F, abs(fftshift(X))); figure(4) f4=200; t=0:1/(10000):20/f; x=2*cos(2*pi*f4*t); subplot(2,1,1) plot(t,x) title('f=200') X = fft(x); F = [-length(X)/2:(length(X)/2)-1]*1000/length(X); subplot(2,1,2) plot(F, abs(fftshift(X)));