如何快速设计一个FIR滤波器(二)通俗易懂,膜拜

一、理想低通滤波器单位脉冲响应是什么样

如何快速设计一个FIR滤波器(一)中,我们介绍了一种简单设计FIR的方法——零极点法。这个方法非常简单,稍加培训,用笔和纸就能完成;当然缺点也很显而易见:零极点设计出的滤波器,只能给出大概的频率响应,对于一些要求较高的系统,显得无能为力。今天我们介绍一种更加严谨的方法。

现在假设我们要设计一个低通滤波器,截止频率为 omega_c ,其理想频率响应可以用如下函数表示:

H_d(e^{jomega})=left{ egin{aligned} 1  ,quad &left| omega 
ight| leq omega_c\ 0, quad& left| omega 
ight|   >omega_c \ end{aligned} 
ight.

则该滤波器的脉冲响应为:

h_d(t)=frac{1}{2pi}int_{-infty}^{+infty}H_d(e^{jomega})e^{jomega t}domega=frac{sin(omega_c t)}{piomega_c}

可见脉冲响应为一个sinc函数。画出来大概张这个样子:

这个函数非常重要哦,建议都自己手动推导一下,非常有意思,我们把幅值最大的那一瓣称之为主瓣,剩下称之为从瓣,是不是很形象呢?注意主瓣的周期是从瓣的2倍哦。

可能你已经看见了,脉冲响应是由无穷多个——不对呀,我们是想设计一个有限脉冲响应的,结果出现了无穷多个响应,这可如何搞?——既然臣妾做不到把所有的脉冲响应都用上(因为我是有限脉冲滤波器FIR啊),那我只能截取有限一部分进行代替了,也就是——加窗。

窗函数这个东西很不好理解,我们就多花点功夫说说这个事。


二、什么是加窗

现在假设有一个函数,表达式为 y(t)=Acos(2pi f t) ,为简单起见,令 A=1 , f=1,很显然,余弦函数的频域包含两个分量 ±f ,即 ±1 Hz。

假如我们现在以10Hz的频率进行采样,显然是满足香浓采样定理的,采样脉冲及其频谱如下图:

采样过程就是拿采样信号和原始信号进行乘积,那采样之后的信号长啥样呢?

上图即为采样后的信号及其频谱。注意:此时采样后的信号是在时域无限扩展的,频谱也一样。但实际的计算中是不可能处理无穷多数量的信号的,那怎么办呢?——截短呗,只观察一部分样本行不?貌似我们也没有更好的选择了。

那怎么截短呢?最简单的方式就是矩形函数了,如下图所示:

矩形函数的方法简单粗暴,只保留一部分(幅值为1),其余全部设为0,然后那矩形函数和采样信号相乘,就得到加窗(矩形窗)之后的采样信号了。矩形信号的频谱就是我们熟悉的sinc函数。加窗之后的采样信号及其频谱如下:

加窗后的离散余弦信号可以看成窗函数(矩形,频谱为sinc函数)与余弦离散信号的乘积。那在频域呢就是sinc函数和余弦频谱的卷积。看一个图就明白加窗后的频谱是怎么来的了:

可以看出,由于加窗行为,频谱中频率成分变复杂了,由于矩形窗函数的频谱是无限宽的,在频谱的最大值处也出现了一定的畸变(略大于1)。这种有限长的离散信号在时域中计算机是可以处理(离散的),但是在频域还是不能,因为还是连续信号。在如何理解离散傅里叶变换及Z变换一文中,我们介绍了将有限长的信号进行周期延拓,就得到一个周期的离散信号,其傅里叶变换也是周期离散的,这样就可以在计算机中处理了。关于离散傅里叶变换,可以参照:

J Pan:如何理解离散傅里叶变换及Z变换​zhuanlan.zhihu.com图标

现在问题就变成了如何进行周期延拓——卷积,什么意思呢?我们现在找一个采样脉冲,时域中周期和加窗后的采样信号的周期一致(都是1),其形貌如下图(左)所示:

我们已经知道,周期延拓采样脉冲的频谱仍未采样脉冲如上图(右)所示。现在就可以做一个有意思的事情:周期延拓在数学上可以看成是一个有限长的信号(长度即为周期)和一个同样周期的无限长采样脉冲的卷积。什么意思呢?看一个图就一目了然了:

感兴趣的可以用卷积的数学定义推导一下看看是不是这样,很简单的。由卷积定理我们知道,卷积和乘积在时域和频域是对偶的,也就是说刚才我们所做的周期延拓(时域卷积)在频域是乘积哦,具体如下图所示。

所以最终得到的经过周期延拓的加窗信号及其频谱如下图所示。可以看出,进过周期延拓,加窗导致的频谱泄漏似乎得到了抑制,加窗对频谱没什么影响。

事实是不是这样呢?——直觉告诉我们应该不是。我们之所以最终得到这个结论,是因为我们选的初始函数太简单了。余弦函数是周期信号,频谱是有限宽度的,如果窗函数的宽度是信号周期的整数倍,则频域采样(时域周期延拓)将正好采在sinc函数的从瓣取值为零的地方,泄漏的频谱没有采到而已,具体如下图所示。

可能心细有童鞋会有疑问,为啥主瓣的最大值会有一定的偏移?——原因就是sinc函数主瓣和从瓣的周期不一样,卷积时主瓣与从半叠加,最大值就偏移了。

实际工程中,我们可能不知道被处理信号的周期,或者压根就不是周期的信号,那会出现什么呢?现在假如我们把窗函数加宽,比如增加到原来的1.4倍,如下图所示:

则可以得到加窗后的信号及频谱为:

加窗并进行周期延拓后的信号及频谱为:

可见,由于窗函数宽度和原信号周期不一致,加窗后的信号进行周期延拓时出现了不连续,频域采样时不像前面那样只采集到sinc为零的地方,频谱出现了较大的泄漏,如下图所示。

实际工程上绝大多数信号都不是周期信号,也不是有限长信号,而是我们只能观察有限长度而已。要想不出现频谱泄漏,需要窗函数长度是被处理信号频谱所有分量的整数倍,这几乎是不能达到的。因此加窗后,频谱泄漏是不可避免的,只能尽量减小。

本节中所用的算例来自离散傅里叶变换DFT基本原理图解 ,感兴趣请前往阅读。

三、窗函数都有哪些

最开始我们介绍了一个理想的低通滤波器,在频域内其幅值响应是一个矩形。这个矩形进行傅里叶逆变换,发现单位冲击响应有无穷多个,我们没办法,用了一个矩形窗进行截短,注意着两个矩形不是同一个事情哦,一个是幅值响应(频域),一个截短函数(时域),只不过样子都呈现矩形。

现在我们都放在频域来看:首先时域矩形函数在频域为sinc函数,时域的矩形函数截短(乘积)在频域对应卷积。也就是说在频域里,脉冲响应截短后的滤波器频域幅值响应是矩形信号和sinc信号的卷积,仔细想想这段话哦!下图展示了卷积之后的信号的样子。

可以想象:假设窗函数宽度无穷大,则对应的sinc函数无穷窄,则卷积之后的函数非常接近理想矩形;反之,若窗函数很窄,在sinc函数就会变得很宽,则卷积之后的函数与理想矩形差别就比较大,这就是传说中的不确定原理哦,具体可参照:

J Pan:如何理解不确定性原理—不确定or测不准?​zhuanlan.zhihu.com图标

那什么样的窗函数比较好呢?——在频域看就是要能量相对集中,也就是旁瓣要低。因为出现泄漏的主要原因就是窗函数的频谱是无限长的,与信号的频谱卷积时,主瓣与从瓣会出现叠加,因此从瓣的能量越小,影响就越小。在时域看就是加窗后的函数在进行周期延拓时尽量不要有非连续,因为非连续就需要很多分量来模拟,就造成了频谱泄露。那什么样函数的频谱主瓣能量较集中呢?——主要有以下是几种常见的窗函数:矩形窗、汉宁窗、汉明窗、凯撒窗等。

这些窗函数对应的频谱分别为:

现在我们拿汉宁窗(Hanning)来举个例子。还是原来的余弦函数y(t)=Acos(2pi f t) , A=1 , f=1 。假设窗函数的宽度 wT=1.4 ,Hanning窗窗函数及其频谱为:

可见其频谱中主瓣能量与矩形窗比,集中了很多。余弦信号采样、加窗之后的信号及其频谱为:

由于Hanning窗从零开始,结束的时候也是零,因此加窗之后的信号开始和结束也都是零,连续性好了,从频谱看,泄漏也少了。这个从周期延拓后(DTFT)的信号与频谱看的更清楚。


四、如何基于窗函数法设计FIR

这篇文章中,我们主要说怎么较为精确的设计一个FIR,我们花了大量篇幅说了窗函数,这是因为窗函数在实际的工程应用中很多,比如基于窗函数的FIR就是FIR设计中最广泛运用的一个方法。接下来,我们要介绍一下什么是用窗函数设计FIR。

在文章开始的时候,我们已经介绍了,对于理想低通滤波器,其脉冲响应为sinc函数,且有无穷多个,我们必须加窗才能处理,加窗就会有一定的频谱泄漏,因此实际设计出来的滤波器和理想滤波器还是有差别的。

基于窗函数的FIR设计就是通过选择合适常函数来找到满足要求的滤波器,一般步骤是这样的:

1)确定频域的响应函数 H_d(e^{jomega }) ,低通、高通、带通或者其他;

