雙二階濾波器之MATLAB設計及C語言實現

搬運自我的CSDN https://blog.csdn.net/u013213111/article/details/90529164

參考：
雙二階濾波器

本文中的例子和代碼放在Github

First，什么是雙二階濾波器？wiki上是這么說的：二階、遞歸、線性，含有兩個極點和兩個零點，“雙二階”的名字來源於它的傳遞函數是兩個二次多項式的比值。

In signal processing, a digital biquad filter is a second order recursive linear filter, containing two poles and two zeros. "Biquad" is an abbreviation of "biquadratic", which refers to the fact that in the Z domain, its transfer function is the ratio of two quadratic functions: H(z)=(b₀+b₁z⁻¹+b₂z⁻²)/(a₀+a₁z⁻¹+a₂z⁻²) The coefficients are often normalized such that a₀ = 1: H(z)=(b₀+b₁z⁻¹+b₂z⁻²)/(1+a₁z⁻¹+a₂z⁻²) High-order IIR filters can be highly sensitive to quantization of their coefficients, and can easily become unstable.

歸一化傳遞函數寫成這樣：

\[H(z)= \frac{b₀+b₁z⁻¹+b₂z⁻²}{1+a₁z⁻¹+a₂z⁻²} \]

用MATLAB的Filter Designer來設計一個：400Hz帶通IIR，需要用4個Sections來實現，默認給出的濾波器結構是Direct-Form II。

在菜單欄的Analysis中選擇Filter Coeffients就能看到濾波器系數了：

Numerator，分子，也就是傳遞函數中的b項們，從上到下依次為\(b_0\)，\(b_1\)和\(b_2\)。
Denominator，分母，也就是傳遞函數中的a項，從上到下依次為\(a_0\)，\(a_1\)和\(a_2\)，其中\(a_0\)總是為1。
Gain，增益。

用數組來存放濾波器系數：

//8-order IIR filter with 4 sections
const int sections = 4;

//nominator
const float b[4][3] = {
	{ 1, -1.984454421, 1 },
	{ 1, -1.999405318, 1 },
	{ 1, -1.993167556, 1 },
	{ 1, -1.998644244, 1 }
};

//denominator
const float a[4][3] = {
	{ 1, -1.984532740, 0.9884094716 },
	{ 1, -1.988571923, 0.9909378613 },
	{ 1, -1.991214225, 0.9962624248 },
	{ 1, -1.995917854, 0.9977478940 }
};

const float gain[4] = { 0.583805364, 0.583805364, 0.170388576, 0.170388576 };

以Direct-Form II為例來寫個實現代碼，計算過程是：

\[y[n]=b_{0}w[n]+b_{1}w[n-1]+b_{2}w[n-2] \]

其中$$w[n]=x[n]-a_{1}w[n-1]-a_{2}w[n-2]$$
代碼如下：
用一個數組來存放濾波器的中間狀態量\(w[n-1]\)和\(w[n-2]\)：

float w[sections][2]; //filter states

在濾波計算之前，初始化濾波器：

for (int i = 0; i < sections; i++) {
		w[i][0] = 0; //w[n-1]
		w[i][1] = 0; //w[n-2]
}

正式開始計算：pcmIn[i]是原始的輸入信號，輸入section 1，section 1的輸出則作為section 2的輸入，以此類推，由於這個濾波器由4個sections構成，因此要循環4次。
注意輸出要乘gain。

y[0] = pcmIn[i];
		for (j = 0; j < sections; j++) {
			tmp[j] = y[j] - a[j][1] * w[j][0] - a[j][2] * w[j][1]; //calculate w[n]
			y[j+1] = tmp[j] + b[j][1] * w[j][0] + b[j][2] * w[j][1]; //calculate the j-th section filter output y[n]
			w[j][1] = w[j][0]; //move w[n-1] -> w[n-2]
			w[j][0] = tmp[j]; //move w[n] -> w[n-1]
			y[j+1] = gain[j] * y[j+1]; //multiply with gain
		}
			
out = y[j];

如果需要得到例如PCM16的輸出，那么再對out進行限幅，that's all.