MATLAB cơ bản 6: FILTERING PERIODIC SIGNALS

[Lu ROm]

Các câu lệnh và nguồn dữ liệu của MATLAB:
+ Một số câu lệnh dùng trong lap :
- fft : biến đổi Fourier thuận (tức là miền thời gian sang miền tần số)
- iff : 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).
- tffshift: 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).
3.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 = 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.
4.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+1N/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 * (maabs(x))/maabs(y)));

Ví dụ 1:
Code:

clc;
clear all;
close all;
[x,Fs]=wavread('castanets44m.wav'); % load file am thanh
X=fft(x); %bien doi Fourier khong nam o trung tam
N=length(x); %Xac dinh kich thuoc cua x
a = 500*2*pi; %tan so cat
w1 = (-N/2+1:(N/2)); % Vector tan so trung tam
w = w1.*Fs/N; % LAY MOT TAN SO MAU
H = a./(a + 1i*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));
sound(x,Fs); % am thanh goc
sound(y,Fs); % am thanh sau khi qua bo loc thong thap
subplot(2,1,1);
plot(w,abs(fftshift(X))) % abs lay bien do cua so phuc, dich tan so ve trung tam
title('Tin hieu ban dau');
subplot(2,1,2);
plot(w,abs(fftshift(Y)))
title(Tin hieu loc');

Mô phỏng:

Ví dụ 2
Code:

clc;
clear all;
close all;
x = repmat([zeros(1,99) 1], 1, 5);
X = fft(x); %bien doi Fourier khong nam o trung tam
N = length(x); %Xac dinh kich thuoc cua x
a = 20; %tan so cat
w = (-N/2+1:(N/2)); % Vector tan so trung tam
H = a./(a + 1i*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)); %bien doi Fourier nguoc 
subplot(4,1,1)
plot(w,abs(fftshift(x)))
title('Tin hieu ban dau trong mien thoi gian');
subplot(4,1,2)
plot(w,abs(fftshift(X)))
title('Tin hieu ban dau trong mien tan so');
subplot(4,1,3)
plot(w,abs(fftshift(Y)))
title('Tin hieu loc trong mien tan so');
subplot(4,1,4)
plot(w,abs(fftshift(y)))
title('Tin hieu loc trong mien thoi gian');

5,High -pass Filtering :
- 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.

- Ví dụ 3:
Code:

%lap 3.3
clc;
clear all;
close all;
[x,Fs] = wavread('castanets44m.wav'); % load file am thanh
X = fft(x); %bien doi Fourier khong nam o trung tam
N = length(x); %Xac dinh kich thuoc cua x
a = 2000*2*pi; %tan so cat
w0 =(-N/2+1:(N/2)); % Vector tan so trung tam
w = w0.*Fs/N; % LAY MOT TAN SO MAU
H = 1-(a./(a + (1i*w)));  %H nam o trung tam
Hshift = fftshift(H);
Y = X .*Hshift' ;  % loc tin hieu, chuyen vi de nhan 2 ma tran
y = real(ifft(Y));
sound(x,Fs) % am thanh goc
sound(y,Fs) %am thanh sau khi qua bo loc thong thap
subplot(2,1,1);
plot(w,abs(fftshift(X))) % neu so phuc thi lay do lon cua no
% abs lay bien do cua so phuc, dich tan so ve trung tam
title('Tin hieu ban dau');
subplot(2,1,2);
plot(w,abs(fftshift(Y)))
title('Tin hieu loc');

Mô phỏng:

- Ví dụ 4:
Code:

