IIR濾波器軟體實現（Matlab+C++）

標籤：數組   ...   技術   變數   fopen   c++編程   #define   col   理論   

使用C++來寫一個IIR濾波器

我們首先要在MATLAB中設計一個IIR濾波器，並產生一個標頭檔，這個標頭檔中反映了IIR濾波器的頻率響應特性

理論支援

IIR濾波叫做遞迴濾波器，它是一種具有反饋的濾波器。當階數較大時一般採取多個二階節濾波進行串聯，這樣可以提高系統穩定性。

一個二階節係數規律：

可以寫出第K個二階節的差分方程

N個二階節的級連接構如所示：

根據二階節圖，把前一級的輸出作為後一級的輸入，就可以通過軟體實現IIR數字濾波的功能。

 

使用Matlab產生標頭檔

首先開啟MATLAB中Filter Design & Analysis Tool

 

 

 

 

 

 

 

 

 

 

 

 

 

 

這裡我們先設計一個低通濾波器

Fs代表採樣頻率，採樣頻率必須大於原訊號最高頻率的兩倍，

否則會產生頻譜混疊。

Fpass為通帶頻率，Fstop為阻帶截止頻率

這些參數設定好就可以點擊Design Filter

 

產生的是一個二階節濾波組合，一共有31階，也就是多個二階濾波器的組合

接著在Target選項中產生C Header File

 Numerator為分子係數數組的命名，Numerator length為分子係數數組的長度，

 Denominator為分母。

對產生標頭檔進行分析

以下以Fpass為10K，Fstop為12K的低通濾波器舉栵

在使用標頭檔前需要根據情況將Matlab的資料類型轉換為C++支援的資料類型，這裡我們使用double類型

 

在分析標頭檔前先看下Matlab提供的第一節濾波參數

以第一個二階節的資料舉例：

• Numerator: 1  2  1

• Denominator: 1  -0.55930961405944157  0.92579835996619642

• Gain：0.34162218647668868

Numerator為分子的係數，分別為b0，b1，b2

Denominator為分母的係數，分別為a0，a1，a2.

Gain為各節的增益，此項為了穩定各節，穩定訊號大小

 

接著對照標頭檔，以下為標頭檔主要部分的一段截取：

#define MWSPT_NSEC 41int NL[MWSPT_NSEC] = { 1,3,1,3,1,3,1,3,1,3,1,3,1,3,1,3,1,3,1,3,1,3,1,3,1,3,1,3,1,3,1,3,1,3,1,3,1,3,1,2,1 };double NUM[MWSPT_NSEC][3] = {  {     0.3416221864767,                 0,                 0   },  {                   1,                 2,                 1   },  {     0.3180955154747,                 0,                 0   },  {                   1,                 2,                 1   },......

MWSPT_NSEC為濾波器階數，具體的節數在標頭檔開頭的注釋中

NL[MWSPT_NSEC]這個數組定義了NUM[MWSPT_NSEC][3]數組每一行的有用資料個數（可以不用）

在NUM[MWSPT_NSEC][3]數組（分子參數）奇數行第一項都為增益項，偶數行為3個係數，分別為b0,b1,b2。

由此可以找出規律，定義K為目前所在的階數，p為數組的首指標，則，每一節的增益項為（p+6*K）,第一個係數為（p+3+6*K），

第二個係數為（p+3+6*K+1），第三個係數為（p+3+6*K+2）。

C++編程實現

在軟體設計的過程中，每個二階節的延遲變數只取 和 ， 作為中間變數在過程中直接賦給 。這是因為對於下一個輸入資料n+1的延遲變數即為上一個輸入資料的 和 ，採用這種方式進行設計，可以節省寄存器的空間。

為了提高處理速度，程式中需要使用指標進行參數傳遞，特別注意二維數組的首地址傳遞方式為&a[0][0]->double* a。

 

濾波函數

double iir(double *a, double *b,double* w, double xin, int N_IIR)

{

    int k;

double temp = xin;

 

for (k = 0; k<N_IIR; k++)

{

   

    *(w+k*3) = temp - *(a + 3+6 * k  + 1) *(*(w + k * 3+1)) - *(a + 3 + 6 * k + 2) *(*(w + k * 3+2));

         //這裡temp為本二階節的輸入，也是上一個二階節的輸出

         temp = *(b + 3 + 6 * k )* (*(w + k * 3)) + *(b + 3 + 6 * k + 1) * (*(w + k * 3+1)) + *(b + 3+6 * k + 2)* (*(w + k * 3+2));

//這裡temp為本二階節的輸出，也是下一個二階節的輸入

 

        

         *(w + k * 3 + 2) = *(w + k * 3 + 1);

         *(w + k * 3 + 1) = *(w + k * 3);

        

         temp = temp*(*(b + 6 * k));//放大倍數，穩定訊號

}

return temp;

}

實際測試

測試Fpass為10K，Fstop為12K的低通濾波器

在程式中輸入三個頻率為2K，11K，20K的訊號，理論上2k完全通過，11k部分衰減，20K完全濾除。

為原訊號，為濾波後訊號。

實際測試發現符合設計要求，而且在過渡帶訊號也基本完全衰減。

測試用C++程式void main(){    const int N = 100;    int i,j;    double xn[N];    double w[20][3];    double yn[N];    for (i = 0; i < 20; i++)//初始化    {        for (j = 0; j < 3; j++)            w[i][j] = 0;    }    for (i = 0; i < N; i++)    {        xn[i] = sin(2 * 3.1416 * 20 / 50 * i)+ sin(2 * 3.1416 * 2 / 50 * i)+ sin(2 * 3.1416 * 11 / 50 * i);        yn[i]=iir(&DEN[0][0], &NUM[0][0], &w[0][0],xn[i], 20);        }    ofstream SaveFile_a("xn.txt");    for (i = 0; i<N; i++)        SaveFile_a << " " << xn[i] << endl;    SaveFile_a.close();    ofstream SaveFile_b("yn.txt");    for (i = 0; i<N; i++)        SaveFile_b << " " << yn[i] << endl;    SaveFile_a.close();}分析用Matlab程式xn1=fopen(‘xn.txt‘,‘r‘);[xn,count]=fscanf(xn1,‘%f‘);fclose(xn1);N = length(xn);%求取抽樣點數xn_f = fft(xn);%對訊號進行傅裡葉變換xn_f=abs(xn_f(1:N/2));f = 50000/N*(0:N/2-1);subplot(211);stem(f,abs(xn_f));xlabel(‘Frequency / (s)‘);ylabel(‘Amplitude‘);title(‘原訊號頻譜‘);grid;yn1=fopen(‘yn.txt‘,‘r‘);[yn,count]=fscanf(yn1,‘%f‘);fclose(yn1);yn_f = fft(yn);%對訊號進行傅裡葉變換yn_f=abs(yn_f(1:N/2));subplot(212);stem(f,abs(yn_f));xlabel(‘Frequency / (s)‘);ylabel(‘Amplitude‘);title(‘濾波後訊號頻譜‘);grid;

 

 

 

IIR濾波器軟體實現（Matlab+C++）