programs.txt

Guilde To DSP 讲DSP(数字信号处理)原理的好书, 练习英语阅读能力的好机会. ~_^ seabird Nov 13,1

1350     JM1% = J%-1
1360     FOR I% = JM1% TO NM1% STEP LE%     'Loop for each butterfly
1370       IP% = I%+LE2%
1380       TR = REX[IP%]*UR - IMX[IP%]*UI   'Butterfly calculation
1390       TI = REX[IP%]*UI + IMX[IP%]*UR
1400       REX[IP%] = REX[I%]-TR
1410       IMX[IP%] = IMX[I%]-TI
1420       REX[I%]  = REX[I%]+TR
1430       IMX[I%]  = IMX[I%]+TI   
1440     NEXT I%
1450     TR = UR
1460     UR = TR*SR - UI*SI
1470     UI = TR*SI + UI*SR
1480   NEXT J%
1490 NEXT L%
1500 '
1510 RETURN


TABLE 12-5
2000 'INVERSE FAST FOURIER TRANSFORM SUBROUTINE
2010 'Upon entry, N% contains the number of points in the IDFT, REX[ ] & IMX[]
2020 'contain the real & imaginary parts of the complex frequency domain. 
2030 'Upon return, REX[ ] and IMX[ ] contain the complex time domain.
2040 'All signals run from 0 to N%-1.
2050 '
2060 FOR K% = 0 TO N%-1           'Change the sign of IMX[ ]
2070   IMX[K%] = -IMX[K%]
2080 NEXT K%
2090 '
2100 GOSUB 1000                   'Calculate forward FFT  (Table 12-3)
2110 '
2120 FOR I% = 0 TO N%-1           'Divide the time domain by N% and
2130   REX[I%] =  REX[I%]/N% 'change the sign of IMX[ ]
2140   IMX[I%] = -IMX[I%]/N%
2150 NEXT I%
2160 '
2170 RETURN


TABLE 12-6
4000 'INVERSE FFT FOR REAL SIGNALS
4010 'Upon entry, N% contains the number of points in the IDFT, REX[ ] and
4020 'IMX[ ] contain the real & imaginary parts of the frequency domain 
4030 'running from index 0 to N%/2.  The remaining samples in REX[ ] and IMX[] 
4040 'are ignored. Upon return, REX[ ] contains the real time domain, IMX[ ]
4045 'contains zeros.
4050 '
4060 '
4070 FOR K% = (N%/2+1) TO (N%-1)      'Make frequency domain symmetrical
4080   REX[K%] =  REX[N%-K%]          '(as in Table 12-1)
4090   IMX[K%] = -IMX[N%-K%]
4100 NEXT K%
4110 '
4120 FOR K% = 0 TO N%-1               'Add real and imaginary parts together
4130   REX[K%] =  REX[K%]+IMX[K%]
4140 NEXT K%
4150 '
4160 GOSUB 3000                       'Calculate forward real DFT (TABLE 12-6)
4170 '
4180 FOR I% = 0 TO N%-1               'Add real and imaginary parts together 
4190   REX[I%] = (REX[I%]+IMX[I%])/N% 'and divide the time domain by N%
4200   IMX[I%] = 0
4210 NEXT I%
4220 '
4230 RETURN