%lap 3.4
clc;
clear all;
close all;
[x,Fs] = wavread('mixed.wav'); % load file am thanh
X = fft(x); %bien doi Fourier khong nam o trung tam
N = length(x); %Xac dinh kich thuoc cua x
% Thong cao lay tin hieu ken
a3 = 2000; %tan so cat
w3 =(-N/2+1:(N/2)); % Vector tan so trung tam
w4 = w3.*Fs/N; % LAY MOT TAN SO MAU
H = 1-(a3./(a3 + (1i*w4)));  %H nam o trung tâm
Hshift3 = fftshift(H);
a5 = 5000; %tan so cat
w5 =(-N/2+1:(N/2)); % Vector tan so trung tam
w6 = w5.*Fs/N; % LAY MOT TAN SO MAU
H5 = 1-(a5./(a5 + (1i*w6)));  %H nam o trung tâm
Hshift5 = fftshift(H5);
a = 6000; %tan so cat
w0 =(-N/2+1:(N/2)); % Vector tan so trung tam
w = w0.*Fs/N; % LAY MOT TAN SO MAU
H2 = 1-(a./(a + (1i*w)));
Hshift = fftshift(H2);
Y = X .*Hshift5' .* Hshift3' .*Hshift' ;  % loc tin hieu, chuyen vi de nhan 2 ma tran
y = real(ifft(Y));
y = y*(max(abs(x))/ max(abs(y))); % bang tin hieu ban dau
%thong thap lay tin hieu bass
a1 = 100; %tan so cat
w1 = (-N/2+1:(N/2)); % Vector tan so trung tam
w2 = w1.*Fs/N; % LAY MOT TAN SO MAU
H1 = a1./(a1 + 1i*w2);  %H nam o trung tam
Hshift1 = fftshift(H1); %H khong nam o trung tam
a10 = 50; %tan so cat
w10 = (-N/2+1:(N/2)); % Vector tan so trung tam
w20 = w10.*Fs/N; % LAY MOT TAN SO MAU
H10 = a10./(a10 + 1i*w20);  %H nam o trung tam
Hshift10 = fftshift(H10); %H khong nam o trung tam
a12 = 5; %tan so cat
w12 = (-N/2+1:(N/2)); % Vector tan so trung tam
w22 = w12.*Fs/N; % LAY MOT TAN SO MAU
H12 = a10./(a10 + 1i*w22);  %H nam o trung tam
Hshift12 = fftshift(H12); %H khong nam o trung tam
Y1 = X .*Hshift1' .*Hshift10' .*Hshift12';  % loc tin hieu
y1 = real(ifft(Y1));
y1 = y1*(max(abs(x))/ max(abs(y1))); % bang tin hieu ban dau
sound(x,Fs) % am thanh goc
sound(y1,Fs) %am thanh sau khi qua bo loc thong thap
sound(y,Fs) %am thanh sau khi qua bo loc thong cao
subplot(3,1,1);
plot(w5,abs(fftshift(X))) % neu so phuc thi lay do lon cua no
title('tin hieu ban dau');
subplot(3,1,2);
plot(w5,abs(fftshift(Y)))
title('tin hieu khi loc thong cao');
subplot(3,1,3);
plot(w2,abs(fftshift(Y1)))
title(' tin hieu khi loc thong thap');

Mô phỏng:

Ví dụ 5
Code:

%lap 3.5
clc;
clear all;
%load file co san
load rainstick;
load shake;
mixedsig = shake + 10*rainstick
%bieu dien tin hieu sau tron,
x = mixedsig;
a = 1000*pi;
Fs = 8000;
X = fft(x);
N = length(x);
w0 =(-N/2+1:(N/2)); % Vector tan so trung tam
w = w0.*Fs/N;
%LOC THONG CAO
H = 1-(a./(a + 1i*w));  %H nam o trung tam
Hshift = fftshift(H);
Y = X .* Hshift';
y = real(ifft(Y));
% LOC THONG THAP
a1 = 500;
w2 =(-N/2+1:(N/2)); % Vector tan so trung tam
w1 = w2.*Fs/N;
H1 = a1./(a1 + 1i*w1);  %H nam o trung tam
Hshift1 = fftshift(H1); %H khong nam o trung tam
Y1 = X .*Hshift1' ;  % loc tin hieu
y1 = real(ifft(Y1));
%bieu do
subplot (3,1,1)
plot(w,abs(fftshift(X))) % neu so phuc thi lay do lon cua no,% abs lay bien do cua so phuc, dich tan so ve trung tam
title('TIN HIEU BAN DAU');
subplot (3,1,2)
plot(w,abs(fftshift(Y)))
title('TIN HIEU LOC THONG CAO');
subplot(3,1,3)
plot(w,abs(fftshift(Y1)))
title('TIN HIEU LOC THONG THAP');

Mô phỏng:

++ Nhận xét:
Sau khi nghe lại đoạn nhạc file WAV thì ta nhận thấy:
- Đối với bộ lọc thông thấp, tín hiệu phổ ra có biên độ giảm, mạch lọc chọn các tín hiệu có tần số thấp, âm thanh nghe êm tai.
- Với bộ lọc thông cao, tín hiệu phổ ra có biên độ tăng, mạch lọc chọn các tín hiệu có tần số cao, âm thanh nghe trong và vang.
- Muốn lọc tín hiệu hiệu quả tên nên dùng nhiều bộ lọc cho tín hiệu đó.