STM32实现IIR滤波器，可用matlab生成的头文件
放假实在无聊，即将到来的高三非常恐怖，先偷闲一把。matlab的fdatool是好东西，不过很多人不知道该怎么使用它生成的C头文件。
趁着放假有时间，摸索了几天，终于搞定。希望阿莫给条裤子。
该程序已经用于心电采集实验
导联aVF，带宽1-25Hz
http://cache.amobbs.com/bbs_upload782111/files_31/ourdev_569832.JPG
实验过程中图片 (原文件名:DSCF6003.JPG)
http://cache.amobbs.com/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.8641357422, 0, 0
},
{
1, -2, 1
},
{
0.9949035645, 0, 0
},
{
1, -1.999938965, 1
},
{
0.9985351563, 0, 0
},
{
1, -1.99987793, 1
},
{
0.9996337891, 0, 0
},
{
1, -1.99987793, 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.938049316, 0.9401855469
},
{
1, 0, 0
},
{
1, -1.989501953, 0.9900512695
},
{
1, 0, 0
},
{
1, -1.996887207, 0.9971923828
},
{
1, 0, 0
},
{
1, -1.999084473, 0.9993286133
},
{
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.8641357422, 0, 0
},
{
1, -2, 1
},
{
0.9949035645, 0, 0
},
{
1,-1.999938965, 1
},
{
0.9985351563, 0, 0
},
{
1, -1.99987793, 1
},
{
0.9996337891, 0, 0
},
{
1, -1.99987793, 1
},
{
1, 0, 0
}
};
const float IIR_A = {
{
1, 0, 0
},
{
1,-1.938049316,0.9401855469
},
{
1, 0, 0
},
{
1,-1.989501953,0.9900512695
},
{
1, 0, 0
},
{
1,-1.996887207,0.9971923828
},
{
1, 0, 0
},
{
1,-1.999084473,0.9993286133
},
{
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"
static float y;
static float x;
int16 iir_filter(int16 in)
{
uint16 i;
x = in;
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 = y;
}
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);
#endif
使用方法:
首先写好iir_coefs.h，然后调用iir_filter.c对数据流进行滤波
一个伪代码例子：
while(运行中)
{
保存到SD卡(iir_filter(读取ADC采样值()));
}
这个函数比STM32 DSP库中的函数要好很多，DSP库中的2个IIR滤波函数都不能连续处理数据流。
记得在开始滤波之前重置滤波器
iir_reset(); 沙发佩服高中生就会这这些了 但愿学习没有拉下 很不简单啊，楼主高中能搞明白什么是IIR滤波器，是个人才 高中生…… 我不如楼主。 沙发佩服高中生就会这这些了 但愿学习没有拉下 值得学习 高中生！ mark 回复【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采样进行滤波效果非常好，能够有效减弱特定的干扰 楼主是个MIT的料！ 楼主方法用的挺前沿，在掌握点基础就MIT了~ 呵呵,看到你的帖子不错!! mark 楼主真不错！ 你厉害高中生都懂这些了,高三不要读了,直接申请进大学吧. 可以直接读研了。 回复【20楼】showgu
可以直接读研了。
-----------------------------------------------------------------------
同上！ 楼主真是个人才，直接读研吧 人才啊，想办法出国吧 现在的高中生真实不得了啊，真实很羡慕。这么年轻就发现了自己的兴趣方向。 不错 建议置酷贴 记得我高二研究的是linux，
matlib只看了看语法，
高三考完了才用matlab，siclab,g什么的弄了弄fft，架了个数据库分析股票， MARK 条件所致,走的早和走的远没什么必然联系 mark 茶具啊~ 顶 mark 试了下，楼主的程序有问题，
编译通不过的
我上传个经过修改整理的
把滤波器常数直接放到了iir.h
iirourdev_600415YQIMLR.rar(文件大小:2K) (原文件名:iir.rar)
http://cache.amobbs.com/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 = y;
}
y是数组吧
x,IIR_B也是数组.
这里面全是数组相乘相加，我用arduino(gcc-avr编译器)编译报错，
根据数据结构这里应该是
x,IIR_B，y
我只是把系数全部复制到iir.h中，不再include原来导出的.h
另外刚才在atmega168上试用这个处理个29阶的，可能是溢出了，芯片资源要求还是比较多的
另外，这个库是你自己写的还是其他地方找的？找了很久都没找到c的算法。你有没有fir的库？ mark mark 大囧。。。严重失误。。。我复制错了文件
这个是正确代码，也包含了直接型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"
static float y;
static float x;
//IIR_NSEC阶直接型II IIR滤波器
//IIR_NSEC个二阶biquad串联
int16 iir_filter(int16 in)
{
uint16 i;
x = in;
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 = y;
}
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
这个库完全是由我自己完成的 http://cache.amobbs.com/bbs_upload782111/files_35/ourdev_600575I82CLW.jpg
新的效果图 (原文件名:DSCF6773.jpg) 晕～
你改了之后多出来的if是干什么的？好像意思和原来一样的，提高代码效率？
const float IIR_B = {
{
0.02256447077, 0, 0
},
{
1, 2, 1
},
{
0.01625524461, 0, 0
},
{
1, 2, 1
},
{
0.03617158905, 0, 0
},
{
1, 2, 1
},
{
0.02432162315, 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 count;
int64 sum=0;
for(count=0;count<taps;count++)
{
sum += input*coefs;
//其实是一样的,在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_div;
if(sum>32767) sum=32767;
if(sum<-32768) sum=-32768;
return sum;
}
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好像比较复杂 Mark 我没找到导出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)什么是窗函数，有什么特点.
以上不是考试题目，考题比这个要难。但如果是真正把理论和实践同时掌握好的工程师，一定不难回答。 mark armk 还真的不太清楚，我尝试着回答一下吧：
(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 有限长离散信号(视为周期信号) 周期性离散频谱
注：个人理解. mark mark 回复【62楼】mitchell
-----------------------------------------------------------------------
汗，看来没仔细学过离散系统确实不行。。。求大虾联系方式，以后多聊聊，谢谢了。 mark 回复【62楼】mitchell
-----------------------------------------------------------------------
麻省大神求拜师学艺。。。 大神啊，都是。 回复【65楼】warmonkey
-----------------------------------------------------------------------
谦虚了~
不敢当, 我QQ: 799842867, 随时交流
回复【67楼】reloaded 电子浪人
-----------------------------------------------------------------------
您太客气了，mitchell是我中文名字的音译 回复【楼主位】warmonkey
-----------------------------------------------------------------------
“这个函数比STM32 DSP库中的函数要好很多，DSP库中的2个IIR滤波函数都不能连续处理数据流。 ”
可以吧？这个biquad 代码写的相当不错，可以随意剪成1个sector或者增加更多。 神人~mark！ 神人 从楼主的这只手来看不像高中生的... 现在matlab 主推基于模型的设计 m3内核的芯片也开始支持，第三方软件rapidstm32 回复【77楼】tiandewen
-----------------------------------------------------------------------
这个不错，下来看看。 mark mark MARK 佩服，研究一下 mark mark mark...................... mark mark MARK 服得很，不得不顶 回复【楼主位】warmonkey
-----------------------------------------------------------------------
楼主莫怪啊，我想问问你是用什么采集心跳信号的，用的换能器？ 楼主真是个人才 MARK，楼主强人 mark! 回复【89楼】chao8828276
心电信号的实际波行我不太清楚，不过我也搞过fir，我没搞过算法，我也是用matlab产生的系数，但是32阶的fir来处理已经效果很好了，我当时仿真了一下，32阶处理掉50hz以及以上的频率信号是没问题的，当然截止频率要很低
-----------------------------------------------------------------------
我说的是高通滤波 楼主，请教个问题
#define MWSPT_NSEC 9
这个是什么意思
那个数组为什么定义成 MWSPT_NSEC行3列？这个有点没看懂 mark！这个是用的纯数字滤波么？是不是直接把信号放大，然后AD转换？ mark 过来佩服一下楼主 自愧不如 mark，自勉 。。。全部看完了。。。
启示：继续潜水
页:
[1]
2