2)确定频域的响应函数 H_d(e^{jomega }) 的傅里叶逆变换,找到连续脉冲响应函数 h_d(t) ;

3)对脉冲响应函数 h_d(t) 按照一定的采样频率进行采样,获得离散信号 h_d[m] ;

4)选择合适的窗函数,对离散信号 h_d[m] 加窗,获得有限长度的脉冲响应 h[m] ;

当然实际实践中,我们只需要知道这个过程,具体细节不需要我们考虑,因为MATLAB提供了一个强大的函数,帮我们都做好了——fir1。

fir1的基本语法如下:h = fir1(n,Wn,ftype,window)

其中n表示滤波器的阶数(order),Wn表示归一化后的滤波器的截止频率,可表示成 [f_{low} quad f_{high}] 的形式。举个例子,比如一个带状滤波器,采样频率 f_s 是200Hz,两个截止频率分别为10Hz和50Hz,则 [f_{low} quad f_{high}]=[0.1 quad 0.5] ,即截止频率除以 frac{1}{2}f_s 。ftype表示滤波器类型,可以是lowpass, highpass, bandpass, bandstop, or multiband filter。window表示窗函数类型,前面我们列到的窗函数都可以选择。h为所设计的滤波器的系数。

举个简单的例子:

h = fir1(48,[0.35 0.65]);
freqz(b,1,512)

