数字滤波器设计之一:巴特沃斯(Butterworth)滤波器
我们常说的经典滤波器是根据傅里叶分析和变换设计出来的,只允许一定频率范围内的信号成分正常通过,而阻止另一部分频率成分通过。按照最佳逼近特性或者滤波通带特性分类,主要为巴特沃斯滤波器(Butterworth)、切比雪夫滤波器(Chebyshev)、贝塞尔滤波器(Bessel)和椭圆滤波器(Elliptic)四种。每种MATLAB都有相应的函数,用起来也比较方便,但是却缺少C/C++的程序,于是自己仔细研究了每种滤波器的特性和原理,并且部分滤波器实现了C语言的代码化,接下来的时间会对这些滤波器的原理和C语言的实现进行介绍。
该系列均以低通滤波器为原型来介绍,其他类型的滤波器可以由低通滤波器通过频率变换转换得到,这里不过多介绍。低通滤波器的主要性能指标有以下几个:通带截止频率fp、阻带截止频率fs、通带衰减 ( Ap)、阻带衰减 ( As) 以及归一化频率时需要用到的-3dB的转折频率fc。
1. Butterworth滤波器原理
Butterworth滤波器因其在通带内的幅值特性具有最大平坦的特性而闻名,是四种经典滤波器中最简单的,巴特沃斯滤波器只需要两个参数表征,滤波器的阶数N和-3dB处的截止频率。其幅度平方函数为:
N是滤波器的阶数,从幅度平方函数可以看出,N阶滤波器有2N个极点,而且这2N个极点均布在一个圆上,圆的半径为,称之为Butterworth圆,Butterworth滤波器系统是一个线性系统,要使其稳定,其极点必须位于S平面的左半平面,所以取左半平面内的N个极点作为滤波器的极点,滤波器就是稳定的了,求出极点之后,计算模拟滤波器的系数as、bs,然后通过双线性变换(不懂得自行查书)由模拟域到数字域,求出系数az和bz 。最后通过差分方程就可以计算滤波结果了。
2. C语言实现
A.求阶数
公式为:
代码:
N = ceil(0.5*( log10 (( pow (10, Stopband_attenuation/10) - 1)/( pow (10, Passband_attenuation/10) - 1)) / log10 (Stopband/Passband) ));
B.求极点
公式:
代码:
for(k = 0;k <= ((2*N)-1) ; k++){if(Cutoff*cos((2*k+N-1)*(pi/(2*N))) < 0.0){ poles[count].x = -Cutoff*cos((2*k+N-1)*(pi/(2*N)));poles[count].y = -Cutoff*sin((2*k+N-1)*(pi/(2*N))); count++;if (count == N) break;}}
C.计算模拟滤波器系数
公式:
代码(这部分需要自己推导,多算几步就找到规律了):
Res[0].x = poles[0].x; Res[0].y = poles[0].y;Res[1].x = 1; Res[1].y= 0;for(count_1 = 0;count_1 < N-1;count_1++)//N个极点相乘次数{for(count = 0;count <= count_1 + 2;count++){if(0 == count){Res_Save[count] = ComplexMul( Res[count], poles[count_1+1] ); }else if((count_1 + 2) == count){Res_Save[count].x += Res[count - 1].x;Res_Save[count].y += Res[count - 1].y; } else {Res_Save[count] = ComplexMul( Res[count], poles[count_1+1] ); Res_Save[count].x += Res[count - 1].x;Res_Save[count].y += Res[count - 1].y; }}for(count = 0;count <= N;count++)//Res[i]=a(i),i越大次数越高{Res[count].x = Res_Save[count].x; Res[count].y = Res_Save[count].y;*(a + N - count) = Res[count].x;} }*(b+N) = *(a+N);
D.双线性变换
用下式进行替换 中的
变量,得到
,然后类似上面的计算方法计算Z域的系数az和bz,其中T为采样周期,但是因为在计算中会被约去,所以简化计算,这里取1。
公式:
代码:
int Count = 0,Count_1 = 0,Count_2 = 0,Count_Z = 0;double *Res, *Res_Save;Res = new double[N+1]();Res_Save = new double[N+1](); memset(Res, 0, sizeof(double)*(N+1));memset(Res_Save, 0, sizeof(double)*(N+1));for(Count_Z = 0;Count_Z <= N;Count_Z++){*(az+Count_Z) = 0;*(bz+Count_Z) = 0;}for(Count = 0;Count<=N;Count++){ for(Count_Z = 0;Count_Z <= N;Count_Z++){Res[Count_Z] = 0;Res_Save[Count_Z] = 0; }Res_Save [0] = 1;for(Count_1 = 0; Count_1 < N-Count;Count_1++)//计算(1-Z^-1)^N-Count的系数,{ //Res_Save[]=Z^-1多项式的系数,从常数项开始for(Count_2 = 0; Count_2 <= Count_1+1;Count_2++){if(Count_2 == 0) {Res[Count_2] += Res_Save[Count_2];} else if((Count_2 == (Count_1+1))&&(Count_1 != 0)) {Res[Count_2] += -Res_Save[Count_2 - 1]; } else {Res[Count_2] += Res_Save[Count_2] - Res_Save[Count_2 - 1];} }for(Count_Z = 0;Count_Z<= N;Count_Z++){Res_Save[Count_Z] = Res[Count_Z] ;Res[Count_Z] = 0; } }for(Count_1 = (N-Count); Count_1 < N;Count_1++)//计算(1-Z^-1)^N-Count*(1+Z^-1)^Count的系数,{ //Res_Save[]=Z^-1多项式的系数,从常数项开始for(Count_2 = 0; Count_2 <= Count_1+1;Count_2++){if(Count_2 == 0) {Res[Count_2] += Res_Save[Count_2];} else if((Count_2 == (Count_1+1))&&(Count_1 != 0)) {Res[Count_2] += Res_Save[Count_2 - 1];} else {Res[Count_2] += Res_Save[Count_2] + Res_Save[Count_2 - 1];} }for(Count_Z = 0;Count_Z<= N;Count_Z++){Res_Save[Count_Z] = Res[Count_Z] ;Res[Count_Z] = 0; }}for(Count_Z = 0;Count_Z<= N;Count_Z++){*(az+Count_Z) += pow(2,N-Count) * (*(as+Count)) * Res_Save[Count_Z];*(bz+Count_Z) += (*(bs+Count)) * Res_Save[Count_Z]; } }//最外层for循环for(Count_Z = N;Count_Z >= 0;Count_Z--){*(bz+Count_Z) = (*(bz+Count_Z))/(*(az+0));*(az+Count_Z) = (*(az+Count_Z))/(*(az+0));}
E.由由差分方程计算滤波结果
采用滤波器直接II型结构,可以减少一半的中间缓存内存,具体代码如下:
double FiltButter(double *pdAz, //滤波器参数表1double *pdBz, //滤波器参数表2int nABLen, //参数序列的长度double dDataIn,//输入数据double *pdBuf) //数据缓冲区
{int i;int nALen;int nBLen;int nBufLen;double dOut;if( nABLen<1 )return 0.0;//根据参数,自动求取序列有效长度nALen = nABLen;for( i=nABLen-1; i; --i ){if( *(pdAz+i) != 0.0 )//从最后一个系数判断是否为0{nALen = i+1; break;}}//printf("%lf ", nALen);if( i==0 ) nALen = 0;nBLen = nABLen;for( i=nABLen-1; i; --i ){if( *(pdBz+i) != 0.0 ){nBLen = i+1;break;}}//printf("%lf ", nBLen);if( i==0 ) nBLen = 0;//计算缓冲区有效长度nBufLen = nALen;if( nALen < nBLen)nBufLen = nBLen;//滤波: 与系数a卷乘dOut = ( *pdAz ) * dDataIn; // a(0) * x(i) for( i=1; i<nALen; i++) // a(i) * w(n-i),i=1toN{dOut -= *(pdAz+i) * *(pdBuf + (nBufLen - 1) - i);} //卷乘结果保存为缓冲序列的最后一个*(pdBuf + nBufLen - 1) = dOut;//滤波: 与系数b卷乘dOut = 0.0;for( i=0; i<nBLen; i++) // b(i) * w(n-i){ dOut += *(pdBz+i) * *(pdBuf + (nBufLen - 1) - i);}//丢弃缓冲序列中最早的一个数, 最后一个数清零for( i=0; i<nBufLen-1; i++){*(pdBuf + i) = *(pdBuf + i + 1);}*(pdBuf + nBufLen - 1) = 0;//返回输出值return dOut;
}
VC6.0开发环境滤波结果如下:
至此就完成了Butterworth滤波器的设计过程,完整代码请自行下载。
参考资料:http://blog.csdn.net/zhoufan900428/article/details/9069475
源码地址:http://download.csdn.net/detail/zhwzhaowei/9830114
数字滤波器设计之一:巴特沃斯(Butterworth)滤波器相关推荐
- java巴特沃斯滤波器编程_巴特沃斯(Butterworth)滤波器 (1)
下面深入浅出讲一下Butterworth原理及其代码编写. 1. 首先考虑一个归一化的低通滤波器(截止频率是1),其幅度公式如下: 当n->∞时,得到一个理想的低通滤波反馈: ω<1时,增 ...
- 巴特沃斯(Butterworth)滤波器 (2) - 双线性变换
这里接着上篇讲一下双线性变换Bilinear Transformation,它实现了模拟信号(连续域)与数字信号(离散域)之间的转换. 双线性变换公式如下: 反推可得到: 因此可以根据连续域传递函数推 ...
- 使用matlab设计IIR巴特沃斯低通滤波器
目的和要求 设计IIR巴特沃斯低通滤波器 将滤波器用于加噪声信号的处理 改变参数指标比较不同 实验结果与分析 (1)滤波器频率特性如图 输入信号特性如图: 经过滤波器之后信号特性如图: (2)改变指标 ...
- 【老生谈算法】matlab实现巴特沃斯IIR滤波器程序设计源码
matlab巴特沃斯IIR滤波器程序设计 1.文档下载: 本算法已经整理成文档如下,有需要的朋友可以点击进行下载 序号 文档(点击下载) 本项目文档 [老生谈算法]matlab巴特沃斯IIR滤波器程序 ...
- IIR滤波器设计代码(巴特沃斯+脉冲响应不变法/双线性变换法) Matlab代码
引言 说实话我感觉自己滤波器不算学到位了,一般来说我是需要把整个过程都弄得非常清楚,但是这个模拟滤波器设计是真的麻烦,至少我现在不确定以后会从事DSP相关的内容,就没有对细节考量,但或许也没谁会对它去 ...
- python巴特沃斯滤波器_butterworth-实时数字巴特沃斯IIR滤波器-Kei Imada
作者:Kei Imada ### 作者邮箱:kimada1@swarthmore.edu ### 首页:https://github.com/keikun555/butter ### 文档:None ...
- 常用模拟低通滤波器的设计~经典 IIR 滤波器之巴特沃斯滤波器
目录 常用模拟低通滤波器的设计--巴特沃斯(Butterworth)滤波器 1.确定系统函数的极点 2.巴特沃斯(Butterworth)滤波器 2.1.buttap 函数 2.2.buttord 函 ...
- 巴特沃斯低通滤波器归一化matlab,基于MATLAB设计巴特沃斯低通滤波器
<基于MATLAB设计巴特沃斯低通滤波器>由会员分享,可在线阅读,更多相关<基于MATLAB设计巴特沃斯低通滤波器(6页珍藏版)>请在人人文库网上搜索. 1.基于MATLAB设 ...
- [DSP] Butterworth (巴特沃斯)数字滤波器设计参考
Butterworth (巴特沃斯)滤波器设计参考 -- By Water 在嵌入式音频产品开发过程中经常会到LPF(Low Pass Filter 低通滤波器)和HPF(High Pass Filt ...
- 算法学习 - 模拟滤波器(巴特沃斯、贝塞尔、切比雪夫、椭圆)及IIR滤波器设计
摘要: 在进行信号处理算法的嵌入式代码实现时,经常需要设计滤波器来对信号进行预处理或者后处理.常用的有IIR滤波器和FIR滤波器.虽然IIR滤波器在稳定性和相位响应上比FIR滤波器较差,但是在相同滤波 ...
最新文章
- Linux上oracle的安装
- [deviceone开发]-doSpace应用源码开源
- 前端编码规范之JavaScript
- ios沙箱模式开启_iOS我眼中的沙盒机制
- Fibonacci (hdu1568)数学公式
- 为什么会出现NoSQL数据库
- PHP实现网页截图?还可以实现登录截图!
- 将一个文件夹下的MP4文件合并为一个
- 苹果硅之后:PC市场会走向何方?
- 来自微信团队的 6 个开源项目
- 牵一发动全身【Nhibernate基本映射】
- 笔记本onenote绘画快捷键_OneNote 超全快捷键
- python自动化测试面试题大全带答案_Python自动化测试笔试面试题精选
- python携程酒店评论_python爬取携程景点评论信息
- LCD1602显示总结
- R语言 编写循环语句
- 【装机心得】关于系统启动U盘的那些事(下)
- Mac OS Catalina (10.15)下编译Redis Desktop Manager(RDM)
- 井柏然自己的字体,手写语录合集
- ieframe.dll修复方法