之前搞了一些數字訊號處理演算法程式設計（OC），一直沒來得及整理，現在整理一下，包括FFT、巴特沃斯濾波器（高通、低通、帶通、帶阻）、資料差值（線性、sinc、三次樣條*）、資料壓縮（等距、平均、峰值檢測）和模仿matlab的STFT功能（spectrogram函式三維繪圖）。

注：我比較喜歡使用matlab（也喜歡自己修改matlab的原始碼做測試，所以重灌了好幾次了，囧。。。），有部分參考matlab的演算法進行（或者直接翻譯為OC），演算法的執行結果經常會跟matlab運算結果進行比較，演算法只做學習用，轉載請註明。

另外可能會有不足或者理解偏差的地方，路過的高人請不吝賜教。

今天來說一下數字濾波器的程式碼實現（IIR）。

---------------------------------------------------------------------------------------

一、基本概念

1.數字濾波器（DF）和模擬濾波器（AF）一樣，都是用來濾波的，它將訊號的某些頻率（頻段）的訊號加以放大，而將另外一些頻率（頻段）的訊號加以抑制。藉助A/D、D/A轉換器，數字濾波器可以處理模擬訊號也可以輸出模擬訊號。

2.數字濾波器有4種表示方法

①線性差分方程：也就是它的濾波作用的基本構成事數值運算的部件：相加器、相乘器、延時器。而模擬濾波器的基本部件則是電感器、電容器、電阻器及有源器件。

②系統函式。將a0歸一化為a0=1，即用z^(-1)或z的有理分式表示系統函式。則有

（IIR數字濾波器）

③單位抽樣響應h(n)拉氏變換得到

④線性訊號流圖。系統函式如下

其直接II型結構圖如下

3.分類

（1）按衝激響應分為無限長衝激響應（IIR）數字濾波器和有限長衝激響應（FIR）數字濾波器。

IIR：從2②中的IIR系統函式可以看出，它是一個z^(-1)的有理分式，包括有輸出到輸入的反饋網路結構，此係統分母多項式A(z)決定了反饋網路，同時確定了有限z平面的極點，而分子多項式B(z)決定了正饋網路，同時確定了有限z平面的零點。

FIR：FIR濾波器系統函式可表示為z^(-1)的多項式，即IIR系統函式中ai=0（分母為1）,可以寫成如下形式

其中h(n)時系統的單位衝激響應，顯然h(i)=bi（當ai=0時），H(z)在有限z平面只有零點，如果是因果系統則全部極點都在z=0處。系統不存在反饋網路。

（2）按濾波器幅度響應可分為低通、高通、帶通、帶阻、全通等。按照奈奎斯特抽樣定理，訊號最高頻率fh只能限於fh≤fs/2（fs為臭氧頻率），即wh≤π。與模擬濾波器不同，數字濾波器頻率響應是以2π為週期的周期函式，這些理想的幅度響應特性（如下圖所示），即“突變”型幅度響應會造成無限長的非因果的單位衝激響應，是不可實現的，只能用可實現的實際濾波來逼近它。

（3）按相位響應分類：線性相位、非線性相位。如果要求嚴格的線性相位，則必須使用FIR線性相位濾波器。

（4）按特殊要求分類：最小相位滯後濾波器、梳狀濾波器、陷波器、全通濾波器、諧振器，甚至包括波形產生器等。採用零極點的適當配置方法，可以得到這些濾波器。

二、IIR數字濾波器的技術指標

以數字低通濾波器為例（見下圖），指標包括：通帶截止頻率wp，阻帶截止頻率wst， 通帶波紋δ1（Rp(dB)）， 阻帶波紋δ2（As(dB)） 。

通帶允許的最大衰減Rp(dB)以及阻帶應達到的最小衰減As(dB)範圍如下（已歸一化為H(e^(jw))=1）：

注：

①高通與低通的指標一致；

②帶通的技術指標通帶截止頻率分通帶上截止頻率wp2和通帶下截止頻率wp1，阻帶截止頻率分上阻帶截止頻率wst2和下阻帶截止頻率wst1；

③帶阻的技術指標通帶截止頻率分為上通帶截止頻率wp2和下通帶截止頻率wp1，阻帶截止頻率分阻帶上截止頻率wst2和阻帶下截止頻率wst1。

三、濾波器的設計步驟

濾波器的設計就是要找到一個滿足技術指標要求的可實現的因果穩定的數字濾波器來逼近理想的濾波器幅度特性。

1.濾波器設計思路

