Fortran语言设计切比雪夫滤波器的特性分析

作者：佚名点击：135次评论：4次标签： Fortran 滤波器

类型：休闲益智大小：113M语言：中文评分：10.0
标签：

立即下载

一、实验背景
数字滤波器是DSP中非常重要的组成部分。滤波器有两种用途：分离混合的信号，复原失真的信号。在基频分析器中，伴奏是对人声的一种噪音干扰，因此需要尽可能减少噪音的影响。由于人的音高在300Hz~2000Hz之间，因此可以构造一个带通滤波器去除这个频段以外的信号。在常见的滤波器中，由于切比雪夫滤波器运算性能较好，且工作在频域上，因此最适合应用在基频分析器中。

二、实验目的
分析切比雪夫滤波器的特性，并完成切比雪夫系数的设计。

三、实验环境和工具
实验平台：Windows
开发语言：Fortran

四、常见滤波器的分类和切比雪夫滤波器的选择

根据应用场合和实现方式，常见滤波器可以如表1所示来进行分类：

	卷积	递归
时域	移动平均滤波器	单极点滤波器
频域	加窗sinc滤波器	切比雪夫滤波器
任意频响	FIR任意频响滤波器	迭代设计滤波器

表1 常见滤波器的分类

根据表1，由于基频的分析提取是工作在频域上的，因此滤波器主要在加窗sinc滤波器和切比雪夫滤波器两者间进行选择。根据系统的需求，为了实现快速获取基频分析结果，选择的标准主要取决于两者在速度上的优劣。

对加窗sinc和切比雪夫滤波器的速度的比较，即是对卷积滤波器（FIR）和递归滤波器（IIR）的比较。图1的曲线展示了加窗sinc滤波器和等效的6极点递归滤波器的相应执行时间。其中，标准卷积和快速傅里叶卷积是两种算法下构造的加窗sinc滤波器。加窗sinc滤波器的运行时间在低频和高频段上升，这是因为为了跟上递归滤波器在这些效率上更好的性能，必须增加滤波器内核的长度。总的来说，IIR滤波器要比相应的FIR滤波器快1个数量级。

图1 FIR滤波器和IIR滤波器的速度对比

五、切比雪夫滤波器的特性分析

切比雪夫滤波器（又译车比雪夫滤波器）是在通带或阻带上频率响应幅度等波纹波动的滤波器。在通带波动的为“I型切比雪夫滤波器”，在阻带波动的为“II型切比雪夫滤波器”。

1. 切比雪夫的频率响应特性

切比雪夫的频率响应如图2所示：

（a）低通频率响应

（b）低通频率响应（dB）

（c）高通频率响应

（d）高通频率响应（dB）

图2 切比雪夫频率响应

如图1所示，图a和图b是0.5%纹波的低通切比雪夫滤波器的频率响应，图c和图d是相应的高通滤波器的响应。

2. 阶跃响应的过冲

切比雪夫滤波器的阶跃响应有5%~30%的过冲，且随着极点数的增加而增大。图3是切比雪夫滤波器阶跃响应的两个例子。

图3 切比雪夫阶跃响应

3. 稳定性

递归滤波器运行很快，但在性能上受限。如果是低通高通滤波器，极点实际的数量将略微取决于纹波。单精度的近似数字如表2所示。

截止效率	0.02	0.05	0.10	0.25	0.40	045	0.48
最大极点数	4	6	10	20	10	6	4

当接近这一极限时，滤波器的性能将变差。阶跃响应将出现更大的过冲，阻带衰减将变差，频率响应将出现过多纹波。如果滤波器的设计指针太高，或系数中出现了错误，输出可能振荡直至有溢出发生。

为了增加极点的最大数目，以提高运算结果的准确性，可以采用分阶段实现滤波器，在本系统的基频分析模块中，我们将滤波器分为9级。如果每一级产生n个极点，则该滤波器将会产生9*n个极点的滤波器。

六、切比雪夫滤波器的设计

递归滤波器是由差分方程描述的。

（1）

