MATLAB窗函數法實現(xiàn)FIR的高通帶通和低通濾波器的程序要點

MATLAB課程設計報告學 院：地球物理與石油資源學院 班 級： 測井（基）11001 姓 名： 大牛啊啊啊 學 號： 班內編號： 指導教師： 陳義群 完成日期： 2013年6月3日 一、 題目 FIR濾波器的窗函數設計法及性能比較1. FIR濾波器簡介數字濾波器是一種用來過濾時間離散信號的數字系統(tǒng)，通過對抽樣數據進行數學處理來達到頻域濾波的目的。根據其單位沖激響應函數的時域特性可分為兩類：無限沖激響應（IIR）濾波器和有限沖激響應（FIR）濾波器。與IIR濾波器相比，F(xiàn)IR濾波器的主要特點為：a. 線性相位；b.非遞歸運算。2. FIR濾波器的設計FIR濾波器的設計方法主要有三種：a.窗函數設計法；b.頻率抽樣發(fā)；c.最小平法抽樣法；這里我主要討論在MATLAB環(huán)境下通過調用信號分析與處理工具箱的幾類窗函數來設計濾波器并分析與比較其性能。窗函數法設計FIR濾波器的一般步驟如下：a. 根據實際問題確定要設計的濾波器類型；b. 根據給定的技術指標，確定期望濾波器的理想頻率特性；c. 求期望濾波器的單位脈沖響應；d. 求數字濾波器的單位脈沖響應；e. 應用。常用的窗函數有 4. 常用窗函數的參數5. FIR濾波器的MATLAB實現(xiàn)方式在MATLAB信號分析與處理工具箱中提供了大量FIR窗函數的設計函數，本次用到主要有以下幾種：hanning(N) hanning窗函數的調用hamming(N) hamming窗函數的調用blackman(N) blackman窗函數的調用kaiser(n+1,beta) kaiser窗函數的調用kaiserord 計算kaiser窗函數的相關參數freqz 求取頻率響應filter 對信號進行濾波的函數6. 實驗具體步驟本次實驗分別通過調用hanning ,hamming ,Blackman,kaiser窗函數，給以相同的技術參數，來設計低通，帶通，高通濾波器，用上述窗函數的選擇標準來比較各種窗函數的優(yōu)劣，并給以一個簡諧波進行濾波處理，比較濾波前后的效果。達到綜合比較的效果。二、源代碼 1.利用hanning hamming blackman kaiser窗,設計一個低通FIRfunction lowpassfilterclc;clear all;Fs=100;%采樣頻率fp=20;%通帶截止頻率fs=30;%阻帶起始頻率wp=2*pi*fp/Fs;%將模擬通帶截止頻率轉換為數字濾波器頻率ws=2*pi*fs/Fs;%將模擬阻帶起始頻率轉換為數字濾波器頻率wn=(wp+ws)/2/pi;%標準化的截止頻率響應Bt=ws-wp;N0=ceil(6.2*pi/Bt);%濾波器長度N=N0+mod(N0+1,2);window1=hanning(N);%使用hanning窗函數window2=hamming(N);%使用hamming窗函數window3=blackman(N);%使用blackman窗函數n,Wn,beta,ftype=kaiserord(20 25,1 0,0.01 0.01,100);window4=kaiser(n+1,beta);%使用kaiser窗函數%設計加窗函數fir1b1=fir1(N-1,wn,window1);b2=fir1(N-1,wn,window2);b3=fir1(N-1,wn,window3);b4=fir1(n,Wn/pi,window4 ,'noscale');%求取頻率響應H1,W1=freqz(b1,1,512,2);H2,W2=freqz(b2,1,512,2);H3,W3=freqz(b3,1,512,2);H4,W4=freqz(b4,1,512,2);figure(1);subplot(2,2,1),plot(W1,20*log10(abs(H1);%繪制頻率響應圖形axis(0,1,-100,100);title('低通hanning窗的頻率響應圖形');xlabel('頻率（Hz)');ylabel('幅值');subplot(2,2,2),plot(W2,20*log10(abs(H2);%繪制頻率響應圖形axis(0,1,-100,100);title('低通hamming窗的頻率響應圖形');xlabel('頻率（Hz)');ylabel('幅值');subplot(2,2,3),plot(W3,20*log10(abs(H3);%繪制頻率響應圖形axis(0,1,-100,100);title('低通blackman窗的頻率響應圖形');xlabel('頻率（Hz)');ylabel('幅值');subplot(2,2,4),plot(W4,20*log10(abs(H4);%繪制頻率響應圖形axis(0,1,-100,100);title('低通kaiser窗的頻率響應圖形');xlabel('頻率（Hz)');ylabel('幅值');T=1/Fs;L=100;%信號長度t=(0:L-1)*T;%定義時間范圍和步長y=sin(2*pi*5*t)+5*sin(2*pi*15*t)+8*sin(2*pi*40*t);%濾波前的圖形NFFT = 2nextpow2(L); % Next power of 2 from length of yY = fft(y,NFFT)/L;%將時域信號變換到頻域f = Fs/2*linspace(0,1,NFFT/2+1);%頻域采樣figure(2);plot(f,2*abs(Y(1:NFFT/2+1);xlabel('frequency/Hz');ylabel('Amuplitude') ;%濾波前頻譜title('濾波前的頻譜');%濾波后頻譜%采用hanning窗濾波器yy1=filter(b1,1,y);%調用濾波函數YY1=fft(yy1,NFFT)/L;%進行傅里葉變換，下同。f1=Fs/2*linspace(0,1,NFFT/2+1);figure(3);subplot(2,2,1),plot(f1,2*abs(YY1(1:NFFT/2+1) ;xlabel('frequency/Hz');ylabel('Amuplitude');title('hanning窗的濾波效果');%采用hammning窗濾波器yy2=filter(b2,1,y);YY2=fft(yy2,NFFT)/L;f1=Fs/2*linspace(0,1,NFFT/2+1);subplot(2,2,2),plot(f1,2*abs(YY2(1:NFFT/2+1) ;xlabel('frequency/Hz');ylabel('Amuplitude');title('hamming窗的濾波效果');%采用blackman窗濾波器yy3=filter(b3,1,y);YY3=fft(yy3,NFFT)/L;f1=Fs/2*linspace(0,1,NFFT/2+1);subplot(2,2,3), plot(f1,2*abs(YY3(1:NFFT/2+1) ;xlabel('frequency/Hz');ylabel('Amuplitude');title('blackman窗的濾波效果');%采用kaiser窗濾波器yy4=filter(b4,1,y);YY4=fft(yy4,NFFT)/L;f1=Fs/2*linspace(0,1,NFFT/2+1);subplot(2,2,4),plot(f1,2*abs(YY4(1:NFFT/2+1) ;xlabel('frequency/Hz');ylabel('Amuplitude');xlabel('frequency/Hz');ylabel('Amuplitude');title('kaiser窗函數濾波效果');%濾波前后的信號的時域對比figure(4);plot(y);xlabel('時間/s');ylabel('振幅');title('濾波前振幅特性');figure(5);subplot(2,2,1),plot(yy1);xlabel('時間/s');ylabel('振幅');title('hanning窗函數濾波振幅特性');subplot(2,2,2),plot(yy2);xlabel('時間/s');ylabel('振幅');title('hamming窗函數濾波振幅特性');subplot(2,2,3),plot(yy3);xlabel('時間/s');ylabel('振幅');title('blackman窗函數濾波振幅特性');subplot(2,2,4),plot(yy4);xlabel('時間/s');ylabel('振幅');title('kaiser窗函數濾波振幅特性');%濾波前后的信號的相位對比figure(6);plot(angle(Y);xlabel('時間/s');ylabel('相位');title('濾波前的相位特性');figure(7);subplot(2,2,1),plot(angle(YY1);xlabel('時間/s');ylabel('相位');title('hanning窗函數濾波相位特性');subplot(2,2,2),plot(angle(YY2);xlabel('時間/s');ylabel('相位');title('hamming窗函數濾波相位特性');subplot(2,2,3),plot(angle(YY3);xlabel('時間/s');ylabel('相位');title('blackman窗函數濾波相位特性');subplot(2,2,4),plot(angle(YY4);xlabel('時間/s');ylabel('相位');title('kaiser窗函數濾波相位特性');2.設計一個hanning hamming blackman kaiser窗函數bandpass_FIR%設計一個hanning hamming blackman kaiser窗函數bandpass_FIRfunction bandpassfilterFs=100;%采樣頻率fp1=15;%通帶下限截止頻率fp2=20;%通帶上限截止頻率fs1=10;fs2=25;wp1=2*pi*fp1/Fs;%將通帶下限截止頻率轉換為數字濾波器頻率wp2=2*pi*fp2/Fs;%將通帶上限截止頻率轉換為數字濾波器頻率ws1=2*pi*fs1/Fs;%將通帶下限截止頻率轉換為數字濾波器頻率ws2=2*pi*fs2/Fs;%將通帶上限截止頻率轉換為數字濾波器頻率Bt=wp1-ws1;N0=ceil(6.2*pi/Bt);N=N0+mod(N0+1,2);wn=(wp1+ws1)/2/pi,(wp2+ws2)/2/pi;window1=hanning(N);%使用hanning窗函數window2=hamming(N);%使用hamming窗函數window3=blackman(N);%使用blackman窗函數%設過渡帶寬度為5Hzn,Wn,beta,ftype=kaiserord(10 15 20 25,0 1 0,0.01 0.01 0.01,100);%求階數n以及參數betawindow4=kaiser(n+1,beta);%使用kaiser窗函數%設計加窗函數fir1b1=fir1(N-1,wn,window1);b2=fir1(N-1,wn,window2);b3=fir1(N-1,wn,window3);b4=fir1(n,Wn,window4,'noscale');%求取頻率響應H1,W1=freqz(b1,1,512,2);H2,W2=freqz(b2,1,512,2);H3,W3=freqz(b3,1,512,2);H4,W4=freqz(b4,1,512,2);figure(1);subplot(2,2,1),plot(W1,20*log10(abs(H1);%繪制頻率響應圖形axis(0,1,-100,100);title('帶通hanning窗的頻率響應圖形');xlabel('頻率（Hz)');ylabel('幅值');subplot(2,2,2),plot(W2,20*log10(abs(H2);%繪制頻率響應圖形axis(0,1,-100,100);title('帶通hamming窗的頻率響應圖形');xlabel('頻率（Hz)');ylabel('幅值');subplot(2,2,3),plot(W3,20*log10(abs(H3);%繪制頻率響應圖形axis(0,1,-100,100);title('帶通blackman窗的頻率響應圖形');xlabel('頻率（Hz)');ylabel('幅值');subplot(2,2,4),plot(W4,20*log10(abs(H4);%繪制頻率響應圖形axis(0,1,-100,100);title('帶通kaiser窗的頻率響應圖形');xlabel('頻率（Hz)');ylabel('幅值');T=1/Fs;L=100;%信號長度t=(0:L-1)*T;%定義時間范圍和步長y=sin(2*pi*5*t)+5*sin(2*pi*15*t)+8*sin(2*pi*40*t);%濾波前的圖形NFFT = 2nextpow2(L); % Next power of 2 from length of yY = fft(y,NFFT)/L;%將時域信號變換到頻域f = Fs/2*linspace(0,1,NFFT/2+1);%頻域采樣figure(2);plot(f,2*abs(Y(1:NFFT/2+1);xlabel('frequency/Hz');ylabel('Amuplitude') ;%濾波前頻譜title('濾波前的頻譜');%濾波后頻譜%采用hanning窗濾波器yy1=filter(b1,1,y);%調用濾波函數YY1=fft(yy1,NFFT)/L;%進行傅里葉變換，下同。f1=Fs/2*linspace(0,1,NFFT/2+1);figure(3);subplot(2,2,1),plot(f1,2*abs(YY1(1:NFFT/2+1) ;xlabel('frequency/Hz');ylabel('Amuplitude');title('hanning窗的濾波效果');%采用hammning窗濾波器yy2=filter(b2,1,y);YY2=fft(yy2,NFFT)/L;f1=Fs/2*linspace(0,1,NFFT/2+1);subplot(2,2,2),plot(f1,2*abs(YY2(1:NFFT/2+1) ;xlabel('frequency/Hz');ylabel('Amuplitude');title('hamming窗的濾波效果');%采用blackman窗濾波器yy3=filter(b3,1,y);YY3=fft(yy3,NFFT)/L;f1=Fs/2*linspace(0,1,NFFT/2+1);subplot(2,2,3), plot(f1,2*abs(YY3(1:NFFT/2+1) ;xlabel('frequency/Hz');ylabel('Amuplitude');title('blackman窗的濾波效果');%采用kaiser窗濾波器yy4=filter(b4,1,y);YY4=fft(yy4,NFFT)/L;f1=Fs/2*linspace(0,1,NFFT/2+1);subplot(2,2,4),plot(f1,2*abs(YY4(1:NFFT/2+1) ;xlabel('frequency/Hz');ylabel('Amuplitude');xlabel('frequency/Hz');ylabel('Amuplitude');title('kaiser窗函數濾波效果');%濾波前后的信號的時域對比figure(4);plot(y);xlabel('時間/s');ylabel('振幅');title('濾波前振幅特性');figure(5);subplot(2,2,1),plot(yy1);xlabel('時間/s');ylabel('振幅');title('hanning窗函數濾波振幅特性');subplot(2,2,2),plot(yy2);xlabel('時間/s');ylabel('振幅');title('hamming窗函數濾波振幅特性');subplot(2,2,3),plot(yy3);xlabel('時間/s');ylabel('振幅');title('blackman窗函數濾波振幅特性');subplot(2,2,4),plot(yy4);xlabel('時間/s');ylabel('振幅');title('kaiser窗函數濾波振幅特性');%濾波前后的信號的相位對比figure(6);plot(angle(Y);xlabel('時間/s');ylabel('相位');title('濾波前的相位特性');figure(7);subplot(2,2,1),plot(angle(YY1);xlabel('時間/s');ylabel('相位');title('hanning窗函數濾波相位特性');subplot(2,2,2),plot(angle(YY2);xlabel('時間/s');ylabel('相位');title('hamming窗函數濾波相位特性');subplot(2,2,3),plot(angle(YY3);xlabel('時間/s');ylabel('相位');title('blackman窗函數濾波相位特性');subplot(2,2,4),plot(angle(YY4);xlabel('時間/s');ylabel('相位');title('kaiser窗函數濾波相位特性');3.分別設計hanning hamming blackman kaiser窗函數highpass_FIRfunction highpassfilterclc;clear all;Fs=100;%采樣頻率fs=35;%高通阻帶模擬截止頻率fp=40;%高通通帶模擬起始頻率ws=2*pi*fs/Fs;wp=2*pi*fp/Fs;wn=(wp+ws)/2/pi;Bt=wp-ws;N0=ceil(55*pi/Bt);N=N0+mod(N0+1,2);%調用窗函數window1=hanning(N);window2=hamming(N);window3=blackman(N);n,Wn,beta,ftype=kaiserord(35,40,0 1,0.01 0.01,100);window4=kaiser(n+1,beta);%設計加窗函數fir1b1=fir1(N-1,wn,'high',window1);b2=fir1(N-1,wn,'high',window2);b3=fir1(N-1,wn,'high',window3);b4=fir1(n,Wn,'high',window4 ,'noscale');%求取頻率響應H1,W1=freqz(b1,1,512,2);H2,W2=freqz(b2,1,512,2);H3,W3=freqz(b3,1,512,2);H4,W4=freqz(b4,1,512,2);figure(1);subplot(2,2,1),plot(W1,20*log10(abs(H1);%繪制頻率響應圖形axis(0,1,-100,100);title('高通hanning窗的頻率響應圖形');xlabel('頻率（Hz)');ylabel('幅值');subplot(2,2,2),plot(W2,20*log10(abs(H2);%繪制頻率響應圖形axis(0,1,-100,100);title('高通hamming窗的頻率響應圖形');xlabel('頻率（Hz)');ylabel('幅值');subplot(2,2,3),plot(W3,20*log10(abs(H3);%繪制頻率響應圖形axis(0,1,-100,100);title('高通blackman窗的頻率響應圖形');xlabel('頻率（Hz)');ylabel('幅值');subplot(2,2,4),plot(W4,20*log10(abs(H4);%繪制頻率響應圖形axis(0,1,-100,100);title(' 高通kaiser窗的頻率響應圖形');xlabel('頻率（Hz)');ylabel('幅值');T=1/Fs;L=100;%信號長度t=(0:L-1)*T;%定義時間范圍和步長y=sin(2*pi*5*t)+5*sin(2*pi*15*t)+8*sin(2*pi*40*t);%濾波前的圖形NFFT = 2nextpow2(L); % Next power of 2 from length of yY = fft(y,NFFT)/L;%將時域信號變換到頻域f = Fs/2*linspace(0,1,NFFT/2+1);%頻域采樣figure(2);plot(f,2*abs(Y(1:NFFT/2+1);xlabel('frequency/Hz');ylabel('Amuplitude') ;%濾波前頻譜title('濾波前的頻譜');%濾波后頻譜%采用hanning窗濾波器yy1=filter(b1,1,y);%調用濾波函數YY1=fft(yy1,NFFT)/L;%進行傅里葉變換，下同。f1=Fs/2*linspace(0,1,NFFT/2+1);figure(3);subplot(2,2,1),plot(f1,2*abs(YY1(1:NFFT/2+1) ;xlabel('frequency/Hz');ylabel('Amuplitude');title('hanning窗的濾波效果');%采用hammning窗濾波器yy2=filter(b2,1,y);YY2=fft(yy2,NFFT)/L;f1=Fs/2*linspace(0,1,NFFT/2+1);subplot(2,2,2),plot(f1,2*abs(YY2(1:NFFT/2+1) ;xlabel('frequency/Hz');ylabel('Amuplitude');title('hamming窗的濾波效果');%采用blackman窗濾波器yy3=filter(b3,1,y);YY3=fft(yy3,NFFT)/L;f1=Fs/2*linspace(0,1,NFFT/2+1);subplot(2,2,3), plot(f1,2*abs(YY3(1:NFFT/2+1) ;xlabel('frequency/Hz');ylabel('Amuplitude');title('blackman窗的濾波效果');%采用kaiser窗濾波器yy4=filter(b4,1,y);YY4=fft(yy4,NFFT)/L;f1=Fs/2*linspace(0,1,NFFT/2+1);subplot(2,2,4),plot(f1,2*abs(YY4(1:NFFT/2+1) ;xlabel('frequency/Hz');ylabel('Amuplitude');xlabel('frequency/Hz');ylabel('Amuplitude');title('kaiser窗函數濾波效果');%濾波前后的信號的時域對比figure(4);plot(y);xlabel('時間/s');ylabel('振幅');title('濾波前振幅特性');figure(5);subplot(2,2,1),plot(yy1);xlabel('時間/s');ylabel('振幅');title('hanning窗函數濾波振幅特性');subplot(2,2,2),plot(yy2);xlabel('時間/s');ylabel('振幅');title('hamming窗函數濾波振幅特性');subplot(2,2,3),plot(yy3);xlabel('時間/s');ylabel('振幅');title('blackman窗函數濾波振幅特性');subplot(2,2,4),plot(yy4);xlabel('時間/s');ylabel('振幅');title('kaiser窗函數濾波振幅特性');%濾波前后的信號的相位對比figure(6);plot(angle(Y);xlabel('時間/s');ylabel('相位');title('濾波前的相位特性');figure(7);subplot(2,2,1),plot(angle(YY1);xlabel('時間/s');ylabel('相位');title('hanning窗函數濾波相位特性');subplot(2,2,2),plot(angle(YY2);xlabel('時間/s');ylabel('相位');title('hamming窗函數濾波相位特性');subplot(2,2,3),plot(angle(YY3);xlabel('時間/s');ylabel('相位');title('blackman窗函數濾波相位特性');subplot(2,2,4),plot(angle(YY4);xlabel('時間/s');ylabel('相位');title('kaiser窗函數濾波相位特性');三、運行結果1.給定的簡諧信號： 圖一：輸入簡諧信號濾波前的頻譜圖二：輸入簡諧信號濾波前的振幅圖三：輸入簡諧信號濾波前的相位2.低通濾波器的設計低通濾波器的技術指標：采樣頻率 Fs=100Hz；通帶截止頻率 fp=20Hz； 阻帶起始頻率fs=30HzHanning Hamming Blackman Kaiser采用相同的技術指標。以下即是四個窗函數的頻響圖及對簡諧信號濾波后的效果圖。 圖四 不同低通窗函數低通濾波器的歸一化頻響圖 圖五 不同窗函數低通濾波器對信號的濾波后頻率域效果圖六 不同窗函數低通濾波器對信號的濾波后時間域效果圖七 不同窗函數低通濾波器對信號的濾波后相位變化由以上濾波后頻率，相位，振幅變化以觀察到：hanning窗 hamming窗及blackman窗的濾波效果基本相當，但三者相比：hamming窗的過渡帶衰減最快，blackman窗旁瓣幅度最小。而kaiser窗只有5Hz信號，15Hz信號被截斷，與設計要求有出入。3. 帶通濾波器的設計采樣頻率 Fs=100Hz；阻帶截止頻率1：fs1=10Hz；通帶起始頻率1：fp1=15；通帶截止頻率2：fp2=20；阻帶截止頻率2：fs2=25Hz；Hanning Hamming Blackman Kaiser采用相同的技術指標。以下即是四個窗函數的頻響圖及對簡諧信號濾波后的效果圖。 圖八 不同帶通窗函數低通濾波器的歸一化頻響圖圖九 不同窗函數帶通濾波器對信號的濾波后頻率域效果圖十 不同窗函數帶通濾波器對信號的濾波后時間域振幅效果 圖十一 不同窗函數帶通濾波器對信號的濾波后相位變化由濾波后的相位，振幅及頻率譜可看出：四個窗函數所設計的帶通濾波器的濾波效果基本相當，但hanning窗對振幅譜有失真。對比圖八可以看出：Blackman窗函數過渡帶較窄，旁瓣幅度較小。因此，同等情況下，blackman窗函數設計帶通濾波器效果較好。4. 高通濾波器的設計高通濾波器的技術指標：采樣頻率 Fs=100Hz；阻帶截止頻率 fs=35Hz；通帶起始頻率fp=40HzHanning Hamming Blackman Kaiser采用相同的技術指標。以下即是四個窗函數的頻響圖及對簡諧信號濾波后的效果圖。 圖十二 不同窗函數高通濾波器的歸一化頻響圖 圖十三 不同窗函數高通濾波器對信號的濾波后頻率域效果圖十四 不同窗函數高通濾波器對信號的濾波后時間域振幅效果圖十五 不同窗函數高通濾波器對信號的濾波后相位變化由上面濾波后振幅，相位及頻率譜可看出：hanning Blackman窗函數的振幅及頻率特性有失真。hamming窗的過渡帶較 Kaiser窗函數要陡，旁瓣幅度要小。但從頻率域與振幅特性上看，kaiser效果較好。因此，設計高通濾波器kaiser窗濾波效果比較好。 四、總結體會信號的分析與處理對我們學習地球物理的學生來說，是一個必須要掌握的知識，它對提高重，磁，電，地震信號信噪比，分辨率極其有用，同時也能用于對測井信號的分析與處理。這也是我選這個題目的原因。結合已學的信號分析與處理課程知識調用MATLAB信號處理的相關函數，通過窗函數設計法實現(xiàn)FIR濾波器的設計。一方面使我更加深刻的理解了窗函數法設計FIR濾波器的精髓，另一方面也讓我認識到MATLAB功能的強大性，特別是它簡潔的語句，突出的運算能力，豐富的函數庫以及十分強大的繪圖功能。不過，在做這個題目的過程中，我也遇到的一些困難，比如：kaiserord函數的調用方式，F(xiàn)FT傅立葉變換的實現(xiàn)，以及一些語法等細節(jié)方面的問題。通過查閱相關資料，調用MATLAB中的HELP功能，咨詢老師等，最后，問題基本得到了解決?？偟膩頃r，程序仍存在一些運算速度，簡潔性及穩(wěn)定性方面的問題。這也是我以后需要進一步學習和加強的地方。