IIR滤波器分阶计算

此链接的不太好不过可以看看，下边的挺好
生成的滤波器系数文件如下：

/*

• Filter Coefficients (C Source) generated by the Filter Design and Analysis Tool

• Generated by MATLAB® 7.6 and the Signal Processing Toolbox 6.9.

• Generated on: 03-Dec-2010 10:41:03

*/

/*

• Discrete-Time IIR Filter (real)

• ---

• Filter Structure : Direct-Form II, Second-Order Sections

• Number of Sections : 3

• Stable : Yes

• Linear Phase : No

*/

/* General type conversion for MATLAB generated C-code */

#include “tmwtypes.h”

/*

• Expected path to tmwtypes.h

• D:/MATLAB/R2008a/extern/include/tmwtypes.h

*/

#define MWSPT_NSEC 7

const int NL[MWSPT_NSEC] = { 1,3,1,3,1,3,1 };

const real64_T NUM[MWSPT_NSEC][3] = {

{

0.09065504059673,                 0,                 0

},

{

               1,  -0.1323122149853,   0.9999674089086

},

{

               1,                 0,                 0

},

{

               1,   0.1670351201308,   0.9999889247428

},

{

               1,                 0,                 0

},

{

               1,    1.417032671609,   0.9978019623105

},

{

               1,                 0,                 0

}

};

const int DL[MWSPT_NSEC] = { 1,3,1,3,1,3,1 };

const real64_T DEN[MWSPT_NSEC][3] = {

{

               1,                 0,                 0

},

{

               1,    -0.40321085647,   0.7334254033056

},

{

               1,                 0,                 0

},

{

               1,  -0.6868636040216,   0.2670185171768

},

{

               1,                 0,                 0

},

{

               1,  -0.2880720042256,   0.9480010462991

},

{

               1,                 0,                 0

}

};

上述系数是3个2阶节IIR结构的级联。可以转换为我们熟悉的b/a的形式如下：

double a[3][3] = {

 {1,  -0.6868636040216,   0.2670185171768},

 {1,    -0.40321085647,   0.7334254033056},

 {1,  -0.2880720042256,   0.9480010462991}

};

double b[3][3] = {

 {1,   1.417032671609,   0.9978019623105},

 {1,   0.1670351201308,   0.9999889247428},

 {1,  -0.1323122149853,   0.9999674089086},

};

注意上面系数文件中还有一个增益：

double gain = 0.09065504059673；

这个增益最好在第一级实现以后加入运算。这样可以减小误差，保证数据 动态范围不被溢出。尤其是在定点计算的时候尤为如此。
2阶节IIR滤波的直接实现
一个2阶节结构是下面这样一个表达式：

实现上面这个表达式需要4个过去的历史值，把它定义在结构体
typedef struct tag_ IIR_State_2order

{

 float y2;

 float y1;

 float x1;

 float x0;

} IIR_State_2order;

调用下面函数之前需要把上述结构体所有值初始化为零。滤波按一帧一帧数据进行。

#define ONE_FRAME_SAMPLE_SIZE 1024

void cy_signal_filter_by_iir(signed short* pcmIn, IIR_State_2order* filter_state, float a[], float b[], signed short* pcmOut)

{

int i;

float x2;

float tmp;

for ( i = 0; i < ONE_FRAME_SAMPLE_SIZE; i++ )

{

  x2 = filter_state->x1;

  filter_state->x1 = filter_state->x0;

  filter_state->x0 = pcmIn[i];

  tmp = ( float )(b[0]* filter_state->x0 + b[1]* filter_state->x1 +

b[2] * x2 - a[1] * filter_state->y1 - a[2] * filter_state->y2);

   if(tmp >= 32767)

   {

       tmp = 32767;

   }

   if(tmp <= -32768)

   {

      tmp = -32768;

   }

  pcmOut[i] = (signed short)tmp;

  filter_state->y2 = filter_state->y1;

  filter_state->y1 = tmp;

}

}

有一个简单的技巧可以把上面的计算简化，使得历史状态数由4减少为2。定义下面的表达式：

结构体定义如下：

typedef struct tag_ IIR_State_2order

{

 float st1;

 float st2;

} IIR_State_2order;

void cy_signal_filter_by_iir(signed short* pcmIn, IIR_State_2order* filter_state, float a[], float b[], signed short* pcmOut)

{

int i;

float st;

float Tmp_fl;

for ( i = 0; i < ONE_FRAME_SAMPLE_SIZE; i++ )

{

    st = (float)(pcmIn[i] -  a[1] * filter_state->st1  - a[2] * filter_state-> st2);

    Tmp_fl = b[0] *  st + b[1] * filter_state->st1 + b[2] * filter_state->st2;

    filter_state->st2= filter_state->st1;

    filter_state->st1 = st;

    if(Tmp_fl >= 32767.0)

    {

        Tmp_fl = 32767;

    }

    if(Tmp_fl <= -32768)

    {

       Tmp_fl = -32768;

    }



     pcmOut[i] = (signed short)Tmp_fl;

}

}

6阶节IIR滤波的实现
有个上面的基础，我们来实现上面设计的6阶IIR滤波器。

6阶节分解为3个2阶节级联实现。每个2阶节需要2个历史状态，总共需要6个历史状态。结构体定义如下：

typedef struct tag_IIR_State_3Order

{

 double w01;

 double w02;

 double w11;

 double w12;

 double w21;

 double w22;

}IIR_State_6order;

代码中数组a,b,还有gain的定义见第一部分。

void cy_signal_filter_by_6th_iir(signed short* pcmIn, IIR_State_6order* filter_state, int sample_size)

{

 double x1, x2, x3, Tmp_f00, Tmp_f10, Tmp_f20;

 int i;

 double Tmp_pcm;

 for (i = 0; i < sample_size; i++)

 {

     Tmp_pcm = pcmIn[i];

     Tmp_f00 = Tmp_pcm - a[0][1] * filter_state->w01 - a[0][2] * filter_state->w02;

     x1 = Tmp_f00 + b[0][1] * filter_state->w01 + b[0][2] * filter_state->w02;

     filter_state->w02 = filter_state->w01;

     filter_state->w01 = Tmp_f00;



     x1 = gain * x1;



     Tmp_f10 = x1 - a[1][1] * filter_state->w11 - a[1][2] * filter_state->w12;

     x2 = Tmp_f10 + b[1][1] * filter_state->w11 + b[1][2] * filter_state->w12;

     filter_state->w12 = filter_state->w11;

     filter_state->w11 = Tmp_f10;





     Tmp_f20 = x2 - a[2][1] * filter_state->w21 - a[2][2] * filter_state->w22;

     x3 = Tmp_f20 * b[2][0] + b[2][1] * filter_state->w21 + b[2][2] * filter_state->w22;



     filter_state->w22 = filter_state->w21;

     filter_state->w21 = Tmp_f20;



     if (x3 >= 32767)

     {

          x3 = 32767;

     }

     if (x3 <= -32768)

     {

          x3 = -32768;

     }

     pcmIn[i] = (signed short)x3;



 }

}
最后看下滤波的效果：

滤波之后的频谱：