MATLAB窗函数法实现FIR的高通,带通和低通滤波器的程序

MATLAB窗函数法实现FIR的高通,带通和低通滤波器的程序

MATLAB 学院：地球物理与石油资源学院班级：姓名：学号：班内编号：指导教师：完成日期：测井11001大牛啊啊啊陈义群2013年6月3日课程设计报告一、题目FIR滤波器的窗函数设计法及性能比较 1. FIR滤波器简介数字滤波器是一种用来过滤时间离散信号的数字系统，通过对抽样数据进行数学处理来达到频域滤波的目的。根据其单位冲激响应函数的时域特性可分为两类：无限冲激响应滤波器和有限冲激响应滤波器。与IIR滤波器相比，FIR滤波器的主要特点为： a. 线性相位；b.非递归运算。 2. FIR 滤波器的设计FIR滤波器的设计方法主要有三种：a.窗函数设计法；b.频率

抽样发；c.最小平法抽样法；这里我主要讨论在MATLAB环境下通过调用信号分析与处理工具箱的几类窗函数来设计滤波器并分析与比较其性能。窗函数法设计FIR滤波器的一般步骤如下： a. 根据实际问题确定要设计的滤波器类型； b. 根据给定的技术指标，确定期望滤波器的理想频率特性；c. 求期望滤波器的单位脉冲响应；d. 求数字滤波器的单位脉冲响应； e. 应用。常用的窗函数有(1)Hanningwindoww(n)?[?((2)Hammingw indoww(n)?[?((3)Balckmanwindoww(n)?[ ?((4)KaiserwindowI0{?1?[2n/(N?1)]2}w(n )?RN(n)I0(?)式中I0(x)是零阶Bessel函数，可定义为()2?n4?n)?()]RN(n)N?1N?1()2?n)]RN(n)N ?1() ?nN?1)]RN(n)() (x/2)m2I0(x)?1??m!m?1? 当x?0时与矩形窗一致；当x?时与海明窗结果相同；当x?时与布莱克曼窗结果相同。3．窗函数的选择标准 1. 较低的旁瓣

幅度，尤其是第一旁瓣； 2. 旁瓣幅度要下降得快，以利于增加阻带衰减；3. 主瓣宽度要窄，这样滤波器过渡带较窄。

4. 常用窗函数的参数

5. FIR 滤波器的MATLAB实现方式在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窗函数，给以相同的技术参数，来设计低通，带通，高通滤波器，用上述窗函数的选择标准来比较各种窗函数的优劣，并给以一个简谐波进行滤波处理，比较滤波前后的效果。达到综合比较的效果。

二、源代码 1.利用hanning hamming blackman kaiser窗,设计一个低通FIR function lowpassfilter clc; 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(*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],[ ],100); window4=kaiser(n+1,beta);%使用kaiser窗函数%设计加窗函数fir1 b1=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 = 2

extpow2(L); % Next power of 2 from length of y Y = 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(‘freq uency/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(‘Am uplitude’); 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(‘Am uplitude’); 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(‘Am uplitude’);

xlabel(‘frequency/Hz’);ylabel(‘Amuplitud e’); 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_FIR function bandpassfilter Fs=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(*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 窗函数%设过渡带宽度为5Hz [n,Wn,beta,ftype]=kaiserord([10 15 20 25],[0 1 0],[ ],100);%求阶数n以及参数beta window4=kaiser(n+1,beta);%使用kaiser窗函数%设计加窗函数fir1 b1=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 = 2

extpow2(L); % Next power of 2 from length of y Y = 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(‘freq uency/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(‘Am uplitude’); 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(‘Am uplitude’); 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(‘Am uplitude’);

xlabel(‘frequency/Hz’);ylabel(‘Amuplitud e’); 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(‘hann ing窗函数滤波相位特性’); 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_FIR function highpassfilter clc; 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],[ ],100); window4=kaiser(n+1,beta); %设计加窗函数fir1 b1=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 = 2

extpow2(L); % Next power of 2 from length of y Y = fft(y,NFFT)/L;%将时域信号变换到频域 f = Fs/2*linspace(0,1,NFFT/2+1);%频域采样figure(2);

plot(f,2*abs(Y(1:NFF T/2+1)));xlabel(‘freq uency/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(‘Am uplitude’); 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/