手把手教你设计FIR数字滤波器

手把手教你设计FIR数字滤波器1.滤波器的时域、频域、s域以及离散化首先，我们要搞清楚一个概念就是滤波器，其实所谓的滤波器就是一个传递函数，它可是通过改变不同频段上信号的幅值来实现滤波，在知道这一点的前提下，下面的讲述就容易了很多。这里我们假设滤波器的时域传递函数（连续）s域为()Hs，时域为()ht；原始信号s域为()Xs，时域为()xt；滤波器的输出s域为()Ys，时域为()yt，如图1所示。（如果分不清楚s域和时域的童鞋，我觉得你就不要反省了，这个领域不适合你！）我们知道传递函数之间的关系是相乘，而时域的关系是卷积。那么我们就有了下面的两个关系式。()()()YsHsXs（1）()()*()ythtxt（其中的‘*’表示卷积）（2）H(s)或者h(t)X(s)或者x(t)Y(s)或者y(t)图1滤波器的传递函数当然，我们知道在数字滤波器中，当然不可能存在连续函数的，所以我们要对连续函数进行离散化，其实就是一个采样的过程，假设采样频率为Sf，这里我们就需要提到一个定理，就是香农采样定理（采样频率一定要大于传递函数截止频率的两倍，否则就会发生高次谐波的混迭，如果这个道理你不懂的话，建议你去恶补一下信号与系统，这里我只解释一点，便于理解，采样就是以时间间隔为ST的脉冲采样，那么传递函数的频域就变成了周期为Sf的周期函数，当然如果不理解，也不耽误对下面讲解的理解），所以信号就离散成了()SxnT、()ShnT、()SynT，其中ST为采样时间间隔，为了方便下面内容的讲述，这里我们做一个频率的归一化，我们取1ST，那么采样频率就变成了1Sf，就有了下面的离散时域表示方法。()()()kynhkxnk（3）值得指出的是再这样种1ST的归一化过程中，我们认为2Sf为传递函数系统的最大频率，这里就是11[,]22，不理解的请努力理解。2.FIR滤波器的线性相位特性这里我们讨论两个概念：FIR滤波器的线性相位特性、传递函数的群延时。这里给出一个定理：对于()hn，(0,1)nN，N为奇数，如果()hn关于12N呈奇对称或者偶对称，那么FIR滤波器的频域相位具有线性特性，区别在于当()hn关于12N呈奇对称时，初始时刻相位为2，当()hn关于12N呈奇对称时，初始时刻相位为0.我们知道控制系统中的延时传递函数为ste，对应于频域就是jte，而对于离散之后归一化的结果（延时dn个步长）djne（其中[,]）。而这个延时系统其实就是一个线性相位系统，angle(djne)=dn,可见，延时系统具有线性相位特性。3.理想滤波器特性首先，我们讲述一下理想滤波器，众所周知，最简单的理想滤波器就是矩形窗滤波器，其相位角为0，幅值在频域上是一个矩形窗函数，如图2所示。图2理想滤波器的时域和频域特性图2中，截止频率c为3，其频域表达式为：1()0ccH（4）其中时域方程可以通过连续函数的傅里叶反变换获得，如下公式所示：1()()2jthtHed（5）带入可以求解可得：sin()()cthtt（6）再进行离散化以后可得sin()()cnhnn（7）4.FIR滤波器设计过程首先，针对于实际中存在的系统，n为有限的实数，假设n的个数为N，这就相当于在时域中用一个窗函数()wn与滤波器的时域信号(),[,]hnn相乘获得。其中，112()102NnwnNn（8）我们知道当0t时，()0ht，所以我们需要将这个时域的窗函数向右平移12N个单位，就变成了如下式子。10n1()0elseRNwn（9）这个时域的平移之后的窗函数经过傅里叶变换之后的表达式如下所示：1()2sin(/2)sin(/2)NjRNWe（10）如果有疑问请自行推导，我也是犯懒了，直接从书上抄过来的！接下来的事情就简单了很多，因为我们在第三部分中已经讲述了理想的矩形窗滤波器（频域，用于滤波）的传递函数（公式4所示），只是当时没有加入相位信息，这里由于时域的窗函数（截断数据使得N为有限数据）为延时12N个单位，所以频域的窗函数所对应的时域信号也应该延时12N个单位，由于滤波器的延时特性可知就是窗函数的频谱特性为线性，故可以得到其传递函数如下所示（这里有一个时域的窗函数，有一个频域的窗函数，要分清楚，搞混了的童鞋请面壁！），12e()0NjccH（11）有了这两个函数的时域和频域特性之后，因为截断数据用的窗函数是在时域上，所以就是具有延时的矩形窗滤波器的时域函数()hn与截断窗函数的时域函数()Rwn时域相乘，这里有一个对应关系，需要大家知道，就是时域相乘与频域卷积的对应关系、频域相乘与时域卷积的对应关系，请大家注意其中的差别，不要弄错。12121212()*()()()1()()()*()2ftftFFftftFF（12）于是FIR滤波器的频域特性可以得到如下式所示：121()()*()2sin[()]122sin()2CCRNjFIRHWNed（13）为了方便大家的理解，这里我做了一张图来便于大家的理解，如图3所示：图3FIR滤波器的时域和频域构造方式图3中红字是时域的窗函数（截断函数），绿字是频域的窗函数（理想滤波器函数）。5.FIR滤波器设计举例我们来举一个例子来方便大家理解，首先是滤波器的阶数为N=81（滤波器的阶数80），截止频率c=200HZ，采样频率S=2000HZ。我们可以将设计的滤波器的特性展示如图4.为了方便大家交流，我把我写的程序贴到最后，欢迎大家交流。图4滤波器的特性本文纯属手打，感谢大家对我的支持，本人现在正在创办一个matlab交流平台，会定期发布一些有创意的小程序和matlab的一些使用技巧，欢迎大家关注，微信公众号如下所示：clear;clc;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%Â˲¨Æ÷µÄʱÓòN=81;%Â˲¨Æ÷µÄ½×ÊýNΪÆæÊýn=0:N-1;wc=200*2*pi;%Â˲¨Æ÷µÄ½ØֹƵÂÊ£¨½ÇƵÂÊ£©ws=2000*2*pi;%ÐźŵIJÉÑùƵÂÊ£¨½ÇƵÂÊ£©forii=1:length(n)ifn(ii)==(N-1)/2hn(ii)=wc/ws*2*pi/pi;elsehn(ii)=sin(wc/ws*2*pi*(n(ii)-(N-1)/2))./(pi*(n(ii)-(N-1)/2));%N½×Â˲¨Æ÷µÄÀëÉ¢Âö³åÏìÓ¦£¨ÀëɢʱÓò£©endend%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%Â˲¨Æ÷µÄƵÓòw=0:2*pi/ws*2*pi:pi;fori=1:length(w)zeta=-wc/ws*2*pi:pi*2/ws*2*pi:wc/ws*2*pi;foriii=1:length(zeta)ifzeta(iii)==w(i)t_trap(iii)=N;elset_trap(iii)=sin(N/2*(w(i)-zeta(iii)))./(sin((w(i)-zeta(iii))/2));endendHw_amp(i)=trapz(zeta,t_trap)/2/pi;endHw_ang=-(N-1)/2*w;figuresubplot(5,1,1)plot(n,hn)subplot(5,1,2)plot(w/pi*(ws/2/pi)/2,Hw_amp)subplot(5,1,3)plot(w/pi*(ws/2/pi)/2,Hw_ang)subplot(5,1,4)plot(w/pi*(ws/2/pi)/2,20*log(abs(Hw_amp))/log(10))%subplot(5,1,5)plot(w/pi*(ws/2/pi)/2,Hw_ang)支持原创，转载请注明出处！作者：leo_silverbullet