对于fir滤波器，已经在前面的文章中记录了（https://blog.csdn.net/suiji2442/article/details/112394026POWER-Z仿制DIY&关于MATLAB中滤波器设计工具的使用心得记录），其设计和实现都非常简单。如果在嵌入式系统中可以满足且有必要实时iir运算，那么使用iir滤波器相对fir滤波器可以在使用更小的阶数的情况下实现更好的效果。实验证明，可能20阶的iir效果堪比500阶左右的fir滤波器效果。

首先放出iir的matlab仿真代码：

%本程序为直接2型iir滤波器的实现
%当阶数较高或fc与fs相差多个数量级时，建议使用基本二阶结的级联形式，避免计算出来的b和a系数非常小或者非常大，因运算中有限字长导致运算出错NN=8000;
NNN=20+1;%在这里写iir滤波器的阶数+1t=0:1/NN:1;%采样1秒的信号
%y= 10+10*cos(2*pi*5*t-pi*30/180)+5*cos(2*pi*50*t-pi*30/180);
%y= 10+10*cos(2*pi*50*t-pi*30/180)+10*cos(2*pi*200*t-pi*30/180);
%y= 10+10*cos(2*pi*5*t-pi*30/180)+10*cos(2*pi*50*t-pi*30/180)+2*cos(2*pi*200*t-pi*30/180);
y= 10+2*cos(2*pi*5*t-pi*30/180)+2*cos(2*pi*30*t-pi*30/180)+10*cos(2*pi*100*t-pi*30/180)+100*cos(2*pi*1000*t-pi*30/180);
N=length(t); %样点个数
figure(1);
plot(t,y);
fs=NN;%采样频率
df=fs/(N-1);%分辨率
f=(0:N-1)*df;%其中每点的频率
Y=fft(y(1:N))/N*2;%真实的幅值
figure(2)
plot(f(1:floor(N/2)),abs(Y(1:floor(N/2))));[b,a]=sos2tf(SOS,G);pre_x = zeros(1,NNN);% 前几个时刻的输入值
pre_y = zeros(1,NNN);% 前几个时刻的输出值
res_iir = zeros(1,N);for i=1:Npre_x(1) = y(i);% 差分方程% res_iir(i) = b(1)*pre_x(1)+b(2)*pre_x(2)+b(3)*pre_x(3)+b(4)*pre_x(4)-a(2)pre_y(2)-a(3)pre_y(3)-a(4)pre_y(4);bx = 0;ay = 0;for k=1:NNNbx = bx + b(k) * pre_x(k);endfor k=2:NNNay = ay + a(k) * pre_y(k);  end  res_iir(i) = bx - ay;pre_y(1) = res_iir(i);for j = NNN : -1 : 2  %更新历史数据pre_x(j) = pre_x(j-1);pre_y(j) = pre_y(j-1);end
endfigure(3);
plot(t,res_iir);
fs=NN;%采样频率
df=fs/(N-1);%分辨率
f=(0:N-1)*df;%其中每点的频率
Y2=fft(res_iir(1:N))/N*2;%真实的幅值
figure(4)
plot(f(1:floor(N/2)),abs(Y2(1:floor(N/2))));

上面代码中，需要使用fdatool生成对应阶数的iir系数，并将系数导出到matlab的工作空间，即导出SOS矩阵和G数组。上述代码是20阶IIR滤波器的实现代码。

上面的代码是直接结构的IIR实现，首先放出一个出问题的情况：
其噪声信号为：y= 10+2cos(2pi5t-pi30/180)+2cos(2pi30t-pi30/180)+10cos(2pi100t-pi30/180)+100cos(2pi1000t-pi30/180);
使用fdatool设计20阶iir，采样频率8000Hz，低通截止频率10Hz（与8000Hz差别很大）。
代码执行效果如下：

图1 实验1的原始含噪声信号的时域及频域波形

图2 实验1的20阶iir滤波器滤波之后信号的时域及频域波形

可以从上图看到，滤波结果明显出错。
那么问题在哪呢？问题就是数字滤波器的运算中有限字长效应。可以看到，将fdatool生成的10个基本二阶结的参数转换为直接形式，得到下面的b和a系数



