如何快速設計一個FIR濾波器（一）真是通俗易懂，膜拜

時間 2020-05-07

標籤如何快速設計一個 fir 濾波器真是通俗易懂膜拜简体版

原文原文鏈接

在工做中，咱們最佩服的一羣人就是那種只用「紙和筆」就能把問題說清楚甚至解決的人，這須要超強的理論基礎以及模型抽象能力——一言不合就上公式，簡單、粗暴、有效。函數

今天，咱們也來裝一裝X，看看如何經過簡單「寫寫畫畫」來設計一個FIR濾波器。測試

濾波器的概念相比你們都很熟悉了，通常按照頻率特性能夠分爲低通、高通、帶通及帶阻濾波器，這是從輸出特性來講的。spa

設計常規濾波器的時候，咱們通常採用另一種分類，FIR（Finite Impulse Response）和IIR（Infinite Impulse Response）filter，即有限脈衝響應濾波器和無限脈衝響應濾波器。前面的文章中，咱們已經介紹了，理想脈衝信號，其傅里葉變換恆爲1，也就時包含了全部頻率份量，是一個理想的測試信號，可以激發出全部單位頻率份量的響應，所以理想脈衝信號的響應，就表明了系統的特性。.net

濾波器也能夠當作一個系統，若是用一個理想脈衝信號激勵，就會有輸出，咱們把輸出個數有限的稱爲有限脈衝響應濾波器（FIR）；輸出無限多的稱爲無限脈衝響應濾波器（IIR）。今天先說一下FIR。設計

1、Z傳遞函數的零點和極點表明什麼

咱們知道域的零極點表徵了系統的響應特性：極點表明了系統的模態，零點表明了系統能屏蔽的模態，具體參看文章3d

J Pan：如何入門自動控制理論zhuanlan.zhihu.comcode

在「如何理解離散傅里葉變換及z變換」一文中，orm

J Pan：如何理解離散傅里葉變換及Z變換zhuanlan.zhihu.comci

咱們介紹了域和域的關係： $z=e^{sT}$ ， $s=\sigma+j\omega$ ，因此get

當 $\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')