设计四种IIR数字滤波器VB核心程序

1、四种 IIR 数字滤波器的核心程序 本文给出双线性变换法设计四种IIR数字滤波器的核心程序。程序用 VB写成。这四种滤波器是1 巴特沃思滤波器 2. 切比雪夫1型滤波器3. 切比雪夫2型滤波器 4. 椭圆滤波器图1、2示出以上几种滤波器的程序框图。 图 1 图 2 使用双线性变换法的 Butterworth 型 IIR 数字滤波器设计程序 形参说明如下 : PbType - 输入整型量,滤波器通带类型 : PbType = 0 : 低通滤波器; PbType = 1 : 高通滤波器; PbType = 2 : 带通滤波器; PbType = 3 : 带阻滤波器. fp1 - 输入双精度量, 低

2、通或高通滤波器的通带边界频率( Hz ); 带通或带阻滤波器的通带低端边 界频率( Hz ). fp2 - 输入双精度量, 带通或带阻滤波器的通带低端边界频率( Hz ). Apass -输入双精度量, 通带衰减( dB ). fs1 - 输入双精度量, 低通或高通滤波器的阻带边界频率( Hz ); 带通或带阻滤波器的阻带高端边 界频率( Hz ). fs2 - 输入双精度量, 带通或带阻滤波器的阻带高端边界频率( Hz ). Astop - 输入双精度量, 阻带衰减( dB ). fsamp - 输入双精度量, 采样频率( Hz ). points - 输入整型量, 幅频特性计算点数. or

3、d - 输入整型量, 滤波器阶数. NumSec( ) - 输出双精度量, 转移函数二阶节的分子多项式系数二维数组. 元素 NumSec( k, i ) 中, k : 二阶节序号; i : 多项式系数, i = 0 相应于常数项. DenSec( ) - 输出双精度量 转移函数二阶节的分母多项式系数二维数组. 元素 DenSec( k, i ) 中, k : 二阶节序号; i : 多项式系数, i = 0 相应于常数项. NumSec_Z( ) - 输出双精度量 系统函数二阶节的分子多项式系数二维数组. 元素 NumSec_Z( k, i ) 中, k : 二阶节序号; i : 多项式系数,

4、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 Int

5、eger, 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 通

6、带、阻带边界频率 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

7、_t) End If If PbType = 1 Then 高通滤波器; wpass = 2# * Pi * fpass / fsamp: wstop = 2# * Pi * fstop / fsamp 通带、阻带边界频率 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) +

8、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)

9、 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) or

10、d = Fix(orde) + 1 omk0 = omika0(omikaP, epass, ord) 调用 Fz_BP 子程序，将带通模拟滤波器的转移函数变量 s 映射为带通数字滤波器的系统函数变量 z Call Fz_BP(fp1, fp2, fsamp, F1(), F2(), ord_t) End If If PbType = 3 Then 带阻滤波器; wp1 = 2# * Pi * fp1 / fsamp: wp2 = 2# * Pi * fp2 / fsamp 通带上下边界频率 ws1 = 2# * Pi * fs1 / fsamp: ws2 = 2# * Pi * fs2 /

