[導讀]：在嵌入式系統中經常需要采集模擬信號，采集模擬信號的信號鏈中難免引入干擾，那么如何濾除干擾呢？今天就來個一步一步描述如何設計部署一個IIR濾波器到你的系統，
何為IIR濾波器？
無限沖激回應（IIR： Infinite Impulse Response）是一種適用于許多線性時不變系統的屬性，這些系統的特征是具有一個沖激回應h（t），該沖激回應h（t）不會在特定點上完全變為零，而是無限期地持續， 這與有限沖激回應（FIR： Finite Impulse Response）系統形成對比，在有限沖激回應（FIR）系統中，對于某個有限T，在時間t> T時，沖激回應確實恰好變為零， 線性時不變系統的常見示例是大多數電子和數字濾波器， 具有此屬性的系統稱為IIR系統或IIR濾波器，
這是常見的教科書式數學嚴謹定義，很多人看到這一下就蒙了，能說人話嗎？
線性時不變系統理論俗稱LTI系統理論，源自應用數學，直接在核磁共振頻譜學、地震學、電路、信號處理和控制理論等技術領域運用，它研究的是線性、非時變系統對任意輸入信號的回應，雖然這些系統的軌跡通常會隨時間變化（例如聲學波形）來測量和跟蹤，但是應用到影像處理和場論時，LTI系統在空間維度上也有軌跡，因此，這些系統也被稱為線性非時變平移，在最一般的范圍理論給出此理論，在離散（即采樣）系統中對應的術語是線性非時變平移系統，由電阻、電容、電感組成的電路是LTI系統的一個很好的例子，比如一個運放系統在一定頻帶范圍內滿足信號的時域疊加，輸入一個100Hz和200Hz正弦信號，輸出頻率是這兩種信號的線性疊加，
用數學對LTI系統描述：
線性：輸入\(x_1(t)\),產生回應\(y_1(t)\),而輸入\(x_2(t)\),產生回應\(y_2(t)\),那么縮放和加和輸入\(a_1x_1(t)+a_2x_2(t)\),產生縮放、加和的回應\(a_1y_1(t)+a_2y_2(t)\),其中\(a_1\)和\(a_2\)是標量，對于任意的有：輸入\(\sum_{k=0}^{N}a_kx_k(t)\),產生回應\(\sum_{k=0}^{N}a_ky_k(t)\)
時不變性：指如果將系統的輸入信號延遲\(\tau\)秒，那么得到的輸出回應也相應延時\(\tau\)秒，用數學描述，也即如果輸入\(x_1(t)\),產生回應\(y_1(t)\),而輸入\(x_1(t+\tau)\)，產生回應\(y_1(t+\tau)\)，這么描述還是不易懂，來個圖，有圖有真相：
假定一個信號放大電路對100Hz正弦信號放大2倍：
則輸出為：
而對200Hz的正弦信號，其放大倍數為1.7倍，（做過運放電路設計的朋友應該有經驗，在其同頻帶其放大倍數往往并不平坦，也即幅頻回應在頻帶內不平坦，這是比較常見的），
也即輸入為：
回應為:
那么如果輸入100Hz和200Hz的時域疊加信號，則其輸入為：
則其回應為：
由這些圖可看出，輸入信號的形狀保持不變，輸出為對應輸入的線性時域疊加，
對于時不變，就不用圖描述了，在一個真實電路中，如果輸入延遲一定時間，則回應對應延遲相同時間輸出，
上面這么多文字只是為了描述在什么場合可以使用IIR濾波器對信號進行數字濾波，總結而言，就是在線性時不變系統中適用，換言之，在大多數電路系統中我們都可以嘗試采用IIR濾波器進行數字濾波，
那么究竟什么是IIR濾波器呢？從數字信號處理的書籍中我們能看到這樣的Z變換信號流圖：
x(n-1)表示延遲一拍,在數字系統中表示對于輸入信號而言,即為上一次采樣值,對于輸出而言,即為上一次的輸出值，
在時域中對于上述流圖，用時域描述即為：
如果熟悉Z變換，則Z變換傳遞函式為：
上述數字濾波器，如果從編程的角度來看，x(n-1),表示上一次的信號，可能是來自ADC的上次采樣，而y(n-1)則為上一次濾波器的輸出值，對應就比較好理解x(n-N)就表示前第n次輸入樣本信號，而y(n-M)則為前第M次濾波器的輸出，
說了這么多，只是為了更好的理解概念，只有概念理解正確，才能用爭取， 概念理解這對工程師而言，非常之重要，
如何設計呢？
打開fdatool
MATLAB提供了非常容易使用的FDATool幫助我們設計數字濾波器，真正精彩的地方開始了，讓我們拭目以待究竟如何一步一步設計并實施一個IIR濾波器，
首先打開MATLAB,在命令列中敲fdatool
彈出表單就是fdatool了，如下：
在設計具體,有幾個相關概念需要澄清：
Fs：采樣率，單位為Hz,真實部署在系統中，請務必確保樣本是按照恒定采樣率進行采樣，否則將得不到想要的效果，
Fpass: 通頻帶，單位為Hz,即系統中期望通過的最高頻率，
Fstop: 截至頻率，即幅頻回應的-3dB處的頻率，這個如不理解，請自行查閱相關書籍，
分貝dB: 這是一個無單位反應輸出與輸入倍數的一個術語，電學中分貝與放大倍數的轉換關系為：
- A(V)(dB)=20lg(Vo/Vi)；電壓增益,Vo 為輸出電壓，Vi為輸入電壓
- A(I)(dB)=20lg(Io/Ii)；電流增益,Io 為輸出電流，Ii為輸入電流
- A(p)(dB)=10lg(Po/Pi)；功率增益，Po 為輸出功率，Pi為輸入功率
濾波器型別：這里有Butterworth(巴特沃斯)、Chebyshev Type I,Chebyshev Type II、(切比雪夫)、Elipic 等可選，
- 巴特沃斯 Butterworth，也被稱作最大平坦濾波器，巴特沃斯濾波器的特點是通頻帶內的頻率回應曲線最大限度平坦，沒有紋波，
- 切比雪夫 Chebyshev，是在通帶或阻帶上頻率回應幅度等波紋波動的濾波器，切比雪夫濾波器在過渡帶比巴特沃斯濾波器的衰減快，但頻率回應的幅頻特性不如后者平坦，
- 橢圓 Elliptic，橢圓濾波器是在通帶和阻帶等波紋的一種濾波器，
- ...這里就不一一介紹了，有興趣可以去查信號處理書籍，
就其特點，這里對其中幾種略作介紹：
- 巴特沃斯具有最平坦的通帶，
- 橢圓濾波器衰減最快，但是通帶、阻帶都有波紋，
- 切比雪夫濾波器衰減比巴特沃斯快，但比橢圓濾波器慢，波紋區域可選擇，
假設我們需要設計一個IIR濾波器，采樣率為32000Hz, 有用信號頻率在2000Hz內，設計IIR濾波器對信號進行數字濾波，這里為節省算力，我們指定濾波器的階數，也即傳遞函式中N/M中的最大值，一般而言N大于M，
這里指定階數為8階，型別指定為巴特沃斯型IIR濾波器，輸入階數8階，采樣率32000Hz,截至頻率設定為2000Hz然后點擊Design Filter如下圖所示：
其相頻回應曲線如下：
除此之外，我們還可以將幅頻與相頻曲線放在一個頻率坐標上去看設計結果：
匯出濾波器引數，這里我們選擇，
選擇ASCII碼
然后就得到了一個檔案，保存2kHz_LPF.fcf，檔案名隨你喜歡，
檔案內容如下：
% Generated by MATLAB(R) 8.4 and the Signal Processing Toolbox 6.22.
% Generated on: 27-Mar-2020 21:27:06
% Coefficient Format: Decimal
% Discrete-Time IIR Filter (real)
% -------------------------------
% Filter Structure : Direct-Form II, Second-Order Sections
% Number of Sections : 4
% Stable : Yes
% Linear Phase : No
SOS Matrix:
1 2 1 1 -1.7193929141691948 0.8610574795347461
1 2 1 1 -1.5237898734101736 0.64933827386370635
1 2 1 1 -1.4017399331200424 0.51723237044751591
1 2 1 1 -1.3435020629061745 0.45419615396638446
Scale Values:
0.035416141341387819
0.031387100113383172
0.028873109331868367
0.027673522765052503
至此設計作業就結束了，馬上進入濾波器的部署測驗階段，
部署測驗濾波器
到這里，沒有經驗的朋友可能會說，這么一堆引數我該咋用呢？
需要自己去寫前面描述的計算公式嗎？當然你也可以這么做，這里就不寫了,ARM的CMSIS庫已經幫大家設計好了種類繁多的數字信號處理函式實作了，而且經過了測驗，這里直接拿來用即可，有興趣自己寫也不難，只要理解Z傳遞函式概念內涵，非常容易實作，這里我們采用32位浮點實作函式：arm_biquad_cascade_df1_f32，
該函式位于：CMSIS\DSP\Source\FilteringFunctions\arm_biquad_cascade_df1_init_f32.c中
以及CMSIS\DSP\Source\FilteringFunctions\arm_biquad_cascade_df1_f32.c
我們來看一看這個函式：
arm_biquad_cascade_df1_init_f32.c：
/*
*作用 :初始化濾波器
*S :指向浮點SOS級聯結構的實體，
*numStages:濾波器中二階SOS的數量
*pCoeffs :濾波器引數指標,引數按下列順序存盤
* {b10, b11, b12, a11, a12, b20, b21, b22, a21, a22, ...}
*pState :歷史狀態緩沖區指標
*/
void arm_biquad_cascade_df1_init_f32(
arm_biquad_casd_df1_inst_f32 * S,
uint8_t numStages,
const float32_t * pCoeffs,
float32_t * pState)
{
/* Assign filter stages */
S->numStages = numStages;
/* Assign coefficient pointer */
S->pCoeffs = pCoeffs;
/* Clear state buffer and size is always 4 * numStages */
memset(pState, 0, (4U * (uint32_t) numStages) * sizeof(float32_t));
/* Assign state pointer */
S->pState = pState;
}
arm_math.h 定義了須用到的結構體，對于本例相關的結構體為arm_biquad_casd_df1_inst_f32
typedef struct
{
unsigned int numStages; /*2階節的個數，應為2*numStages. */
float *pState; /*狀態系數陣列指標，陣列長度為4*numStages*/
float *pCoeffs; /*系數陣列指標， 陣列的長度為5*numStages.*/
} arm_biquad_casd_df1_inst_f32;
濾波器具體濾波函式為arm_biquad_cascade_df1_f32
/**
* *S ：指向浮點Biquad級聯結構的實體.
* *pSrc ：指向輸入資料塊，
* *pDst ：指向輸出資料塊，
* blockSize：每次呼叫要處理的樣本數，
* 回傳值 ：無.
*/
void arm_biquad_cascade_df1_f32(
const arm_biquad_casd_df1_inst_f32 * S,
float * pSrc,
float * pDst,
unsigned int blockSize)
{
float *pIn = pSrc; /*源指標 */
float *pOut = pDst; /*目的指標 */
float *pState = S->pState; /*狀態指標 */
float *pCoeffs = S->pCoeffs; /*引數指標 */
float acc; /*累加器 */
float b0, b1, b2, a1, a2; /*濾波器引數 */
float Xn1, Xn2, Yn1, Yn2; /*濾波器狀態變數*/
float Xn; /*臨時輸入 */
unsigned int sample, stage = S->numStages; /*回圈計數 */
do
{
/* Reading the coefficients */
b0 = *pCoeffs++;
b1 = *pCoeffs++;
b2 = *pCoeffs++;
a1 = *pCoeffs++;
a2 = *pCoeffs++;
Xn1 = pState[0];
Xn2 = pState[1];
Yn1 = pState[2];
Yn2 = pState[3];
sample = blockSize >> 2u;
while(sample > 0u)
{
/* 讀第一個輸入 */
Xn = *pIn++;
/* acc = b0 * x[n] + b1 * x[n-1] + b2 * x[n-2] + a1 * y[n-1] + a2 * y[n-2] */
Yn2 = (b0 * Xn) + (b1 * Xn1) + (b2 * Xn2) + (a1 * Yn1) + (a2 * Yn2);
/* Store the result in the accumulator in the destination buffer. */
*pOut++ = Yn2;
/* 每次計算輸出后，狀態都應更新. */
/* 狀態應更新為: */
/* Xn2 = Xn1 */
/* Xn1 = Xn */
/* Yn2 = Yn1 */
/* Yn1 = acc */
/* Read the second input */
Xn2 = *pIn++;
/* acc = b0 * x[n] + b1 * x[n-1] + b2 * x[n-2] + a1 * y[n-1] + a2 * y[n-2] */
Yn1 = (b0 * Xn2) + (b1 * Xn) + (b2 * Xn1) + (a1 * Yn2) + (a2 * Yn1);
/* 將結果存盤在目標緩沖區的累加器中. */
*pOut++ = Yn1;
/* 每次計算輸出后，狀態都應更新. */
/* 狀態應更新為: */
/* Xn2 = Xn1 */
/* Xn1 = Xn */
/* Yn2 = Yn1 */
/* Yn1 = acc */
/*讀第三個輸入 */
Xn1 = *pIn++;
/* acc = b0 * x[n] + b1 * x[n-1] + b2 * x[n-2] + a1 * y[n-1] + a2 * y[n-2] */
Yn2 = (b0 * Xn1) + (b1 * Xn2) + (b2 * Xn) + (a1 * Yn1) + (a2 * Yn2);
/* 將結果存盤在目標緩沖區的累加器中. */
*pOut++ = Yn2;
/* 每次計算輸出后，狀態都應更新. */
/* 狀態應更新為: */
/* Xn2 = Xn1 */
/* Xn1 = Xn */
/* Yn2 = Yn1 */
/* Yn1 = acc */
/* 讀第四個輸入 */
Xn = *pIn++;
/* acc = b0 * x[n] + b1 * x[n-1] + b2 * x[n-2] + a1 * y[n-1] + a2 * y[n-2] */
Yn1 = (b0 * Xn) + (b1 * Xn1) + (b2 * Xn2) + (a1 * Yn2) + (a2 * Yn1);
/* 將結果存盤在目標緩沖區的累加器中. */
*pOut++ = Yn1;
/* 每次計算輸出后，狀態都應更新. */
/* 狀態應更新為: */
/* Xn2 = Xn1 */
/* Xn1 = Xn */
/* Yn2 = Yn1 */
/* Yn1 = acc */
Xn2 = Xn1;
Xn1 = Xn;
/* 遞級訓圈計數器 */
sample--;
}
/* 如果blockSize不是4的倍數，
*請在此處計算任何剩余的輸出樣本，
*不使用回圈展開. */
sample = blockSize & 0x3u;
while(sample > 0u)
{
/* 讀取輸入 */
Xn = *pIn++;
/* acc = b0 * x[n] + b1 * x[n-1] + b2 * x[n-2] + a1 * y[n-1] + a2 * y[n-2] */
acc = (b0 * Xn) + (b1 * Xn1) + (b2 * Xn2) + (a1 * Yn1) + (a2 * Yn2);
/* 將結果存盤在目標緩沖區的累加器中. */
*pOut++ = acc;
/* 每次計算輸出后，狀態都應更新， */
/* 狀態應更新為: */
/* Xn2 = Xn1 */
/* Xn1 = Xn */
/* Yn2 = Yn1 */
/* Yn1 = acc */
Xn2 = Xn1;
Xn1 = Xn;
Yn2 = Yn1;
Yn1 = acc;
/* d遞級訓圈計數器 */
sample--;
}
/* 將更新后的狀態變數存盤回pState陣列中 */
*pState++ = Xn1;
*pState++ = Xn2;
*pState++ = Yn1;
*pState++ = Yn2;
/*第一階段從輸入緩沖區到輸出緩沖區. */
/*隨后的numStages在輸出緩沖區中就地發生*/
pIn = pDst;
/* 重置輸出指標 */
pOut = pDst;
/* 遞級訓圈計數器 */
stage--;
} while(stage > 0u);
}
開始測驗：
#include <stdio.h>
#include <math.h>
/*
SOS Matrix:
1 2 1 1 -1.7193929141691948 0.8610574795347461
1 2 1 1 -1.5237898734101736 0.64933827386370635
1 2 1 1 -1.4017399331200424 0.51723237044751591
1 2 1 1 -1.3435020629061745 0.45419615396638446
Scale Values:
0.035416141341387819
0.031387100113383172
0.028873109331868367
0.027673522765052503
做如下轉換：
1.縮放
[1 2 1] * 0.035416141341387819
[1 2 1] * 0.031387100113383172
[1 2 1] * 0.028873109331868367
[1 2 1] * 0.027673522765052503
得到：
[0.035416141341387819 2*0.035416141341387819 0.035416141341387819]
[0.031387100113383172 2*0.031387100113383172 0.031387100113383172]
[0.028873109331868367 2*0.028873109331868367 0.028873109331868367]
[0.027673522765052503 2*0.027673522765052503 0.027673522765052503]
2.舍掉第四列引數
3.將后兩列分別乘以-1，即：
0.035416141341387819 2*0.035416141341387819 0.035416141341387819 -1.7193929141691948 0.8610574795347461
0.031387100113383172 2*0.031387100113383172 0.031387100113383172 -1.5237898734101736 0.64933827386370635
0.028873109331868367 2*0.028873109331868367 0.028873109331868367 -1.4017399331200424 0.51723237044751591
0.027673522765052503 2*0.027673522765052503 0.027673522765052503 -1.3435020629061745 0.45419615396638446
這樣就得到了濾波器系陣列了
*/
#define IIR_SECTION 4 /*見前面設計輸出為4個SOS塊*/
static float iir_state[4*IIR_SECTION];/*歷史狀態緩沖區 */
const float iir_coeffs[5*IIR_SECTION]={
0.035416141341387819,2*0.035416141341387819,0.035416141341387819,1.7193929141691948,-0.8610574795347461, 0.031387100113383172,2*0.031387100113383172,0.031387100113383172,1.5237898734101736,-0.64933827386370635, 0.028873109331868367,2*0.028873109331868367,0.028873109331868367,1.4017399331200424,-0.51723237044751591, 0.027673522765052503,2*0.027673522765052503,0.027673522765052503,1.3435020629061745,-0.45419615396638446
};
static arm_biquad_casd_df1_inst_f32 S;
/*假定采樣512個點*/
#define BUF_SIZE 512
#define PI 3.1415926
#define SAMPLE_RATE 32000 /*32000Hz*/
int main()
{
float raw[BUF_SIZE];
float raw_4k[BUF_SIZE];
float raw_out[BUF_SIZE];
float raw_noise[BUF_SIZE];
float raw_noise_out[BUF_SIZE];
arm_biquad_casd_df1_inst_f32 S;
FILE *pFile=fopen("./simulation.csv","wt+");
if(pFile==NULL)
{
printf("file opened failed");
return -1;
}
for(int i=0;i<BUF_SIZE;i++)
{
/*模擬9000Hz輸入信號 */
raw[i] = 0.5*1024.0/3*sin(2*PI*800*i/32000.0f)+rand()%50;
raw_4k[i] = 0.5*1024.0/3*sin(2*PI*4000*i/32000.0f);
/*模擬9000Hz +11000Hz疊加輸入*/
raw_noise[i] = raw[i] + raw_4k[i];
}
arm_biquad_cascade_df1_init_f32(&S, IIR_SECTION, (float *)&iir_coeffs[0], (float *)&iir_state[0]);
arm_biquad_cascade_df1_f32(&S, raw, raw_out, BUF_SIZE);
for(int i=0;i<BUF_SIZE;i++)
{
fprintf(pFile,"%f,",raw[i]);
}
fprintf(pFile,"\n");
for(int i=0;i<BUF_SIZE;i++)
{
fprintf(pFile,"%f,",raw_4k[i]);
}
fprintf(pFile,"\n");
for(int i=0;i<BUF_SIZE;i++)
{
fprintf(pFile,"%f,",raw_out[i]);
}
/*重新初始化*/
arm_biquad_cascade_df1_init_f32(&S, IIR_SECTION, (float *)&iir_coeffs[0], (float *)&iir_state[0]);
arm_biquad_cascade_df1_f32(&S, raw_noise, raw_noise_out, BUF_SIZE);
fprintf(pFile,"\n");
for(int i=0;i<BUF_SIZE;i++)
{
fprintf(pFile,"%f,",raw_noise[i]);
}
fprintf(pFile,"\n");
for(int i=0;i<BUF_SIZE;i++)
{
fprintf(pFile,"%f,",raw_noise_out[i]);
}
fclose(pFile);
return 0;
}
利用csv檔案，將模擬資料存盤，直接用excel打開，將行資料生成曲線圖如下：
- 第一幅圖，為800Hz信號混入隨機噪聲的波形
- 第二幅圖，為4000Hz信號，對假定系統為無用干擾信號
- 第三幅圖， 為800Hz 混入隨機噪聲過濾后，已經很好的還原有用信號頻率
- 第四幅圖， 為800Hz信號混入隨機噪聲，同時疊加4000Hz干擾的波形，對系統而言，從時域中，明顯可見，有用信號已經完全扭曲
- 第五幅圖，為800Hz信號混入隨機噪聲，同時疊加4000Hz干擾的輸入，經過該低通濾波器后的波形，與第三幅圖基本一樣，已經非常好的濾除了干擾信號，
總結：
-
IIR濾波器在線性時不變系統中可以很好的解決工程中一般噪聲問題
-
如果需要設計帶通、高通濾波器其步驟基本類似，只是濾波器的引數以及SOS塊個數可能不一樣而已
-
需要提醒的時，IIR的相頻回應不線性，如果系統對相頻回應有嚴格要求，就需要采用其他的數字濾波器拓撲形式了
-
實際應用中，如果階數不高時，現在算力強勁的單片機或者DSP以及可以直接使用浮點處理，
-
如果對處理速度有嚴格的實時要求，需要在極短時間進行濾波處理，可以考慮降低階數，或采用定點IIR濾波演算法實作，也或者將文中函式進行匯編級優化，
文章出自微信公眾號：嵌入式客堆疊，著作權所有，嚴禁商用，關注公眾號，獲取更新更多內容
轉載請註明出處，本文鏈接：https://www.uj5u.com/houduan/27304.html
標籤：C
上一篇：手把手教系列之移動平均濾波器實作
下一篇：推箱子游戲-C語言