IIR数字滤波器设计的核心程序

' ord ------------ 输入整型量, 滤波器阶数.

' H01 --------------- 输出双精度量,归一化转移函数的常数因子.

' NumSec( ) ------ 输出双精度量, 转移函数二阶节的分子多项式系数二维数组. ' 元素 NumSec( k, i ) 中, ' k : 二阶节序号;

' i : 多项式系数, i = 0 相应于常数项.

' DenSec( ) ------- 输出双精度量 转移函数二阶节的分母多项式系数二维数组. ' 元素 DenSec( k, i ) 中, ' k : 二阶节序号;

' i : 多项式系数, i = 0 相应于常数项.

' H0２ -------------- 输出双精度量,去归一化的转移函数和系统函数的常数因子. ' NumSec_Z( ) ------ 输出双精度量 系统函数二阶节的分子多项式系数二维数组. ' 元素 NumSec_Z( k, i ) 中, ' k : 二阶节序号;

' i : 多项式系数, i = 0 相应于常数项.

' DenSec_Z( ) ------ 输出双精度量 系统函数二阶节的分母多项式系数二维数组. ' 元素 DenSec_Z( k, i ) 中, ' k : 二阶节序号;

' i : 多项式系数, i = 0 相应于常数项. ' AR( ) ----------------- 输出双精度量,滤波器的幅频特性数组. '

Sub Elliptic(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, H01 As Double, NumSec1() As Double, DenSec1() As Double, H02 As Double, NumSec2() As Double, DenSec2() As Double, NumSec_Z() As Double, DenSec_Z() As Double, AR() As Double) Dim i%, j%, start%, k%, NR%, Flag% Dim F1#, F2#, angle#

Dim ratio(0 To 10) As Double ''''''''''''''''''''''''''' ' 低通滤波器

If PbType = 0 Then

wpass = 2# * Pi * fpass / fsamp wstop = 2# * Pi * fstop / fsamp omikaP = Tan(wpass / 2#) omikaS = Tan(wstop / 2#)

' 调用 Fz_LP 子程序，将低通模拟滤波器的转移函数变量 s 映射为低通数字滤波器的系统函数变量 z Call Fz_LP(Numz(), Denz(), order_z) 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#)

' 调用 Fz_HP 子程序，将高通模拟滤波器的转移函数变量 s 映射高通数字滤波器的系统函数变量 z Call Fz_HP(Numz(), Denz(), order_z) End If

''''''''''''''''''''''''''' ' 带通滤波器

If PbType = 2 Then

wp1 = 2# * Pi * fp1 / fsamp

13

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

' 调用 Fz_BP 子程序，将带通模拟滤波器的转移函数变量 s 映射为带通数字滤波器的系统函数变量 z Call Fz_BP(fp1, fp2, fsamp, Numz(), Denz(), order_z) 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 / fsamp Ci = BpC(wp1, wp2)

omikaP = Abs(Sin(wp1) / (Cos(wp1) - Ci))

omikaS1 = Abs(Sin(ws1) / (Cos(ws1) - Ci)): omikaS2 = Abs(Sin(ws2) / (Cos(ws2) - Ci)) If omikaS1 <= omikaS2 Then omikaS = omikaS1 Else

omikaS = omikaS2 End If

