文档章节

IIR数字滤波器实现

abcijkxyz
 abcijkxyz
发布于 2016/11/22 16:46
字数 1232
阅读 28
收藏 0

题目:16k采样率音频数据下采样到8k采样率

求解方案分析:直接每隔一个取一个采样值,这样就可以得到8k采样率的数据。但是这样明显会有问题。按照采样率变换理论,首先应该通过一个低通滤波器,滤掉[pi/2, pi]这个区间上的频率,以防止下采样造成的频率混叠。这个低通滤波器在很多书上都用FIR滤波去实现,并且可以用FIR滤波的多相结构去实现。这样滤波和下采样过程可以互换位置。即先下采样再进行多相FIR滤波。

在嵌入式设备上FIR滤波会占用较长的时间,为此,我们可以采用IIR滤波器来做。滤波器的设计采用MATLABFDATool工具。相关的参数如下:

 

 

采用最小P阶范数设计IIR滤波器的幅频特性较好。滤波器阶数设定为6阶,需要滤除4k8k的频率段,才能保证无混叠失真。实际由于滤波器的特性没法做到理想的状态,选择滤波器截止频率为3800hz36003800为过渡带宽。其它选项采用默认设置。设计的滤波器幅度响应如下图:

 

生成的滤波器系数文件如下:

/*

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

 *

 * Generated by MATLAB(R) 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

  }

};

 

上述系数是32阶节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滤波的实现

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

6阶节分解为32阶节级联实现。每个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;

 

     }

}

 最后看下滤波的效果:

 

 

滤波之后的频谱:

滤波效果不错,下面可以进行下采样了。

 

本文转载自:http://www.cnblogs.com/celerychen/archive/2010/12/03/3588220.html

共有 人打赏支持
abcijkxyz
粉丝 61
博文 6196
码字总数 1876
作品 0
深圳
项目经理
【工具使用系列】关于 MATLAB 数字滤波器设计,你需要知道的事

如何设计滤波器 MATLAB 实现 FIR滤波器设计 窗设计技术 频率采样设计技术 最优等波动设计技术 IIR滤波器设计 模拟滤波器原型的特征 低通滤波器设计 频带变化 应用实例 系数调整 LMS 算法 系统...

AllenMoore
02/03
0
0
用Matlab的FDAtool生成IIR滤波器参数

首先我们要明白相关的概念。 数字滤波器设计采用角频率,如何与实际信号频率对应? 角频率w,采样频率fs ,实际信号频率f的转换关系为: W = 2pi f / fs 采样频率的角频率为 2 *pi. 数字滤波器...

linmin418
2017/05/16
0
0
如何消除IIR滤波器的过冲效应呢?

例如一个5Hz的IIR高通滤波器,在滤除小于5Hz频率成分时,会造成过冲效应,就是如一段幅度全部大于0,频率在4Hz左右的正弦信号,经过5Hz的IIR高通滤波后,信号的波峰部分基本被滤掉,但是由于...

中中中
2013/02/24
320
4
SP++3.0已发布,欢迎大家使用(同心协力,共创开源)

SP++ (Signal Processing in C++) 是一个关于信号处理与数值计算的开源C++程序库,该库提供了信号处理与数值计算中常用算法的C++实现。SP++中所有算法都以C++类模板方法实现,以头文件形式组...

张明
2011/02/12
0
55
SP++3.0 发布,欢迎大家使用

消息来自 Jerry 的博客: SP++ (Signal Processing in C++) 是一个关于信号处理与数值计算的开源C++程序库,该库提供了信号处理与数值计算中常用算法的C++实现。SP++中所有算法都以C++类模板...

红薯
2011/02/12
4.5K
4

没有更多内容

加载失败,请刷新页面

加载更多

MySQL autocommit探究

-- sessionA:tx_isolation=REPEATABLE-READmysql> select connection_id();+-----------------+| connection_id() |+-----------------+| 28 |+-----------------+......

安小乐
7分钟前
4
0
c++多线程锁 Mutex  自动判断死锁

c++多线程锁可以使用absl::Mutex std::mutex这两种,下面是demo代码。 使用absl:Mutex的时候打印: [mutex.cc : 1338] RAW: Cycle: [mutex.cc : 1352] RAW: mutex@0x683b68 stack: @ 0x43856......

青黑
26分钟前
1
0
Blockathon2018(成都站)比赛落幕,留给我们这些区块链应用思考

9月14日,HiBlock区块链社区主办的第二届Blockathon在成都菁融国际广场成功举行,30名参赛者分为5支队伍在48小时内完成区块链项目的创意、开发及路演,经过紧张的开发及现场评选,最终币托(...

HiBlock
31分钟前
0
0
71.告警系统主脚本 配置文件 监控项目

20.20 告警系统主脚本(main.sh) 20.21 告警系统配置文件 20.22 告警系统监控项目 20.20 告警系统主脚本(main.sh): ~1.约定:把以后所有的shell脚本放在/usr/local/sbin下,也方便我们查...

王鑫linux
38分钟前
0
0
装饰者模式

装饰者模式 Q:何为装饰模式? ()地给一个对象添加一些额外的(),并且()时,并不影响原对象。扩展功能来说,装饰器模式相比生成子类更为灵活。 Q:使用场景? 1.想要在不影响其他对象的情况下...

阿元
58分钟前
0
0

没有更多内容

加载失败,请刷新页面

加载更多

返回顶部
顶部