可以看到，其中b系数的值非常非常小，基本都是10的-几十次方，而a系数的值相对正常。那么，会导致什么问题呢？即便是使用双精度浮点进行运算，也会在运算的过程中出现运算字长不足从而使得运算出错的现象，或者iir滤波器的极点出现在单位圆上，出现极限环振荡现象，或者iir滤波器的极点直接出了单位圆，出现不稳定现象，上面的现象便是iir滤波器不稳定，不能正常工作。

那么有什么解决办法呢？办法有三种：
一、降低iir滤波器的阶数
因为在级联的过程中，每一个基本二阶结是级联相乘的，所以当iir的阶数较高时计算出来的直接结构的b、a参数很病态（非常小或者非常大），某些情况下即使使用双精度浮点计算，也是远远不够的，会出现稳定性问题。
但是当阶数较少的时候就没问题了，例如，还是上面的输入信号，将iir的阶数改为2阶iir，即一个基本二阶结，得到如下的效果（注意将matlab代码中的NNN改成2+1）：

图3 实验2的原始含噪声信号的时域及频域波形

图4 实验2的2阶iir滤波器滤波之后信号的时域及频域波形

可以看到，上面的iir输出波形，已经明显将高频的30Hz、100Hz、1000Hz的高频噪声滤除，只剩下直流分量10，和5Hz的交流低频分量，滤波效果还可以接受，试想，这只是一个2阶的iir呀，效果真好了，使用得当足以秒杀fir滤波器的效果。

但是上面的2阶iir效果还不是很好，那么再试一下4阶的iir吧！见下图：

图5 实验3的原始含噪声信号的时域及频域波形

图6 实验3的4阶iir滤波器滤波之后信号的时域及频域波形

看到没，iir滤波器的效果就是这么棒，仅用了一个4阶的iir，上面已经完全把高频分量滤除了，5Hz的原始信号是多么光滑。更加印证了iir滤波器的厉害之处。
但是，小心，万事都有其两面性。不要沉醉在iir滤波器使用较小的阶数便可以获得较好的结果的优势上，想想，我们为什么要对2阶和4阶iir进行实验？
因为高阶的iir直接实现形式不稳定！！！
此时再看，同样的滤波器参数，使用4阶的时候，计算出来的b，a参数：


可以看到，其实b参数已经特别小了，达到了10的-10次方数量级，这样的精度在计算过程中其实也比较难了，也可能会出现稳定性问题，但是相对前面20阶的iir，已经好很多了，因此上面的滤波结果也可以正常运行。

可以看到，降低iir滤波器的阶数，可以帮助稳定iir滤波器。
那么还有没有其它办法呢？
有！

二、使用直接形式的iir滤波器时，fc与fs不要相差多个数量级

上面的实验中，fs是8000Hz，而低通滤波器的截止频率仅为10Hz，相差800倍，太大了，因此使用直接形式的iir计算出来的b和a系数可能会很病态。
下面进行实验：
fs为8000Hz，而iir低通滤波器的截止频率为1000Hz，iir的阶数为20阶，使用直接形式进行计算，看能不能稳定。
此时输入信号改为：y= 10+2cos(2pi10t-pi30/180)+2cos(2pi1500t-pi30/180)+10cos(2pi2000t-pi30/180)+100cos(2pi3000t-pi30/180);

实验结果如下：

图7 实验4的原始含噪声信号的时域及频域波形

图8 实验4的20阶iir滤波器滤波之后信号的时域及频域波形

如上图，可以看到，上面的滤波效果也非常好，将1000Hz以上的高频分量已经滤除。而且即便使用20阶的iir的直接形式，其输出也能基本稳定，原因就是不让fs和fc差别特别大即可。
此时，计算得到的直接形式的iir滤波器的b，a参数如下：


可以看到，只要让fs和fc别差别很大，那么计算得到的b，a参数就不会很病态，此时iir滤波器就能正常运行，当然，前提是使用双精度浮点进行运算是这样的，如果使用单精度浮点进行运算，也不一定能稳定。。。。。。这里顺带再夸一句fir滤波器！！！yyds

上面，也是可以稳定的一种方法，但是，第一中方法没法做高阶的iir，第二种方法又不得不改变iir的截止频率fc和fs之间的关系。。。。显然，在某些情况下，这两种方法都不能被同时接受，那么，还有没有更好的办法呢？