' 调用 Fz_BS 子程序，将带阻模拟滤波器的转移函数变量 s 映射为带阻数字滤波器的系统函数变量 z Call Fz_BS(fp1, fp2, fsamp, Numz(), Denz(), order_z) End If ''''''''''''''''

F1 = omikaP / (2# * Pi): F2 = omikaS / (2# * Pi) ' 求低通模拟椭圆滤波器的归一化转移函数

Call HS_E(F1, F2, Apass, Astop, ord, Flag, H01, NumSec1(), DenSec1())

If ord Mod 2 <> 0 Then

NR = (ord - 1) / 2: start = 0 Else

NR = ord / 2: start = 1 End If

' 将所求得的低通模拟滤波器的归一化转移函数去归一化

Call DeNormHs(Flag, NR, omk0, H01, NumSec1(), DenSec1(), H02, NumSec2(), DenSec2())

' 系统函数由一阶、二阶节级联而成，k 是节序号。此处根据模拟滤波器的一、二阶节求数字滤波器的一、 ‘ 二阶节的分子分母多项式系数 For k = start To NR

Nums(0) = NumSec2(k, 0): Nums(1) = NumSec2(k, 1): Nums(2) = NumSec2(k, 2) Dens(0) = DenSec2(k, 0): Dens(1) = DenSec2(k, 1): Dens(2) = DenSec2(k, 2) If k = 0 Then order_s = 1 Else

order_s = 2

14

End If

' 根据模拟滤波器的一、二阶节求数字滤波器的一、二阶节的分子分母多项式系数 Call BSF2(Nums(), Dens(), order_s, Numz(), Denz(), order_z, num(), Den(), order) ' 将二维数组 NumSec_Z() 和 DenSec_Z() 清零. For i = 0 To 3

NumSec_Z(k, i) = 0 DenSec_Z(k, i) = 0 Next i

' 将系统函数各个二阶节的分子分母多项式的系数分别记入二维数组 NumSec_Z() 和 DenSec_Z(). For j = 0 To order

NumSec_Z(k, j) = num(j) ' 第 k 个二阶节的分子多项式的第 j 个系数 DenSec_Z(k, j) = Den(j) ' 第 k 个二阶节的分母多项式的第 j 个系数 Next j Next k

For i = 0 To points - 1

For k = start To NR ' k 是二阶节序号

angle = i * 2# * Pi / points ' 在单位圆上,第 i 点的幅值为 1, 相角为 angle。此即第 i 个数字频

‘ 率。.

' 给定数字频率，求第 k 个二阶节的频率响应

Call Find_ratio(order, k, angle, NumSec_Z(), DenSec_Z(), ratio(k)) 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) * H02

' 此处已求出第 i 点频率响应 Next i End Sub

' 求 Butterworth 滤波器转移函数的二阶节分子分母多项式的系数组 '

Sub HS_B(n As Integer, i As Integer, omikaP As Double, ep As Double, Nums() As Double, Dens() As Double)

Dim omk0#, arg#

omk0 = omika0(omikaP, ep, n) arg = thetas(n, i)

Nums(0) = 1#: Nums(1) = 0#: Nums(2) = 0#

Dens(0) = 1#: Dens(1) = (-2# / omk0) * Cos(arg): Dens(2) = (1# / omk0) ^ 2 End Sub ‘

' 求切比雪夫1型模拟滤波器转移函数的二阶节分子分母多项式的系数组 '

Sub HS_Ch1(n As Integer, i As Integer, omikaP As Double, ep As Double, Nums() As Double, Dens() As Double, NumSec() As Double, DenSec() As Double) Dim j%

Dim phase#, A#, omk0#, omki#

15

'

phase = thetas(n, i)

A = (1# / n) * arcsinh(1# / ep)

omk0 = omikaP * sinh(A): omki = omikaP * Sin(phase)

Nums(0) = omk0 ^ 2 + omki ^ 2: Nums(1) = 0#: Nums(2) = 0#

Dens(0) = Nums(0): Dens(1) = -2# * omk0 * Cos(phase): Dens(2) = 1# For j = 0 To 2

NumSec(i, j) = Nums(j): DenSec(i, j) = Dens(j) Next j End Sub

' 求切比雪夫2型模拟滤波器转移函数的二阶节分子分母多项式的系数组

'

Sub HS_Ch2(n As Integer, i As Integer, omikaS As Double, eS As Double, num() As Double, Den() As Double)

Dim phase#, A#, omk0#, omki#

phase = thetas(n, i)

A = (1# / n) * arcsinh(eS) omk0 = omikaS / sinh(A) omki = omikaS / Sin(phase)

num(0) = 1#: num(1) = 0: num(2) = omki ^ (-2)

Den(0) = 1#: Den(1) = -2# * (omk0 ^ (-1)) * Cos(phase): Den(2) = omk0 ^ (-2) + omki ^ (-2) End Sub '

' 求归一化椭圆模拟低通滤波器的转移函数 '

Sub HS_E(fpass As Double, fstop As Double, Apass As Double, Astop As Double, ord As Integer, Flag As Integer, H0 As Double, NumSec() As Double, DenSec() As Double) Dim i%, j%, k%, start%, M%, order%, NR%, sign%

Dim angle#, temp#, Au#, PinAu#, Omi#, Omi2#, Vi#, SigVi#, Denor#, ep#, ep2#, B2#, omkC#, D#, q# Dim Q0#, Q04#, Q05#, Q09#, Q013#, Q14#, PiN#, PiN2#, TR#, TN#, Rk#, RkPr# Dim arg#, Arg1#, Arg2#, Sum1#, Sum2#, S0#, SO2#, w#, omk02# Dim ratio(0 To 10) As Double

omk0 = 2# * Pi * Sqr(fpass * fstop) omk02 = omk0 ^ 2

omikaP = 2# * Pi * fpass / omk0 omikaS = 1 / omikaP sign = 1 Sum2 = 0#

ep2 = 10# ^ (0.1 * Apass) ep = Sqr(ep2)

B2 = 10# ^ (0.1 * Astop) D = (B2 - 1#) / (ep2 - 1#) Rk = omikaP ^ 2

RkPr = (1# - Rk * Rk) ^ 0.25

Q0 = 0.5 * (1# - RkPr) / (1# + RkPr) Q04 = Q0 ^ 4 Q05 = Q0 ^ 5 Q09 = Q04 * Q05

16