本人本科渣渣一個，前兩天導師讓我設計一個數字濾波器。由於本人基本沒有數字訊號處理基礎，於是只能依靠百度和matlab，折騰了半天總算是摸索明白了。百度上有一些文章不靠譜，很容易誤導別人，故在此發一篇部落格。

濾波器設計目標：設計一個1Hz截止頻率的2階低通巴特沃斯數字濾波器，並轉化成C語言函式。（國標裡提的要求）

濾波器指標：指標：截止頻率Fc = 1Hz，階數N=2，低通巴特沃斯濾波器，取樣頻率Fs = 15Hz。

一、Matlab計算濾波器係數

Matlab計算巴特沃斯低通濾波器係數過程如下：

①根據給定的通帶截止頻率、通帶截止增益、阻帶截止頻率、阻帶截止增益，利用buttord函式計算巴特沃斯濾波器所需的最小階數和截止頻率。（由於要求中已給出截止頻率，故這步省去）

（圖是網上扒拉的，指標與本設計不符）

②根據上述計算得到的階數，利用buttap函式計算出巴特沃斯濾波器原型。

③利用lp2lp函式，將原型濾波器轉換成目標截止頻率的濾波器。

④利用脈衝響應不變法（impinvar函式）或是雙線性變換法（bilinear函式）將模擬濾波器轉換為數字濾波器。數字濾波器形式為z域有理函式，分子分母系數即為濾波器係數。

（從左到右：濾波器原型、模擬濾波器、數字濾波器）

設計過程matlab原始碼如下：

Fs = 15;        %取樣頻率
Nn = 12800;
N = 2;          %階數
Wc = 1*2*pi;    %截止頻率
[z,p,k] = buttap(N);  %計算巴特沃斯濾波器原型
[Bap,Aap] = zp2tf(z,p,k);  %轉換成多項式模式
[b,a] = lp2lp(Bap,Aap,Wc);  %根據截止頻率計算模擬巴特沃斯濾波器係數
[bz,az] = impinvar(b,a,Fs);  %用脈衝響應不變法離散化
figure(1)
[H,W] = freqz(bz,az,Nn,Fs);  %繪製頻率特性曲線
subplot(2,1,1)
plot(W,20*log10(abs(H)));
grid on;
subplot(2,1,2)
plot(W,180/pi*unwrap(angle(H)));
grid on;

二、Matlab計算驗證

先在Matlab中驗證濾波函式。先編寫帶噪聲的輸入函式，然後經過濾波器函式後，觀察濾波效果。其中濾波器函式寫法為：

Filter函式為Matlab自帶函式，其演算法為：

其中，a即為z域傳遞函式的分母系數，b為分子係數。例如本應用中：

則演算法為az(1)*y(k) = bz(1)*x(k) + bz(2)*x(k-1) – az(2)*y(k-1) – az(3)*y(k-2)

 Matlab中得到的結果如下（訊號頻率0.1Hz，噪聲頻率6Hz）：

三、C語言函式編寫與驗證

C函式原始碼如下：

const float bz[2] = {0,0.128580115806658};  //分子
const float az[3] = {1,-1.42252474659021,0.553007125840971};  //分母

float Data_Output[DATA_LENTH];  //輸出資料

float* but_filter(unsigned int len, float* x)  //len為輸入資料陣列長度，x為輸入資料陣列指標
{

unsigned int i = 2;
static float init[2] = {0,0};  //初值，一開始設為0

  if(len<2)                   //如果長度小於2，直接返回
return Data_Output;

Data_Output[0] = init[0];   //賦初值
Data_Output[1] = init[1];
for(i = 2;i < len;i++)
{
Data_Output[i] = bz[0]*x[i] + bz[1]*x[i-1] - az[1]*Data_Output[i-1] - az[2]*Data_Output[i-2];
/*演算法為a1*y(k) = b1*x(k) + b2*x(k-1) - a(2)*y(k-1) - a(3)*y(k-2)*/
/*由於a1 = 1，故不做除法*/
}
init[0] = Data_Output[len-2];  //考慮到會被連續呼叫，此次的終值作為下次的初值
init[1] = Data_Output[len-1];

return Data_Output; 

}