如何快速设计一个FIR滤波器（一）真是通俗易懂，膜拜

时间 2020-07-26 标签如何快速设计一个 fir 滤波器一真是通俗易懂膜拜

在工做中，咱们最佩服的一群人就是那种只用“纸和笔”就能把问题说清楚甚至解决的人，这须要超强的理论基础以及模型抽象能力——一言不合就上公式，简单、粗暴、有效。html

今天，咱们也来装一装X，看看如何经过简单“写写画画”来设计一个FIR滤波器。函数

滤波器的概念相比你们都很熟悉了，通常按照频率特性能够分为低通、高通、带通及带阻滤波器，这是从输出特性来讲的。post

设计常规滤波器的时候，咱们通常采用另一种分类，FIR（Finite Impulse Response）和IIR（Infinite Impulse Response）filter，即有限脉冲响应滤波器和无限脉冲响应滤波器。前面的文章中，咱们已经介绍了，理想脉冲信号，其傅里叶变换恒为1，也就时包含了全部频率份量，是一个理想的测试信号，可以激发出全部单位频率份量的响应，所以理想脉冲信号的响应，就表明了系统的特性。测试

滤波器也能够当作一个系统，若是用一个理想脉冲信号激励，就会有输出，咱们把输出个数有限的称为有限脉冲响应滤波器（FIR）；输出无限多的称为无限脉冲响应滤波器（IIR）。今天先说一下FIR。spa

1、Z传递函数的零点和极点表明什么

咱们知道域的零极点表征了系统的响应特性：极点表明了系统的模态，零点表明了系统能屏蔽的模态，具体参看文章设计

J Pan：如何入门自动控制理论zhuanlan.zhihu.com3d

在“如何理解离散傅里叶变换及z变换”一文中，code

J Pan：如何理解离散傅里叶变换及Z变换zhuanlan.zhihu.comorm

咱们介绍了域和域的关系： $z=e^{sT}$ ， $s=\sigma+j\omega$ ，因此htm

当 $\sigma=0$ 时， $s=j\omega$ ，对应的是域的虚轴，而此时 $z=e^{j\omega T}$ 对应的是单位圆，也就是说变换将域的虚轴映射成域的单位圆。

当 $\sigma>0$ 时， $s=\sigma +j\omega$ ，对应的是域的正半轴，而此时 $z=e^{\sigma}e^{j\omega T}$ ，因为 $e^{\sigma}>1$ ，也就是说此时变换将域正半轴映射到了域的单位圆外部。

当 $\sigma<0$ 时， $s=\sigma +j\omega$ ，对应的是域的负半轴，而此时 $z=e^{\sigma}e^{j\omega T}$ ，因为 $e^{\sigma}<1$ ，也就是说此时变换将域负半轴映射到了域的单位圆内部。

继续扩展， $z=e^{sT}=e^{j\omega/f_s}=e^{j2\pi f/f_s }$ ,很显然：

当时 ;
当 $f=\frac{1}{4}f_s$ 时 $z=e^{j\pi/2}=j$ ;
当 $f=\frac{1}{2}f_s$ 时 $z=e^{j\pi}=-1$ ;
当 $f=-\frac{1}{4}f_s$ 时 $z=e^{-j\pi/2}=-j$ ;

咱们知道在域上，虚轴上不一样的点对应不一样的频率，而域上单位圆与域虚轴对应，可见，域单位圆上不一样的点，表明了不一样的频率。

很容易获得，对于域的传递函数的零极点，也有和域零极点相似的结论：

规律1：若是在单位圆上有零点，则在零点所对应的频率上幅值响应为零；
规律2：对于不在单位圆上的零点，在单位圆上离零点最近的点对应的频率上幅值响应最小。
规律3：对于在单位圆内部的极点，在单位圆上离极点最近的点对应的频率上幅值响应最大。
规律4：若是极点和零点重合，对系统的频率响应没有影响。

在“如何理解离散傅里叶变换及z变换”一文中，咱们还介绍了若是一个信号的频谱以下：

频谱中最大的频率为 $f_{max}$ ，用一个周期为狄拉克梳状函数进行采采样后的频谱为原频谱的周期延拓，延拓的周期为采样周期，示意图以下：