其中，x[]和y[]分别是输入和输出信号，a和b是递归系数。可以写成下面的形式：

（2）

滤波器的设计关键就在于a和b系数的获取。

递归系数可以使用下面的Fortran程序来获得：

1: 100 'CHEBYSHEV FILTER- RECURSION COEFFICIENT CALCULATION
2:
3: 110 '
4:
5: 120 'INITIALIZE VARIABLES
6:
7: 130 DIM A[22] 'holds the "a" coefficients upon program completion
8:
9: 140 DIM B[22] 'holds the "b" coefficients upon program completion
10:
11: 150 DIM TA[22] 'internal use for combining stages
12:
13: 160 DIM TB[22] 'internal use for combining stages
14:
15: 170 '
16:
17: 180 FOR I% = 0 TO 22
18:
19: 190 A[I%] = 0
20:
21: 200 B[I%] = 0
22:
23: 210 NEXT I%
24:
25: 220 '
26:
27: 230 A[2] = 1
28:
29: 240 B[2] = 1
30:
31: 250 PI = 3.14159265
32:
33: 260 'ENTER THE FOUR FILTER PARAMETERS
34:
35: 270 INPUT "Enter cutoff frequency (0 to .5): ", FC
36:
37: 280 INPUT "Enter 0 for LP, 1 for HP filter: ", LH
38:
39: 290 INPUT "Enter percent ripple (0 to 29): ", PR
40:
41: 300 INPUT "Enter number of poles (2,4,...20): ", NP
42:
43: 310 '
44:
45: 320 FOR P% = 1 TO NP/2 'LOOP FOR EACH POLE-PAIR
46:
47: 330 '
48:
49: 340 GOSUB 1000 'The subroutine in TABLE 20-5
50:
51: 350 '
52:
53: 360 FOR I% = 0 TO 22 'Add coefficients to the cascade
54:
55: 370 TA[I%] = A[I%]
56:
57: 380 TB[I%] = B[I%]
58:
59: 390 NEXT I%
60:
61: 400 '
62:
63: 410 FOR I% = 2 TO 22
64:
65: 420 A[I%] = A0*TA[I%] + A1*TA[I%-1] + A2*TA[I%-2]
66:
67: 430 B[I%] = TB[I%] - B1*TB[I%-1] - B2*TB[I%-2]
68:
69: 440 NEXT I%
70:
71: 450 '
72:
73: 460 NEXT P%
74:
75: 470 '
76:
77: 480 B[2] = 0 'Finish combining coefficients
78:
79: 490 FOR I% = 0 TO 20
80:
81: 500 A[I%] = A[I%+2]
82:
83: 510 B[I%] = -B[I%+2]
84:
85: 520 NEXT I%
86:
87: 530 '
88:
89: 540 SA = 0 'NORMALIZE THE GAIN
90:
91: 550 SB = 0
92:
93: 560 FOR I% = 0 TO 20
94:
95: 570 IF LH = 0 THEN SA = SA + A[I%]
96:
97: 580 IF LH = 0 THEN SB = SB + B[I%]
98:
99: 590 IF LH = 1 THEN SA = SA + A[I%] * (-1)^I%
100:
101: 600 IF LH = 1 THEN SB = SB + B[I%] * (-1)^I%
102:
103: 610 NEXT I%
104:
105: 620 '
106:
107: 630 GAIN = SA / (1 - SB)
108:
109: 640 '
110:
111: 650 FOR I% = 0 TO 20
112:
113: 660 A[I%] = A[I%] / GAIN
114:
115: 670 NEXT I%
116:
117: 680 ' 'The final recursion coefficients are in A[ ] and B[ ]
118:
119: 690 END

第270~300行，4个参数被输入到程序中。截止频率FC，被表达为抽样频率的分数形式，因此必须在0到0.5的范围内。变量LH，设置为1时对应高通滤波器，设置为0时对应低通滤波器。变量PR的输入值要在0到29之间，表示滤波器频率响应的纹波在0%到29%之间。滤波器的极点数通过变量NP输入，且必须是2到20见的一个偶数。程序运行完成时，a和b系数存放在数组A[]和B[]中（a0 = A[0]，a[1] = A[1]，等等）。

