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

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

濾波器的概念相比大家都很熟悉了，一般按照頻率特性可以分為低通、高通、帶通及帶阻濾波器，這是從輸出特性來說的。

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

濾波器也可以看成一個系統，如果用一個理想脈衝信號激勵，就會有輸出，我們把輸出個數有限的稱為有限脈衝響應濾波器（FIR）；輸出無限多的稱為無限脈衝響應濾波器（IIR）。今天先說一下FIR。

一、Z傳遞函數的零點和極點代表什麼

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

J Pan：如何入門自動控制理論?

zhuanlan.zhihu.com圖標

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

J Pan：如何理解離散傅里葉變換及Z變換?

zhuanlan.zhihu.com圖標

我們介紹了 z 域和s 域的關係：z=e^{sT} ， s=sigma+jomega ，所以

當 sigma=0 時，s=jomega ，對應的是 s 域的虛軸，而此時 z=e^{jomega T} 對應的是單位圓，也就是說z 變換將 s 域的虛軸映射成 z 域的單位圓。

當 sigma>0 時，s=sigma +jomega ，對應的是 s 域的正半軸，而此時 z=e^{sigma}e^{jomega T} ，由於 e^{sigma}>1，也就是說此時z 變換將s 域正半軸映射到了z 域的單位圓外部。

當 sigma<0 時，s=sigma +jomega ，對應的是 s 域的負半軸，而此時 z=e^{sigma}e^{jomega T} ，由於 e^{sigma}<1，也就是說此時z 變換將s 域負半軸映射到了z 域的單位圓內部。

繼續擴展， z=e^{sT}=e^{jomega/f_s}=e^{j2pi f/f_s } ,很顯然：

  • 當 f=0 時 z=1 ;
  • 當 f=frac{1}{4}f_s 時 z=e^{jpi/2}=j ;
  • 當 f=frac{1}{2}f_s 時 z=e^{jpi}=-1 ;
  • 當 f=-frac{1}{4}f_s 時 z=e^{-jpi/2}=-j ;

我們知道在 s 域上，虛軸上不同的點對應不同的頻率，而 z 域上單位圓與 s 域虛軸對應，可見， z 域單位圓上不同的點，代表了不同的頻率。

很容易得到，對於 z 域的傳遞函數的零極點，也有和 s 域零極點類似的結論：

  • 規律1：如果在單位圓上有零點，則在零點所對應的頻率上幅值響應為零；
  • 規律2：對於不在單位圓上的零點，在單位圓上離零點最近的點對應的頻率上幅值響應最小。
  • 規律3：對於在單位圓內部的極點，在單位圓上離極點最近的點對應的頻率上幅值響應最大。
  • 規律4：如果極點和零點重合，對系統的頻率響應沒有影響。

在「如何理解離散傅里葉變換及z變換」一文中，我們還介紹了如果一個信號的頻譜如下：

頻譜中最大的頻率為 f_{max} ，用一個周期為 f_s 狄拉克梳狀函數進行采採樣後的頻譜為原頻譜的周期延拓，延拓的周期為採樣周期 f_s ，示意圖如下：

也就是說，採樣之後的頻譜是一個周期函數，我們把 left[ 0 quad f_s ight]稱為主值區間，其中 left[ frac{1}{2}f_squad f_s ight] 是 left[ -frac{1}{2}f_squad 0 ight] 延拓一個 f_s 周期得來的，與 left[ 0quad frac{1}{2}f_s ight] 完全對稱，因此我們一般只考慮 left[ 0quad frac{1}{2}f_s ight] 區間，也就是半個單位圓的區域。

二、零、極點分布如何影響頻率響應

我們先看個簡單的例子，熟悉一下套路。

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

對於這個系統，在 z=0 有一個極點，在 z=1 時有一個零點。零、極點分布如下：

其中 o 表示零點， imes 表示極點。我們來粗略分析一下這個系統的響應會是什麼樣。從 z=1 也就是單位圓上角度為零（也是頻率為零）的點開始，此處z=1有一個零點，根據規律1，顯然在頻率為零時系統響應為零，順著單位圓沿逆時針方向旋轉，我們離零點越來越遠，零點的影響也越來越小，因此幅值響應會逐漸增大。當我們到達 z=-1 ,也就是頻率為 1/2f_s 時，此時離零點最遠，因此響應會達到一個最大值，當頻率繼續增大時，由於離零點又開始接近了，幅值響應又開始變小。

細心的童鞋可能發現了另外一個端倪，你剛分析了零點，可系統明明還有一個極點啊！——沒錯，為你的細心點贊——我們仔細觀察，發現極點正好位於圓心位置，也就是說所有頻率段離極點的距離都一樣，因此可以認為都沒影響。

用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) 為系統輸出， X(z) 為系統輸入，則 Y(z)=H(z)X(z)=(1-z^{-1})X(z)=X(z)-z^{-1}X(z)

