Thiết kế bộ lọc FIR thông thấp phương pháp cửa sổ (matlab code)
Thiết kế bộ lọc số FIR sử dụng cửa sổ hamming với các thông số cho dưới đây :
wp = 0.2 pi, Rp = 0.25 dB
ws = 0.3 pi, As = 50 dB.
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)))
% vẽ
subplot(1,1,1)
subplot(2,2,1); stem(n,hd); title('Đáp ứng xung lý tưởng')
axis([0 M-1 -0.1 0.3]); xlabel('n'); ylabel('hd(n)')
subplot(2,2,2); stem(n,w_ham);title('Cửa sổ Hamming')
axis([0 M-1 0 1.1]); xlabel('n'); ylabel('w(n)')
subplot(2,2,3); stem(n,h);title('Đáp ứng xung thực tế')
axis([0 M-1 -0.1 0.3]); xlabel('n'); ylabel('h(n)')
subplot(2,2,4); plot(w/pi,db);title('Đáp ứng biên độ theo dB');grid
axis([0 1 -100 10]); xlabel('frequency in pi units'); ylabel('Decibels')
set(gca,'XTickMode','manual','XTick',[0,0.2,0.3,1])
set(gca,'YTickMode','manual','YTick',[-50,0])
set(gca,'YTickLabelMode','manual','YTickLabels',['50';' 0'])
Dùng cửa sổ Kaiser
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)))
% vẽ
subplot(1,1,1)
subplot(2,2,1); stem(n,hd); title('Đáp ứng xung lý tưởng')
axis([0 M-1 -0.1 0.3]); xlabel('n'); ylabel('hd(n)')
subplot(2,2,2); stem(n,w_ham);title('Cửa sổ Hamming')
axis([0 M-1 0 1.1]); xlabel('n'); ylabel('w(n)')
subplot(2,2,3); stem(n,h);title('Đáp ứng xung thực tế')
axis([0 M-1 -0.1 0.3]); xlabel('n'); ylabel('h(n)')
subplot(2,2,4); plot(w/pi,db);title('Đáp ứng biên độ theo dB');grid
axis([0 1 -100 10]); xlabel('frequency in pi units'); ylabel('Decibels')
set(gca,'XTickMode','manual','XTick',[0,0.2,0.3,1])
set(gca,'YTickMode','manual','YTick',[-50,0])
set(gca,'YTickLabelMode','manual','YTickLabels',['50';' 0'])
wp = 0.2*pi; ws = 0.3*pi; As = 50;
tr_width = ws - wp;
M = ceil((As-7.95)/(14.36*tr_width/(2*pi))+1) + 1
n=[0:1:M-1];
beta = 0.1102*(As-8.7)
wc = (ws+wp)/2;
hd = ideal_lp(wc,M);
w_kai = (kaiser(M,beta))';
h = hd .* w_kai;
[db,mag,pha,grd,w] = freqz_m(h,[1]);
delta_w = 2*pi/1000;
As = -round(max(db(ws/delta_w+1:1:501)))
% Vẽ
subplot(1,1,1)
subplot(2,2,1); stem(n,hd); title('Đáp ứng xung lý tưởng')
axis([0 M-1 -0.1 0.3]); xlabel('n'); ylabel('hd(n)')
subplot(2,2,2); stem(n,w_kai);title('Cửa sổ Kaiser')
axis([0 M-1 0 1.1]); xlabel('n'); ylabel('w(n)')
subplot(2,2,3); stem(n,h);title('Đáp ứng xung thực tế')
axis([0 M-1 -0.1 0.3]); xlabel('n'); ylabel('h(n)')
subplot(2,2,4);plot(w/pi,db);title('Đáp ứng biên độ theo dB');grid
axis([0 1 -100 10]); xlabel('frequency in pi units'); ylabel('Decibels')
set(gca,'XTickMode','manual','XTick',[0,0.2,0.3,1])
set(gca,'YTickMode','manual','YTick',[-50,0])
set(gca,'YTickLabelMode','manual','YTickLabels',['50';' 0'])
tr_width = ws - wp;
M = ceil((As-7.95)/(14.36*tr_width/(2*pi))+1) + 1
n=[0:1:M-1];
beta = 0.1102*(As-8.7)
wc = (ws+wp)/2;
hd = ideal_lp(wc,M);
w_kai = (kaiser(M,beta))';
h = hd .* w_kai;
[db,mag,pha,grd,w] = freqz_m(h,[1]);
delta_w = 2*pi/1000;
As = -round(max(db(ws/delta_w+1:1:501)))
% Vẽ
subplot(1,1,1)
subplot(2,2,1); stem(n,hd); title('Đáp ứng xung lý tưởng')
axis([0 M-1 -0.1 0.3]); xlabel('n'); ylabel('hd(n)')
subplot(2,2,2); stem(n,w_kai);title('Cửa sổ Kaiser')
axis([0 M-1 0 1.1]); xlabel('n'); ylabel('w(n)')
subplot(2,2,3); stem(n,h);title('Đáp ứng xung thực tế')
axis([0 M-1 -0.1 0.3]); xlabel('n'); ylabel('h(n)')
subplot(2,2,4);plot(w/pi,db);title('Đáp ứng biên độ theo dB');grid
axis([0 1 -100 10]); xlabel('frequency in pi units'); ylabel('Decibels')
set(gca,'XTickMode','manual','XTick',[0,0.2,0.3,1])
set(gca,'YTickMode','manual','YTick',[-50,0])
set(gca,'YTickLabelMode','manual','YTickLabels',['50';' 0'])
e viết code giống vậy nhưng matlap vẫn báo lỗi là
Trả lờiXóaUndefined function 'ideal_lp' for input arguments of type 'double'.
Error in kaiser2 (line 11)
hd= ideal_lp(wc,M);
là sao vậy add?
bạn xem thử đã thêm hàm ideal_lp() vào chưa?
Xóathêm hàm ideal_lp() tức là như nào anh ơi
Xóabạn tạo function ideal_lp() như http://ratdongian.blogspot.com/2016/07/thiet-ke-bo-loc-fir-code-matlab.html, sau đó include vào code và sử dụng nha.
Xóacho em hỏi 3 dòng lệnh dưới này delta_w là gì sao lại có công thức delta_w = 2*pi/1000;
Trả lờiXóatại sao tính Rp và As lại có thể tính được khi lấy min và max ở mảng db vậy ạ..
delta_w = 2*pi/1000;
Rp = -(min(db(1:1:wp/delta_w+1)))
As = -round(max(db(ws/delta_w+1:1:501)))
hihi. Lâu quá Anh không đụng tới code matlab nên anh không nhớ sao lại như vậy nữa.
XóaNếu em biết lý do thì em post lại trên này để chia sẻ cho mọi người nhé.
Cảm ơn em!
subplot(1,1,1); plot(w/pi,db);title('dap ung bien do theo dB');grid
Trả lờiXóaaxis([0 1 -100 10]); xlabel('frequency in pi units'); ylabel('Decibels')
Error using stem (line 31)
X must be same length as Y.
Error in barlett (line 25)
subplot(1,1,1); stem(n,h);title('dap ung xung thuc te')
lỗi này thì sửa ntn ạ?