劳动人民（此处指代研究人员=打工人）的智慧是无穷的！！！

三、使用级联结构实现高阶的iir滤波器
一般情况下，使用的都是多个基本二阶结形式进行级联从而实现高阶的iir滤波器，这样的话，便可以让每个基本二阶结的系数不那么病态了，可以避免因为计算字长的限制而导致的iir滤波器不稳定情况。
下面便是基本二阶结结构实现的高阶iir形式：


下面是基本二阶结级联形式的代码（需要SOS和G矩阵，与上面相同）：

%本程序为多个基本二阶结形式的iir滤波器的实现
%当阶数较高或fc与fs相差多个数量级时，建议使用基本二阶结的级联形式，避免计算出来的b和a系数非常小或者非常大，因运算中有限字长导致运算出错NN=8000;
NNN=20+1;%在这里写iir滤波器的阶数+1
Sections = floor(NNN/2);%得到iir滤波器的级联形式基本二阶结的级数t=0:1/NN:1;%采样1秒的信号
%y= 10+10*cos(2*pi*5*t-pi*30/180)+5*cos(2*pi*50*t-pi*30/180);
%y= 10+10*cos(2*pi*50*t-pi*30/180)+10*cos(2*pi*200*t-pi*30/180);
%y= 10+10*cos(2*pi*5*t-pi*30/180)+10*cos(2*pi*50*t-pi*30/180)+2*cos(2*pi*200*t-pi*30/180);
y= 1+0.5*cos(2*pi*5*t-pi*30/180)+2*cos(2*pi*15*t-pi*30/180)+1*cos(2*pi*50*t-pi*30/180)+5*cos(2*pi*2000*t-pi*30/180);
N=length(t); %样点个数
figure(1);
plot(t,y);
fs=NN;%采样频率
df=fs/(N-1);%分辨率
f=(0:N-1)*df;%其中每点的频率
Y=fft(y(1:N))/N*2;%真实的幅值
figure(2)
plot(f(1:floor(N/2)),abs(Y(1:floor(N/2))));%[b,a]=sos2tf(SOS,G);%不需要再计算直接形式的b和a的系数了%首先需要先定义存放数据的空间
data = zeros(Sections,6);%有3个维度。第一个维度存放基本二阶结的信息，第二个维度是当前和历史x,y信息
%data(k,1)表示第k级的x
%data(k,2)表示第k级的x（n-1）
%data(k,3)表示第k级的x（n-2）
%data(k,4)表示第k级的y
%data(k,5)表示第k级的y（n-1）
%data(k,6)表示第k级的y（n-2）res_iir = zeros(1,N);%存放滤波的结果for i=1:Ndata(1,1) = y(i);%将当前点的信息放进来for k=1:Sectionsdata(k,4) =  G(k) * (SOS(k,1)*data(k,1)+SOS(k,2)*data(k,2)+SOS(k,3)*data(k,3)) -SOS(k,5)*data(k,5)-SOS(k,6)*data(k,6);data(k,6) = data(k,5);data(k,5) = data(k,4);data(k,3) = data(k,2);data(k,2) = data(k,1);if(k~=Sections)%如果不是计算到了最后一级，就将当前的输出送给下一级输入data(k+1,1) = data(k,4);elseres_iir(i) = data(k,4);%如果是计算到了最后一级，就将当前得到结果输出endend
endfigure(3);
plot(t,res_iir);
fs=NN;%采样频率
df=fs/(N-1);%分辨率
f=(0:N-1)*df;%其中每点的频率
Y2=fft(res_iir(1:N))/N*2;%真实的幅值
figure(4)
plot(f(1:floor(N/2)),abs(Y2(1:floor(N/2))));

设计一个20阶的iir滤波器，同样，采样频率为8000Hz，低通截止频率为10Hz。信号为：y= 1+0.5cos(2pi5t-pi30/180)+2cos(2pi15t-pi30/180)+1cos(2pi50t-pi30/180)+5cos(2pi2000t-pi30/180);效果如下：

图9 实验5的原始含噪声信号的时域及频域波形

图10 实验5的20阶iir级联滤波器滤波之后信号的时域及频域波形

