数字信号处理公式变程序（四）——巴特沃斯滤波器（中）

上一篇写了巴特沃斯滤波器设计的所有理论知识，这篇文章继续~

 

注：可能会有不足或者理解偏差的地方，路过的高人请不吝赐教。

 

闲谈少絮，开始！

 

========================================================

说明：所有流程图都是伪流程图，认真你就输了，嘿嘿。。。

具备了所有理论知识之后，编程变得简单粗暴。上一篇文章已经确定了开发的流程和思路，即：

①根据给定的滤波器的指标通带截止频率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处理结果对比。