本文采用間接法設計數字濾波器（先設計模擬低通濾波器在通過雙線性變換法得到數字低通、高通、帶通、帶阻濾波器），設計模型如下圖所示。引數說明：fp——通帶截止頻率，fs——阻帶截止頻率，rp通帶最大衰減，rs阻帶最小衰減，N——濾波器階數，fc——3dB截止頻率，sa——原型濾波器分母多項式的係數，sb原型濾波器分子多項式的係數，za——數字濾波器系統函式分母多項式的係數，zb——數字濾波器系統函式分子多項式的係數，filter——濾波方法（將訊號代入，返回濾波後的訊號）。

即，數字濾波器的設計及濾波思路圖如下圖：

2.模擬濾波器設計

（1）模擬濾波器的設計步驟

①給定模擬濾波器技術指標Ωp，Ωst，(1-δ1)（或Rp dB），δ2（或As dB）；

②計算濾波器所需階數N；

注：本文設計的濾波器階數N不能大於10，因為階數越高計算出的值越不準確。

③查表確定歸一化低通濾波器系統函式Has(s)；

④將Has(s)轉換為所需型別的低通濾波器系統函式Ha(s)。

（2）具體實現

①確定濾波器的階數N，公式如下所示，計算得到的結果向上取整。注意此處的Ωst和Ωp為預畸後的值（預畸公式見上思路設計圖中公式）

②確定3dB頻率，公式如下，當選Ωc=(Ωcp+Ωcs)/2時通帶阻帶衰減皆可超過要求（檢視matlab原始碼發現採用Ω=Ωcs，因此本文采取Ω=Ωcs）

③根據N的值，查表得歸一化低通濾波器系統函式Has(s)的分母多項式係數，公式如下：

d0一般由Ω=0時|Han(j0)|=1（增益為1）來確定，由於a0=1，一次d0=a0=1。

附歸一化巴特沃斯濾波器分母多項式的係數表如下：

N    a0    a1    a2    a3    a4    a5    a6    a7    a8    a9
1    1    1
2    1    1.4142136
3    1    2    2
4    1    2.6131259    3.4142136    2.6131259
5    1    3.236068    5.236068    5.236068    3.236068
6    1    3.8637033    7.4641016    9.1416202    7.4641016    3.8637033
7    1    4.4939592    10.0978347    14.5917939    14.5917939    10.0978347    4.4939592
8    1    5.1258309    13.1370712    21.846151    25.6883559    21.846151    13.1370712    5.1258309
9    1    5.7587705    16.5817187    31.1634375    41.9863857    41.9863857    31.1634375    16.5817187    5.7587705
10    1    6.3924532    20.4317291    42.8020611    64.8823963    74.2334292    64.8823963    42.8020611    20.4317291    6.3924532

④轉換成所需型別的低通濾波器系統函式Ha(s)，去歸一化另s=s/Ωc。

3.模擬濾波器的數字化知識

模擬濾波器的數字化方法有很多種，包括衝激響應不變法（脈衝響應不變法）、階躍響應不變法和雙線性變換法等。本文采用雙線性變換法實現。

（1）基本思路

雙線性變換法是使數字濾波器的頻率響應與模擬濾波器的頻率響應相似的一種變換，它使得Ω和ω之間是單值對映關係可以避免頻率響應的混疊失真。

下圖表示了變換思路，即將s平面整個變換到一箇中介s1平面的一個窄帶Ω：-π/T→π/T之中，然後經過z=e^(s1*T)的變換，將s1平面對映到z平面。從s1到z平面的變換是單值變換，從而使整個變換過程（s到z）成為單值的變換。

（2）變換公式

低通、高通、帶通、帶阻的變換關係式都在本節第一模組的“數字濾波器的設計及濾波思路圖”中展現了。

4.濾波方法

如前所述，濾波過程就是解常係數線性差分方程的過程，形式如：

其中，x(n)序列為濾波前的訊號序列，ak, bm為H(z)系統函式分母與分子的系統陣列，求出的y(n)即為濾波後的訊號序列。

注：x(n)與y(n)的長度要相等，且a0=1。

公式的化簡工程如下：

預設條件，當k<0時，x(k), y(k)都為0。例如n=0時，y(0)=b0*x0+b1*x(-1)+...+bM*x(-M)-a1*y(-1)-...-aN*y(-N)=b0*x(0)。經過迭代，可以求出y(n)序列的所有值。

========================================================

OK，理論到此結束，下一節上程式碼！
轉：https://blog.csdn.net/shengzhadon/article/details/46784509