FIR数字滤波器的Matlab实现.docx
FIR数字滤波器的Matlab实现第7章 FIR数字滤波器的Matlab实现 7.1 实验目的 l 学习用窗函数法设计FIR数字滤波器的原理及其设计步骤； l 学习编写数字滤波器的设计程序的方法，并能进行正确编程； l 根据给定的滤波器指标，给出设计步骤。 7.2 实验原理及实例分析 7.2.1 FIR低通数字滤波器的设计原理 如果系统的冲激响应hd(n)为已知，则系统的输入/输出关系为： y(n)=x(n)×hd(n) 对于低通滤波器，只要设计出低通滤波器的冲激响应函数，就可以由上式得到系统的输出了。假设所希望的数字滤波器的频率响应为Hd(ejw)，它是频域的周期函数，周期为2p，那么它与Hd(ejw)相对应的傅立叶系数为 hd(n)=12pòp-pHd(ejw)ejnwdw 以hd(n)为冲激响应的数字滤波器将具有频域响Hd(ejw)。 但是将hd(n)作为滤波器脉冲响应有两个问题： 它是无限长的，与FIP滤波器脉冲响应有限长这一前提不一致 它是非因果的，hd(n)¹0,n<0 对此，要采取以下的措施， 将hd(n)截短 将其往右平移， 由此得到h2(n)的实际频域响应Hd(e)=jwåh(n)e2n=0N-1jnw，与理想频域响应Hd(ejw)相近，但不完全一致。理论证明上述现象是对hd(n)进行简单截短处理的必然结果，一般称为吉布斯现象，为尽可能的减少吉布斯现象，应对hd(n)进行加窗截取，即以h(n)=hd(n)×WN(n)作为FIR滤波器的系数。 常用的窗函数有矩形窗、海明窗和布莱克曼窗等。 7.2.2 用窗函数法设计FIR滤波器 Matlab设计FIR滤波器有多种方法和对应的函数，见表71。 表71 matlab设计FIR滤波器的方法和函数 方法 窗方法 多带方法 最小二乘法 任意响应法 余弦法 描述 使用窗函数和逆傅立叶变换实现 包含子带频率域 使用最小二乘法将整个频率域上的错误几率压缩到最小 使用任意响应，包括非线性相位以及复滤波器 使用三角函数的低通响应 函数 fir1,fir2,kaiserord等 firls，remez等 fircls ，fircls1等 cremez等 firrcos等 窗函数方法不仅在数字滤波器的设计中占有重要的地位，同时可以用于功率谱的估计，从根本上讲，使用窗函数的目的就是消除由无限序列的截短而引起的Gibbs现象所带来的影响。 窗函数设计线性相位FIR滤波器步骤如下： 确定数字滤波器的性能要求，临界频率wk，滤波器单位脉冲响应长度N 根据性能要求，合理选择单位脉冲响应h(n)的奇偶对称性，从而确定理想频率响应Hd(ejw)的幅频特性和相频特性 求理想单位脉冲响应hd(n)，在实际计算中，可对Hd(ejw)采样，并对其求IDFT的hM(n)，用hM(n)代替hd(n) 选择适当的窗函数w，根据h(n)=hd(n)×WN(n)求所需设计的FIR滤波器单位脉冲响应 求Hd(ejw)，分析其幅频特性，若不满足要求，可适当改变窗函数形式或长度N，重复上述设计过程，以得到满意的结果 FIR1.m 调用格式： a) b=fir1(N,wn) b) b=fir1(N,wn,hign) N为滤波器的阶次，wn是通带截止频率，其值在01之间，1对应抽样频率的一半；b是设计好的滤波器系数h(n)。 对于格式a)，若wn是标量，则用来设计低通滤波器，若wn是1×L的向量，则用来设计带通(带阻)滤波器。其格式为： b=fir1(N,wn,DC-1)或b=fir1(N,wn,DC-0),其中，前者保证第一个带为通带，后者保证第一个带为阻带。 b）用来设计高通滤波器。 用于产生窗函数的matlab文件有如下8个： bartlett.m(巴特利特窗) blackman.m(布莱克曼窗) rectangle.m(矩形窗) hamming.m hanning.m triang.m(三角窗) chebwin.m Kaiser.m 设计一个34阶的高通滤波器，截止频率为0.48p，使用具有30dB波纹的chebyshev窗。 解：源程序如下： b=fir1(34,0.48,'high',chebwin(35,30); freqz(b,1,512) 其响应波形如图7-1所示。 图71 带通FIR滤波器 设计一个24阶FIR带通滤波器，通带频率0.35p£w£0.65p。 解：源程序如下： wn=0.35 0.65; N=24; b=fir1(2*N,wn,'DC-0');% DC-0表示第一个带为阻带，这样设计出来的滤波器为带通滤波器。 freqz(b,1,512) 图72 带通FIR滤波器 设计具有下面指标的低通FIR滤波器 wp=0.2p,Rp=0.25dB,ws=0.3p,AS=50dB由于其最小阻带衰减为50dB，因此可以选择hamming窗来实现这个滤波器，因为它具有较小的过渡带。 解：MATLAB源程序为 %数字滤波器指标 wp=0.2*pi; ws=0.3*pi; tr_width=ws-wp; %通带宽度 N=ceil(6.6*pi/tr_width)+1; %滤波器长度 if rem(N,2)=0 N=N+1; end Nw=N; %滤波器长度为奇数 n=0:1:Nw-1; wc=(ws+wp)/2; %截止频率 alpha=(Nw-1)/2; m=n-alpha+0.00001; %如果不加0.00001，当n= alpha 时，m为0，下式分母为0，不能执行 hd=sin(wc*m)./(pi*m); %理想滤波器的冲激响应 % 生成hamming窗 w_ham=(hamming(Nw)' % 频域图像的绘制 h=hd.*w_ham; % 滤波器实际冲激响应 freqz (h,1) figure(2); subplot(2,2,1),stem(n,hd);title('理想脉冲响应') axis(0 Nw-1 -0.3 0.3);xlabel('n');ylabel('hd(n)') subplot(2,2,2),stem(n,w_ham);title('hamming窗') axis(0 Nw-1 -0.3 1.2);xlabel('n');ylabel('w(n)') subplot(2,2,3),stem(n,h);title('实际脉冲响应') axis(0 Nw-1 -0.3 0.3);xlabel('n');ylabel('h(n)') 其响应波形如图7-2所示。 图73 hamming窗函数设计FIR滤波器 设带通滤波器的指标为 wp1=0.35p,wp2=0.8p,ws1=0.2p,ws2=0.65p,As=60dB选择Blackman窗来实现这个滤波器。 解：MATLAB源程序为 % 数字滤波器指标 ws1=0.2*pi; wp1=0.35*pi; ws2=0.65*pi; wp2=0.8*pi; As=60; tr_width=min(wp1-ws1),(wp2-ws2); % 过渡带带宽 N=ceil(11*pi/tr_width)+1; if rem(N,2)=0 N=N+1; end Nw=N; n=0:1:Nw-1; wc1=(ws1+wp1)/2; wc2=(ws2+wp2)/2; alpha=(Nw-1)/2; m=n-alpha+0.00001; %如果不加0.00001，当n= alpha 时，m为0，下式分母为0，不能执行 hd=sin(wc2*m)./(pi*m)- sin(wc1*m)./(pi*m); %理想滤波器的冲激响应 % 生成blackman窗 w_bla=(blackman(M)' h=hd.*w_bla; % 频域图像的绘制 freqz(h,1) figure(2); subplot(2,2,1),stem(n,hd);title('idael impulse response') axis(0 Nw-1 -0.4 0.5);xlabel('n');ylabel('hd(n)') subplot(2,2,2),stem(n,w_bla);title('blackman window') axis(0 Nw-1 0 1.1);xlabel('n');ylabel('w(n)') subplot(2,2,3),stem(n,h);title('actual impulse response') axis(0 Nw-1 -0.4 0.5);xlabel('n');ylabel('h(n)') 其响应波形如图7-4所示。 图74 blackman窗函数设计FIR滤波器 7.3 编程练习 1. 设计一个48阶FIR带阻滤波器，通带频率0£w£0.35p,0.65p£w£p。 2.用窗函数设计一个线性相位FIR高通滤波器，并满足性能指标：通带边界频率wp=0.66p,阻带边界频率ws=0.5p，阻带衰减不小于40db。 提示：参考例73，4；首先根据阻带衰减，选择合适的窗函数，见课本P342，表7-3。再根据P345，式7-74写出理想滤波器的冲激响应hd 。
