逆向iir滤波器的实现

这是一个新的算法,来自一个文档。不同于最广泛应用的block process,还有iir尾部序列抵消(cancellation)的策略。

首先看一个单极点的iir滤波器,其z函数为:

 (1)

其中|c|<1。然后根据高中学的等比序列之和的公式,展开记得到后面的一长串加法,再一变就成了一系列的乘法了,这就说明,对于这个iir滤波器,可以通过延迟(z^(-n)),衰减(c^n),以及旁路(bypass)操作形成的单元滤波器级联而成。每一级滤波器的延迟时间为n=2^N,而且幅度相应衰减。图1画出了对应的框图:


图 1

由于iir滤波器的冲击响应为无限长度,因此要准确表达iir的精度需要考虑无限多个单元级联项。但是实际上,可以在某个给定的阈值水平上进行iir截断,使得其变成一个FIR。因为每一个单元衰减是指数级的,截断误差可以做到任意小的水平。既然成了FIR滤波器,一切都好说了。FIR 的逆滤波器实现很直接:只要将冲击响应序列反转然后加上一个总的延迟,使得逆滤波器为因果滤波器即可。图1中的IIR结构被截断后,逆滤波器结构为:


图2 的结构对应的逆滤波器的传输函数为式(2)。式(2)提供了一个稳定的算法,运用一系列的延迟和衰减因子实现了iir 的逆滤波器。

文章后面还有对共轭极点的传输函数的逆滤波实现。

这个方式的优点是,滞后量比block processing的方式少,结构稳定,性能优异。不过说实话,我还没能完全体会这个算法的优点和精髓,sigh。。。以后有机会再来研究吧。

