Digital Signal Processing. The Scientist and Engineers Guide. Steven W. Smith
.pdf
Chapter 20Chebyshev Filters |
339 |
Stability
The main limitation of digital filters carried out by convolution is execution time. It is possible to achieve nearly any filter response, provided you are willing to wait for the result. Recursive filters are just the opposite. They run like lightning; however, they are limited in performance. For example, consider a 6 pole, 0.5% ripple, low-pass filter with a 0.01 cutoff frequency. The recursion coefficients for this filter can be obtained from Table 20-1:
a0= 1.391351E-10
a1= 8.348109E-10 b1= 5.883343E+00 a2= 2.087027E-09 b2= -1.442798E+01 a3= 2.782703E-09 b3= 1.887786E+01 a4= 2.087027E-09 b4= -1.389914E+01 a5= 8.348109E-10 b5= 5.459909E+00 a6= 1.391351E-10 b6= -8.939932E-01
Look carefully at these coefficients. The "b" coefficients have an absolute value of about ten. Using single precision, the round-off noise on each of these numbers is about one ten-millionth of the value, i.e., 10&6 . Now look at the "a" coefficients, with a value of about 10&9 . Something is obviously wrong here. The contribution from the input signal (via the "a" coefficients) will be 1000 times smaller than the noise from the previously calculated output signal (via the "b" coefficients). This filter won't work! In short, round-off noise limits the number of poles that can be used in a filter. The actual number will depend slightly on the ripple and if it is a high or low-pass filter. The approximate numbers for single precision are:
TABLE 20-3 |
Cutoff frequency |
0.02 |
0.05 |
0.10 |
0.25 |
0.40 |
0.45 |
0.48 |
The maximum number of |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
poles for single precision. |
Maximum poles |
4 |
6 |
10 |
20 |
10 |
6 |
4 |
|
The filter's performance will start to degrade as this limit is approached; the step response will show more overshoot, the stopband attenuation will be poor, and the frequency response will have excessive ripple. If the filter is pushed too far, or there is an error in the coefficients, the output will probably oscillate until an overflow occurs.
There are two ways of extending the maximum number of poles that can be used. First, use double precision. This requires using double precision in the coefficient calculation as well (including the value for pi ).
The second method is to implement the filter in stages. For example, a six pole filter starts out as a cascade of three stages of two poles each. The program in Table 20-4 combines these three stages into a single set of recursion coefficients for easier programming. However, the filter is more stable if carried out as the original three separate stages. This requires knowing the "a" and "b" coefficients for each of the stages. These can
340 |
The Scientist and Engineer's Guide to Digital Signal Processing |
||||
100 'CHEBYSHEV FILTERRECURSION COEFFICIENT CALCULATION |
|||||
110 |
' |
|
|
|
|
120 |
|
|
|
|
'INITIALIZE VARIABLES |
130 |
DIM A[22] |
|
|
|
'holds the "a" coefficients upon program completion |
140 |
DIM B[22] |
|
|
|
'holds the "b" coefficients upon program completion |
150 |
DIM TA[22] |
|
|
|
'internal use for combining stages |
160 |
DIM TB[22] |
|
|
|
'internal use for combining stages |
170 |
' |
|
|
|
|
180 FOR I% = 0 TO 22 |
|
|
|
||
190 |
A[I%] = 0 |
|
|
|
|
200 |
B[I%] = 0 |
|
|
|
|
210 NEXT I% |
|
|
|
|
|
220 |
' |
|
|
|
|
230 |
A[2] = 1 |
|
|
|
|
240 |
B[2] = 1 |
|
|
|
|
250 |
PI = 3.14159265 |
|
|
|
|
260 |
|
|
|
|
'ENTER THE FOUR FILTER PARAMETERS |
270 |
INPUT "Enter cutoff frequency |
(0 to .5): ", FC |
|||
280 |
INPUT "Enter |
0 for LP, 1 for HP filter: |
", LH |
||
290 |
INPUT "Enter percent ripple |
(0 to 29): |
", PR |
||
300 |
INPUT "Enter number of poles (2,4,...20): |
", NP |
|||
310 |
' |
|
|
|
|
320 FOR P% = 1 TO NP/2 |
|
|
'LOOP FOR EACH POLE-PAIR |
||
330 |
' |
|
|
|
|
340 |
GOSUB 1000 |
|
|
'The subroutine in TABLE 20-5 |
|
350 |
' |
|
|
|
|
360 |
FOR I% = 0 TO 22 |
|
|
'Add coefficients to the cascade |
|
370 |
TA[I%] = A[I%] |
|
|
|
|
380 |
TB[I%] = B[I%] |
|
|
|
|
390 |
NEXT I% |
|
|
|
|
400 |
' |
|
|
|
|
410 |
FOR I% = 2 TO 22 |
|
|
|
|
420 |
A[I%] = A0*TA[I%] + A1*TA[I%-1] + A2*TA[I%-2] |
||||
430 |
B[I%] = |
TB[I%] - B1*TB[I%-1] - B2*TB[I%-2] |
|||
440 |
NEXT I% |
|
|
|
|
450 |
' |
|
|
|
|
460 NEXT P% |
|
|
|
|
|
470 |
' |
|
|
|
|
480 |
B[2] = 0 |
|
|
|
'Finish combining coefficients |
490 FOR I% = 0 TO 20 |
|
|
|
||
500 |
A[I%] = A[I%+2] |
|
|
|
|
510 |
B[I%] = -B[I%+2] |
|
|
|
|
520 NEXT I% |
|
|
|
|
|
530 |
' |
|
|
|
|
540 SA = 0 |
|
|
|
'NORMALIZE THE GAIN |
|
550 |
SB = 0 |
|
|
|
|
560 FOR I% = 0 TO 20 |
|
|
|
||
570 |
IF LH = 0 THEN SA |
= SA |
+ A[I%] |
|
|
580 |
IF LH = 0 THEN SB |
= SB |
+ B[I%] |
|
|
590 |
IF LH = 1 THEN SA |
= SA |
+ A[I%] * (-1)^I% |
||
600 |
IF LH = 1 THEN SB |
= SB |
+ B[I%] * (-1)^I% |
||
610 NEXT I% |
|
|
|
|
|
620 |
' |
|
|
|
|
630 |
GAIN = SA / (1 - SB) |
|
|
|
|
640 |
' |
|
|
|
|
650 FOR I% = 0 TO 20 |
|
|
|
||
660 |
A[I%] = A[I%] / GAIN |
|
|
||
670 NEXT I% |
|
|
|
|
|
680 |
' |
|
|
|
'The final recursion coefficients are in A[ ] and B[ ] |
690 END |
|
|
|
|
|
TABLE 20-4
|
|
Chapter 20Chebyshev Filters |
341 |
|
1000 'THIS SUBROUTINE IS CALLED FROM TABLE 20-4, LINE 340 |
|
|||
1010 ' |
|
|
|
|
1020 |
' |
Variables entering subroutine: |
PI, FC, LH, PR, HP, P% |
|
1030 |
' |
Variables exiting subroutine: |
A0, A1, A2, B1, B2 |
|
1040 |
' |
Variables used internally: |
RP, IP, ES, VX, KX, T, W, M, D, K, |
|
1050 ' |
|
X0, X1, X2, Y1, Y2 |
|
|
1060 ' |
|
|
|
|
1070 ' |
|
'Calculate the pole location on the unit circle |
|
|
1080 |
RP = -COS(PI/(NP*2) + (P%-1) * PI/NP) |
|
||
1090 |
IP = SIN(PI/(NP*2) + (P%-1) * PI/NP) |
|
|
|
1100 ' |
|
|
|
|
1110 ' |
|
'Warp from a circle to an ellipse |
|
|
1120 IF PR = 0 THEN GOTO 1210 |
|
|
||
1130 |
ES = SQR( (100 / (100-PR))^2 -1 ) |
|
|
|
1140 |
VX = (1/NP) * LOG( (1/ES) + SQR( (1/ES^2) + 1) ) |
|
||
1150 |
KX = (1/NP) * LOG( (1/ES) + SQR( (1/ES^2) - 1) ) |
|
||
1160 |
KX = (EXP(KX) + EXP(-KX))/2 |
|
|
|
1170 |
RP = RP * ( (EXP(VX) - EXP(-VX) ) /2 ) / KX |
|
||
1180 |
IP = IP * ( (EXP(VX) + EXP(-VX) ) /2 ) / KX |
|
||
1190 ' |
|
|
|
|
1200 ' |
|
's-domain to z-domain conversion |
|
|
1210 |
T = 2 * TAN(1/2) |
|
|
|
1220 |
W = 2*PI*FC |
|
|
|
1230 |
M = RP^2 + IP^2 |
|
|
|
1240 |
D = 4 - 4*RP*T + M*T^2 |
|
|
|
1250 |
X0 = T^2/D |
|
|
|
1260 |
X1 = 2*T^2/D |
|
|
|
1270 |
X2 = T^2/D |
|
|
|
1280 |
Y1 = (8 - 2*M*T^2)/D |
|
|
|
1290 |
Y2 = (-4 - 4*RP*T - M*T^2)/D |
|
|
|
1300 ' |
|
|
|
|
1310 ' |
|
'LP TO LP, or LP TO HP transform |
|
|
1320 |
IF LH = 1 THEN K = -COS(W/2 + 1/2) / COS(W/2 - 1/2) |
|
||
1330 |
IF LH = 0 THEN K = SIN(1/2 - W/2) / SIN(1/2 + W/2) |
|
||
1340 |
D = 1 + Y1*K - Y2*K^2 |
|
|
|
1350 |
A0 = (X0 - X1*K + X2*K^2)/D |
|
|
|
1360 |
A1 = (-2*X0*K + X1 + X1*K^2 - 2*X2*K)/D |
|
||
1370 |
A2 = (X0*K^2 - X1*K + X2)/D |
|
|
|
1380 |
B1 = (2*K + Y1 + Y1*K^2 - 2*Y2*K)/D |
|
|
|
1390 |
B2 = (-(K^2) - Y1*K + Y2)/D |
|
|
|
1400 IF LH = 1 THEN A1 = -A1
1410 IF LH = 1 THEN B1 = -B1 1420 '
1430 RETURN
TABLE 20-5
TABLE 20-4 and 20-5
Program to calculate the "a" and "b" coefficients for Chebyshev recursive filters. In lines 270-300, four parameters are entered into the program. The cutoff frequency, FC, is expressed as a fraction of the sampling frequency, and therefore must be in the range: 0 to 0.5. The variable, LH, is set to a value of one for a high-pass filter, and zero for a low-pass filter. The value entered for PR must be in the range of 0 to 29, corresponding to 0 to 29% ripple in the filter's frequency response. The number of poles in the filter, entered in the variable NP, must be an even integer between 2 and 20. At the completion of the program, the "a" and "b" coefficients are stored in the arrays A[ ] and B[ ] (a0 = A[0], a1 = A[1], etc.). TABLE 20-5 is a subroutine called from line 340 of the main program. Six variables are passed to this subroutine, and five variables are returned. Table 20-6 (next page) contains two sets of data to help debug this subroutine. The functions: COS and SIN, use radians, not degrees. The function: LOG is the natural (base e) logarithm. Declaring all floating point variables (including the value of B) to be double precision will allow more poles to be used. Tables 20-1 and 20-2 were generated with this program and can be used to test for proper operation. Chapter 33 describes the mathematical operation of this program.
342 |
The Scientist and Engineer's Guide to Digital Signal Processing |
|||||
|
DATA SET 1 |
DATA SET 2 |
||||
|
Enter the subroutine with these values: |
|
|
|
||
|
FC |
= |
0.1 |
FC |
= |
0.1 |
|
LH |
= |
0 |
LH |
= |
1 |
|
PR |
= |
0 |
PR |
= |
10 |
|
NP |
= |
4 |
NP |
= |
4 |
|
P% |
= |
1 |
P% |
= |
2 |
|
PI |
= |
3.141592 |
PI |
= |
3.141592 |
These values should be present at line 1200:
RP |
= -0.923879 |
RP |
= -0.136178 |
||
IP |
= |
0.382683 |
IP |
= |
0.933223 |
ES |
= |
not used |
ES |
= |
0.484322 |
VX |
= |
not used |
VX |
= |
0.368054 |
KX |
= |
not used |
KX |
= |
1.057802 |
These values should be present at line 1310:
T |
= |
1.092605 |
T |
= |
1.092605 |
W |
= |
0.628318 |
W |
= |
0.628318 |
M |
= |
1.000000 |
M |
= |
0.889450 |
D |
= |
9.231528 |
D |
= |
5.656972 |
X0 |
= |
0.129316 |
X0 |
= |
0.211029 |
X1 |
= |
0.258632 |
X1 |
= |
0.422058 |
X2 |
= |
0.129316 |
X2 |
= |
0.211029 |
Y1 |
= |
0.607963 |
Y1 |
= |
1.038784 |
Y2 |
= -0.125227 |
Y2 |
= -0.789584 |
||
These values should be return to the main program:
A0 |
= |
0.061885 |
A0 |
= |
0.922919 |
A1 |
= |
0.123770 |
A1 |
= -1.845840 |
|
A2 |
= |
0.061885 |
A2 |
= |
0.922919 |
B1 |
= |
1.048600 |
B1 |
= |
1.446913 |
B2 |
= |
-0.296140 |
B2 |
= |
-0.836653 |
TABLE 20-6
Debugging data. This table contains two sets of data for debugging the subroutine listed in Table 20-5.
be obtained from the program in Table 20-4. The subroutine in Table 20-5 is called once for each stage in the cascade. For example, it is called three times for a six pole filter. At the completion of the subroutine, five variables are return to the main program: A0, A1, A2, B1, & B2. These are the recursion coefficients for the two pole stage being worked on, and can be used to implement the filter in stages.
344 |
The Scientist and Engineer's Guide to Digital Signal Processing |
capacitors. Even if a Butterworth response is designed (i.e., 0% ripple), filters of this complexity will have a residue ripple of, perhaps, 1%. On the other hand, the flatness of digital filters is primarily limited by round-off error, making them hundreds of times flatter than their analog counterparts. Score one point for the digital filter.
Next, look at the frequency response on a log scale, as shown in (c) and (d). Again, the digital filter is clearly the victor in both roll-off and stopband attenuation. Even if the analog performance is improved by adding additional stages, it still can't compare to the digital filter. For instance, imagine that you need to improve these two parameters by a factor of 100. This can be done with simple modifications to the windowed-sinc, but is virtually impossible for the analog circuit. Score two more for the digital filter.
The step response of the two filters is shown in (e) and (f). The digital filter's step response is symmetrical between the lower and upper portions of the step, i.e., it has a linear phase. The analog filter's step response is not symmetrical, i.e., it has a nonlinear phase. One more point for the digital filter. Lastly, the analog filter overshoots about 20% on one side of the step. The digital filter overshoots about 10%, but on both sides of the step. Since both are bad, no points are awarded.
In spite of this beating, there are still many applications where analog filters should, or must, be used. This is not related to the actual performance of the filter (i.e., what goes in and what comes out), but to the general advantages that analog circuits have over digital techniques. The first advantage is speed: digital is slow; analog is fast. For example, a personal computer can only filter data at about 10,000 samples per second, using FFT convolution. Even simple op amps can operate at 100 kHz to 1 MHz, 10 to 100 times as fast as the digital system!
The second inherent advantage of analog over digital is dynamic range. This comes in two flavors. Amplitude dynamic range is the ratio between the largest signal that can be passed through a system, and the inherent noise of the system. For instance, a 12 bit ADC has a saturation level of 4095, and an rms quantization noise of 0.29 digital numbers, for a dynamic range of about 14000. In comparison, a standard op amp has a saturation voltage of about 20 volts and an internal noise of about 2 microvolts, for a dynamic range of about ten million. Just as before, a simple op amp devastates the digital system.
The other flavor is frequency dynamic range. For example, it is easy to design an op amp circuit to simultaneously handle frequencies between 0.01 Hz and 100 kHz (seven decades). When this is tried with a digital system, the computer becomes swamped with data. For instance, sampling at 200 kHz, it takes 20 million points to capture one complete cycle at 0.01 Hz. You may have noticed that the frequency response of digital filters is almost always plotted on a linear frequency scale, while analog filters are usually displayed with a logarithmic frequency. This is because digital filters need
