1. 程式人生 > >在matlab中使用不同的窗函式構造FIR數字低通濾波器

在matlab中使用不同的窗函式構造FIR數字低通濾波器

目的:
1.通過窗函式法設計FIR濾波器
2.使用窗函式法設計FIR濾波器,瞭解窗函式的形式和長度對濾波器效能的影響。

原理:
ex5.1

實現環境:
Matlab

理想低通濾波器設計:

function my_output=ideallp(wc,N)
%Ideal Lowpass filter computation
%------------------------------------
%[hd]=ideal_lp(wc,M)
% hd=ideal impulse response between 0 to M-1
% wc=cutoff frequency in radians
% M=length of the ideal filter % alpha=(N-1)/2; n=0:1:(N-1); m=n-alpha+eps; my_output=sin(wc*m)./(pi*m); end
  1. Boxcar窗
wp=0.2*pi;                          
ws=0.3*pi;                         
deltaw=ws-wp;                     
N=ceil(1.8*pi/deltaw);              
wc=(wp+ws)/2;                       
hd
=ideallp(wc,N); wd1=boxcar(N)'; h1=hd.*wd1; [H1,w]=freqz(h1,1); plot(w/pi,20*log10(abs(H1))); title('矩形窗');

得到頻率響應曲線如下:

2.Kaiser窗
wp=0.2*pi;
ws=0.5*pi;
As=50;                              
deltaf=(ws-wp)/(2*pi);
N
=ceil((As-7.95)/(14.36*deltaf))+1; % beta=0.1102*(As-8.7); wdkai=kaiser(N,beta)'; wc=(ws+wp)/2; hd=ideallp(wc,N); h=hd.*wdkai; [H,w]=freqz(h,1); plot(w/pi,20*log10(abs(H))); title('凱澤窗');

得到頻率響應影象如下:

3.triang窗
wp=0.2*pi;                          %
ws=0.3*pi;
As=50;
deltaw=ws-wp;
N=ceil(6.1*pi/deltaw);
wdtri=triang(N)';
wc=(wp+ws)/2;
hd=ideallp(wc,N);
h=hd.*wdtri;
[H,w]=freqz(h,1);
plot(w/pi,20*log10(abs(H)));
title('Triang');

得到頻率響應影象如下:

4.hanning窗 wp=0.2*pi; % ws=0.3*pi; As=50; deltaw=ws-wp; N=ceil(6.2*pi/deltaw); wdhann=hanning(N)’; wc=(wp+ws)/2; hd=ideallp(wc,N); h=hd.*wdhann; [H,w]=freqz(h,1); plot(w/pi,20*log10(abs(H))); title(‘Hanning’); 得到頻率響應影象如下: 5.hamming窗
wp=0.2*pi;                          %
ws=0.3*pi;
As=50;
deltaw=ws-wp;
N=ceil(6.6*pi/deltaw);
wdhamm=hamming(N)';
wc=(wp+ws)/2;
hd=ideallp(wc,N);
h=hd.*wdhamm;
[H,w]=freqz(h,1);
plot(w/pi,20*log10(abs(H)));
title('Hamming');

得到頻率響應影象如下:

6.blackman窗
wp=0.2*pi;                          %
ws=0.3*pi;
As=50;
deltaw=ws-wp;
N=ceil(6.1*pi/deltaw);
wdblack=blackman(N)';
wc=(wp+ws)/2;
hd=ideallp(wc,N);
h=hd.*wdblack;
[H,w]=freqz(h,1);
plot(w/pi,20*log10(abs(H)));
title('Blackman');

得到頻率響應影象如下:


一、 主瓣的寬度由N決定,旁瓣的衰減由窗函式的形狀決定。
二、 Kaiser的不同在於可以同時調整主瓣寬度和旁瓣衰減,這是其他窗函式不具備的。

三、 實際濾波的應用
設定初始訊號為10Hz的餘弦與40Hz餘弦疊加而成的訊號,分別採用以上濾波方法濾波。

  1. boxcar
%初始訊號
t=0:0.01:2;
y=cos(2*pi*10*t)+cos(2*pi*40*t);
subplot(2,1,1);
plot(t,y);
title('initial signal');
yfft=fftshift(fft(y));
w=-50:0.5:50;
subplot(2,1,2);
plot(w,abs(yfft));
axis([0 50 0 100]);
title('frequence');
 %濾波器內容