这是一个48阶的带通滤波器,归一化的截止频率为[0.35 0.65],其幅值响应和相位响应如下图所示。

可见,采用fir1函数设计FIR滤波器非常方便。

除了fir1函数,fir2函数也有时候会用到,fir2函数是什么意思呢?——我们称之为基于频率采样法的设计方法,具体做法就是:把滤波器期望的幅值响应用一组分立的点 (f_i,A(f_i)) quad i=1,2,3,...M 来表示, A(f_i) 表示频率为 f_i 时的期望幅值响应,然后在不同的点之间进行线性插值,得到完成的期望幅值响应,然后进行傅里叶逆变换并用Hamming(汉明窗)截短来获得滤波器的系数。

基本语法为:h = fir2(n,f,m)

n表示滤波器阶数,f表示分离点的矢量,m表示分立点对应的幅值响应矢量,举个例子就明白了。现在要设计一个低通滤波器,幅值响应如下图所示:

显然四个特征点的坐标分别为(0,1), (0.2,1), (0.2,0), (1,0) ,于是我们可以获得分立点的矢量为f=[0 0.2 0.2 1] ,对应幅值响应矢量为m=[1 1 0 0]。

f=[0 0.2 0.2 1] ;
m=[1 1 0 0];
h = fir2(48,f,m);
freqz(b,1)


五、如何基于响应最优法设计FIR

前面我们说了两种基于窗函数的方法,主要是用窗函数对无限长的脉冲进行截短,获得近似的频域响应。除了窗函数法,还有另外的一种解决方案也比较常用,我们称之为“最优法”,主要思路就是找到一组脉冲响应,让它的频域响应 H(e^{jomega}) 与期望的滤波器的频域响应 H_d(e^{jomega}) 尽可能的一致,主要通过两种方法来实现,一个是最小二乘法,另一个是切比雪夫法。

(一)最小二乘法

与频率采样法近似,期望的频率响应用一组分立的点 (f_i,A(f_i)) quad i=1,2,3,...M 来表示,在不同的点之间进行插值,然后来寻找滤波器的系数 {h} 时期满足

sum_{i=1}^{N}W_ileft[ H(f_i)-H_d(f_i) 
ight]^2
ightarrow min

其中 W_i 为不同频率下的权重。

MATLAB指令为:h=firls(n,f,m)

n为滤波器阶数,f表示分离点的矢量,m表示分立点对应的幅值响应矢量。仍以前面说的低通滤波器为例:

f=[0 0.2 0.2 1] ;
m=[1 1 0 0];
h = firls(48,f,m);
freqz(b,1)

(二)切比雪夫法(又叫帕克斯-麦克劳林法)

与最小二乘法不同(方差最小),切比雪夫法采用的方案是最大误差最小,即:

maxleft| W_ileft( H(f_i)-H_d(f_i) 
ight) 
ight|
ightarrow min

MATLAB指令为:h=firpm(n,f,m)

n、f、m定义和前面一致,不同的是频率响应在同一频率下必须去不同的值。前面的例子中当频率为0.2(归一化后)时,频率响应可直接从1降到0,对于切比雪夫法,规定不能这样做,我们将第二个0.2改为0,21。

f=[0 0.2 0.21 1] ;
m=[1 1 0 0];
h = firpm(48,f,m);
freqz(b,1)


六、如何利用MATLAB设计FIR滤波器

要想快速设计一个FIR滤波器,MATLAB可以说是一个最有利的帮手,当你知道较多滤波器知识和相关函数时,可以直接调用函数进行设计,灵活且方便。当你不知道时,也没有关系,打开MATLAB,在命令行输入:filterDesigner或者fdatool(老版MATLAB),你就能看到如下界面:

通过勾勾选选就能设计滤波器哦!

如果你又进步了一点,已经设计出一个滤波器,相看看性能怎么样,也简单,打开MATLAB,输入:fvtool(h),h就是你设计的滤波器的系数,你就会看到各种你想要的滤波器特性:幅值响应、相位响应、脉冲响应、零极点等等。

转载在知乎——膜拜能把知识理解的怎么透彻

原文地址:https://www.cnblogs.com/smallqing/p/10239966.html