数字信号处理公式变程序（四）——巴特沃斯滤波器（上）之前搞了⼀些数字信号处理算法编程（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 a0a1a2a3a4a5a6a7a8a9
111
21 1.4142136
3122
41 2.6131259 3.4142136 2.6131259
51 3.236068 5.236068 5.236068 3.236068
61 3.86370337.46410169.14162027.4641016 3.8637033
71 4.493959210.097834714.591793914.591793910.0978347 4.4939592
81 5.125830913.137071221.84615125.688355921.84615113.1370712 5.1258309
91 5.758770516.581718731.163437541.986385741.986385731.163437516.5817187 5.7587705
101 6.392453220.431729142.802061164.882396374.233429264.882396342.802061120.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，理论到此结束，下⼀节上代码！
豫公网安备41160202000603
站长QQ:729038198
关于我们
投诉建议