TABLE 12-7
3000 'FFT FOR REAL SIGNALS
3010 'Upon entry, N% contains the number of points in the DFT, REX[ ] contains
3020 'the real input signal, while values in IMX[ ] are ignored.  Upon return,
3030 'REX[ ] & IMX[ ] contain the DFT output. All signals run from 0 to N%-1.
3040 '
3050 NH% = N%/2-1                   'Separate even and odd points
3060 FOR I% = 0 TO NH%
3070   REX(I%) = REX(2*I%)
3080   IMX(I%) = REX(2*I%+1)
3090 NEXT I%
3100 '
3110 N% = N%/2                      'Calculate N%/2 point FFT
3120 GOSUB 1000                     '(GOSUB 1000 is the FFT in Table 12-3)
3130 N% = N%*2
3140 '
3150 NM1% = N%-1                    'Even/odd frequency domain decomposition
3160 ND2% = N%/2
3170 N4% = N%/4-1                                
3180 FOR I% = 1 TO N4%
3190   IM%  = ND2%-I%
3200   IP2% = I%+ND2%
3210   IPM% = IM%+ND2%
3220   REX(IP2%) =  (IMX(I%) + IMX(IM%))/2
3230   REX(IPM%) =  REX(IP2%)
3240   IMX(IP2%) = -(REX(I%) - REX(IM%))/2
3250   IMX(IPM%) = -IMX(IP2%)
3260   REX(I%)   =  (REX(I%) + REX(IM%))/2
3270   REX(IM%)  =  REX(I%)
3280   IMX(I%)   =  (IMX(I%) - IMX(IM%))/2
3290   IMX(IM%)  = -IMX(I%)
3300 NEXT I%
3310 REX(N%*3/4) = IMX(N%/4)
3320 REX(ND2%) = IMX(0)
3330 IMX(N%*3/4) = 0
3340 IMX(ND2%) = 0
3350 IMX(N%/4) = 0
3360 IMX(0) = 0
3370 '
3380 PI = 3.14159265                'Complete the last FFT stage
3390 L% = CINT(LOG(N%)/LOG(2))                   
3400 LE% = CINT(2^L%)
3410 LE2% = LE%/2
3420 UR = 1
3430 UI = 0
3440 SR =  COS(PI/LE2%)
3450 SI = -SIN(PI/LE2%)
3460 FOR J% = 1 TO LE2%
3470   JM1% = J%-1
3480   FOR I% = JM1% TO NM1% STEP LE%
3490     IP% = I%+LE2%
3500     TR = REX[IP%]*UR - IMX[IP%]*UI
3510     TI = REX[IP%]*UI + IMX[IP%]*UR
3520     REX[IP%] = REX[I%]-TR
3530     IMX[IP%] = IMX[I%]-TI
3540     REX[I%]  = REX[I%]+TR
3550     IMX[I%]  = IMX[I%]+TI
3560   NEXT I%
3570   TR = UR
3580   UR = TR*SR - UI*SI
3590   UI = TR*SI + UI*SR
3600 NEXT J%
3610 RETURN


TABLE 15-1
100 'MOVING AVERAGE FILTER
110 'This program filters 5000 samples with a 101 point moving
120 'average filter, resulting in 4900 samples of filtered data.
130 '
140 DIM X[4999]                    'X[ ] holds the input signal
150 DIM Y[4999]                    'Y[ ] holds the output signal
160 '
170 GOSUB XXXX                     'Mythical subroutine to load X[ ]
180 '
190 FOR I% = 50 TO 4949            'Loop for each point in the output signal
200   Y[I%] = 0                    'Zero, so it can be used as an accumulator 
210   FOR J% = -50 TO 50           'Calculate the summation
220     Y[I%] = Y[I%] + X(I%+J%]
230   NEXT J%
240   Y[I%] = Y[I%]/101            'Complete the average by dividing
250 NEXT I%
260 '
270 END


