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

手把手教你设计FIR 数字滤波器
1. 滤波器的时域、频域、s 域以及离散化
首先，我们要搞清楚一个概念就是滤波器，其实所谓的滤波器就是一个传递函数，它可是通过改变不同频段上信号的幅值来实现滤波，在知道这一点的前提下，下面的讲述就容易了很多。

这里我们假设滤波器的时域传递函数（连续）s 域为()H s ，时域为()h t ；原始信号s 域为()X s ，时域为()x t ；滤波器的输出s 域为()Y s ，时域为()y t ，如图1所示。

（如果分不清楚s 域和时域的童鞋，我觉得你就不要反省了，这个领域不适合你！）我们知道传递函数之间的关系是相乘，而时域的关系是卷积。

那么我们就有了下面的两个关系式。

()()()Y s H s X s = （1） ()()*()y t h t x t =（其中的‘*’表示卷积） （2）
图1 滤波器的传递函数
当然，我们知道在数字滤波器中，当然不可能存在连续函数的，所以我们要对连续函数进行离散化，其实就是一个采样的过程，假设采样频率为S f ，这里我们就需要提到一个定理，就是香农采样定理（采样频率一定要大于传递函数截止频率的两倍，否则就会发生高次谐波的混迭，如果这个道理你不懂的话，建议你去恶补一下信号与系统，这里我只解释一点，便于理解，采样就是以时间间隔为S T 的脉冲采样，那么传递函数的频域就变成了周期为S f 的周期函数，当然如果不理解，也不耽误对下面讲解的理解），所以信号就离散成了()S x nT 、()S h nT 、()S y nT ，其中S T 为采样时间间隔，为了方便下面内容的讲述，这里我们做一个频率的归一化，我们取1S T =，那么采样频率就变成了1S f =，就有了下面的离散时域表示方法。

()()()k y n h k x n k +∞=-∞=
-∑ （3）
值得指出的是再这样种1S T =的归一化过程中，我们认为
2S f 为传递函数系统的最
大频率，这里就是11[,]22
-，不理解的请努力理解。

2. FIR 滤波器的线性相位特性
这里我们讨论两个概念：FIR 滤波器的线性相位特性、传递函数的群延时。

这里给出一个定理：对于()h n ，(0,1)n N ∈-，N 为奇数，如果()h n 关于12
N -呈奇对称或者偶对称，那么FIR 滤波器的频域相位具有线性特性，区别在于当()h n 关于12N -呈奇对称时，初始时刻相位为2
π，当()h n 关于12N -呈奇对称时，初始时刻相位为0.
我们知道控制系统中的延时传递函数为st e -，对应于频域就是j t e ω-，而对于离散之后归一化的结果（延时d n 个步长）d j n e ω-（其中[,]ωππ∈-）。

而这个延时系统其实就是一个线性相位系统，angle(d j n e ω-)=d n ω-,可见，延时系统具有线性相位特性。

3. 理想滤波器特性
首先，我们讲述一下理想滤波器，众所周知，最简单的理想滤波器就是矩形窗滤波器，其相位角为0，幅值在频域上是一个矩形窗函数，如图2所示。

图2 理想滤波器的时域和频域特性
图2中，截止频率c ω为3
π，其频域表达式为： 1 ()0 c c H ωωωωωπ
⎧≤⎪=⎨<<⎪⎩ （4）
其中时域方程可以通过连续函数的傅里叶反变换获得，如下公式所示：
1()()2j t h t H e d πωπωωπ-=⎰ （5）
带入可以求解可得：
sin()()c t h t t
ωπ= （6） 再进行离散化以后可得
sin()()c n h n n
ωπ= （7） 4. FIR 滤波器设计过程
首先，针对于实际中存在的系统，n 为有限的实数，假设n 的个数为N ，这就相当于在时域中用一个窗函数()w n 与滤波器的时域信号(),[,]h n n ∈-∞∞相乘获得。

其中，11 2()1
0 2N n w n N n -⎧≤⎪⎪=⎨-⎪>⎪⎩
（8）
我们知道当0t
<时，()0h t =，所以我们需要将这个时域的窗函数向右平移12
N -个单位，就变成了如下式子。

1 0n 1()0 else
R N w n ≤≤-⎧=⎨⎩ （9） 这个时域的平移之后的窗函数经过傅里叶变换之后的表达式如下所示： 1()2sin(/2)sin(/2)
N j R N W e ωωω--= （10） 如果有疑问请自行推导，我也是犯懒了，直接从书上抄过来的！
接下来的事情就简单了很多，因为我们在第三部分中已经讲述了理想的矩形窗滤波器（频域，用于滤波）的传递函数（公式4所示），只是当时没有加入相
位信息，这里由于时域的窗函数（截断数据使得N 为有限数据）为延时12
N -个单位，所以频域的窗函数所对应的时域信号也应该延时12
N -个单位，由于滤波器的延时特性可知就是窗函数的频谱特性为线性，故可以得到其传递函数如下所示（这里有一个时域的窗函数，有一个频域的窗函数，要分清楚，搞混了的童鞋请面壁！），
12e ()0 N
j c c H ωωωωωωπ
--⎧⎪≤=⎨⎪<<⎩ （11）
有了这两个函数的时域和频域特性之后，因为截断数据用的窗函数是在时域上，所以就是具有延时的矩形窗滤波器的时域函数()h n 与截断窗函数的时域函数()R w n 时域相乘，这里有一个对应关系，需要大家知道，就是时域相乘与频域卷积的对应关系、频域相乘与时域卷积的对应关系，请大家注意其中的差别，不要弄错。

12121212()*()()()
1()()()*()2f t f t F F f t f t F F ωωωωπ
↔↔ （12） 于是FIR 滤波器的频域特性可以得到如下式所示：
121()()*()2sin[()]122sin()2C C R N j FIR H W N e d ωωωωωωπ
ωθθωθ---=
-=-⎰ （13） 为了方便大家的理解，这里我做了一张图来便于大家的理解，如图3所示：
图3 FIR 滤波器的时域和频域构造方式
图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ÉÑùƵÂÊ£¨½ÇƵÂÊ£©
for ii=1:length(n)
if n(ii)==(N-1)/2
hn(ii)=wc/ws*2*pi/pi;
else
hn(ii)=sin(wc/ws*2*pi*(n(ii)-(N-1)/2))./(pi*(n(ii)-(N-1)/2));%N½×Â˲¨Æ÷µÄÀëÉ¢Âö³åÏìÓ¦£¨ÀëɢʱÓò£©
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%Â˲¨Æ÷µÄƵÓò
w=0:2*pi/ws*2*pi:pi;
for i=1:length(w)
zeta=-wc/ws*2*pi:pi*2/ws*2*pi:wc/ws*2*pi;
for iii=1:length(zeta)
if zeta(iii)==w(i)
t_trap(iii)=N;
else
t_trap(iii)=sin(N/2*(w(i)-zeta(iii)))./(sin((w(i)-zeta(iii))/2));
end
end
Hw_amp(i)=trapz(zeta,t_trap)/2/pi;
end
Hw_ang=-(N-1)/2*w;
figure
subplot(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。