巴特沃兹滤波器的vb.net程序。 ' 使用双线性变换法的 Butterworth 型 IIR 数字滤波器设计程序 ' ' 形参说明如下 : ' ' PbType ----------- 输入整型量 ,滤波器通带类型 : ' PbType = 0 : 低通滤波器 ; ' PbType = 1 : 高通滤波器 ; 2 ' PbType = 2 : 带通滤波器 ; ' PbType = 3 : 带阻滤波器 . ' fp1 ----------- 输入双精度量 , 低通或高通滤波器的通带边界频率 ( Hz ); 带通或带阻滤波器的通带低端边 ‘ 界频率 ( Hz ). ' fp2 ----------- 输入双精度量 , 带通或带阻滤波器的通带低端边界频率 ( Hz ). ' Apass ----------- 输入双精度量 , 通带衰减 ( dB ). ' fs1 ----------- 输入双精度量 , 低通或高通滤波器的阻带边界频率 ( Hz ); 带通或带阻滤波器的阻带高端边 ‘ 界频率 ( Hz ). ' fs2 ----------- 输入双精度量 , 带通或带阻滤波器的阻带高端边界频率 ( Hz ). ' Astop ----------- 输入双精度量 , 阻带衰减 ( dB ). ' fsamp ----------- 输入双精度量 , 采样频率 ( Hz ). ' points ----------- 输入整型量 , 幅频特性计算点数 . ' ord ----------- 输入整型量 , 滤波器阶数 . ' NumSec( ) -------- 输出双精度量 , 转移函数二阶节的分子多项式系数二维数组 . ' 元素 NumSec( k, i ) 中 , ' k : 二阶节序号 ; ' i : 多项式系数 , i = 0 相应于常数项 . ' DenSec( ) -------- 输出双精度量 转移函数二阶节的分母多项式系数二维数组 . ' 元素 DenSec( k, i ) 中 , ' k : 二阶节序号 ; ' i : 多项式系数 , i = 0 相应于常数项 . ' NumSec_Z( ) ------ 输出双精度量 系统函数二阶节的分子多项式系数二维数组 . ' 元素 NumSec_Z( k, i ) 中 , ' k : 二阶节序号 ; ' i : 多项式系数 , i = 0 相应于常数项 . ' DenSec_Z( ) ------ 输出双精度量 系统函数二阶节的分母多项式系数二维数组 . ' 元素 DenSec_Z( k, i ) 中, ' k : 二阶节序号 ; ' i : 多项式系数 , i = 0 相应于常数项 . ' AR( ) ------------ 输出双精度量 ,滤波器的幅频特性数组 . ' Sub Butterworth(PbType As Integer, fp1 As Double, fp2 As Double, Apass As Double, fs1 As Double, fs2 As Double, Astop As Double, fsamp As Double, points As Integer, ord As Integer, NumSec() As Double, DenSec() As Double, NumSec_Z() As Double, DenSec_Z() As Double, AR() As Double) Dim i%, j%, k%, ord_t% Dim angle#, emp1#, temp2#, temp3# Dim ratio(0 To 50) As Double '''''''''''''''''''' If PbType = 0 Then ' 低通滤波器 ; wpass = 2# * Pi * fpass / fsamp: wstop = 2# * Pi * fstop / fsamp ' 通带、阻带边界频率 omikaP = Tan(wpass / 2#): omikaS = Tan(wstop / 2#) epass = epson(Apass): estop = epson(Astop) ' 根据对幅频特性的技术要求 ,计算模拟滤波器的阶数 orde = Ne_B(estop, epass, omikaS, omikaP) ord = Fix(orde) + 1 omk0 = omika0(omikaP, epass, ord) ' 调用 Fz_LP 子程序, 将低通模拟滤波器的转移函数变量 s 映射为低通数字滤波器的系统函数变量 z Call Fz_LP(F1(), F2(), ord_t) ' End If '''''''''''''''''''' If PbType = 1 Then ' 高通滤波器 ; wpass = 2# * Pi * fpass / fsamp: wstop = 2# * Pi * fstop / fsamp ' 通带、阻带边界频率 3 omikaP = 1# / Tan(wpass / 2#): omikaS = 1# / Tan(wstop / 2#) epass = epson(Apass): estop = epson(Astop) ' 根据对幅频特性的技术要求 ,计算模拟滤波器的阶数 orde = Ne_B(estop, epass, omikaS, omikaP) ord = Fix(orde) + 1 omk0 = omika0(omikaP, epass, ord) ' 调用 Fz_HP 子程序,将高通模拟滤波器的转移函数变量 s 映射为高通数字滤波器的系统函数变量 z Call Fz_HP(F1(), F2(), ord_t) End If '''''''''''''''''''' If PbType = 2 Then ' 带通滤波器 ; wp1 = 2# * Pi * fp1 / fsamp: wp2 = 2# * Pi * fp2 / fsamp ' 通带上下边界频率 ws1 = 2# * Pi * fs1 / fsamp: ws2 = 2# * Pi * fs2 / fsamp ' 阻带上下边界频率 Ci = BpC(wp1, wp2) omikaP = Abs((Ci - Cos(wp2)) / Sin(wp2)) omikaS1 = Abs((Ci - Cos(ws1)) / Sin(ws1)) omikaS2 = Abs((Ci - Cos(ws2)) / Sin(ws2)) If omikaS1 <= omikaS2 Then omikaS = omikaS1 Else omikaS = omikaS2 End If epass = epson(Apass): estop = epson(Astop) ' 根据对幅频特性的技术要求 ,计算模拟滤波器的阶数 orde = Ne_B(estop, epass, omikaS, omikaP) ord = Fix(orde) + 1 omk0 = omika0(omikaP, epass, ord) ' 调用 Fz_BP 子程序,将带通模拟滤波器的转移函数变量 s 映射为带通数字滤波器的系统函数变量 z Call Fz_BP(fp1, fp2, fsamp, F1(), F2(), ord_t) End If
评论 2
成就一亿技术人!
拼手气红包6.0元
还能输入1000个字符
 
红包 添加红包
表情包 插入表情
 条评论被折叠 查看
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值