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

之前搞了一些数字信号处理算法编程（OC），一直没来得及整理，现在整理一下，包括FFT、巴特沃斯滤波器（高通、低通、带通、带阻）、数据差值（线性、sinc、三次样条*）、数据压缩（等距、平均、峰值检测）和模仿matlab的STFT功能（spectrogram函数三维绘图）。

 

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

 

今天来说一下数字滤波器的代码实现（IIR）。

---------------------------------------------------------------------------------------

一、基本概念

1.数字滤波器（DF）和模拟滤波器（AF）一样，都是用来滤波的，它将信号的某些频率（频段）的信号加以放大，而将另外一些频率（频段）的信号加以抑制。借助A/D、D/A转换器，数字滤波器可以处理模拟信号也可以输出模拟信号。

2.数字滤波器有4种表示方法

①线性差分方程：也就是它的滤波作用的基本构成事数值运算的部件：相加器、相乘器、延时器。而模拟滤波器的基本部件则是电感器、电容器、电阻器及有源器件。

②系统函数。将a0归一化为a0=1，即用z^(-1)或z的有理分式表示系统函数。则有

（IIR数字滤波器）

③单位抽样响应h(n)拉氏变换得到

④线性信号流图。系统函数如下

其直接II型结构图如下

3.分类

（1）按冲激响应分为无限长冲激响应（IIR）数字滤波器和有限长冲激响应（FIR）数字滤波器。

IIR：从2②中的IIR系统函数可以看出，它是一个z^(-1)的有理分式，包括有输出到输入的反馈网络结构，此系统分母多项式A(z)决定了反馈网络，同时确定了有限z平面的极点，而分子多项式B(z)决定了正馈网络，同时确定了有限z平面的零点。

FIR：FIR滤波器系统函数可表示为z^(-1)的多项式，即IIR系统函数中ai=0（分母为1）,可以写成如下形式

其中h(n)时系统的单位冲激响应，显然h(i)=bi（当ai=0时），H(z)在有限z平面只有零点，如果是因果系统则全部极点都在z=0处。系统不存在反馈网络。

（2）按滤波器幅度响应可分为低通、高通、带通、带阻、全通等。按照奈奎斯特抽样定理，信号最高频率fh只能限于fh≤fs/2（fs为臭氧频率），即wh≤π。与模拟滤波器不同，数字滤波器频率响应是以2π为周期的周期函数，这些理想的幅度响应特性（如下图所示），即“突变”型幅度响应会造成无限长的非因果的单位冲激响应，是不可实现的，只能用可实现的实际滤波来逼近它。

（3）按相位响应分类：线性相位、非线性相位。如果要求严格的线性相位，则必须使用FIR线性相位滤波器。

（4）按特殊要求分类：最小相位滞后滤波器、梳状滤波器、陷波器、全通滤波器、谐振器，甚至包括波形产生器等。采用零极点的适当配置方法，可以得到这些滤波器。

二、IIR数字滤波器的技术指标

以数字低通滤波器为例（见下图），指标包括：通带截止频率wp，阻带截止频率wst， 通带波纹δ1（Rp(dB)）， 阻带波纹δ2（As(dB)） 。

通带允许的最大衰减Rp(dB)以及阻带应达到的最小衰减As(dB)范围如下（已归一化为H(e^(jw))=1）：

 

注：

①高通与低通的指标一致；

②带通的技术指标通带截止频率分通带上截止频率wp2和通带下截止频率wp1，阻带截止频率分上阻带截止频率wst2和下阻带截止频率wst1；

③带阻的技术指标通带截止频率分为上通带截止频率wp2和下通带截止频率wp1，阻带截止频率分阻带上截止频率wst2和阻带下截止频率wst1。

 

三、滤波器的设计步骤

滤波器的设计就是要找到一个满足技术指标要求的可实现的因果稳定的数字滤波器来逼近理想的滤波器幅度特性。

1.滤波器设计思路