wp=0.2*pi;                          
ws=0.3*pi;                         
deltaw=ws-wp;                      
N=ceil(1.8*pi/deltaw);            
wc=(wp+ws)/2;                       
hd=ideallp(wc,N);                 
wd1=boxcar(N)';                   
h1=hd.*wd1;
final=fftfilt(h1,y,200);
  %濾波之後的訊號
figure;
subplot(2,1,1);
plot(final);
title('processed signal');
ynfft=fftshift(fft(final));
subplot(2,1,2);
plot(w,abs(ynfft));
axis([0 50 0 100]);
title('frequence filted');

執行結果如下:
初始訊號如下


處理後訊號如下:


可見40Hz分頻成分明顯降低,而10Hz分頻保留較好。
同理,更改程式碼中的濾波器部分分別得到下面幾個濾波結果對比
2. Kaiser

t=0:0.01:2;
y=cos(2*pi*10*t)+cos(2*pi*40*t);
w=-50:0.5:50;

wp=0.2*pi;
ws=0.5*pi;
As=50;                              %
deltaf=(ws-wp)/(2*pi);
N=ceil((As-7.95)/(14.36*deltaf))+1; %
beta=0.1102*(As-8.7); 
wdkai=kaiser(N,beta)';
wc=(ws+wp)/2;
hd=ideallp(wc,N);
h=hd.*wdkai;

final=fftfilt(h1,y,200);
subplot(2,1,1);
plot(final);
title('processed signal');
ynfft=fftshift(fft(final));
subplot(2,1,2);
plot(w,abs(ynfft));
axis([0 50 0 100]);
title('frequence filted');
3. Triang
t=0:0.01:2;
y=cos(2*pi*10*t)+cos(2*pi*40*t);
w=-50:0.5:50;

wp=0.2*pi;                          %
ws=0.3*pi;
As=50;
deltaw=ws-wp;
N=ceil(6.1*pi/deltaw);
wdtri=triang(N)';
wc=(wp+ws)/2;
hd=ideallp(wc,N);
h=hd.*wdtri;

final=fftfilt(h1,y,200);
subplot(2,1,1);
plot(final);
title('processed signal');
ynfft=fftshift(fft(final));
subplot(2,1,2);
plot(w,abs(ynfft));
axis([0 50 0 100]);
title('frequence filted');

  1. Hanning
t=0:0.01:2;
y=cos(2*pi*10*t)+cos(2*pi*40*t);
w=-50:0.5:50;

wp=0.2*pi;                          %
ws=0.3*pi;
As=50;
deltaw=ws-wp;
N=ceil(6.2*pi/deltaw);
wdhann=hanning(N)';
wc=(wp+ws)/2;
hd=ideallp(wc,N);
h=hd.*wdhann;

final=fftfilt(h1,y,200);
subplot(2,1,1);
plot(final);
title('processed signal');
ynfft=fftshift(fft(final));
subplot(2,1,2);
plot(w,abs(ynfft));
axis([0 50 0 100]);
title('frequence filted');
5. Hamming
t=0:0.01:2;
y=cos(2*pi*10*t)+cos(2*pi*40*t);
w=-50:0.5:50;

wp=0.2*pi;                          %
ws=0.3*pi;
As=50;
deltaw=ws-wp;
N=ceil(6.6*pi/deltaw);
wdhamm=hamming(N)';
wc=(wp+ws)/2;
hd=ideallp(wc,N);
h=hd.*wdhamm;

final=fftfilt(h1,y,200);
subplot(2,1,1);
plot(final);
title('processed signal');
ynfft=fftshift(fft(final));
subplot(2,1,2);
plot(w,abs(ynfft));
axis([0 50 0 100]);
title('frequence filted');
6. Blackman
t=0:0.01:2;
y=cos(2*pi*10*t)+cos(2*pi*40*t);
w=-50:0.5:50;

wp=0.2*pi;                          %
ws=0.3*pi;
As=50;
deltaw=ws-wp;
N=ceil(6.1*pi/deltaw);
wdblack=blackman(N)';
wc=(wp+ws)/2;
hd=ideallp(wc,N);
h=hd.*wdblack;

final=fftfilt(h1,y,200);
subplot(2,1,1);
plot(final);
title('processed signal');
ynfft=fftshift(fft(final));
subplot(2,1,2);
plot(w,abs(ynfft));
axis([0 50 0 100]);
title('frequence filted');

不同濾波方式衰減過程不同,但由於此次訊號只有兩個分頻,體現的不太明顯。