進行逆變換就可以得到：

y[n]=x[n]-x[n-1]

這本質就是一個差分，對應連續系統的微分，我們知道微分對應的是傳遞函數是 s ，穩態時為 s=jomega ,這顯然是一個高通濾波器，與前面的分析是一致的。

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

很容易看出系統的零極點圖如下：

顯然，零點跑到了 1/2f_s 處，因此，系統的頻響會先減小，到 1/2f_s 處達到最小值，然後又增加，具體頻響如下圖，這本質上是一個低通濾波器。

%%
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)

很容易得到時域的表達式為： y[n]=x[n]+x[n-1]

這本質就是一個離散求和，對應連續系統的積分，我們知道微分對應的是傳遞函數是 1/s ，穩態時為 1/s=1/{(jomega)} ,這顯然是一個低通濾波器，與前面的分析是一致的。

Ex3：假如我們在0到 frac{1}{2}f_s 之間放置一個零點，那會不會是一個帶阻濾波器呢？比如我們想在頻率在 frac{3f_s}{8} 這個點的系統頻率響應為零。

頻率frac{3f_s}{8} 所在點對應的相角為 frac{3pi}{4} ，由第一部分可知，頻率響應在 left[ 0quad frac{1}{2}f_s ight] 與 left[ frac{1}{2}f_squad f_s ight] 之間具有對稱性，因此上述系統在相角為 -frac{3pi}{4} 處也有一個零點。

轉化成傳遞函數就是：

H(z)= left( z-e^{jfrac{3pi}{4}} ight)left( z-e^{-jfrac{3pi}{4}} ight)

展開可以獲得：

H(z)=z^2-2zcos(frac{3pi}{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) 為系統輸出， X(z) 為系統輸入，則

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：H(z)=z-0.5

前面，我們把零點和極點都放在了單位圓上，那能不能放在其他位置呢？——單位圓外面是不行的，因為外面對應著s域的正半軸，系統是不穩定；內部呢？我不妨把零點先放在x軸上試試，放在 z=0.5 這個點上。

粗略分析，當 z=1 時（對應頻率為零）離零點最近，此時頻率響應應該最小，但不為零。當 z=-1 時（對應 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：H(z)=z^6-1

這個傳遞函數有點意思了，它有6個根——都是複數哦！

我們可以將上述方程寫成如下格式： z^6=e^{j2pi n} quad n=0,1,2,3,4,5

所以解為： z=e^{jfrac{2pi 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}

對於該函數，其零點位於 z^6-1=0 處，6個零點均勻分布在單位圓上。

另外，在 z=1 處還有一個極點，與該處的零點重合，零極點如下圖所示。

可見，在 z=1 處由於零極點重合，並未對系統產生影響，也就是說，如果想消除某零點給系統帶來的影響，我們可以再該位置同時也放置一個極點；反之亦然。

觀察一下與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)

三、如何簡單粗暴的設計一個FIR濾波器

前面套路差不多說完了——有高通、低通、帶阻濾波器，好像還沒有帶通濾波器，下面我們就拿帶通濾波器來練練手。

要求如下：在125Hz時頻率響應達最大值，採樣頻率 f_s=1000Hz 。

由於125Hz是採樣頻率1000Hz的1/8，我們先均勻布置8個零點在單位圓上，這樣就能保證有一個零點是位於125Hz處的。

我們知道，零點位置時頻率響應最小的點，我們現在是想要在125Hz（相角 jfrac{pi}{4} ）處響應最大，貌似有矛盾——沒關係，我們可以再放個極點，將這個零點抵消掉（規律4）！

由於對稱性呢，在-125Hz（相角 -jfrac{pi}{4} ）也要放置一個極點。最終零、極點分布如下圖：

由Ex5可知，均勻分布的8個點，對應的傳遞函數的分子為 z^8-1 ，兩個極點對應傳遞函數的分母為 left(z-e^{jfrac{pi}{4}} ight)left(z-e^{-jfrac{pi}{4}} ight) ，所以總的傳遞函數為：

H(z)=frac{z^8-1}{left(z-e^{jfrac{pi}{4}} ight)left(z-e^{-jfrac{pi}{4}} ight)}

化簡一下：

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]

細心地童鞋可能注意到： x[n+6] 是未來的數啊，這顯然不太合理，那怎麼辦呢？——還記得Ex1嗎？我們可以再圓點處加極點啊，當其他零極點都在單位圓上時不影響頻率響應，如下圖：

則傳遞函數變成了：

H(z)=frac{z^8-1}{z^6left(z-e^{jfrac{pi}{4}} ight)left(z-e^{-jfrac{pi}{4}} ight)}

最終時域關係變為了：

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)

怎麼樣，讀到這的時候你是不是已經開始躍躍欲試了？那就開始拿起筆和紙，試試吧！


推薦閱讀：
相关文章