本文采用间接法设计数字滤波器（先设计模拟低通滤波器在通过双线性变换法得到数字低通、高通、带通、带阻滤波器），设计模型如下图所示。参数说明：fp——通带截止频率，fs——阻带截止频率，rp通带最大衰减，rs阻带最小衰减，N——滤波器阶数，fc——3dB截止频率，sa——原型滤波器分母多项式的系数，sb原型滤波器分子多项式的系数，za——数字滤波器系统函数分母多项式的系数，zb——数字滤波器系统函数分子多项式的系数，filter——滤波方法（将信号代入，返回滤波后的信号）。

 

即，数字滤波器的设计及滤波思路图如下图：

2.模拟滤波器设计

（1）模拟滤波器的设计步骤

①给定模拟滤波器技术指标Ωp，Ωst，(1-δ1)（或Rp dB），δ2（或As dB）；

②计算滤波器所需阶数N；

注：本文设计的滤波器阶数N不能大于10，因为阶数越高计算出的值越不准确。

③查表确定归一化低通滤波器系统函数Has(s)；

④将Has(s)转换为所需类型的低通滤波器系统函数Ha(s)。

（2）具体实现

①确定滤波器的阶数N，公式如下所示，计算得到的结果向上取整。注意此处的Ωst和Ωp为预畸后的值（预畸公式见上思路设计图中公式）

②确定3dB频率，公式如下，当选Ωc=(Ωcp+Ωcs)/2时通带阻带衰减皆可超过要求（查看matlab源码发现采用Ω=Ωcs，因此本文采取Ω=Ωcs）

③根据N的值，查表得归一化低通滤波器系统函数Has(s)的分母多项式系数，公式如下：

d0一般由Ω=0时|Han(j0)|=1（增益为1）来确定，由于a0=1，一次d0=a0=1。

附归一化巴特沃斯滤波器分母多项式的系数表如下：

 

N	a0	a1	a2	a3	a4	a5	a6	a7	a8	a9
1	1	1
2	1	1.4142136
3	1	2	2
4	1	2.6131259	3.4142136	2.6131259
5	1	3.236068	5.236068	5.236068	3.236068
6	1	3.8637033	7.4641016	9.1416202	7.4641016	3.8637033
7	1	4.4939592	10.0978347	14.5917939	14.5917939	10.0978347	4.4939592
8	1	5.1258309	13.1370712	21.846151	25.6883559	21.846151	13.1370712	5.1258309
9	1	5.7587705	16.5817187	31.1634375	41.9863857	41.9863857	31.1634375	16.5817187	5.7587705
10	1	6.3924532	20.4317291	42.8020611	64.8823963	74.2334292	64.8823963	42.8020611	20.4317291	6.3924532

 

④转换成所需类型的低通滤波器系统函数Ha(s)，去归一化另s=s/Ωc。

 

3.模拟滤波器的数字化知识

模拟滤波器的数字化方法有很多种，包括冲激响应不变法（脉冲响应不变法）、阶跃响应不变法和双线性变换法等。本文采用双线性变换法实现。

（1）基本思路

双线性变换法是使数字滤波器的频率响应与模拟滤波器的频率响应相似的一种变换，它使得Ω和ω之间是单值映射关系可以避免频率响应的混叠失真。

下图表示了变换思路，即将s平面整个变换到一个中介s1平面的一个窄带Ω：-π/T→π/T之中，然后经过z=e^(s1*T)的变换，将s1平面映射到z平面。从s1到z平面的变换是单值变换，从而使整个变换过程（s到z）成为单值的变换。

（2）变换公式

低通、高通、带通、带阻的变换关系式都在本节第一模块的“数字滤波器的设计及滤波思路图”中展现了。

4.滤波方法

如前所述，滤波过程就是解常系数线性差分方程的过程，形式如：

其中，x(n)序列为滤波前的信号序列，ak, bm为H(z)系统函数分母与分子的系统数组，求出的y(n)即为滤波后的信号序列。

注：x(n)与y(n)的长度要相等，且a0=1。

公式的化简工程如下：

默认条件，当k<0时，x(k), y(k)都为0。例如n=0时，y(0)=b0*x0+b1*x(-1)+...+bM*x(-M)-a1*y(-1)-...-aN*y(-N)=b0*x(0)。经过迭代，可以求出y(n)序列的所有值。

 

 

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

OK，理论到此结束，下一节上代码！