淺談頻域數字濾波器

01-28

1.系統帶寬與雜訊

大家都知道反饋是閉環控制系統的命根子，然而在實際應用中，反饋信號總是伴隨著雜訊。搞控制的人都對卡爾曼濾波這種時域濾波器耳熟能詳，反倒是對頻域濾波器了解可能不是太好。其實在現實中，對於有效信號和雜訊的頻率特性區別明顯的時候，頻域濾波器可以是一個好的選擇。

例如下圖：

期望的信號是頻率為5Hz的正弦波，但是混雜了雜訊，我們來看一下它的頻譜分析：

原來它在15Hz和20Hz的頻率上還有著明顯的雜訊，那麼我們在這個信號之後，再設計一個類似於這樣的低通濾波器：

不就可以把高頻的雜訊給過濾了嗎？

在控制系統的反饋中，一般會結合被控對象的帶寬來設置反饋信號的濾波器截止頻率，高於系統帶寬的信號往往是無意義的雜訊，比如我以前寫的寫過的：

控制研究中的『帶寬』怎樣理解？

以及現實中如何標定系統的帶寬

如何用實驗方法確定系統的頻域響應 - 知乎專欄

一般來說，宏觀機械系統的帶寬在50Hz以下，微機電系統可以達到100KHz，電路系統的帶寬可達100MHz甚至10GHz以上。

2.從滑均濾波到FIR濾波器

問題是上面的濾波器該如何數字實現呢？凡是有過工程經驗的人都有過一個樸素的概念：多采幾個樣，求平均值，也就是常說的「滑動平均濾波」

比如上圖，就是根據三個歷史值和當前值，加起來求一個平均來當做當前濾波的結果。

也就是

$y(n)=(u(n)+u(n-1)+u(n-2)+u(n-3))/4$

用c語言偽代碼可以寫一個簡單的程序

double AveragingFilter(double dataN)

{

static double dataHistory[3]=0;

double result=0;

result=(dataN+dataHistory[0]+dataHistory[1]+dataHistory[2])/4;

dataHistory[2]=dataHistory[1];

dataHistory[1]=dataHistory[0];

dataHistory[0]=dataN;

return result;

}

相信這是很多人都會做的，但是我們換一種角度來看它，對上述差分方程做Z變換

$frac{Y(z)}{U(z)}=(1+z^{-1}+z^{-2}+z^{-3})/4$

我們可以看到它的頻域響應（假定採樣頻率/執行頻率為1KHz）：

滑動平均濾波器實質上是一個低通濾波器。

我們再來看更一般N階的滑動平均濾波器的情況：

$frac{Y(z)}{U(z)}=(1+z^{-1}+z^{-2}+z^{-3}+.....+z^{-n-1})/n$

它的截止頻率為Fs/n。其中Fs為採樣頻率。再來看一個n=10的例子：

用信號處理的語言來說，這叫做有限衝擊響應（Finite impulse response，FIR）濾波器，即系統輸出和過去的有限個的輸入有關係。從控制人的角度，也可以理解為這是一個純零點系統，沒有極點，再換句話說，這個濾波器的系統沒有反饋。

但是滑動平均濾波器並不能算一個好的低通濾波器，我們能否通過改變它的零點，來設計更好的低通濾波器呢？

可以。在matlab命令行輸入fadtool，我們就可以方便的設計低通濾波器了。

我們想要lowpass filter，FIR類型，再設定好採樣頻率，截止頻率，我們就可以設計出一個FIR濾波器了。再點擊Analysis-Coefficients，就可以看到這個濾波器的係數

知道了z傳遞函數的係數，編程實現可參照

如何將連續域的控制器進行離散化實現 - 知乎專欄

其中的道理是一樣的，在fadtool里也可以將它的係數通過c語言頭文件的形式導出來，方便編程。

FIR由於只有前饋通道，沒有反饋，可以不擔心穩定性問題，這使得它在處理高速信號的時候具有優勢。

但是，FIR是非常消耗硬體資源的。

上面那個濾波器要101階！手寫程序寫的都得累了吧？單片機算101次浮點乘法和加法需要多久？我還手寫過一個21階的FIR濾波器，用FPGA並行實現的：

消耗了我非常多的乘法器和加法器。而且，截止頻率和採樣頻率差的越遠，階數越高，再看一個例子：

502階！我那麼貴的開發板，一共才840個乘法器和加法器。這時候，是該談談IIR濾波器了。

3.反饋系統-IIR濾波器

其實最簡單的IIR濾波器我們搞控制的一點也不陌生，一階慣性環節就是

$G(s)=frac{1}{Ts+1}$

我們看它的bode圖（T=0.1），

就是截止頻率w約等於1/T的低通濾波器么把它離散化實現了不就完了。但是再仔細看看它的相位響應，在控制系統裡面，你確定這麼大的相位延遲它不會失穩？

這時候就需要更高階的設計方法了，濾波器的差分方程/時域表達式形式就變成了：

$y(n)=(a_0u(n)+a_0u(n-1)+...+a_mu(m))+(b_1y(n-1)+b_2y(n-2)+...+b_ny(n))$

在此形式下，看起來系統當前時刻的輸出不僅和之前時刻的輸入有關，而且還和歷史輸出有關，看起來當前的結果是迭代了無數過去的歷史狀態，所以叫做無窮脈衝響應(infinite impulse response，IIR) 濾波器。

即引入了反饋的形式，寫成Z函數

$frac{Y(z)}{U(z)}=frac{a_0+a_1z^{-1}+..+a_mz^{-m}}{1+b_1z^{-1}+..+a_bz^{-n}}$

同樣可用fdatool來幫助設計，記得把設計方法換成IIR

階數只剩下15階了，消耗的資源大為減少。但是IIR是反饋系統，在高速的信號處理中，即使是使用並行實現，由於計算延遲和有限字長來造成的穩定性仍然是非常有挑戰性的。