TABLE 15-2
100 'MOVING AVERAGE FILTER IMPLEMENTED BY RECURSION
110 'This program filters 5000 samples with a 101 point moving
120 'average filter, resulting in 4900 samples of filtered data.
130 'A double precision accumulator is used to prevent round-off drift.
140 '
150 DIM X[4999]              'X[ ] holds the input signal
160 DIM Y[4999]              'Y[ ] holds the output signal
170 DEFDBL  ACC              'Define the variable ACC to be double precision
180 '
190 GOSUB XXXX               'Mythical subroutine to load X[ ]
200 '
210 ACC = 0                  'Find Y[50] by averaging points X[0] to X[100]
220 FOR I% = 0 TO 100        
230   ACC = ACC + X[I%]
240 NEXT I%
250 Y[[50] = ACC/101
260 '                        'Recursive moving average filter (Eq. 15-3)
270 FOR I% = 51 TO 4949
280   ACC = ACC + X[I%+50] - X[I%-51]
290   Y[I%] = ACC
300 NEXT I%                       
310 '
320 END


TABLE 16-1
100 'LOW-PASS WINDOWED-SINC FILTER 
110 'This program filters 5000 samples with a 101 point windowed-sinc filter, 
120 'resulting in 4900 samples of filtered data.
130 '
140 DIM X[4999]             'X[ ] holds the input signal
150 DIM Y[4999]             'Y[ ] holds the output signal
160 DIM H[100]              'H[ ] holds the filter kernel
170 '
180 PI = 3.14159265
190 FC = .14                'Set the cutoff frequency (between 0 and 0.5)
200 M% = 100                'Set filter length (101 points)
210 '
220 GOSUB XXXX              'Mythical subroutine to load X[ ]
230 '
240 '                       'Calculate the low-pass filter kernel via Eq. 16-4
250 FOR I% = 0 TO 100
260    IF (I%-M%/2) = 0 THEN H[I%] = 2*PI*FC
270    IF (I%-M%/2) <> 0  THEN H[I%] = SIN(2*PI*FC * (I%-M%/2)) / (I%-M%/2)
280    H[I%] = H[I%] * (0.54 - 0.46*COS(2*PI*I%/M%) )
290 NEXT I%
300 '
310 SUM = 0                 'Normalize the low-pass filter kernel for 
320 FOR I% = 0 TO 100       'unity gain at DC
330    SUM = SUM + H[I%]
340 NEXT I%
350 '
360 FOR I% = 0 TO 100
370    H[I%] = H[I%] / SUM
380 NEXT I%
390 '               
400 FOR J% = 100 TO 4999    'Convolve the input signal & filter kernel
410    Y[J%] = 0                  
420    FOR I% = 0 TO 100
430       Y[J%] = Y[J%] + X[J%-I%] * H[I%]
440    NEXT I%
450 NEXT J%
460 '
470 END


TABLE 16-2
100 'BAND-PASS WINDOWED-SINC FILTER
110 'This program calculates an 801 point band-pass filter kernel
120 '
130 DIM A[800]           'A[ ]  workspace for the lower cutoff
140 DIM B[800]           'B[ ]  workspace for the upper cutoff
150 DIM H[800]           'H[ ]  holds the final filter kernel
160 '
170 PI = 3.1415926
180 M% = 800             'Set filter kernel length (801 points)
190 '
200 '                    'Calculate the first low-pass filter kernel via Eq.
210 FC = 0.196           '16-4, with a cutoff frequency of 0.196, store in A[]
220 FOR I% = 0 TO 800
230    IF (I%-M%/2) = 0 THEN A[I%] = 2*PI*FC
240    IF (I%-M%/2) <> 0  THEN A[I%] = SIN(2*PI*FC * (I%-M%/2)) / (I%-M%/2)
250    A[I%] = A[I%] * (0.42 - 0.5*COS(2*PI*I%/M%) + 0.08*COS(4*PI*I%/M%))     
260 NEXT I%
270 '
280 SUM = 0              'Normalize the first low-pass filter kernel for
290 FOR I% = 0 TO 800    'unity gain at DC
300    SUM = SUM + A[I%]
310 NEXT I%
320 '
330 FOR I% = 0 TO 800
340    A[I%] = A[I%] / SUM
350 NEXT I%
360 '                    'Calculate the second low-pass filter kernel via Eq.
370 FC = 0.204           '16-4, with a cutoff frequency of 0.204, store in B[]
380 FOR I% = 0 TO 800
390    IF (I%-M%/2) = 0 THEN B[I%] = 2*PI*FC
400    IF (I%-M%/2) <> 0  THEN B[I%] = SIN(2*PI*FC * (I%-M%/2)) / (I%-M%/2)
410    B[I%] = B[I%] * (0.42 - 0.5*COS(2*PI*I%/M%) + 0.08*COS(4*PI*I%/M%))     
420 NEXT I%
430 '
440 SUM = 0              'Normalize the second low-pass filter kernel for
450 FOR I% = 0 TO 800    'unity gain at DC
460    SUM = SUM + B[I%]
470 NEXT I%
480 '
490 FOR I% = 0 TO 800
500    B[I%] = B[I%] / SUM
510 NEXT I%
520 '               
530 FOR I% = 0 TO 800     'Change the low-pass filter kernel in B[ ] into a 
540   B[I%] = - B[I%]     'high-pass filter kernel using spectral inversion 
550 NEXT I%               '(as in Fig. 14-5)
560 B[400] =  B[400] + 1     
570 ' 
580 '
590 FOR I% = 0 TO 800         'Add the low-pass filter kernel in A[ ], to the 
600    H[I%] = A[I%] + B[I%]  'high-passfilter kernel in B[ ], to form a band
610 NEXT I%                   'reject filter kernel, stored in H[ ] (Fig.14-8)
620 '
630 FOR I% = 0 TO 800     'Change the band-reject filter kernel into a 
640    H[I%] = -H[I%]     'band-passfilter kernel by using spectral inversion
650 NEXT I%                            
660 H[400] =  H[400] + 1
670 '                     'The band-pass filter kernel now resides in H[ ]
680 END


TABLE 17-1
100 'CUSTOM FILTER DESIGN 
110 'This program converts an aliased 1024 point impulse response into an M+1
120 'point filter kernel (such as Fig. 17-1b being converted into Fig. 17-1c) 
130 '
140 DIM REX[1023]           'REX[ ] holds the signal being converted 
150 DIM T[1023]             'T[ ] is a temporary storage buffer
160 '
170 PI = 3.14159265
180 M% = 40                 'Set filter kernel length (41 total points)
190 '                                  
200 GOSUB XXXX              'Mythical sub to load REX[ ] with impulse response
210 '
220 FOR I% = 0 TO 1023      'Shift (rotate) the signal M/2 points to the right
230    INDEX% = I% + M%/2
240    IF INDEX% > 1023 THEN INDEX% = INDEX%-1024
250    T[INDEX%] = REX[I%]
260 NEXT I%
270 '
280 FOR I% = 0 TO 1023       
290    REX[I%] = T[I%]
300 NEXT I%
310 '                       'Truncate and window the signal
320 FOR I% = 0 TO 1023
330   IF I% <= M% THEN REX[I%] = REX[I%] * (0.54 - 0.46 * COS(2*PI*I%/M%))
340   IF I% > M%  THEN REX[I%] = 0
350 NEXT I%
360 '                      'The filter kernel now resides in REX[0] to REX[40]
370 END


TABLE 18-1
100 'FFT CONVOLUTION
110 'This program convolves a 10 million point signal with a 400 point filter
120 'kernel. The input signal is broken into 16000 segments, each with 625 
125 'points.  1024 point FFTs are used.
130 '
130 '                    'INITIALIZE THE ARRAYS
140 DIM XX[1023]         'the time domain signal (for the FFT)
150 DIM REX[512]         'real part of the frequency domain (for the FFT)
160 DIM IMX[512]         'imaginary part of the frequency domain (for the FFT)
170 DIM REFR[512]        'real part of the filter's frequency response
180 DIM IMFR[512]        'imaginary part of the filter's frequency response
190 DIM OLAP[398]        'holds overlapping samples from segment to segment
200 '
210 FOR I% = 0 TO 398    'zero the array holding the overlapping samples
220   OLAP[I%] = 0
230 NEXT I%
240 '
250 '                    'FIND & STORE THE FILTER'S FREQUENCY RESPONSE
260 GOSUB XXXX           'Mythical sub to load the filter kernel into XX[ ]
270 GOSUB XXXX           'Mythical FFT subroutine:  XX[ ] -->  REX[ ] & IMX[ ]
280 FOR F% = 0 TO 512    'Save the frequency response in REFR[ ] & IMFR[ ]
290   REFR[F%] = REX[F%]
300   IMFR[F%] = IMX[F%]
310 NEXT F%
320 '
330 '                    'PROCESS EACH OF THE 16000 SEGMENTS
340 FOR SEGMENT% = 0 TO 15999
350   '
360   GOSUB XXXX         'Mythical sub to load next input segment into XX[ ]
370   GOSUB XXXX         'Mythical FFT sub:  XX[ ]  -->  REX[ ] & IMX[ ]
380   '                 
390   FOR F% = 0 TO 512  'Multiply freq. spectrum by the freq. response
400     TEMP      = REX[F%]*REFR[F%] - IMX[F%]*IMFR[F%]
410     IMX[F%]  = REX[F%]*IMFR[F%] + IMX[F%]*REFR[F%]
420     REX[F%]  = TEMP
430   NEXT F%
440   '
450   GOSUB XXXX         'Mythical IFFT sub:  REX[ ] & IMX[ ]  -->  XX[ ]
460   '                         
470   FOR I% = 0 TO 398  'Add the last segment's overlap to this segment
480     XX[I%] = XX[I%] + OLAP[I%]
490   NEXT I%
500   '                
510   FOR I% = 625 TO 1023   'Save samples that will overlap the next segment
520     OLAP[I%-625] = XX[I%]
530   NEXT I%
540   '
550   GOSUB XXXX         'Mythical sub to output the 625 samples stored 
560   '                  'in XX[0] to XX[624]
570   '
580 NEXT SEGMENT%
590 '
600 GOSUB XXXX           'Mythical sub to output all 399 samples in OLAP[ ]
610 END


TABLE 19-1
100 'RECURSIVE FILTER
110 '
120 DIM X[499]           'holds the input signal
130 DIM Y[499]           'holds the filtered output signal
140 '
150 GOSUB XXXX           'Mythical subroutine to calculate the recursion 
160 '                    'coefficients: A0, A1, A2, B1, B2
170 '
180 GOSUB XXXX           'Mythical subroutine to load X[ ] with the input data
190 '
200 FOR I% = 2 TO 499
210    Y[I%] = A0*X[I%] + A1*X[I%-1] + A2*X[I%-2] + B1*Y[I%-1] + B2*Y[I%-2]
220 NEXT I%
230 '
240 END


TABLE 20-4
100 'CHEBYSHEV FILTER-  RECURSION 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]