文档章节

IIR数字滤波器实现

abcijkxyz
 abcijkxyz
发布于 2016/11/22 16:46
字数 1232
阅读 28
收藏 0
点赞 0
评论 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
粉丝 60
博文 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
Android 怎么实现变声呀   如:男声变女声,男声变机器人声,女声变童声这样的用户体验效果.

实现魔音的思路: 1)获取通话音频流.2)把获取到的音频流通过某种算法,某种数字处理技术,改变原有音频流.3)将改变后的音频流传出去,让通话对方听到改变后的音频流,感觉到说话方声音的变化. ①:...

wang.shuai
2011/07/27
13.1K
4
MIMO声信号处理

1 引言 1.1 MIMO声信号处理 一方面,我们面临着更加复杂的声环境。麦克风接收到的信号,目标语音混杂着以下四种干扰: 1 噪声。 2 回声。回声的产生是由于扬声器和麦克风的耦合。回声的存在会...

chenxiao60
2016/07/07
88
0
重磅!阿里开源自研语音识别模型DFSMN,准确率高达96.04%

阿里开源语音识别模型DFSMN 在近期举行的云栖大会武汉峰会上,装有DFSMN语音识别模型的“AI收银员”在与真人店员的PK中,在嘈杂环境下准确识别了用户的语音点单,在短短49秒内点了34杯咖啡。...

技术小能手
06/08
0
0
有限长滤波器在血压测量中的运用

通过ADC采集到的位于袖袋内的压力传感器的电压值,我们可以换算得到对应的压力值。然后根据血压算法,找到对应的收缩压和舒张压。系统实现的核心就是滤波和寻找算法。在ADC采样之前加入一个硬...

兔之
2014/03/12
0
0
重磅!MaxCompute助力阿里开源自研语音识别模型DFSMN,准确率高达96.04%

阿里开源语音识别模型DFSMN 在近期举行的云栖大会武汉峰会上,装有DFSMN语音识别模型的“AI收银员”在与真人店员的PK中,在嘈杂环境下准确识别了用户的语音点单,在短短49秒内点了34杯咖啡。...

隐林
06/22
0
0

没有更多内容

加载失败,请刷新页面

加载更多

下一页

about git flow

  昨天元芳做了git分支管理规范的分享,为了拓展大家关于git分支的认知,这里我特意再分享这两个关于git flow的链接,大家可以看一下。 Git 工作流程 Git分支管理策略   git flow本质上是...

qwfys
今天
2
0
Linux系统日志文件

/var/log/messages linux系统总日志 /etc/logrotate.conf 日志切割配置文件 参考https://my.oschina.net/u/2000675/blog/908189 dmesg命令 dmesg’命令显示linux内核的环形缓冲区信息,我们可...

chencheng-linux
今天
1
0
MacOS下给树莓派安装Raspbian系统

下载镜像 前往 树莓派官网 下载镜像。 点击 最新版Raspbian 下载最新版镜像。 下载后请,通过 访达 双击解压,或通过 unzip 命令解压。 检查下载的文件 ls -lh -rw-r--r-- 1 dingdayu s...

dingdayu
今天
1
0
spring boot使用通用mapper(tk.mapper) ,id自增和回显等问题

最近项目使用到tk.mapper设置id自增,数据库是mysql。在使用通用mapper主键生成过程中有一些问题,在总结一下。 1、UUID生成方式-字符串主键 在主键上增加注解 @Id @GeneratedValue...

北岩
今天
2
0
告警系统邮件引擎、运行告警系统

告警系统邮件引擎 cd mail vim mail.py #!/usr/bin/env python#-*- coding: UTF-8 -*-import os,sysreload(sys)sys.setdefaultencoding('utf8')import getoptimport smtplibfr......

Zhouliang6
今天
1
0
Java工具类—随机数

Java中常用的生成随机数有Math.random()方法及java.util.Random类.但他们生成的随机数都是伪随机的. Math.radom()方法 在jdk1.8的Math类中可以看到,Math.random()方法实际上就是调用Random类...

PrivateO2
今天
3
0
关于java内存模型、并发编程的好文

Java并发编程:volatile关键字解析    volatile这个关键字可能很多朋友都听说过,或许也都用过。在Java 5之前,它是一个备受争议的关键字,因为在程序中使用它往往会导致出人意料的结果。在...

DannyCoder
昨天
1
0
dubbo @Reference retries 重试次数 一个坑

在代码一中设置 成retries=0,也就是调用超时不用重试,结果DEBUG的时候总是重试,不是0吗,0就不用重试啊。为什么还是调用了多次呢? 结果在网上看到 这篇文章才明白 https://www.cnblogs....

奋斗的小牛
昨天
2
0
数据结构与算法3

要抓紧喽~~~~~~~放羊的孩纸回来喽 LowArray类和LowArrayApp类 程序将一个普通的Java数组封装在LowArray类中。类中的数组隐藏了起来,它是私有的,所以只有类自己的方法才能访问他。 LowArray...

沉迷于编程的小菜菜
昨天
1
0
spring boot应用测试框架介绍

一、spring boot应用测试存在的问题 官方提供的测试框架spring-boot-test-starter,虽然提供了很多功能(junit、spring test、assertj、hamcrest、mockito、jsonassert、jsonpath),但是在数...

yangjianzhou
昨天
8
0

没有更多内容

加载失败,请刷新页面

加载更多

下一页

返回顶部
顶部