上面的程序只能用于获取高通或低通滤波器的系数，对于带通滤波器，则需要再进一步进行构造。在DSP系统中，带通滤波器可通过低通滤波器和高通滤波器级联而形成，如图4所示。

图4 级联构成的带通滤波器

图4所示的带通滤波器由一个低通滤波器和一个高通滤波器串联而成。设两个滤波器的系数分别为a0、a1、a2、b1、b2和A0、A1、A2、B1、B2，则整个系统的传递函数H[z]由两级传递函数相乘得到

（3）

进行多项式相乘和合并同类项后得到下面的形式

（4）

公式（4）符合公式（2）的形式，因此我们可以直接提取实现串联系统的递归系数

上面的计算过程可以交由计算机程序来实现。如下所示：

1: 100 'COMBINING RECURSION COEFFICIENTS OF CASCADE AND PARALLEL STAGES
2:
3: 110 '
4:
5: 120 ' 'INITIALIZE VARIABLES
6:
7: 130 DIM A1[8], B1[8] 'a and b coefficients for system 1, one of the stages
8:
9: 140 DIM A2[8], B2[8] 'a and b coefficients for system 2, one of the stages
10:
11: 150 DIM A3[16], B3[16] 'a and b coefficients for system 3, the combined system
12:
13: 160 '
14:
15: 170 'Indicate cascade or parallel combination
16:
17: 180 INPUT "Enter 0 for cascade, 1 for parallel: ", CP%
18:
19: 190 '
20:
21: 200 GOSUB XXXX 'Mythical subroutine to load: A1[ ], B1[ ], A2[ ], B2[ ]
22:
23: 210 '
24:
25: 220 FOR I% = 0 TO 8 'Convert the recursion coefficients into transfer functions
26:
27: 230 B2[I%] = -B2[I%]
28:
29: 240 B1[I%] = -B1[I%]
30:
31: 250 NEXT I%
32:
33: 260 B1[0] = 1
34:
35: 270 B2[0] = 1
36:
37: 280 '
38:
39: 290 FOR I% = 0 TO 16 'Multiply the polynomials by convolving
40:
41: 300 A3[I%] = 0
42:
43: 310 B3[I%] = 0
44:
45: 320 FOR J% = 0 TO 8
46:
47: 330 IF I%-J% < 0 OR I%-J% > 8 THEN GOTO 370
48:
49: 340 IF CP% = 0 THEN A3[I%] = A3[I%] + A1[J%] * A2[I%-J%]
50:
51: 350 IF CP% = 1 THEN A3[I%] = A3[I%] + A1[J%] * B2[I%-J%] + A2[J%] * B1[I%-J%]
52:
53: 360 B3[I%] = B3[I%] + B1[J%] * B2[I%-J%]
54:
55: 370 NEXT J%
56:
57: 380 NEXT I%
58:
59: 390 '
60:
61: 400 FOR I% = 0 TO 16 'Convert the transfer function into recursion coefficients.
62:
63: 410 B3[I%] = -B3[I%]
64:
65: 420 NEXT I%
66:
67: 430 B3[0] = 0
68:
69: 440 ' 'The recursion coefficients of the combined system now
70:
71: 450 END 'reside in A3[ ] & B3[ ]

在本系统中，选用44100的采样频率，高通频率为30Hz，低通频率为2Khz，为了获得更快滚降特性，选用4%的纹波。使用上面的程序获得a、b系数如下：

	A	B
0	9.684457e-07	0
1	6.779120e-06	-7.616874e+00
2	1.936891e-05	2.634041e+01
3	2.711648e-05	-5.424242e+01
4	1.355824e-05	7.325881e+01
5	-1.355824e-05	-6.726091e+01
6	-2.711648e-05	4.196424e+01
7	-1.936891e-05	-1.715111e+01
8	-6.779120e-06	4.166027e+00
9	-9.684457e-07	-4.581673e-01