可以见到上图，使用级联的模式的情况下，便可以使用高阶的iir滤波器对采样的数据进行滤波，即使fs和fc相差非常大也不会带来稳定性问题。

另外，再进行一个小实验，试一下20阶和50阶的iir带阻滤波器（陷波器）效果怎么样，其中信号中含有幅度非常巨大的49、50、51Hz干扰，输入信号如下：y= 1+2cos(2pi5t-pi30/180)+100cos(2pi49t-pi30/180)+100cos(2pi50t)+100cos(2pi51t-pi30/180)+0.5cos(2pi70t-pi30/180)+5cos(2pi2000t-pi*30/180);

实验中，因为iir阶数变大，所以将数据时间长度调大到了5秒的数据长度。
陷波器阻带为45Hz~55Hz。采样频率8000Hz。

下面是实验结果：

图11 实验6的原始含噪声信号的时域及频域波形

图12 实验6的20阶iir级联带阻滤波器滤波之后信号的时域及频域波形

可以看到仅使用20阶iir滤波效果就非常好了。

下面开始50阶iir滤波器的验证：

图13 实验7的原始含噪声信号的时域及频域波形

图14 实验7的50阶iir级联带阻滤波器滤波之后信号的时域及频域波形

可以看到上面20阶和50阶iir滤波器的对比，发现50阶iir滤波器需要更长的时间才能稳定，即其反馈环路更长，然而他们的实际滤波效果并差不多，因此，给我们的启示是：在实际的系统中，应当合理选择iir滤波器的阶数，并不是越高越好。

下面展示一下两个滤波器（20阶和50阶）的滤波后稳定波形的细节放大图：

图15 20阶iir级联带阻滤波器滤波之后信号的时域波形细节图

图16 50阶iir级联带阻滤波器滤波之后信号的时域波形细节图

对比图15和图16，发现无论是20阶还是50阶iir滤波器稳定之后的滤波效果都非常好，输出的波形里面只包含了幅度为1的直流分量、幅度为2的5Hz低频分量和幅度为0.5的70Hz高频分量。两个图的波形基本一致。然而50阶的iir滤波器需要更长时间稳定。

至此，无论是fir滤波器还是iir滤波器，我们均从其具体实现方法方面进行了一些实验，发现、解决了一些问题，对于里面的代码，可以非常方便地修改为C语言代码应用于嵌入式开发中。

这两篇小文档作为自己进行这些实验的一个小结和心得吧，可以留给以后的自己去看。

但是，这两种滤波器仅当噪声的频率与实际信号频率之间不重叠的时候有用。对于一些随机时间信号噪声，这种经典滤波器是无能为力的，然后就需要一些现代统计最优滤波器去对这些噪声进行滤波，例如维纳滤波器、卡尔曼滤波器、自适应滤波等。这些统计最优滤波器的相关实验等待后面更新吧~
最后，再次感谢本学期“信号处理与数据分析”课程的邱老师，这门课程真的很有用，特别是应用于实践之后。