也就是说，采样以后的频谱是一个周期函数，咱们把 $\left[ 0 \quad f_s \right]$ 称为主值区间，其中 $\left[ \frac{1}{2}f_s\quad f_s \right]$ 是 $\left[ -\frac{1}{2}f_s\quad 0 \right]$ 延拓一个周期得来的，与 $\left[ 0\quad \frac{1}{2}f_s \right]$ 彻底对称，所以咱们通常只考虑 $\left[ 0\quad \frac{1}{2}f_s \right]$ 区间，也就是半个单位圆的区域。

2、零、极点分布如何影响频率响应

咱们先看个简单的例子，熟悉一下套路。

Ex1： $H(z)=1-z^{-1}=\frac{z-1}{z}$

对于这个系统，在有一个极点，在时有一个零点。零、极点分布以下：

其中表示零点， $\times$ 表示极点。咱们来粗略分析一下这个系统的响应会是什么样。从也就是单位圆上角度为零（也是频率为零）的点开始，此处有一个零点，根据规律1，显然在频率为零时系统响应为零，顺着单位圆沿逆时针方向旋转，咱们离零点愈来愈远，零点的影响也愈来愈小，所以幅值响应会逐渐增大。当咱们到达 ,也就是频率为时，此时离零点最远，所以响应会达到一个最大值，当频率继续增大时，因为离零点又开始接近了，幅值响应又开始变小。

细心的童鞋可能发现了另一个端倪，你刚分析了零点，可系统明明还有一个极点啊！——没错，为你的细心点赞——咱们仔细观察，发现极点正好位于圆心位置，也就是说全部频率段离极点的距离都同样，所以能够认为都没影响。

用freqz函数将系统的频响画出来，就长成下图的样子，这也印证了咱们以前的分析，这个系统本质上是一个高通滤波器。

%% clc;clear; fs=1e4;b=[1 -1];a=[0 1]; [h,f] = freqz(b,a,2001,'whole',fs);N=round(0.5*length(h)); plot(f(1:N)/fs,20*log10(abs(h((1:N)))),'linewidth',2,'color','k'); ax = gca;ax.YLim = [-60 10];ax.XTick = 0:.1:0.5; xlabel('Normalized Frequency (f/fs)');ylabel('Magnitude (dB)'); grid on;title('Made by J Pan')

这个系统换作时域是什么样？

$H(z)=1-z^{-1}$

若为系统输出，为系统输入，则 $Y(z)=H(z)X(z)=(1-z^{-1})X(z)=X(z)-z^{-1}X(z)$

进行逆变换就能够获得：

这本质就是一个差分，对应连续系统的微分，咱们知道微分对应的是传递函数是，稳态时为 $s=j\omega$ ,这显然是一个高通滤波器，与前面的分析是一致的。

Ex2： $H(z)=1+z^{-1}=\frac{1+z}{z}$

很容易看出系统的零极点图以下：

显然，零点跑到了处，所以，系统的频响会先减少，到处达到最小值，而后又增长，具体频响以下图，这本质上是一个低通滤波器。

%% clc;clear; fs=1e4;b=[1 1];a=[0 1]; [h,f] = freqz(b,a,2001,'whole',fs);N=round(0.5*length(h)); plot(f(1:N)/fs,20*log10(abs(h((1:N)))),'linewidth',2,'color','k'); ax = gca;ax.YLim = [-60 10];ax.XTick = 0:.1:0.5; xlabel('Normalized Frequency (f/fs)');ylabel('Magnitude (dB)'); grid on;title('Made by J Pan')

很容易获得时域的表达式为：

这本质就是一个离散求和，对应连续系统的积分，咱们知道微分对应的是传递函数是，稳态时为 $1/s=1/{(j\omega)}$ ,这显然是一个低通滤波器，与前面的分析是一致的。

Ex3：假如咱们在0到 $\frac{1}{2}f_s$ 之间放置一个零点，那会不会是一个带阻滤波器呢？好比咱们想在频率在 $\frac{3f_s}{8}$ 这个点的系统频率响应为零。

