第36章FIR滤波器的Matlab设计

安安富富莱莱UUMM440033DDSSPP教教程程SSTTMM3322--VV55开开发发板板系系统统篇篇手手册册22001155年年0011月月1155日日版版本本：：11..00第第11页页共共1188页页第第3366章章FFIIRR滤滤波波器器的的MMaattllaabb设设计计本章节讲解FIR滤波器的Matlab设计。主要是函数fir1和fir2的使用。36.1窗函数36.2fir1函数36.2fir2函数36.4总结3366..11窗窗函函数数在数字信号处理中不可避免地要用到数据截取的问题。例如，在应用DFT的时候，数据x(n)总是有限长的，在滤波器设计中遇到了对理想滤波器抽样响应h(n)的截取问题，在功率谱估计中也要遇到对自相关函数的截取问题。总之，我们在实际工作中所能处理的离散序列总是有限长，把一个长序列变换成有限长的序列不可避免的要用到窗函数。因此，窗函数本身的研究及其应用是信号处理中的一个基本问题。不同的窗函数对信号频谱的影响是不一样的，这主要是因为不同的窗函数，产生泄漏的大小不一样，频率分辨能力也不一样。信号的截断产生了能量泄漏，而用FFT算法计算频谱又产生了栅栏效应，从原理上讲这两种误差都是不能消除的，但是我们可以通过选择不同的窗函数对它们的影响进行抑制。(矩形窗主瓣窄，旁瓣大，频率识别精度最高，幅值识别精度最低；布莱克曼窗主瓣宽，旁瓣小，频率识别精度最低，但幅值识别精度最高)。对于窗函数的选择，应考虑被分析信号的性质与处理要求。如果仅要求精确读出主瓣频率，而不考虑幅值精度，则可选用主瓣宽度比较窄而便于分辨的矩形窗，例如测量物体的自振频率等；如果分析窄带信号，且有较强的干扰噪声，则应选用旁瓣幅度小的窗函数，如汉宁窗、三角窗等；对于随时间按指数衰减的函数，可采用指数窗来提高信噪比。矩形窗：矩形窗属于时间变量的零次幂窗。矩形窗使用最多，习惯上不加窗就是使信号通过了矩形窗。这种窗的优点是主瓣比较集中，缺点是旁瓣较高，并有负旁瓣，导致变换中带进了高频干扰和泄漏，甚至出现负谱现象。三角窗：三角窗亦称费杰（Fejer）窗，是幂窗的一次方形式。与矩形窗比较，主瓣宽约等于矩形窗的两倍，但旁瓣小，而且无负旁瓣。汉宁窗：安安富富莱莱UUMM440033DDSSPP教教程程SSTTMM3322--VV55开开发发板板系系统统篇篇手手册册22001155年年0011月月1155日日版版本本：：11..00第第22页页共共1188页页汉宁窗又称升余弦窗，汉宁窗可以看作是3个矩形时间窗的频谱之和，或者说是3个sinc（t）型函数之和，而括号中的两项相对于第一个谱窗向左、右各移动了π/T，从而使旁瓣互相抵消，消去高频干扰和漏能。可以看出，汉宁窗主瓣加宽并降低，旁瓣则显著减小，从减小泄漏观点出发，汉宁窗优于矩形窗．但汉宁窗主瓣加宽，相当于分析带宽加宽，频率分辨力下降。海明窗：海明窗也是余弦窗的一种，又称改进的升余弦窗。海明窗与汉宁窗都是余弦窗，只是加权系数不同。海明窗加权的系数能使旁瓣达到更小。分析表明，海明窗的第一旁瓣衰减为一42dB．海明窗的频谱也是由3个矩形时窗的频谱合成，但其旁瓣衰减速度为20dB/（10oct），这比汉宁窗衰减速度慢。海明窗与汉宁窗都是很有用的窗函数。高斯窗：三角窗亦称费杰（Fejer）窗，是幂窗的一次方形式。与矩形窗比较，主瓣宽约等于矩形窗的两倍，但旁瓣小，而且无负旁瓣。还有很多其它的窗口这里就不做介绍了，需更详细的了解的话，可以看matlab中help文档中的如下部分（这里只截图了一部分窗口）。或者直接在命令窗口输入wintool可以打开窗口工具：安安富富莱莱UUMM440033DDSSPP教教程程SSTTMM3322--VV55开开发发板板系系统统篇篇手手册册22001155年年0011月月1155日日版版本本：：11..00第第33页页共共1188页页打开后界面如下：3366..22ffiirr11函函数数3366..22..11ffiirr11函函数数介介绍绍函数fir1用来设计标准频率响应的基于窗函数的FIR滤波器，可实现加窗线性相位FIR滤波器设计。语法：b=fir1(n,Wn)b=fir1(n,Wn,'ftype')b=fir1(n,Wn,window)b=fir1(n,Wn,'ftype',window)安安富富莱莱UUMM440033DDSSPP教教程程SSTTMM3322--VV55开开发发板板系系统统篇篇手手册册22001155年年0011月月1155日日版版本本：：11..00第第44页页共共1188页页b=fir1(...,'normalization')其中，n：为了滤波器的阶数；Wn：为滤波器的截止频率；ftype：参数用来决定滤波器的类型，当ftype=high时，可设计高通滤波器，当ftype=stop时，可设计带阻滤波器。Window参数用来指导滤波器采用的窗函数类型。其默认值为汉明（Hamming）窗。使用fir1函数可设计标准的低通，高通，带通和带阻滤波器。滤波器的系数包含在返回值b中，可表示为：b(z)=b(1)+b(2)z-1+……+b(n+1)z-n(1)采用汉明窗设计低通FIR滤波器使用b=fir1(n,Wn)可得到低通滤波器。其中，0≤Wn≤1,Wn=1相当于0.5f。其语法格式为b=fir1(n,Wn)(2)采用汉明窗设计高通FIR滤波器在b=fir1(n,Wn,'ftype')中，当ftype=high时，可设计高通滤波器。其语法格式为b=fir1(n,Wn,'high')(3)采用汉明窗设计带通FIR滤波器在b=fir1(n,Wn)中，当Wn=[W1W2]时，fir1函数可得到带通滤波器，其通带为W1WW2W1和W2分别为通带的下限频率和上限频率。其语法格式为b=fir1(n,[W1W2])(4)采用汉明窗设计带阻FIR滤波器在b=fir1(n,Wn,'ftype')中，当ftype=stop，Wn=[W1W2]时，fir1函数可得到带阻滤波器，其语法格式为b=fir1(n,[W1W2],'stop')(5)采用其他窗口函数设计FIR滤波器使用Window参数，可以用其他窗口函数设计出各种加窗滤波器，Window参数可采用的窗口函数有Boxcar，Hamming，Bartlett，Blackman，Kasier和Chebwin等。其默认时为Hamming窗。例如，采用Bartlett窗设计带阻滤波器，其语法结构为b=fir1(n,[W1W2],'stop',Bartlett[n+1])注意：用fir1函数设计高通和带阻滤波器时，所使用的阶数n应为偶数，当输入的阶数n为奇数时，fir1函数会自动将阶数增加1形成偶数。3366..22..22ffiirr11设设计计低低通通滤滤波波器器实实例例下面我们通过一个实例来讲解fir1的用法。原始信号是由50Hz正弦波和200Hz的正弦波组成，将200Hz的正弦波当做噪声滤掉，下面通过函数fir1设计一组低通滤波器系数，其阶数是30，截止频率为0.25（也就是125Hz）。Matlab运行代码如下：安安富富莱莱UUMM440033DDSSPP教教程程SSTTMM3322--VV55开开发发板板系系统统篇篇手手册册22001155年年0011月月1155日日版版本本：：11..00第第55页页共共1188页页%****************************************************************************************%FIR低通滤波器设计%***************************************************************************************fs=1000;%设置采样频率1kN=1024;%采样点数n=0:N-1;t=0:1/fs:1-1/fs;%时间序列f=n*fs/N;%频率序列Signal_Original=sin(2*pi*50*t);%信号50Hz正弦波Signal_Noise=sin(2*pi*200*t);%噪声200Hz正弦波Mix_Signal=Signal_Original+Signal_Noise;%将信号Signal_Original和Signal_Original合成一个信号进行采样subplot(221);plot(t,Mix_Signal);%绘制信号Mix_Signal的波形xlabel('时间');ylabel('幅值');title('原始信号');gridon;subplot(222);y=fft(Mix_Signal,N);%对信号Mix_Signal做FFTplot(f,abs(y));xlabel('频率/Hz');ylabel('振幅');title('原始信号FFT');gridon;b=fir1(30,0.25);%30阶FIR低通滤波器，截止频率125Hz%y2=filter(b,1,x);y2=filtfilt(b,1,x);%经过FIR滤波器后得到的信号Ps=sum(Signal_Original.^2);%信号的总功率Pu=sum((y2-Signal_Original).^2);%剩余噪声的功率安安富富莱莱UUMM440033DDSSPP教教程程SSTTMM3322--VV55开开发发板板系系统统篇篇手手册册22001155年年0011月月1155日日版版本本：：11..00第第66页页共共1188页页SNR=10*log10(Ps/Pu);%信噪比y3=fft(y2,N);%经过FIR滤波器后得到的信号做FFTsubplot(223);plot(f,abs(y3));xlabel('频率/Hz');ylabel('振幅');title('滤波后信号FFT');gridon;[H,F]=freqz(b,1,512);%通过fir1设计的FIR系统的频率响应subplot(224);plot(F/pi,abs(H));%绘制幅频响应xlabel('归一化频率');title(['Order=',int2str(30),'SNR=',num2str(SNR)]);gridon;Matlab的运行结果如下：从运行结果的FFT和信噪比来看，滤波效果比较明显。00.51-2-1012时间幅值原始信号050010000200400600频率/Hz振幅原始信号FFT050010000200400600频率/Hz振幅滤波后信号FFT00.5100.511.5归一化频率Order=30SNR=25.9258安安富富莱莱UUMM440033DDSSPP教教程程SSTTMM3322--VV55开开发发板板系系统统篇篇手手册册22001155年年0011月月1155日日版版本本：：11..00第第77页页共共1188页页3366..22..33ffiirr11设设计计高高通通滤滤波波器器实实例例下面我们通过一个实例来讲解fir1的高通滤波器的用法。原始信号是由50Hz正弦波和200Hz的正弦波组成，将50Hz的正弦波当做噪声滤掉，下面通过函数fir1设计一组高通滤波器系数，其阶数是30，截止频率为0.25（也就是125Hz）。Matlab运行代码如下：%****************************************************************************************%FIR高通滤波器设计%***************************************************************************************fs=1000;%设置采样频率1kN=1024;%采样点数n=0:N-1;t=0:1/fs:1-1/fs;%时间序列f=n*fs/N;%频率序列Signal_Original=sin(2*pi*2