上一篇写了巴特沃斯滤波器设计的所有理论知识,这篇文章继续~
注:可能会有不足或者理解偏差的地方,路过的高人请不吝赐教。
闲谈少絮,开始!
========================================================
说明:所有流程图都是伪流程图,认真你就输了,嘿嘿。。。
具备了所有理论知识之后,编程变得简单粗暴。上一篇文章已经确定了开发的流程和思路,即:
①根据给定的滤波器的指标通带截止频率wp,阻带截止频率wst, 通带波纹δ1(Rp(dB)), 阻带波纹δ2(As(dB)),和滤波器的类型确定滤波器阶数N和3dB频率wc。
②根据N查表确定归一化滤波器的系数,并带入wc去归一化转化为H(s),求出其分式的系数数组a(N)和b(M);
③根据求得的系数和目标滤波器类型(低通、高通、带通、带阻),带入相应的s=f(z)转化公式(见上一章的思路图),求得数字滤波器的系统函数H(z)的分子分母系数数组。
④将原始信号和H(z)组成常系数线性差分方程,解出滤波后的信号(此为输出结果,可与matlab运行结果进行对比测试)。
代码分为三个模块:计算原型滤波器的参数(阶数和3dB频率)、计算数字滤波器的系统函数H(z)、滤波方法fliter。
一、计算原型滤波器的参数
首先需要输入参数,模仿matlab的实现,将传入参数定为通带截止频率、阻带截止频率、通带最大衰减、阻带最小衰减、抽样频率和滤波器的类型。流程图如下
注:传入的频率为模拟频率,单位为Hz,三种频率的关系如下表。
模拟频率 f (Hz) Ω = 2π f
模拟角频率 Ω (rad/s) => ω = ΩT
数字频率 ω (rad) ω = 2π f T
低通、高通部分求阶数和3dB的流程图如下图所示(频率不满足判断条件返回错误):
带通、带阻部分求阶数和3dB的流程图如下图所示(频率不满足判断条件返回错误):
二、计算数字滤波器的系统函数H(z)
1.总框图
我习惯先说框架(整体的内容),然后分开讲解各部分的具体实现。滤波器设计的总框架流程图如下图所示:
2.各部分的具体实现
下面把总框图中的部分单元(模块)摘出来一一介绍。
(1)首先去归一化的流程图如下:
看上去流程图很复杂的样子,其实代码很简单。。。
由于高通与低通的滤波器设计思路和设计代码近似(仅有变化公式不同),所以将二者代码合并为一个方法,中间用滤波器的类型参数区别。同理,带通和带阻滤波器也做相应的处理。因此在总框图中就出现“低高通处理”和“带通、带阻处理”的模块,下面先介绍低高通处理方法。
(2)高低通处理方法
高低通处理方法的流程图如下图所示:
(3)带通、带阻处理方法
带通阻处理方法的流程图如下图所示:
(4)辅助程序(工具)
在(2)(3)流程图中用到了杨辉三角和多项式相乘的方法,这两个方法有助于多项式的展开时系数的计算,流程图就不贴了,说明如下
①杨辉三角是一个数列,当一个表达式是类似(s+a)(s+a)(s+a)表达式时,其展开结果就是杨辉三角的一行数。当符号为减时,只需要吧a换成(-a)即可。因此,需要传递一下分式的符号作为参数,进行计算。
②当表达式是两个任意的式子时,需要将表达式写成幂连续下降的形式(如a0+a1*s+0*s^2+0*s^3+a4*s^4),当中间的某项没有时,说明其系数为0。
现在已经把数字滤波器的系统函数(其实是系统函数的分子分母的系数)求出来了,但是只有系统函数不能处理信号也没有什么用,也就是说对于滤波还差一步:对信号的处理——fliter方法。
累残了,fliter部分放在下一节吧。。。
=============================================================
OK,具体求滤波器系统函数的代码随笔就此结束,下一节讲fliter方法和与matlab处理结果对比。