频率 $\frac{3f_s}{8}$ 所在点对应的相角为 $\frac{3\pi}{4}$ ，由第一部分可知，频率响应在 $\left[ 0\quad \frac{1}{2}f_s \right]$ 与 $\left[ \frac{1}{2}f_s\quad f_s \right]$ 之间具备对称性，所以上述系统在相角为 $-\frac{3\pi}{4}$ 处也有一个零点。

转化成传递函数就是：

$H(z)= \left( z-e^{j\frac{3\pi}{4}} \right)\left( z-e^{-j\frac{3\pi}{4}} \right)$

展开能够得到：

$H(z)=z^2-2zcos(\frac{3\pi}{4})+1$

即 $H(z)=z^2+\sqrt{2}z+1$

%% clc;clear; fs=1e4;b=[1 sqrt(2) 1];a=[0 0 1]; [h,f] = freqz(b,a,2001,'whole',fs);N=round(0.5*length(h)); plot(f(1:N)/fs,20*log10(abs(h((1:N)))),'linewidth',2,'color','k'); ax = gca;ax.YLim = [-60 10];ax.XTick = 0:.1:0.5; xlabel('Normalized Frequency (f/fs)');ylabel('Magnitude (dB)'); grid on;title('Made by J Pan')

这个系统换作时域是什么样？

$H(z)=z^2+\sqrt{2}z+1$

若为系统输出，为系统输入，则

$Y(z)=H(z)X(z)=(z^2+\sqrt{2}z+1)X(z)$

进行逆变换就能够获得：

$y[n]=x[n+2]+\sqrt{2}x[n+1]+x[n]$

Ex4：

前面，咱们把零点和极点都放在了单位圆上，那能不能放在其余位置呢？——单位圆外面是不行的，由于外面对应着s域的正半轴，系统是不稳定；内部呢？我不妨把零点先放在x轴上试试，放在这个点上。

粗略分析，当时（对应频率为零）离零点最近，此时频率响应应该最小，但不为零。当时（对应 $\frac{1}{2}f_s$ ）离零点最远，响应应该到达最大值。

可见，零极点的位置决定了系统在不一样频率下的响应状况。

%% clc;clear; fs=1e4;b=[1 -0.5];a=[0 1]; [h,f] = freqz(b,a,2001,'whole',fs); N=round(0.5*length(h)); plot(f(1:N)/fs,20*log10(abs(h((1:N)))),'linewidth',2,'color','k'); ax = gca;ax.YLim = [-7 5];ax.XTick = 0:.1:0.5; xlabel('Normalized Frequency (f/fs)'); ylabel('Magnitude (dB)'); grid on;title('Made by J Pan')

Ex5：

这个传递函数有点意思了，它有6个根——都是复数哦！

咱们能够将上述方程写成以下格式： $z^6=e^{j2\pi n} \quad n=0,1,2,3,4,5$

因此解为： $z=e^{j\frac{2\pi n}{6}}$

总共有6个根均布在单位圆上，以下图：

咱们能够画出以下的频率响应，可见其本质是一个多个带阻的滤波器。这种滤波器有啥用呢？咱们知道，市电频率是50Hz，其带来的干扰通常就是50Hz其整数倍谐波100Hz、150Hz，200Hz等，选择这种数字滤波器就能够消除相似于市电50Hz带来的噪声影响。

%% clc;clear; fs=1e4;b=[1 0 0 0 0 0 -1];a=[0 0 0 0 0 0 -1]; [h,f] = freqz(b,a,2001,'whole',fs); N=round(0.5*length(h)); plot(f(1:N)/fs,20*log10(abs(h((1:N)))),'linewidth',2,'color','k'); ax = gca;ax.YLim = [-30 10];ax.XTick = 0:.1:0.5; xlabel('Normalized Frequency (f/fs)'); ylabel('Magnitude (dB)'); grid on;title('Made by J Pan')

Ex6： $H(z)=\frac{z^6-1}{z-1}$

对于该函数，其零点位于处，6个零点均匀分布在单位圆上。

另外，在处还有一个极点，与该处的零点重合，零极点以下图所示。