使用matlab设计iir滤波器并自行编写代码实现iir滤波器（可对应于C语言应用在嵌入式系统中）相关推荐

  1. MATLAB离散卷积的实现_自己编写代码_实现两列数的卷积

    一.实验目的 1.了解熟悉MATLAB中conv.filter函数的使用: 2.回忆卷积,深入了解conv原理,从差分方程和矩阵两个角度重构conv函数. 二.实验原理 卷积是两个变量在某范围内相乘后 ...

  2. jquery 图像滑块_jQuery CSS图像滑块–自行编写代码

    jquery 图像滑块 Today we will learn how to make our own jQuery image slider. This Article is a part of & ...

  3. 嵌入式系统图形用户界面（GUI）的设计与研究

    1 引 言 在工业控制领域里,各种仪器仪表.智能工控设备也广泛采用了嵌入式技术,但由于资源有限, 这些系统一般不希望建立在庞大累赘的.非常消耗系统资源的操作系统和GUI之上,比如Windows或X W ...

  4. 嵌入式系统开发设计---嵌入式系统开发设计

    嵌入式系统设计的主要任务是定义系统的功能.决定系统的架构,并将功能映射到系统实现架构上.这里,系统架构既包括软件系统架构也包括硬件系统架构.一种架构可以映射到各种不同的物理实现,每种实现表示不同的取舍 ...

  5. 编写书籍《C语言嵌入式系统编程修炼之道》序言

    序言        目前,嵌入式系统已经无处不在,遍布于世界的每一个角落.智能家电.手机.PDA.汽车.通信电台等几乎所有的电力.电器与电子产品都包含一个或多个嵌入式系统.有人的地方就有江湖,有电的地 ...

  6. 绝对好文：嵌入式系统的软件架构设计！

    要学嵌入式,关注@我要学嵌入式,嵌入式猛男的加油站. 1. 前言 嵌入式是软件设计领域的一个分支,它自身的诸多特点决定了系统架构师的选择,同时它的一些问题又具有相当的通用性,可以推广到其他的领域. 提 ...

  7. 在linux下进行嵌入式系统设计,一种应用于测控系统的基于Linux的嵌入式系统的设计...

    描述 1.前言 随着网络控制技术的快速发展,工业以太网得到逐步完善,在工业控制领域获得越来越广泛的应用.工业以太网使用了TCP/IP协议,便于联网,并具有高速控制网络的优点.随着32位嵌入式CPU价格 ...

  8. 嵌入式系统开发这六点硬件设计需要细心留意

    嵌入式设计是个庞大的工程,今天就说说硬件电路设计方面的几个注意事项,首先,咱们了解下嵌入式的硬件构架.我们知道,CPU是整个系统的灵魂,所有的外围配置都与其相关联,这也突出了嵌入式设计的一个特点硬件可 ...

  9. 嵌入式系统课程设计题目

    简介:一些嵌入式系统课程设计题目,可以当做对你学习ARM的一个检测. 嵌入式系统课程设计-选题要求及课题 1.嵌入式系统课程设计时长两星期,要求学生分组进行课程设计,每组学生人数为2-3人(可在不超过 ...

  10. c语言嵌入式系统编程软件,C语言嵌入式系统编程软件设计研究论文

    C语言嵌入式系统编程软件设计研究论文 摘要:近年来,C语言编程在嵌入式系统越来越受到广大技术人员的青睐.介绍了C语言系统软件的编程思路,阐述了嵌入式系统编程软件架构的基本知识,包括模块划分.分层架构. ...

最新文章

  1. Swift - 委托（delegate）的介绍，及使用样例
  2. 招聘|腾讯机器人实验室语义视觉方向（实习+社招）
  3. 数据结构 互换二叉树中所有结点的左右子树 C
  4. 内核网络中的GRO、RFS、RPS技术介绍和调优
  5. supersr--时间显示逻辑--NSDate+NSCalendar
  6. Scala中集合类型与java中集合类型转换
  7. JavaScript的正则表达式实现邮箱校验
  8. c语言建立队列（顺序队列、循化队列和链式队列）
  9. java 动态代理深度学习(Proxy,InvocationHandler)
  10. 前端学习（1802）：前端调试之事件伪类练习
  11. Magento教程 22：如何确认订单报表？
  12. (7)Microsoft office Word 2013版本操作入门_常用技巧
  13. 分子排列不同会导致_东华大学《高分子物理》各章选择判断题
  14. 关联分析研究思路及应用：GWASTWAS
  15. Android的jsoup方法,在Android中使用Jsoup
  16. 游戏筑基开发之利用文件函数读出文件数据及处理（反序列化）（C语言）
  17. 快递小哥逆袭自传：用了6年时间做到了IT部门主管
  18. 【深度学习】模型平均误差分析
  19. ArcMap导入数据到ArcSDE报000597或者000224的错误
  20. iOS Quartz2D 渐变图形 CGGradient CGShading

热门文章

  1. SCI论文写作 -- 搜索工具汇总
  2. 月薪30K+的电子工程师应具备什么？
  3. openwrt - transmission
  4. 保护您的眼睛：电脑背景色设置（XP WIN 7）
  5. 数字逻辑：多级门电路
  6. 《RFID技术与应用》试题库（含答案）
  7. （转）技嘉 MA790FXT-UD5P搭配AMD X4 965超频解析
  8. Linux下安装hbase
  9. mw150um 驱动程序win10_mercury无线网卡驱动
  10. 赶紧收藏！不可多得的Instagram运营技巧