11、fsamp 阻带上下边界频率 Ci = BpC(wp1, wp2) omikaP = Abs(Sin(wp2) / (Cos(wp2) - Ci) omikaS1 = Sin(ws1) / (Cos(ws1) - Ci) omikaS2 = Sin(ws2) / (Cos(ws2) - Ci) If Abs(omikaS1) = Abs(omikaS2) Then omikaS = Abs(omikaS1) Else omikaS = Abs(omikaS2) End If epass = epson(Apass): estop = epson(Astop) 根据对幅频特性的技术要求,计算模拟

12、滤波器的阶数 orde = Ne_B(estop, epass, omikaS, omikaP) ord = Fix(orde) + 1 omk0 = omika0(omikaP, epass, ord) 调用 Fz_BS 子程序，将带阻模拟滤波器的转移函数变量 s 映射为带阻数字滤波器的系统函数变量 z Call Fz_BS(fp1, fp2, fsamp, F1(), F2(), ord_t) End If If ord Mod 2 0 Then 滤波器系统函数的阶数为奇数时，级联节的起始序号为 0（序号为 0 的级联节是一阶节，其余为 二阶节） start = 0 Else 滤波器系统函

13、数的阶数为偶数时，级联节的起始序号为 1（级联节都是二阶节，没有一阶节） start = 1 End If NR = ord 2 系统函数由一阶、二阶节级联而成，k 是节序号。此处求各节的分子分母多项式的系数 For k = start To NR If k = 0 Then 求 Butterworth 滤波器转移函数的一阶节分子分母多项式的系数组 Nums(0) = 1#: Nums(1) = 0#: Nums(2) = 0# Dens(0) = 1#: Dens(1) = 1# / omk0: Dens(2) = 0# NumSec(0, 0) = Nums(0): NumSec(0, 1

14、) = Nums(1): NumSec(0, 2) = Nums(2) DenSec(0, 0) = Dens(0): DenSec(0, 1) = Dens(1): DenSec(0, 2) = Dens(2) End If If k 0 Then 求 Butterworth 滤波器转移函数的二阶节分子分母多项式的系数组 Call HS_B(ord, k, omikaP, epass, Nums(), Dens() For j = 0 To 2 NumSec(k, j) = Nums(j): DenSec(k, j) = Dens(j) Next j End If 从模拟滤波器的转移函数求数

15、字滤波器的系统函数 If k = 0 Then Call BSF2(Nums(), Dens(), 1, F1(), F2(), ord_t, Numz(), Denz(), order_z) Else Call BSF2(Nums(), Dens(), 2, F1(), F2(), ord_t, Numz(), Denz(), order_z) End If 求系统函数第 k 节的分子多项式系数组 For j = 0 To order_z NumSec_Z(k, j) = Numz(j) Next j 求系统函数第 k 节的分母多项式系数组 For j = 0 To order_z DenS

16、ec_Z(k, j) = Denz(j) Next j Next k 求滤波器的幅频特性 For i = 0 To points - 1 For k = start To NR angle = i * 2# * Pi / points 在单位圆上,第 i 点的幅值为 1, 相角为 angle. For j = 0 To order Call Find_ratio(ord, k, angle, NumSec_Z(), DenSec_Z(), ratio(k) Call Find_ratio(order, angle, Num(), Den(), ratio(k) If ratio(k) = 0#

17、 Then ratio(k) = 0.000001 End If Next k AR(i) = 1 For k = start To NR AR(i) = AR(i) * ratio(k) Next k Next iEnd Sub 使用双线性变换法的 Chebyshev 1 型 IIR 数字滤波器设计程序 形参说明如下 : PbType - 输入整型量,滤波器通带类型 : PbType = 0 : 低通滤波器; PbType = 1 : 高通滤波器; PbType = 2 : 带通滤波器; PbType = 3 : 带阻滤波器. fp1 - 输入双精度量, 低通或高通滤波器的通带边界频率( H

18、z ); 带通或带阻滤波器的通带低端 边界频率( Hz ). fp2 - 输入双精度量, 带通或带阻滤波器的通带低端边界频率( Hz ). Apass - 输入双精度量, 通带衰减( dB ). fs1 - 输入双精度量, 低通或高通滤波器的阻带边界频率( Hz ); 带通或带阻滤波器的阻带高端 边界频率( Hz ). fs2 - 输入双精度量, 带通或带阻滤波器的阻带高端边界频率( Hz ). Astop - 输入双精度量, 阻带衰减( dB ). fsamp - - 输入双精度量, 采样频率( Hz ). points - 输入整型量, 幅频特性计算点数. ord - 输入整型量, 滤波器

19、阶数. H0 - 输出双精度量,转移函数和系统函数的常数因子. NumSec( ) - 输出双精度量, 转移函数二阶节的分子多项式系数二维数组. 元素 NumSec( k, i ) 中, k : 二阶节序号; i : 多项式系数, i = 0 相应于常数项. DenSec( ) - 输出双精度量 转移函数二阶节的分母多项式系数二维数组. 元素 DenSec( k, i ) 中, k : 二阶节序号; i : 多项式系数, i = 0 相应于常数项. NumSec_Z( ) - 输出双精度量 系统函数二阶节的分子多项式系数二维数组. 元素 NumSec_Z( k, i ) 中, k : 二阶节序

20、号; i : 多项式系数, i = 0 相应于常数项. DenSec_Z( ) - 输出双精度量 系统函数二阶节的分母多项式系数二维数组. 元素 DenSec_Z( k, i ) 中, k : 二阶节序号; i : 多项式系数, i = 0 相应于常数项. AR( ) - 输出双精度量,滤波器的幅频特性数组.Sub Chebyshev1(PbType As Integer, fp1 As Double, fp2 As Double, Apass As Double, fs1 As Double, fs2 As Double, Astop As Double, fsamp As Double,

21、points As Integer, ord As Integer, H0 As Double, NumSec() As Double, DenSec() As Double, NumSec_Z() As Double, DenSec_Z() As Double, AR() As Double) Dim i%, k%, Fg%, ord_t%, order_z0% Dim angle#, A# Dim ratio(0 To 50) As Double If PbType = 0 Then 低通滤波器 wpass = 2# * Pi * fpass / fsamp: wstop = 2# * P

22、i * fstop / fsamp omikaP = Tan(wpass / 2#): omikaS = Tan(wstop / 2#) epass = epson(Apass): estop = epson(Astop) 根据对幅频特性的技术要求,计算模拟滤波器的阶数 orde = Ne_Ch12(omikaP, omikaS, epass, estop) ord = Fix(orde) + 1 A = 1# / epass A = (1# / ord) * Log(A + Sqr(A 2 + 1#) A = sinh(A) omk0 = A * omikaP 调用 Fz_LP 子程序，将低

23、通模拟滤波器的转移函数变量 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 omikaP = 1# / Tan(wpass / 2#): omikaS = 1# / Tan(wstop / 2#) epass = epson(Apass): estop = epson(Astop) 根据对幅频特性的技术要求,计算模拟滤波器的阶数 orde = Ne

24、_Ch12(omikaP, omikaS, epass, estop) ord = Fix(orde) + 1 A = 1# / epass A = (1# / ord) * Log(A + Sqr(A 2 + 1#) A = sinh(A) omk0 = A * omikaP 调用 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 * fp

25、2 / 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): esto

26、p = epson(Astop) 根据对幅频特性的技术要求,计算模拟滤波器的阶数 orde = Ne_Ch12(omikaP, omikaS, epass, estop) ord = Fix(orde) + 1 A = 1# / epass A = (1# / ord) * Log(A + Sqr(A 2 + 1#) A = sinh(A) omk0 = A * omikaP 调用 Fz_BP 子程序，将带通模拟滤波器的转移函数变量 s 映射为带通数字滤波器的系统函数变量 z Call Fz_BP(fp1, fp2, fsamp, F1(), F2(), ord_t) End If If Pb

27、Type = 3 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(Sin(wp2) / (Cos(wp2) - Ci) omikaS1 = Sin(ws1) / (Cos(ws1) - Ci) omikaS2 = Sin(ws2) / (Cos(ws2) - Ci) If Abs(omikaS1) = Abs(omikaS2) Th

28、en omikaS = Abs(omikaS1) Else omikaS = Abs(omikaS2) End If epass = epson(Apass): estop = epson(Astop) 根据对幅频特性的技术要求,计算模拟滤波器的阶数 orde = Ne_Ch12(omikaP, omikaS, epass, estop) ord = Fix(orde) + 1 A = 1# / epass A = (1# / ord) * Log(A + Sqr(A 2 + 1#) A = sinh(A) omk0 = A * omikaP Call Fz_BS(fp1, fp2, fsam

29、p, F1(), F2(), ord_t) End If NR = ord 2 If ord Mod 2 0 Then start = 0 Nums(0) = omk0: Nums(1) = 0#: Nums(2) = 0# Dens(0) = omk0: Dens(1) = 1#: Dens(2) = 0# For i = 0 To 2 NumSec(0, i) = Nums(i): DenSec(0, i) = Dens(i) Next i H0 = 1 Else start = 1 H0 = Sqr(1# / (1# + epass 2) End If 系统函数由一阶、二阶节级联而成，k

30、 是节序号。此处求各节的分子分母多项式的系数 For k = start To NR If k 0 Then Call HS_Ch1(ord, k, omikaP, epass, Nums(), Dens(), NumSec(), DenSec() End If If k = 0 Then Call BSF2(Nums(), Dens(), 1, F1(), F2(), ord_t, Numz(), Denz(), order_z0) order_z = order_z0 Else Call BSF2(Nums(), Dens(), 2, F1(), F2(), ord_t, Numz(),

31、Denz(), order_z) End If For i = 0 To order_z NumSec_Z(k, i) = Numz(i) Next i For i = 0 To order_z DenSec_Z(k, i) = Denz(i) Next i Next k 求滤波器的幅频特性 For i = 0 To points - 1 angle = i * 2# * Pi / points 在单位圆上,第 i 点的幅值为 1, 相角为 angle. For k = start To NR If k = 0 Then Call Find_ratio(order_z0, k, angle,

32、NumSec_Z(), DenSec_Z(), ratio(k) Else Call Find_ratio(order_z, k, angle, NumSec_Z(), DenSec_Z(), ratio(k) End If If ratio(k) = 0# Then ratio(k) = 0.000001 End If Next k AR(i) = 1 For k = start To NR AR(i) = AR(i) * ratio(k) Next k AR(i) = AR(i) * H0 Next iEnd Sub 使用双线性变换法的 Chebyshev 2 型 IIR 数字滤波器设计程

33、序 形参说明如下 : PbType - 输入整型量,滤波器通带类型 : PbType = 0 : 低通滤波器; PbType = 1 : 高通滤波器; PbType = 2 : 带通滤波器; PbType = 3 : 带阻滤波器. fp1 - 输入双精度量, 低通或高通滤波器的通带边界频率( Hz ); 带通或带阻滤波器的通带低端 边界频率( Hz ). fp2 - 输入双精度量, 带通或带阻滤波器的通带低端边界频率( Hz ). Apass - 输入双精度量, 通带衰减( dB ). fs1 - 输入双精度量, 低通或高通滤波器的阻带边界频率( Hz ); 带通或带阻滤波器的阻带高端 边界频

34、率( Hz ). fs2 - 输入双精度量, 带通或带阻滤波器的阻带高端边界频率( Hz ). Astop - 输入双精度量, 阻带衰减( dB ). fsamp - 输入双精度量, 采样频率( Hz ). points - 输入整型量, 幅频特性计算点数. ord - 输入整型量, 滤波器阶数. H0 - 输出双精度量,转移函数和系统函数的常数因子. NumSec( ) - 输出双精度量, 转移函数二阶节的分子多项式系数二维数组. 元素 NumSec( k, i ) 中, k : 二阶节序号; i : 多项式系数, i = 0 相应于常数项. DenSec( ) - 输出双精度量 转移函数二

35、阶节的分母多项式系数二维数组. 元素 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( ) - 输出双精度量,滤波器的幅频特性数组

36、.Sub Chebyshev2(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, H0 As Double, NumSec() As Double, DenSec() As Double, NumSec_Z() As Double, DenSec_Z() As Double, AR() As Double) Dim i

37、%, j%, k%, Fg%, incOrd%, ord_t%, order_z0% Dim angle#, A# Dim ratio(0 To 20) 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_Ch12(omikaP, omikaS, epas

38、s, estop) ord = Fix(orde) + 1 ord = ord + incOrd A = (1# / ord) * Log(estop + Sqr(estop 2 + 1#) omk0 = omikaS / sinh(A) 调用 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 omikaP =

39、 1# / Tan(wpass / 2#): omikaS = 1# / Tan(wstop / 2#) epass = epson(Apass): estop = epson(Astop) orde = Ne_Ch12(omikaP, omikaS, epass, estop) ord = Fix(orde) + 1 ord = ord + incOrd A = (1# / ord) * Log(estop + Sqr(estop 2 + 1#) omk0 = omikaS / sinh(A) 调用 Fz_HP 子程序，将高通模拟滤波器的转移函数变量 s 映射高通数字滤波器的系统函数变量 z

40、 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

41、(ws2) If omikaS1 = omikaS2 Then omikaS = omikaS1 Else omikaS = omikaS2 End If epass = epson(Apass): estop = epson(Astop) orde = Ne_Ch12(omikaP, omikaS, epass, estop) ord = Fix(orde) + 1 ord = ord + incOrd A = (1# / ord) * Log(estop + Sqr(estop 2 + 1#) omk0 = omikaS / sinh(A) 调用 Fz_BP 子程序，将带通模拟滤波器的转移函数变量 s 映射为带通数字滤波器的系统函数变量 z Call Fz_BP(fp1, fp2, fsamp, F1(), F2(), ord_t) End If If PbType = 3 Then wp1 = 2# * Pi * fp1 / fsamp: wp2 = 2# * Pi * fp2 / fsamp ws1