基于MATLAB和DSP的IIR滤波器的设计与仿真_图文_百度文库
两大类热门资源免费畅读
续费一年阅读会员，立省24元！
评价文档：
基于MATLAB和DSP的IIR滤波器的设计与仿真
调​试​器
阅读已结束，如果下载本文需要使用
想免费下载本文？
下载文档到电脑，查找使用更方便
还剩1页未读，继续阅读
你可能喜欢 上传我的文档
 下载
 收藏
该文档贡献者很忙，什么也没留下。
 下载此文档
正在努力加载中...
改进的直接I型IIR数字滤波器的DSP实现
下载积分：100
内容提示：改进的直接I型IIR数字滤波器的DSP实现
文档格式：PDF|
浏览次数：89|
上传日期： 04:19:04|
文档星级：
该用户还上传了这些文档
改进的直接I型IIR数字滤波器的DSP实现
官方公共微信STM32实现IIR滤波器，可用matlab生成的头文件
(amoBBS 阿莫电子论坛) -
STM32实现IIR滤波器，可用matlab生成的头文件
放假实在无聊，即将到来的高三非常恐怖，先偷闲一把。
matlab的fdatool是好东西，不过很多人不知道该怎么使用它生成的C头文件。
趁着放假有时间，摸索了几天，终于搞定。希望阿莫给条裤子。
该程序已经用于心电采集实验
导联aVF，带宽1-25Hz
/bbs_upload782111/files_31/ourdev_569832.JPG
实验过程中图片 (原文件名:DSCF6003.JPG)
/bbs_upload782111/files_31/ourdev_569833.jpg
液晶截图 (原文件名:aVF_LCD.jpg)
不多说，切入正题
这里有个fdatool设计的IIR高通滤波器，采样率400Hz时截止频率1Hz。
设计定型之后，要做些调整。
以下说明中的英文名词有些可能对不上fdatool界面上的原文，请大家意会吧
点击菜单中的Edit-&Convert Structure 选择Direct Form I ，SOS，（必须是Direct Form I，II不行）
一般情况下，按照默认设置，fdatool设计都是由二阶部分串联组成的。
这种结构的滤波器稳定性比一个section的要好很多，其他方面的性能也好些。
如果不是的话，点击Convert to second order sections。
这时，滤波器的结构（structure）应该显示为 Direct Form I，second order sections
选择quantize filter，精度选择single precision floating point （单精度浮点）
之所以不用定点是因为噪声太大，也不容易稳定。
点击菜单中的Targets -& generate c header ，选择export as：single precision floating point （单精度浮点）
填写变量名称时，把NUM改成IIR_B，DEN改成IIR_A,其他不用动，保存为iir_coefs.h
保存好的文件如下：
//一大堆注释
/* General type conversion for MATLAB generated C-code*/
#include &tmwtypes.h&
* Expected path to tmwtypes.h
* C:\Program Files\MATLAB\R2010a\extern\include\tmwtypes.h
* Warning - Filter coefficients were truncated to fit specified data type.
*& &The resulting response may not match generated theoretical response.
*& &Use the Filter Design & Analysis Tool to design accurate
*& &single-precision filter coefficients.
#define MWSPT_NSEC 9
const int NL = { 1,3,1,3,1,3,1,3,1 };
const real32_T IIR_B = {
& &0.,& && && && &0,& && && && &0
& && && && && & 1,& && && && & -2,& && && && &1
& &0.,& && && && &0,& && && && &0
& && && && && & 1,& &-1.,& && && && &1
& &0.,& && && && &0,& && && && &0
& && && && && & 1,& & -1.,& && && && &1
& &0.,& && && && &0,& && && && &0
& && && && && &1,& & -1.,& && && && &1
& && && && && & 1,& && && && &0,& && && && &0
const int DL = { 1,3,1,3,1,3,1,3,1 };
const real32_T IIR_A = {
& && && && && & 1,& && && && &0,& && && && &0
& && && && && & 1,& &-1.,& &0.
& && && && && & 1,& && && && &0,& && && && &0
& && && && && & 1,& &-1.,& &0.
& && && && && & 1,& && && && &0,& && && && &0
& && && && && & 1,& &-1.,& &0.
& && && && && & 1,& && && && &0,& && && && &0
& && && && && & 1,& &-1.,& &0.
& && && && && & 1,& && && && &0,& && && && &0
打开iir_coefs.h把MWSPT_NSEC替换成IIR_NSEC,
NL、DL数组删除掉，real32_T改成float ，
其中有一个#include &twmtypes.h&,不要它了，删掉
改完的文件如下：
#define IIR_NSEC 9
//原来叫做MWSPT_NSEC
const float IIR_B = {
//为什么改为float很明显了吧
& & 0.,& && && && &0,& && && && &0
& && && && && & 1,& && && && &-2,& && && && &1
& & 0.,& && && && &0,& && && && &0
& && && && && & 1,-1.,& && && && &1
& & 0.,& && && && &0,& && && && &0
& && && && && & 1,& & -1.,& && && && &1
& &0.,& && && && &0,& && && && &0
& && && && && & 1,& & -1.,& && && && &1
& && && && && & 1,& && && && &0,& && && && &0
const float IIR_A = {
& && && && && & 1,& && && && &0,& && && && &0
& && && && && & 1,-1..
& && && && && & 1,& && && && &0,& && && && &0
& && && && && & 1,-1..
& && && && && & 1,& && && && &0,& && && && &0
& && && && && & 1,-1..
& && && && && & 1,& && && && &0,& && && && &0
& && && && && & 1,-1..
& && && && && & 1,& && && && &0,& && && && &0
保存文件，然后使用以下代码进行滤波
这段代码是根据Direct Form I 2阶IIR滤波的差分方程编写的
a0*y = b0*x + b1*x + b2*x - a1*y -a2*y;
//iir_filter.c
#include &datatype.h&
#include &iir_filter.h&
#include &iir_coefs.h&
int16 iir_filter(int16 in)
& && &uint16
& && &for(i=0;i&IIR_NSEC;i++)
& && && & y =x*IIR_B+x*IIR_B+x*IIR_B-y*IIR_A-y*IIR_A;
& && && & y /= IIR_A;
& && && & y=y;y=y;
& && && & x=x;x=x;
& && && & x =
& && &if( x&32767)x=32767;
& && &if( x&-32768) x=-32768;
& && &return((int16)x);& &
//复位滤波器
void iir_reset(void)
& & uint16 i,j;
& & for(i=0;i&IIR_NSEC+1;i++)
& && &for(j=0;j&3;j++)
& && && & x=0;
& & for(i=0;i&IIR_NSEC;i++)
& && &for(j=0;j&3;j++)
& && && & y=0;
//iir_filter.h
#ifndef _IIR_FILTER_H__
#define _IIR_FILTER_H__
int16 iir_filter(int16 x);
void iir_reset(void);
首先写好iir_coefs.h，然后调用iir_filter.c对数据流进行滤波
一个伪代码例子：
while(运行中)
保存到SD卡(iir_filter(读取ADC采样值()));
这个函数比STM32 DSP库中的函数要好很多，DSP库中的2个IIR滤波函数都不能连续处理数据流。
记得在开始滤波之前重置滤波器
iir_reset();
沙发佩服高中生就会这这些了 但愿学习没有拉下
franklinjin
很不简单啊，楼主高中能搞明白什么是IIR滤波器，是个人才
高中生……
我不如楼主。
沙发佩服高中生就会这这些了 但愿学习没有拉下
huohuansong
回复【2楼】franklinjin
很不简单啊，楼主高中能搞明白什么是iir滤波器，是个人才
-----------------------------------------------------------------------
是啊，怀疑原理搞懂没。
不过这样做iir stm32够呛
回复【9楼】yzak_juel
回复【2楼】franklinjin
很不简单啊，楼主高中能搞明白什么是iir滤波器，是个人才
-----------------------------------------------------------------------
是啊，怀疑原理搞懂没。
不过这样做iir stm32够呛
-----------------------------------------------------------------------
我并不记得系数怎么计算了。。。只记得各种算法的特性，反正有fdatool用着
IIR基本上靠椭圆函数法，FIR则靠等纹波，这2种算法的各项特性都比较好。
然后根据系数和结构把方程写出来，改成程序就ok了。
如果是传递函数也能变换成系数。
stm32还算好，几百条指令而已，IIR阶数少，吃得消的，而且采样率越高IIR需要的阶数越少，计算量增加不多。
尽量不要用定点系数，本底噪声太大，即使是int32型也很厉害，经常比要滤掉的噪声还大。FIR则可以用定点
楼主你把你做这个的完整资料来一份吧 我佩服高中生。。。
回复【10楼】warmonkey
-----------------------------------------------------------------------
回复【11楼】rayz82
楼主你把你做这个的完整资料来一份吧 我佩服高中生。。。
-----------------------------------------------------------------------
整个工程就不发过来了
帖子里的代码复制出来直接就能用的，作为一个模块。
对AD采样进行滤波效果非常好，能够有效减弱特定的干扰
crazy_stone
楼主是个MIT的料！
real_sugar
楼主方法用的挺前沿，在掌握点基础就MIT了~
ggyyll8683
呵呵,看到你的帖子不错!!
shenzhenlang
楼主真不错！
你厉害高中生都懂这些了,高三不要读了,直接申请进大学吧.
可以直接读研了。
回复【20楼】showgu
可以直接读研了。
-----------------------------------------------------------------------
楼主真是个人才，直接读研吧
人才啊，想办法出国吧
现在的高中生真实不得了啊，真实很羡慕。这么年轻就发现了自己的兴趣方向。
shenzhenlang
建议置酷贴
记得我高二研究的是linux，
matlib只看了看语法，
高三考完了才用matlab，siclab,g什么的弄了弄fft，架了个数据库分析股票，
chenyixing
条件所致,走的早和走的远没什么必然联系
myfriend6042
试了下，楼主的程序有问题，
编译通不过的
我上传个经过修改整理的
把滤波器常数直接放到了iir.h
iirourdev_600415YQIMLR.rar(文件大小:2K) (原文件名:iir.rar)
/bbs_upload782111/files_35/ourdev_600416WCHTCZ.jpg
效果 (原文件名:j.jpg)
LS的仁兄说一下报了什么错？我这个代码在发上来之前刚刚编译通过了。
另外，系数还是要分开放，iir.h文件中放置系数可能会导致错误，必须保证整个工程中iir_coefs.h只被iir.c包含。
for(i=0;i&IIR_NSEC;i++)
& && && & y =x*IIR_B+x*IIR_B+x*IIR_B-y*IIR_A-y*IIR_A;
& && && & y /= IIR_A;
& && && & y=y;y=y;
& && && & x=x;x=x;
& && && & x =
x,IIR_B也是数组.
这里面全是数组相乘相加，我用arduino(gcc-avr编译器)编译报错，
根据数据结构这里应该是
x,IIR_B，y
我只是把系数全部复制到iir.h中，不再include原来导出的.h
另外刚才在atmega168上试用这个处理个29阶的，可能是溢出了，芯片资源要求还是比较多的
另外，这个库是你自己写的还是其他地方找的？找了很久都没找到c的算法。你有没有fir的库？
大囧。。。严重失误。。。我复制错了文件
这个是正确代码，也包含了直接型FIR滤波器。
FIR&IIR Filterourdev_600574LYM1IF.rar(文件大小:2K) (原文件名:filter.rar)
核心部分：
//iir_filter.c
#include &../platform.h&
//#include &iir_coefs.float.flat.h&
//#include &iir_coefs.float.sharp.h&
#include &iir_coefs_pass@2Hz_stop@0.8Hz.h&
#include &iir_filter.h&
//IIR_NSEC阶直接型II IIR滤波器
//IIR_NSEC个二阶biquad串联
int16 iir_filter(int16 in)
& && &uint16
& && &for(i=0;i&IIR_NSEC;i++)
& && &//y = x*IIR_B +x*IIR_B +x*IIR_B-y*IIR_A-y*IIR_A;
& && && & y = 0;
& && && & if(IIR_B == 1) y+=x;
& && && & else if(IIR_B == -1) y-=x;
& && && & else if(IIR_B == -2) y=y-x-x;
& && && & else if(IIR_B == 0);& &
& && && & else y += x*IIR_B;
& && && & if(IIR_B == 1) y+=x;
& && && & else if(IIR_B == -1) y-=x;
& && && & else if(IIR_B == -2) y=y-x-x;
& && && & else if(IIR_B == 0);
& && && & else y += x*IIR_B;
& && && & if(IIR_B == 1) y+=x;
& && && & else if(IIR_B == -1) y-=x;
& && && & else if(IIR_B == -2) y=y-x-x;
& && && & else if(IIR_B == 0);
& && && & else y += x*IIR_B;
& && && & if(IIR_A == 1) y-=y;
& && && & else if(IIR_A == -1) y+=y;
& && && & else if(IIR_A == -2) y=y+y+y;
& && && & else if(IIR_A == 0);
& && && & else y -= y*IIR_A;
& && && & if(IIR_A == 1) y-=y;
& && && & else if(IIR_A == -1) y+=y;
& && && & else if(IIR_A == -2) y=y+y+y;
& && && & else if(IIR_A == 0);
& && && & else y -= y*IIR_A;
& && && & if(IIR_A != 1) y /= IIR_A;
& && && & y=y;y=y;
& && && & x=x;x=x;
& && && & x =
& && &if( x&32767)x=32767;
& && &if( x&-32768) x=-32768;
& && &return((int16)x);& &
//复位滤波器
void iir_reset(void)
& & uint16 i,j;
& & for(i=0;i&IIR_NSEC+1;i++)
& && & for(j=0;j&3;j++)
& && && & x=0;
& & for(i=0;i&IIR_NSEC;i++)
& && & for(j=0;j&3;j++)
& && && & y=0;
如果是给AVR单片机应用，请把float改成int32，并且在fdatool中进行系数舍入和稳定性检验。
由于int32的系数经常通不过稳定性检验（红字显示unstable），我才用float
这个库完全是由我自己完成的
/bbs_upload782111/files_35/ourdev_CLW.jpg
新的效果图 (原文件名:DSCF6773.jpg)
你改了之后多出来的if是干什么的？好像意思和原来一样的，提高代码效率？
const float IIR_B = {
& & 0.,& && && && &0,& && && && &0
& && && && && & 1,& && && && &2,& && && && &1
& & 0.,& && && && &0,& && && && &0
& && && && && & 1,& && && && &2,& && && && &1
& & 0.,& && && && &0,& && && && &0
& && && && && & 1,& && && && &2,& && && && &1
& & 0.,& && && && &0,& && && && &0
& && && && && & 1,& && && && &2,& && && && &1
& && && && && & 1,& && && && &0,& && && && &0
};代码中有else if(IIR_B == -2)，但是我的参数表中有2没有-2，只不过比较了下，输出值和我修改的那个是完全一样
if是提高效率的，为了改善通用性，我直接全部写上了。目前我生成出来的系数记得只有1.0.-1.-2这4个常见的整数
佩服楼主，这个年纪就这么厉害了。
FIR.c好像不太对啊。。
根据fir滤波器的公式y(n)=∑h(m)x(n-m);(m: 0～(N-1))，楼主的程序和这个公式不太一样哦
FIR是正确的
int16 fir_filter(uint16 taps,const int16 *coefs,int32 coefs_div,int16 *input)
& & uint16
& & int64 sum=0;
& & for(count=0;count&count++)
& && &sum += input*
//其实是一样的,在t=n时调用函数,sum == y
// = h(0)*x(n) + h(1)*x(n-1) + ... + h(m)*x(n-m) + ... + h(taps-1)*x(n-taps-1)
//也就是m=，而h(m) == coefs& && &--楼上写的N应该为抽头系数个数
//taps（抽头系数个数）== 阶数N+1
//表达式乱七八糟，千万别吐槽了～
& & sum /= coefs_
& & if(sum&32767) sum=32767;
& & if(sum&-32768) sum=-32768;
http://zh.wikipedia.org/zh-cn/%E6%9C%89%E9%99%90%E8%84%89%E5%86%B2%E5%93%8D%E5%BA%94
我不如楼主。
回复【44楼】warmonkey
-----------------------------------------------------------------------
有空帮研究一下 maximally flat FIR low pass filter的公式？居然没找到哪个和matlab对应的。。
带内最平坦是butterworth吧。。。如果要带内带外都平坦，就得损失截止的锐利程度
一般把带内搞定就可以了
mark一下，真的太佩服楼主了。
有相位失真的，对波形要求不高的情况下，IIR的确很好
不过真没找到如何导出matlab的方程的方法
如果是FIR肯定得几百阶才能做出来，延迟太大。椭圆函数法设计的高通，相位特性较其他方法好。例如这次应用的1Hz低通，在0.8Hz处相位滞后不到10度，也就是不到23ms的延迟，问题不大。
如果把数据保存下来进行处理，可以做成零相移滤波，但这里是实时监护。
现在已经换成了3次函数插值法，高通仅仅是纠正AD的偏移，所以现在已经算是彻底解决了延迟问题
回复【49楼】powerpan
-----------------------------------------------------------------------
好像只能导出参数为头文件
导出C好像比较复杂
我没找到导出C的选项，似乎matlab不支持这个功能。
下次用汇编写个FIR代码，库里面那个限制太多了。
太强了，有点像研究生论文，MARK!
给力的哇，楼上多个人参加探讨就是不一样，可惜我还不行
回复【46楼】powerpan
回复【44楼】warmonkey
-----------------------------------------------------------------------
有空帮研究一下 maximally flat fir low pass filter的公式？居然没找到哪个和matlab对应的。。
-----------------------------------------------------------------------
汗，记错了，IIR的maximally flat 才是butterworth，fir的不清楚。其实这个可以看fdatool的源码。
另外，simulink模型是可以导出为C的，也就是说要先把滤波器导出到simulink，再用simulink变换成C源码
...高中用fir，matlab？神啦！
哦，提醒一下，关于心电图还有一种经典算法，自适应滤波。
回复【10楼】warmonkey
回复【9楼】yzak_juel& &
回复【2楼】franklinjin& &
很不简单啊，楼主高中能搞明白什么是iir滤波器，是个人才
-----------------------------------------------------------------------
是啊，怀疑原理搞懂没。
不过这样做iir stm32够呛
-----------------------------------------------------------------------
我并不记得系数怎么计算了。。。只记得各种算法的特性，反正有fdatool用着
iir基本上靠椭圆函数法，fir则靠等纹波，这2种算法的各项特性都比较好。
然后根据系数和结构把方程写出来，改成程序就ok了。
如果是传递函数也能变换成系数。
stm32还算好，几百条指令而已，iir阶数少，吃得消的，而......
-----------------------------------------------------------------------
佩服楼主的能力，如果将来发展好的话，肯定是个了不起的人才。
不过我还是要泼点冷水：一定要注重基础知识！
你现在用的工具，毕竟是别人的。我怀疑你是否参透了其中最基本的原理。不要光知道怎么做，不知道为什么，只知道皮毛，不知道深入的原理。
我的经历和你相似，从小动手能力很强，而且很早接触电子。大二的时候，我的电设作品，全国一等奖，而且是得分最高的，也就是说是第1名。但同时，我对自己的理论水平不满，我放下手里的东西，研究了一年的理论知识。于是，很自然的觉得自己以前很肤浅。之后的实践，因为有了理论的基础，自然有了质的飞跃。
我同样见证过，理论知识欠缺而实践能力很强的工程师，悲剧的结果，后期的发展非常受限。
说什么，工程师不需要很高的理论水平，纯属悲剧人物发表的悲剧的言论。
看到你在做数学信号方面的课题，留给你几个问题，看看你对DSP是不是真的理解了：
(1)z变换和频率特性的关系;
(2)z变换和稳定性的关系;
(3)s域传递函数和z域的传递函数联系，如何转化;
(4)DFT、FFT、离散信号傅里叶变换、连续信号傅里叶变换，四者的关系及相应的频率和时域特点;
(5)什么是窗函数，有什么特点.
以上不是考试题目，考题比这个要难。但如果是真正把理论和实践同时掌握好的工程师，一定不难回答。
onyourmark66
还真的不太清楚，我尝试着回答一下吧：
(1)z变换和频率特性的关系;
变换为差分方程，计算离散化的白噪声信号激励下，输出进行FFT获得频谱。
或者对于每个频点，计算离散的正弦波激励下，输出的幅度。
(2)z变换和稳定性的关系;
极点都在单位圆内则稳定。
(3)s域传递函数和z域的传递函数联系，如何转化;
s域传函转换为典型环节的组合，再用Z域的典型环节表示。
(4)DFT、FFT、离散信号傅里叶变换、连续信号傅里叶变换，四者的关系及相应的频率和时域特点;
连续信号傅立叶变换把连续信号从时域变换到频域，也就是把输入分解为一系列正弦波的叠加。
DFT（离散信号傅里叶变换）把离散信号从时域变换到频域。
fft是dft的快速算法，利用信号的连续性，简化了计算。
(5)什么是窗函数，有什么特点.
将一段信号“变换”为无限长，不断重复的信号。窗函数会导致信号频谱发生变化。不同窗函数对原始信号的“歪曲”不同。
回复【61楼】warmonkey
了解LZ的情况了
特别提一下(1)、(3)和(4)：
(1)z变换和频率特性的关系;
其实将z = e ^ (jω / fs)代入z变换中，就得到了频幅特性;
(3)s域传递函数和z域的传递函数联系，如何转化;
设计IIR滤波器的时候，简单的可以用脉冲响应不变法来转化;
至于两者的关系：
s和z都是复数算子，在s变化中，e ^ (-st) = e ^ (σ + jω) * t，
z变化中，z ^ (-n) = e ^ (σ + jω) / fs * n.
因此，当σ = 0时，s平面中，正好位于虚轴，& & z平面中，正好位于单位圆上；
& && &当σ & 0时，s平面中，正好位于左半平面，z平面中，正好位于单位圆内；
& && &当σ & 0时，s平面中，正好位于右半平面，z平面中，正好位于单位圆外.
(注意，这里的s平面指的是s平面上从 -fs/2 到 fs/2 之间的部分)
至于z域传递函数和稳定性的关系，对应的用s域来解释就太容易了，那就是e ^ σt这个指数函数的收敛性问题了.
(4)dft、fft、离散信号傅里叶变换、连续信号傅里叶变换，四者的关系及相应的频率和时域特点;
“DFT”和“离散信号傅里叶变换”是两回事：“DFT”是离散化的傅里叶变化，是对有限长的离散信号而言的;“离散信号傅里叶变换”是对无限长的离散信号而言的.
再加上傅里叶级数，这个问题就完整了. 另外，FFT是DFT在K = 2 ^ N时的简单情况，放在一边不做讨论.
四者的关系如下：
& && && && && && && && &输入& && && && && && && && && &输出
傅里叶级数& && && && && & 有限长连续信号(周期信号)& && & 无限长离散频谱& && && && &
连续信号傅里叶变换& && &无限长连续信号& && && && && &无限长连续频谱
离散信号傅里叶变换& && &无限长离散信号& && && && && &周期性连续频谱
DFT& && && && && && && &有限长离散信号(视为周期信号)& &周期性离散频谱
注：个人理解.
myworkmail
myworkmail
回复【62楼】mitchell
-----------------------------------------------------------------------
汗，看来没仔细学过离散系统确实不行。。。求大虾联系方式，以后多聊聊，谢谢了。
回复【62楼】mitchell
-----------------------------------------------------------------------
麻省大神求拜师学艺。。。
大神啊，都是。
回复【65楼】warmonkey
-----------------------------------------------------------------------
不敢当, 我QQ: , 随时交流
回复【67楼】reloaded 电子浪人
-----------------------------------------------------------------------
您太客气了，mitchell是我中文名字的音译
回复【楼主位】warmonkey
-----------------------------------------------------------------------
“这个函数比STM32 DSP库中的函数要好很多，DSP库中的2个IIR滤波函数都不能连续处理数据流。 ”
可以吧？这个biquad 代码写的相当不错，可以随意剪成1个sector或者增加更多。
pengknight
神人~mark！
tangwei039
从楼主的这只手来看不像高中生的...
现在matlab 主推基于模型的设计 m3内核的芯片也开始支持，第三方软件rapidstm32
回复【77楼】tiandewen
-----------------------------------------------------------------------
这个不错，下来看看。
佩服，研究一下
xiezheming
kingboy100
selina1983
mark......................
chao8828276
服得很，不得不顶
wangda6222
回复【楼主位】warmonkey
-----------------------------------------------------------------------
楼主莫怪啊，我想问问你是用什么采集心跳信号的，用的换能器？
楼主真是个人才
MARK，楼主强人
高中？佩服！！！
回复【89楼】chao8828276
心电信号的实际波行我不太清楚，不过我也搞过fir，我没搞过算法，我也是用matlab产生的系数，但是32阶的fir来处理已经效果很好了，我当时仿真了一下，32阶处理掉50hz以及以上的频率信号是没问题的，当然截止频率要很低
-----------------------------------------------------------------------
我说的是高通滤波
楼主，请教个问题
#define MWSPT_NSEC 9
这个是什么意思
那个数组为什么定义成 MWSPT_NSEC行3列？这个有点没看懂
查看完整版本:【原创】【安富莱——DSP教程】第40章 IIR滤波器的实现 - STM32 - 意法半导体STM32/STM8技术社区
后使用快捷导航没有帐号？
查看: 890|回复: 8
【原创】【安富莱——DSP教程】第40章 IIR滤波器的实现
主题帖子积分
金牌会员, 积分 1963, 距离下一级还需 1037 积分
金牌会员, 积分 1963, 距离下一级还需 1037 积分
特别说明：完整45期数字信号处理教程，原创高性能示波器代码全开源地址：
第40章 IIR滤波器的实现
& & 本章节讲解IIR滤波器直接I型的低通，高通，带通和带阻滤波器的实现。& & 40.1 IIR滤波器介绍& & 40.2 Matlab工具箱fdatool生成IIR滤波器系数& & 40.3 IIR低通滤波器设计& & 40.4 IIR高通滤波器设计& & 40.5 IIR带通滤波器设计& & 40.6 IIR带阻滤波器设计& & 40.7 总结
40.1 IIR滤波器介绍& & ARM官方提供的直接I型IIR库支持Q7，Q15，Q31和浮点四种数据类型。其中Q15和Q31提供了基于Cortex-M3和Cortex-M4的快速版本。& & 直接I型IIR滤波器是基于二阶Biquad级联的方式来实现的。每个Biquad由一个二阶的滤波器组成：y[n] = b0 * x[n] + b1 * x[n-1] + b2 * x[n-2] + a1 * y[n-1] + a2 * y[n-2]& & 直接I型算法每个阶段需要5个系数和4个状态变量。
40.1.png (17.25 KB, 下载次数: 0)
10:04 上传
& & 这里有一点要特别的注意，有些滤波器系数生成工具是采用的下面公式实现：y[n] = b0 * x[n] + b1 * x[n-1] + b2 * x[n-2] - a1 * y[n-1] - a2 * y[n-2]& & 比如matlab就是使用上面的公式实现的，所以在使用fdatool工具箱生成的a系数需要取反才能用于直接I型IIR滤波器的函数中。& & 高阶IIR滤波器的实现是采用二阶Biquad级联的方式来实现的。其中参数numStages就是用来做指定二阶Biquad的个数。比如8阶IIR滤波器就可以采用numStages=4个二阶Biquad来实现。
40.2.png (14.57 KB, 下载次数: 0)
10:04 上传
如果要实现9阶IIR滤波器就需要将numStages=5，这时就需要其中一个Biquad配置成一阶滤波器(也就是b2=0，a2=0)。
主题帖子积分
金牌会员, 积分 1963, 距离下一级还需 1037 积分
金牌会员, 积分 1963, 距离下一级还需 1037 积分
40.2 Matlab工具箱fdatool生成IIR滤波器系数& & 前面介绍FIR滤波器的时候，我们讲解了如何使用fdatool生成C头文件，从而获得滤波器系数。这里不能再使用这种方法了，主要是因为通过C头文件获取的滤波器系数需要通过ARM官方的IIR函数调用多次才能获得滤波结果，所以我们这里换另外一种方法。& & 下面我们讲解如何通过fdatool工具箱生成滤波器系数。首先在matlab的命令窗口输入fdatool就能打开这个工具箱：
40.3.png (3.45 KB, 下载次数: 1)
10:12 上传
fdatool界面打开效果如下：
40.4.png (41.4 KB, 下载次数: 0)
10:12 上传
IIR滤波器的低通，高通，带通，带阻滤波的设置会在下面一一讲解，这里说一下设置后相应参数后如何生成滤波器系数。参数设置好以后点击如下按钮：
40.5.png (41.08 KB, 下载次数: 0)
10:12 上传
点击Design Filter之后，注意左上角生成的滤波器结构：
40.6.png (45.78 KB, 下载次数: 0)
10:12 上传
默认生成的IIR滤波器类型是Direct-Form II, Second-Order Sections（直接II型，每个Section是一个二阶滤波器）。这里我们需要将其转换成Direct-Form I, Second-Order Sections，因为本章使用的IIR滤波器函数是Direct-Form I的结构。& & 转换方法，点击Edit-&Convert Structure，界面如下，这里我们选择第一项，并点击OK：
40.7.png (13.27 KB, 下载次数: 0)
10:12 上传
转换好以后再点击File-Export，第一项选择Coefficient File（ASCII）：
40.8.png (13.45 KB, 下载次数: 0)
10:12 上传
第一项选择好以后，第二项选择Decimal：
40.9.png (10.29 KB, 下载次数: 0)
10:12 上传
两个选项都选择好以后，点击Export进行导出，导出后保存即可：
40.10.png (22.25 KB, 下载次数: 0)
10:12 上传
保存后Matlab会自动打开untitled.fcf文件，可以看到生成的系数：%
% Generated by MATLAB(R) 7.14 and the Signal Processing Toolbox 6.17.
%
% Generated on: 30-Dec-:50
%
% Coefficient Format: Decimal
% Discrete-Time IIR Filter (real)& && && && && && && && && &
% -------------------------------& && && && && && && && && &
% Filter Structure& & : Direct-Form II, Second-Order Sections
% Number of Sections&&: 2& && && && && && && && && && && && &
% Stable& && && && &&&: Yes& && && && && && && && && && && &
% Linear Phase& && &&&: No& && && && && && && && && && && &&&
& && && && && && && && && && && && && && && && && && && && &
SOS Matrix:& && && && && && && && && && && && && && && && &&&
1&&2&&1&&1&&-1.3479& &0.95477& && &&&
1&&2&&1&&1&&-0.175& && &&&
& && && && && && && && && && && && && && && && && && && && &
Scale Values:& && && && && && && && && && && && && && && && &
0.15171& && && && && && && && && && && && && && &
0.34613复制代码由于咱们前面选择的是4阶IIR滤波，生成的结果就是由两组二阶IIR滤波系数组成，系数的对应顺序如下：SOS Matrix:& && && && && && && && && && && && && && && && &&&
1& &2& &1& &1& &-1.3479& &0.95477& && &&&
b0&&b1&&b2&&a0& && &&&a1& && && && && && && & a2
1 2& &1& &1& &-0.175& && &&&
b0&&b1&&b2&&a0& && &&&a1& && && && && && && & a2复制代码注意，实际使用ARM官方的IIR函数调用的时候要将a1和a2取反。另外下面两组是每个二阶滤波器的增益，滤波后的结果要乘以这两个增益数值才是实际结果：0.15171& && && && && && && && && && && && && && &
0.34613复制代码实际的滤波系数调用方法，看下面的例子即可。
主题帖子积分
金牌会员, 积分 2324, 距离下一级还需 676 积分
金牌会员, 积分 2324, 距离下一级还需 676 积分
沙发，支持原创
主题帖子积分
金牌会员, 积分 1963, 距离下一级还需 1037 积分
金牌会员, 积分 1963, 距离下一级还需 1037 积分
40.3 IIR低通滤波器设计& & 本章使用的IIR滤波器函数是arm_biquad_cascade_df1_f32。下面使用此函数设计IIR低通，高通，带通和带阻滤波器。40.3.1 函数arm_biquad_cascade_df1_f32说明函数定义如下：& & void arm_biquad_cascade_df1_f32(& && && && &const arm_biquad_casd_df1_inst_f32 * S,& && && && &float32_t * pSrc,& && && && &float32_t * pDst,& && && &&&uint32_t blockSize)参数定义：& && & [in]&&*S& && && &points to an instance of the floating-point Biquad cascade structure.& & & && & [in]&&*pSrc& && &points to the block of input data.& & & && &[out] *pDst& && &points to the block of output data.& & & && &[in]&&blockSize&&number of samples to process per call.& & & && &return& &&&none.& & 注意事项：结构arm_fir_instance_f32的定义如下（在文件arm_math.h文件）：& && &typedef struct& && &{& && &/**& number of 2nd order stages in the filter.
Overall order is 2*numStages. */& && &uint32_t numS& && && & /**& Points to the array of state coefficients.
The array is of length 4*numStages. */& && &float32_t *pS & && &/**& Points to the array of coefficients.&&The array is of length 5*numStages. */& && && & float32_t *pC & && &} arm_biquad_casd_df1_inst_f32;特别注意，参数pState指向的数组大小要是4倍的numStages，pCoeffs指向的数组大小要是5倍的numStages。1. 参数pCoeffs指向滤波因数，滤波因数数组长度为numTaps。但要注意pCoeffs指向的滤波因数应该按照如下的顺序进行排列：& && & {b10, b11, b12, a11, a12, b20, b21, b22, a21, a22, ...} & & 先放第一个二阶Biquad系数，然后放第二个，以此类推。2. pState指向状态变量数组。3. blockSize 这个参数的大小没有特殊要求，用户只需保证大于1且小于等于采样点个数即可。
40.3.2 fdatool获取低通滤波器系数& & 设计一个如下的例子：& & 信号由50Hz正弦波和200Hz正弦波组成，采样率1Kbps，现设计一个巴特沃斯滤波器低通滤波器，采用直接I型，截止频率80Hz，采样400个数据，滤波器阶数设置为4。fadtool的配置如下：
40.11.png (47.89 KB, 下载次数: 0)
10:21 上传
配置好低通滤波器后，具体滤波器系数的生成大家参考本章第二小节的方法即可。
40.3.3 低通滤波器实现 & & 通过工具箱fdatool获得低通滤波器系数后在开发板上运行函数arm_biquad_cascade_df1_f32来测试低通滤波器的效果。#define numStages&&2& && && && && & /* 2阶IIR滤波的个数 */
#define TEST_LENGTH_SAMPLES&&400& & /* 采样点数 */
static float32_t testInput_f32_50Hz_200Hz[TEST_LENGTH_SAMPLES]; /* 采样点 */
static float32_t testOutput[TEST_LENGTH_SAMPLES];& && && && && &/* 滤波后的输出 */
static float32_t IIRStateF32[4*numStages];& && && && && && && & /* 状态缓存，大小numTaps + blockSize - 1*/
& && && && && && && && && && && && && && && && && && && && && && && && && &&&
/* 巴特沃斯低通滤波器系数 80Hz*/& && && && && && && && && && && && && && && && && && && && && && && && && && && && && && && && && && && && && && && && && && && && && &&&
const float32_t IIRCoeffs32LP[5*numStages] = {
1.0f,&&2.0f,&&1.0f, 1.2167f,&&-0.86178f,& && &&&
1.0f,&&2.0f,&&1.0f, 1.2184f,&&-0.55354f& && && && && && && && &&&
};
/*
*********************************************************************************************************
*& & & & 函 数 名: arm_iir_f32_lp
*& & & & 功能说明: 调用函数arm_iir_f32_lp实现低通滤波器
*& & & & 形& & 参：无
*& & & & 返 回 值: 无
*********************************************************************************************************
*/
static void arm_iir_f32_lp(void)
{
uint32_
arm_biquad_casd_df1_inst_f32 S;
float32_t ScaleV
/* 初始化 */
arm_biquad_cascade_df1_init_f32(&S, numStages, (float32_t *)&IIRCoeffs32LP[0], (float32_t
*)&IIRStateF32[0]);
/* IIR滤波 */
& & & & arm_biquad_cascade_df1_f32(&S, testInput_f32_50Hz_200Hz, testOutput, TEST_LENGTH_SAMPLES);
/*放缩系数 */
ScaleValue = 0.161221f * 0.58381f
/* 打印滤波后结果 */
for(i=0; i&TEST_LENGTH_SAMPLES; i++)
{
printf(&%f\r\n&, testOutput[i]*ScaleValue);
}
}复制代码运行如上函数可以通过串口打印出函数arm_biquad_cascade_df1_f32滤波后的波形数据，下面通过Matlab绘制波形来对比Matlab计算的结果和ARM官方库计算的结果。& & 对比前需要先将串口打印出的一组数据加载到Matlab中， arm_biquad_cascade_df1_f32的计算结果起名sampledata，加载方法在前面的教程中已经讲解过，这里不做赘述了。Matlab中运行的代码如下：fs=1000;& && && && &&&%设置采样频率 1K
N=400;& && && && && &%采样点数& && &
n=0:N-1;
t=n/& && && && && & %时间序列
f=n*fs/N;& && && && &&&%频率序列
x1=sin(2*pi*50*t);
x2=sin(2*pi*200*t);& &&&%50Hz和200Hz正弦波
subplot(211);
plot(t, x1);
title('滤波后的理想波形');
subplot(212);
plot(t, sampledata);
title('ARM官方库滤波后的波形');
复制代码Matlab计算结果如下：
40.12.png (10.43 KB, 下载次数: 0)
10:21 上传
从上面的波形对比来看，matlab和函数arm_biquad_cascade_df1_f32计算的结果基本是一致的。为了更好的说明滤波效果，下面从频域的角度来说明这个问题，Matlab上面运行如下代码：fs=1000;& && && && && & %设置采样频率 1K
N=400;& && && && && &&&%采样点数& && &
n=0:N-1;
t=n/& && && && && && &%时间序列
f=n*fs/N;& && && && && & %频率序列
x = sin(2*pi*50*t) + sin(2*pi*200*t);& && &%50Hz和200Hz正弦波合成
&&
subplot(211);
y=fft(x, N);& && && && && & %对信号x做FFT& &
plot(f,abs(y));
xlabel('频率/Hz');
ylabel('振幅');
title('原始信号FFT');
y3=fft(sampledata, N);& & %经过IIR滤波器后得到的信号做FFT
subplot(212);& && && && && && && && && && &
plot(f,abs(y3));
xlabel('频率/Hz');
ylabel('振幅');
title('IIR滤波后信号FFT');
复制代码Matlab计算结果如下：
40.13.png (9.69 KB, 下载次数: 0)
10:21 上传
上面波形变换前的FFT和变换后FFT可以看出，200Hz的正弦波基本被滤除。
主题帖子积分
金牌会员, 积分 1963, 距离下一级还需 1037 积分
金牌会员, 积分 1963, 距离下一级还需 1037 积分
40.4 IIR高通滤波器设计
40.4.1 fdatool获取高通滤波器系数& & 设计一个如下的例子：& & 信号由50Hz正弦波和200Hz正弦波组成，采样率1Kbps，现设计一个巴特沃斯滤波器高通滤波器，采用直接I型，截止频率140Hz，采样400个数据，滤波器阶数设置为4。fadtool的配置如下：
40.14.png (48.67 KB, 下载次数: 0)
10:27 上传
配置好高通滤波器后，具体滤波器系数的生成大家参考本章第二小节的方法即可。
40.4.2 高通滤波器实现 & & 通过工具箱fdatool获得高通滤波器系数后在开发板上运行函数arm_biquad_cascade_df1_f32来测试高通滤波器的效果。#define numStages&&2& && && && && & /* 2阶IIR滤波的个数 */
#define TEST_LENGTH_SAMPLES&&400& & /* 采样点数 */
static float32_t testInput_f32_50Hz_200Hz[TEST_LENGTH_SAMPLES]; /* 采样点 */
static float32_t testOutput[TEST_LENGTH_SAMPLES];& && && && && &/* 滤波后的输出 */
static float32_t IIRStateF32[4*numStages];& && && && && && && & /* 状态缓存，大小numTaps + blockSize - 1*/
& && && && && && && && && && && && && && && && && && && && && && && && && &&&
/* 巴特沃斯高通滤波器系数 140Hz */& && && && && && && && && && && && && && && && && && && && && && && && && && && && && && && && && && && && && && && && && && && && && &&&
const float32_t IIRCoeffs32HP[5*numStages] = {
1.0f,&&-2.0f,&&1.0f,&&0.1518f,& &-0.81642f,& && &
1.0f,&&-2.0f,&&1.0f,&&0.32121f,&&-0.97309f,& && && && && && && && && && &&&
};
/*
*********************************************************************************************************
*& & & & 函 数 名: arm_iir_f32_hp
*& & & & 功能说明: 调用函数arm_iir_f32_hp实现高通滤波器
*& & & & 形& & 参：无
*& & & & 返 回 值: 无
*********************************************************************************************************
*/
static void arm_iir_f32_hp(void)
{
uint32_
arm_biquad_casd_df1_inst_f32 S;
float32_t ScaleV
/* 初始化 */
arm_biquad_cascade_df1_init_f32(&S, numStages, (float32_t *)&IIRCoeffs32HP[0], (float32_t
*)&IIRStateF32[0]);
/* IIR滤波 */
& & & & arm_biquad_cascade_df1_f32(&S, testInput_f32_50Hz_200Hz, testOutput, TEST_LENGTH_SAMPLES);
& &
/*放缩系数 */& &
ScaleValue = 0.99203f * 0.07356f;&&
/* 打印滤波后结果 */
for(i=0; i&TEST_LENGTH_SAMPLES; i++)
{
printf(&%f\r\n&, testOutput[i]*ScaleValue);
}
}复制代码运行如上函数可以通过串口打印出函数arm_biquad_cascade_df1_f32滤波后的波形数据，下面通过Matlab绘制波形来对比Matlab计算的结果和ARM官方库计算的结果。& & 对比前需要先将串口打印出的一组数据加载到Matlab中， arm_biquad_cascade_df1_f32的计算结果起名sampledata，加载方法在前面的教程中已经讲解过，这里不做赘述了。Matlab中运行的代码如下：fs=1000;& && && && &&&%设置采样频率 1K
N=400;& && && && && &%采样点数& && &
n=0:N-1;
t=n/& && && && && & %时间序列
f=n*fs/N;& && && && &&&%频率序列
x1=sin(2*pi*50*t);
x2=sin(2*pi*200*t);& &&&%50Hz和200Hz正弦波
subplot(211);
plot(t, x2);
title('滤波后的理想波形');
subplot(212);
plot(t, sampledata);
title('ARM官方库滤波后的波形');
复制代码Matlab计算结果如下：
40.15.png (9.78 KB, 下载次数: 0)
10:27 上传
从上面的波形对比来看，matlab和函数arm_biquad_cascade_df1_f32计算的结果基本是一致的。为了更好的说明滤波效果，下面从频域的角度来说明这个问题，Matlab上面运行如下代码：fs=1000;& && && && && & %设置采样频率 1K
N=400;& && && && && &&&%采样点数& && &
n=0:N-1;
t=n/& && && && && && &%时间序列
f=n*fs/N;& && && && && & %频率序列
x = sin(2*pi*50*t) + sin(2*pi*200*t);& && &%50Hz和200Hz正弦波合成
&&
subplot(211);
y=fft(x, N);& && && && && & %对信号x做FFT& &
plot(f,abs(y));
xlabel('频率/Hz');
ylabel('振幅');
title('原始信号FFT');
y3=fft(sampledata, N);& & %经过IIR滤波器后得到的信号做FFT
subplot(212);& && && && && && && && && && &
plot(f,abs(y3));
xlabel('频率/Hz');
ylabel('振幅');
title('IIR滤波后信号FFT');
复制代码Matlab计算结果如下：
40.16.png (9.54 KB, 下载次数: 0)
10:27 上传
上面波形变换前的FFT和变换后FFT可以看出，50Hz的正弦波基本被滤除。
主题帖子积分
金牌会员, 积分 1963, 距离下一级还需 1037 积分
金牌会员, 积分 1963, 距离下一级还需 1037 积分
40.5 IIR带通滤波器设计
40.5.1 fdatool获取低通滤波器系数& & 设计一个如下的例子：& & 信号由50Hz正弦波和200Hz正弦波组成，采样率1Kbps，现设计一个巴特沃斯滤波器带通滤波器，采用直接I型，截止频率140Hz和，采样400个数据，滤波器阶数设置为4。fadtool的配置如下：
40.17.png (49.2 KB, 下载次数: 0)
10:32 上传
配置好带通滤波器后，具体滤波器系数的生成大家参考本章第二小节的方法即可。
40.5.2 带通滤波器实现 & & 通过工具箱fdatool获得带通滤波器系数后在开发板上运行函数arm_biquad_cascade_df1_f32来测试带通滤波器的效果。#define numStages&&2& && && && && & /* 2阶IIR滤波的个数 */
#define TEST_LENGTH_SAMPLES&&400& & /* 采样点数 */
static float32_t testInput_f32_50Hz_200Hz[TEST_LENGTH_SAMPLES]; /* 采样点 */
static float32_t testOutput[TEST_LENGTH_SAMPLES];& && && && && &/* 滤波后的输出 */
static float32_t IIRStateF32[4*numStages];& && && && && && && & /* 状态缓存，大小numTaps + blockSize - 1*/
& && && && && && && && && && && && && && && && && && && && && && && && && &&&
/* 巴特沃斯带通滤波器系数140Hz 400Hz*/& && && && && && && && && && && && && && && && && && && && && && && && && && && && && && && && && && && && && && && && && && && && && &&&
const float32_t IIRCoeffs32BP[5*numStages] = {
1.0f,&&0.0f,&&-1.0f,& & -1.1668f,& &-0.53411f,& && &
1.0f,&&0.0f,&&-1.0f,& & 0.04886f,&&-0.68387f& && && && && && && && && && &
/*
*********************************************************************************************************
*& & & & 函 数 名: arm_iir_f32_bp
*& & & & 功能说明: 调用函数arm_iir_f32_hp实现带通滤波器
*& & & & 形& & 参：无
*& & & & 返 回 值: 无
*********************************************************************************************************
*/
static void arm_iir_f32_bp(void)
{
uint32_
arm_biquad_casd_df1_inst_f32 S;
float32_t ScaleV
/* 初始化 */
arm_biquad_cascade_df1_init_f32(&S, numStages, (float32_t *)&IIRCoeffs32BP[0], (float32_t
*)&IIRStateF32[0]);
/* IIR滤波 */
& & & & arm_biquad_cascade_df1_f32(&S, testInput_f32_50Hz_200Hz, testOutput, TEST_LENGTH_SAMPLES);
& &
/*放缩系数 */& &
ScaleValue = 0.77365f * 0.77365f;&&
/* 打印滤波后结果 */
for(i=0; i&TEST_LENGTH_SAMPLES; i++)
{
printf(&%f\r\n&, testOutput[i]*ScaleValue);
}
}复制代码运行如上函数可以通过串口打印出函数arm_biquad_cascade_df1_f32滤波后的波形数据，下面通过Matlab绘制波形来对比Matlab计算的结果和ARM官方库计算的结果。& & 对比前需要先将串口打印出的一组数据加载到Matlab中， arm_biquad_cascade_df1_f32的计算结果起名sampledata，加载方法在前面的教程中已经讲解过，这里不做赘述了。Matlab中运行的代码如下：fs=1000;& && && && &&&%设置采样频率 1K
N=400;& && && && && &%采样点数& && &
n=0:N-1;
t=n/& && && && && & %时间序列
f=n*fs/N;& && && && &&&%频率序列
x1=sin(2*pi*50*t);
x2=sin(2*pi*200*t);& &&&%50Hz和200Hz正弦波
subplot(211);
plot(t, x1);
title('滤波后的理想波形');
subplot(212);
plot(t, sampledata);
title('ARM官方库滤波后的波形');
复制代码Matlab计算结果如下：
40.18.png (9.7 KB, 下载次数: 0)
10:32 上传
从上面的波形对比来看，matlab和函数arm_biquad_cascade_df1_f32计算的结果基本是一致的。为了更好的说明滤波效果，下面从频域的角度来说明这个问题，Matlab上面运行如下代码：fs=1000;& && && && && & %设置采样频率 1K
N=400;& && && && && &&&%采样点数& && &
n=0:N-1;
t=n/& && && && && && &%时间序列
f=n*fs/N;& && && && && & %频率序列
x = sin(2*pi*50*t) + sin(2*pi*200*t);& && &%50Hz和200Hz正弦波合成
&&
subplot(211);
y=fft(x, N);& && && && && & %对信号x做FFT& &
plot(f,abs(y));
xlabel('频率/Hz');
ylabel('振幅');
title('原始信号FFT');
y3=fft(sampledata, N);& & %经过IIR滤波器后得到的信号做FFT
subplot(212);& && && && && && && && && && &
plot(f,abs(y3));
xlabel('频率/Hz');
ylabel('振幅');
title('IIR滤波后信号FFT');
复制代码Matlab计算结果如下：
40.19.png (9.54 KB, 下载次数: 0)
10:32 上传
上面波形变换前的FFT和变换后FFT可以看出，50Hz的正弦波基本被滤除。
主题帖子积分
金牌会员, 积分 2595, 距离下一级还需 405 积分
金牌会员, 积分 2595, 距离下一级还需 405 积分
主题帖子积分
金牌会员, 积分 1963, 距离下一级还需 1037 积分
金牌会员, 积分 1963, 距离下一级还需 1037 积分
40.6 IIR带阻滤波器设计
40.6.1 fdatool获取带阻滤波器系数& & 设计一个如下的例子：& & 信号由50Hz正弦波和200Hz正弦波组成，采样率1Kbps，现设计一个巴特沃斯滤波器带阻滤波器，采用直接I型，截止频率100Hz和325Hz，采样400个数据，滤波器阶数设置为4。fadtool的配置如下：
40.20.png (49.37 KB, 下载次数: 0)
12:45 上传
配置好带阻滤波器后，具体滤波器系数的生成大家参考本章第二小节的方法即可。
40.6.2 带阻滤波器实现 & & 通过工具箱fdatool获得带阻滤波器系数后在开发板上运行函数arm_biquad_cascade_df1_f32来测试带阻滤波器的效果。#define numStages&&2& && && && && & /* 2阶IIR滤波的个数 */
#define TEST_LENGTH_SAMPLES&&400& & /* 采样点数 */
static float32_t testInput_f32_50Hz_200Hz[TEST_LENGTH_SAMPLES]; /* 采样点 */
static float32_t testOutput[TEST_LENGTH_SAMPLES];& && && && && &/* 滤波后的输出 */
static float32_t IIRStateF32[4*numStages];& && && && && && && & /* 状态缓存，大小numTaps + blockSize - 1*/
& && && && && && && && && && && && && && && && && && && && && && && && && &&&
/* 巴特沃斯带阻滤波器系数100Hz 325Hz*/& && && && && && && && && && && && && && && && && && && && && && && && && && && && && && && && && && && && && && && && && && && && && &&&
const float32_t IIRCoeffs32BS[5*numStages] = {
1.0f,&&-0.35005f,&&1.0f,&&1.7746f,& &-0.21391f,
1.0f,&&-0.35005f,&&1.0f,&&-0.41883f, -0.08849f& && && && && && && && &
/*
*********************************************************************************************************
*& & & & 函 数 名: arm_iir_f32_bs
*& & & & 功能说明: 调用函数arm_iir_f32_bs实现带阻滤波器
*& & & & 形& & 参：无
*& & & & 返 回 值: 无
*********************************************************************************************************
*/
static void arm_iir_f32_bs(void)
{
uint32_
arm_biquad_casd_df1_inst_f32 S;
float32_t ScaleV
/* 初始化 */
arm_biquad_cascade_df1_init_f32(&S, numStages, (float32_t *)&IIRCoeffs32BS[0], (float32_t
*)&IIRStateF32[0]);
/* IIR滤波 */
& & & & arm_biquad_cascade_df1_f32(&S, testInput_f32_50Hz_200Hz, testOutput, TEST_LENGTH_SAMPLES);
& &
/*放缩系数 */& &
ScaleValue = 0.78698f * 0.78698f;&&
/* 打印滤波后结果 */
for(i=0; i&TEST_LENGTH_SAMPLES; i++)
{
printf(&%f\r\n&, testOutput[i]*ScaleValue);
}
}复制代码运行如上函数可以通过串口打印出函数arm_biquad_cascade_df1_f32滤波后的波形数据，下面通过Matlab绘制波形来对比Matlab计算的结果和ARM官方库计算的结果。& & 对比前需要先将串口打印出的一组数据加载到Matlab中， arm_biquad_cascade_df1_f32的计算结果起名sampledata，加载方法在前面的教程中已经讲解过，这里不做赘述了。Matlab中运行的代码如下：fs=1000;& && && && &&&%设置采样频率 1K
N=400;& && && && && &%采样点数& && &
n=0:N-1;
t=n/& && && && && & %时间序列
f=n*fs/N;& && && && &&&%频率序列
x1=sin(2*pi*50*t);
x2=sin(2*pi*200*t);& &&&%50Hz和200Hz正弦波
subplot(211);
plot(t, x1);
title('滤波后的理想波形');
subplot(212);
plot(t, sampledata);
title('ARM官方库滤波后的波形');
复制代码Matlab计算结果如下：
40.21.png (10.4 KB, 下载次数: 0)
12:45 上传
从上面的波形对比来看，matlab和函数arm_biquad_cascade_df1_f32计算的结果基本是一致的。为了更好的说明滤波效果，下面从频域的角度来说明这个问题，Matlab上面运行如下代码：fs=1000;& && && && && & %设置采样频率 1K
N=400;& && && && && &&&%采样点数& && &
n=0:N-1;
t=n/& && && && && && &%时间序列
f=n*fs/N;& && && && && & %频率序列
x = sin(2*pi*50*t) + sin(2*pi*200*t);& && &%50Hz和200Hz正弦波合成
&&
subplot(211);
y=fft(x, N);& && && && && & %对信号x做FFT& &
plot(f,abs(y));
xlabel('频率/Hz');
ylabel('振幅');
title('原始信号FFT');
y3=fft(sampledata, N);& & %经过IIR滤波器后得到的信号做FFT
subplot(212);& && && && && && && && && && &
plot(f,abs(y3));
xlabel('频率/Hz');
ylabel('振幅');
title('IIR滤波后信号FFT');
复制代码Matlab计算结果如下：
40.22.png (9.63 KB, 下载次数: 0)
12:45 上传
上面波形变换前的FFT和变换后FFT可以看出，200Hz的正弦波基本被滤除。
主题帖子积分
金牌会员, 积分 1963, 距离下一级还需 1037 积分
金牌会员, 积分 1963, 距离下一级还需 1037 积分
40.7 总结& & 本章节主要讲解了巴特沃斯低通，高通，带通和带阻滤波器的实现，有兴趣的可以使用同样的方法实现切比雪夫滤波器的设计。
Tel: 3-8056
备案号: 苏ICP备号-2
Powered by}