可见，在处因为零极点重合，并未对系统产生影响，也就是说，若是想消除某零点给系统带来的影响，咱们能够再该位置同时也放置一个极点；反之亦然。

观察一下与Ex5的频响的区别，是否是颇有意思？

clc;clear; fs=1e4;b=[1 0 0 0 0 0 -1];a=[0 0 0 0 0 1 -1]; [h,f] = freqz(b,a,2001,'whole',fs); N=round(0.5*length(h)); plot(f(1:N)/fs,20*log10(abs(h((1:N)))),'linewidth',2,'color','k'); ax = gca;ax.YLim = [-30 20];ax.XTick = 0:.1:0.5; xlabel('Normalized Frequency (f/fs)'); ylabel('Magnitude (dB)'); grid on;title('Made by J Pan')

3、如何简单粗暴的设计一个FIR滤波器

前面套路差很少说完了——有高通、低通、带阻滤波器，好像尚未带通滤波器，下面咱们就拿带通滤波器来练练手。

要求以下：在125Hz时频率响应达最大值，采样频率。

因为125Hz是采样频率1000Hz的1/8，咱们先均匀布置8个零点在单位圆上，这样就能保证有一个零点是位于125Hz处的。

咱们知道，零点位置时频率响应最小的点，咱们如今是想要在125Hz（相角 $j\frac{\pi}{4}$ ）处响应最大，貌似有矛盾——不要紧，咱们能够再放个极点，将这个零点抵消掉（规律4）！

因为对称性呢，在-125Hz（相角 $-j\frac{\pi}{4}$ ）也要放置一个极点。最终零、极点分布以下图：

由Ex5可知，均匀分布的8个点，对应的传递函数的分子为，两个极点对应传递函数的分母为 $\left(z-e^{j\frac{pi}{4}} \right)\left(z-e^{-j\frac{pi}{4}} \right)$ ，因此总的传递函数为：

$H(z)=\frac{z^8-1}{\left(z-e^{j\frac{pi}{4}} \right)\left(z-e^{-j\frac{pi}{4}} \right)}$

化简一下：

$H(z)=\frac{Y(z)}{X(z)}=\frac{z^8-1}{z^2-\sqrt{2}z+1}$

这就是咱们要的传递函数了，换成时域是什么样呢？

$z^2Y(z)-\sqrt{2}zY(z)+Y(z)=z^8X(z)-X(z)$

逆变换一下：

$y[n+2]-\sqrt{2}y[n+1]+y[n]=x[n+8]-x[n]$

平移一下：

$y[n]-\sqrt{2}y[n-1]+y[n-2]=x[n+6]-x[n-2]$

即

$y[n]=\sqrt{2}y[n-1]-y[n-2]+x[n+6]-x[n-2]$

细心地童鞋可能注意到：是将来的数啊，这显然不太合理，那怎么办呢？——还记得Ex1吗？咱们能够再圆点处加极点啊，当其余零极点都在单位圆上时不影响频率响应，以下图：

则传递函数变成了：

$H(z)=\frac{z^8-1}{z^6\left(z-e^{j\frac{pi}{4}} \right)\left(z-e^{-j\frac{pi}{4}} \right)}$

最终时域关系变为了：

$y[n]=\sqrt{2}y[n-1]-y[n-2]+x[n]-x[n-8]$

%% clc;clear; fs=1e4;b=[1 0 0 0 0 0 0 0 -1];a=[1 -sqrt(2) 1]; [h,f] = freqz(b,a,2001,'whole',fs); N=round(0.5*length(h)); plot(f(1:N)/fs,20*log10(abs(h((1:N)))),'linewidth',2,'color','k'); ax = gca;ax.YLim = [-30 20];ax.XTick = 0:.1:0.5; xlabel('Normalized Frequency (f/fs)'); ylabel('Magnitude (dB)'); grid on;title('Made by J Pan')

怎么样，读到这的时候你是否是已经开始跃跃欲试了？那就开始拿起笔和纸，试试吧！

转载自知乎——讲的真是通俗易懂，膜拜啊。

转载于:https://www.cnblogs.com/smallqing/p/10239753.html