programs.txt

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

Programs from: "The Scientist and Engineer's Guide to Digital Signal Processing.

TABLE 2-1
100 CALCULATION OF THE MEAN AND STANDARD DEVIATION
110 '
120 DIM X[511]           'The signal is held in X[0] to X[511]
130 N% = 512             'N% is the number of points in the signal
140 '
150 GOSUB XXXX           'Mythical subroutine that loads the signal into X[ ]
160 '
170 MEAN = 0             'Find the mean via Eq. 2-1
180 FOR I% = 0 TO N%-1
190   MEAN = MEAN + X[I%]
200 NEXT I%
210 MEAN = MEAN/N%
220 '
230 VARIANCE = 0         'Find the standard deviation via Eq. 2-2
240 FOR I% = 0 TO N%-1
250   VARIANCE = VARIANCE + ( X[I%] - MEAN )^2
260 NEXT I%
270 VARIANCE = VARIANCE/(N%-1)
280 SD = SQR(VARIANCE)
290 '
300 PRINT  MEAN  SD      'Print the calculated mean and standard deviation
310 '
320 END


TABLE 2-2
100 'MEAN AND STANDARD DEVIATION USING RUNNING STATISTICS
110 '
120 DIM X[511]           'The signal is held in X[0] to X[511]
130 '
140 GOSUB XXXX           'Mythical subroutine that loads the signal into X[ ]
150 '
160 N% = 0               'Zero the three running parameters 
170 SUM = 0
180 SUMSQUARES = 0
190 '
200 FOR I% = 0 TO 511    'Loop through each sample in the signal
210   '
220   N% = N%+1          'Update the three parameters 
230   SUM = SUM + X(I%)
240   SUMSQUARES = SUMSQUARES + X(I%)^2
250   '
260   MEAN = SUM/N% 'Calculate mean and standard deviation via Eq. 2-3
270   VARIANCE = (SUMSQUARES -  SUM^2/N%) / (N%-1)
280   SD = SQR(VARIANCE)
290   '
300   PRINT  MEAN   SD   'Print the running mean and standard deviation
310   '
320 NEXT I%
330 '     
340 END


TABLE 2-3
100 'CALCULATION OF THE HISTOGRAM, MEAN, AND STANDARD DEVIATION
110 '
120 DIM X%[25000]        'X%[0] to X%[25000] holds the signal being processed
130 DIM H%[255]          'H%[0] to H%[255] holds the histogram
140 N% = 25001           'Set the number of points in the signal
150 '
160 FOR I% = 0 TO 255    'Zero, so it can be used as an accumulator
170   H%[I%] = 0
180 NEXT I%
190 '
200 GOSUB XXXX           'Mythical subroutine that loads the signal into X%[ ]
210 '
220 FOR I% = 0 TO 25000  'Calculate the histogram for 25001 points
230   H%[ X%[I%] ] = H%[ X%[I%] ] + 1   
240 NEXT I%
250 '
260 MEAN = 0             'Calculate the mean via Eq. 2-6
270 FOR I% = 0 TO 255
280   MEAN = MEAN + I% * H%[I%]
290 NEXT I%
300 MEAN = MEAN / N%
310 '
320 VARIANCE = 0         'Calculate the standard deviation via Eq. 2-7
330 FOR I% = 0 TO 255
340   VARIANCE = VARIANCE + H[I%] * (I%-MEAN)^2
350 NEXT I%
360 VARIANCE = VARIANCE / (N%-1)
370 SD = SQR(VARIANCE)
380 '
390 PRINT  MEAN  SD      'Print the calculated mean and standard deviation.
400 '
410 END


TABLE 2-4
100 'CALCULATION OF BINNED HISTOGRAM
110 '
120 DIM X[25000]         'X[0] to X[25000] holds the floating point signal, 
130 '                    'with each sample being in the range: 0.0 to 10.0 
140 DIM H%[999]          'H%[0] to H%[999] holds the binned histogram
150 '
160 FOR I% = 0 TO 999    'Zero the binned histogram for use as an accumulator
170   H%[I%] = 0
180 NEXT I%
190 '
200 GOSUB XXXX           'Mythical subroutine that loads the signal into X%[ ]
210 '
220 FOR I% = 0 TO 25000  'Calculate the binned histogram for 25001 points
230   BINNUM% = INT( X[I%] * .01 )
240   H%[ BINNUM%] = H%[ BINNUM%] + 1   
250 NEXT I%
260 '
270 END


TABLE 6-1
100  'CONVOLUTION USING THE INPUT SIDE ALGORITHM
110  '
120 DIM X[80]            'The input signal, 81 points
130 DIM H[30]            'The impulse response, 31 points
140 DIM Y[110]           'The output signal, 111 points
150  '
160 GOSUB XXXX           'Mythical subroutine to load X[ ] and H[ ] 
170  '
180 FOR I% = 0 TO 110    'Zero the output array
190   Y(I%) = 0
200 NEXT I%
210  '
220 FOR I% = 0 TO 80     'Loop for each point in X[ ]
230   FOR J% = 0 TO 30   'Loop for each point in H[ ]
240     Y[I%+J%] = Y[I%+J%] + X[I%]*H[J%]
250   NEXT J%
260 NEXT I%
270  '
280 GOSUB XXXX           'Mythical subroutine to store Y[ ]
290  '
300 END


TABLE 6-2
100  'CONVOLUTION USING THE OUTPUT SIDE ALGORITHM
110  '
120 DIM X[80]            'The input signal, 81 points
130 DIM H[30]            'The impulse response, 31 points
140 DIM Y[110]           'The output signal, 111 points
150  '
160 GOSUB XXXX           'Mythical subroutine to load X[ ] and H[ ] 
170  '                             
180 FOR I% = 0 TO 110    'Loop for each point in Y[ ]
190   Y[I%] = 0          'Zero the sample in the output array 
200   FOR J% = 0 TO 30   'Loop for each point in H[ ]
210     IF (I%-J% < 0)     THEN GOTO 240
220     IF (I%-J% > 80)    THEN GOTO 240
230     Y(I%) = Y(I%) + H(J%) * X(I%-J%)
240   NEXT J%
250 NEXT I%
260  '
270 GOSUB XXXX           'Mythical subroutine to store Y[ ]
280  '
290 END


TABLE 8-1
100 'THE INVERSE DISCRETE FOURIER TRANSFORM  
110 'The time domain signal, held in XX[ ], is calculated from the frequency 
120 'domain signals, held in REX[ ] and IMX[ ].    
130 '
140 DIM XX[511]       'XX[ ]  holds the time domain signal
150 DIM REX[256]      'REX[ ] holds the real part of the frequency domain
160 DIM IMX[256]      'IMX[ ] holds the imaginary part of the frequency domain
170 '
180 PI = 3.14159265   'Set the constant, PI
190 N% = 512          'N% is the number of points in XX[ ]
200 '
210 GOSUB XXXX        'Mythical subroutine to load data into REX[ ] & IMX[ ] 
220 '
230 
240 '                 'Find the cosine and sine wave amplitudes using Eq. 8-3
250 FOR K% = 0 TO 256           
260   REX[K%] =  REX[K%] / (N%/2)    
270   IMX[K%]  = -IMX[K%] / (N%/2)  
280 NEXT K%                               
290 '
300 REX[0] = REX[0] / 2
310 REX[256] = REX[256] / 2
320 '
330 '
340 FOR I% = 0 TO 511  'Zero XX[ ] so it can be used as an accumulator
350   XX[I%] = 0
360 NEXT I%
370 '
380 '              Eq. 8-2  SYNTHESIS METHOD #1.  Loop through each 
390 '              frequency generating the entire length of the sine & cosine 
400 '              waves, and add them to the accumulator signal, XX[ ] 
410 '
420 FOR K% = 0 TO 256      'K% loops through each sample in REX[ ] and IMX[ ]
430   FOR I% = 0 TO 511    'I% loops through each sample in XX[ ]
440     '             
450     XX[I%] = XX[I%] + REX[K%] * COS(2*PI*K%*I%/N%)
460     XX[I%] = XX[I%] + IMX[K%] * SIN(2*PI*K%*I%/N%)
470     '
480   NEXT I%
490 NEXT K%
500 '
510 END 

ALTERNATE CODE FOR LINES 380 to 510
380 'Eq. 8-2  SYNTHESIS METHOD #2.  Loop through each 
390 'sample in the time domain, and sum the corresponding 
400 'samples from each cosine and sine wave 
410 '
420 FOR I% = 0 TO 511       'I% loops through each sample in XX[ ]
430   FOR K% = 0 TO 256     'K% loops through each sample in REX[ ] and IMX[ ]
440     '             
450     XX[I%] = XX[I%] + REX[K%] * COS(2*PI*K%*I%/N%)
460     XX[I%] = XX[I%] + IMX[K%] * SIN(2*PI*K%*I%/N%)
470     '
480   NEXT K%
490 NEXT I%
500 '
510 END


TABLE 8-2
100 'THE DISCRETE FOURIER TRANSFORM 
110 'The frequency domain signals, held in REX[ ] and IMX[ ], are 
120 'calculated from the time domain signal, held in XX[ ].
130 '
140 DIM XX[511]       'XX[ ]  holds the time domain signal
150 DIM REX[256]      'REX[ ] holds the real part of the frequency domain
160 DIM IMX[256]      'IMX[ ] holds the imaginary part of the frequency domain
170 '
180 PI = 3.14159265   'Set the constant, PI
190 N% = 512          'N% is the number of points in XX[ ]
200 '
210 GOSUB XXXX        'Mythical subroutine to load data into XX[ ]      
220 '
230 '
240 FOR K% = 0 TO 256 'Zero REX[ ] & IMX[ ] so they can be accumulators 
250   REX[K%] = 0
260   IMX[K%] = 0
270 NEXT K%
280 '
290 '              'Correlate XX[ ] with the cosine and sine waves, Eq. 8-4
300 '
310 FOR K% = 0 TO 256      'K% loops through each sample in REX[ ] and IMX[ ]
320   FOR I% = 0 TO 511    'I% loops through each sample in XX[ ]
330     '
340     REX[K%] = REX[K%] + X[I%] * COS(2*PI*K%*I%/N%)
350     IMX[K%] = IMX[K%]  -  X[I%] * SIN(2*PI*K%*I%/N%)
360     '
370   NEXT I%
380 NEXT K%
390 '
400 END


TABLE 8-3
100 'RECTANGULAR-TO-POLAR & POLAR-TO-RECTANGULAR CONVERSION
110 '
120 DIM REX[256]    'REX[ ]       holds the real part
130 DIM IMX[256]    'IMX[ ]       holds the imaginary part
140 DIM MAG[256]    'MAG[ ]       holds the magnitude
150 DIM PHASE[256]  'PHASE[ ] holds the phase
160 '
170 PI = 3.14159265
180 '
190 GOSUB XXXX      'Mythical subroutine to load data into REX[ ] and IMX[ ]
200 '
210 '
220 '               'Rectangular-to-polar conversion, Eq. 8-6
230 FOR K% = 0 TO 256
240   MAG[K%] = SQR( REX[K%]^2 + IMX[K%]^2 )     'from Eq. 8-6
250   IF REX[K%] = 0 THEN REX[K%] = 1E-20        'prevent divide by 0 
260   PHASE[K%] = ATN( IMX[K%] / REX[K%] )       'from Eq. 8-6  
270   '                                          'correct the arctan   
280   IF REX[K%] < 0 AND IMX[K%] <   0 THEN PHASE[K%] = PHASE[K%]  - PI
290   IF REX[K%] < 0 AND IMX[K%] >= 0 THEN PHASE[K%] = PHASE[K%] + PI
300 NEXT K%
310 '         
320 '    
330 '               'Polar-to-rectangular conversion, Eq. 8-7
340 FOR K% = 0 TO 256
350   REX[K%] = MAG[K%] * COS( PHASE[K%] )
360   IMX[K%] = MAG[K%] *  SIN( PHASE[K%] )
370 NEXT K%
380 ' 
390 END

 
TABLE 8-4     
100 ' PHASE UNWRAPPING 
110 '
120 DIM PHASE[256]      'PHASE[ ]   holds the original phase
130 DIM UWPHASE[256]    'UWPHASE[ ] holds the unwrapped phase
140 '
150 PI = 3.14159265
160 '
170 GOSUB XXXX          'Mythical subroutine to load data into PHASE[ ]
180 '
190 UWPHASE[0] = 0      'The first point of all phase signals is zero    
200 '
210 '                   'Go through the unwrapping algorithm
220 FOR K% = 1 TO 256
230   C% = CINT( (UWPHASE[K%-1] - PHASE[K%]) / (2 * PI) )
240   UWPHASE[K%] = PHASE[K%] + C%*2*PI
250 NEXT K%
260 '
270 END


TABLE 12-1
6000 'NEGATIVE FREQUENCY GENERATION
6010 'This subroutine creates the complex frequency domain from the real 
6020 'frequency domain. Upon entry to this subroutine, N% contains the number
6030 'of points in the signals, and REX[ ] and IMX[ ] contain the real 
6040 'frequency domain in samples 0 to N%/2. On return, REX[ ] and IMX[ ] 
6045 'contain the complex frequency domain in samples 0 to N%-1.
6050 '
6060 FOR K% = (N%/2+1) TO (N%-1)
6070   REX[K%] =  REX[N%-K%]
6080   IMX[K%] = -IMX[N%-K%]
6090 NEXT K%
6100 '
6110 RETURN


TABLE 12-2
5000 'COMPLEX DFT BY CORRELATION
5010 'Upon entry, N% contains the number of points in the DFT, and
5020 'XR[ ] and XI[ ] contain the real and imaginary parts of the time domain.
5030 'Upon return, REX[ ] and IMX[ ] contain the frequency domain data.  
5040 'All signals run from 0 to N%-1.
5050 '
5060 PI = 3.14159265             'Set constants
5070 '
5080 FOR K% = 0 TO N%-1          'Zero REX[ ] and IMX[ ], so they can be used
5090   REX[K%] = 0               'as accumulators during the correlation
5100   IMX[K%] = 0
5110 NEXT K%
5120 '
5130 FOR K% = 0 TO N%-1          'Loop for each value in frequency domain
5140   FOR I% = 0 TO N%-1        'Correlate with the complex sinusoid, SR & SI
5150 '
5160     SR =  COS(2*PI*K%*I%/N%)      'Calculate complex sinusoid
5170     SI  = -SIN(2*PI*K%*I%/N%)
5180     REX[K%] = REX[K%] + XR[I%]*SR - XI[I%]*SI
5190     IMX[K%] = IMX[K%] + XR[I%]*SI + XI[I%]*SR
5200 '
5210   NEXT I%
5220 NEXT K%
5230 '
5240 RETURN


TABLE 12-3
    SUBROUTINE  FFT(X,M)
    COMPLEX X(4096),U,S,T
    PI=3.14159265
    N=2**M
    DO 20 L=1,M
    LE=2**(M+1-L)
    LE2=LE/2
    U=(1.0,0.0)
    S=CMPLX(COS(PI/FLOAT(LE2)),-SIN(PI/FLOAT(LE2)))
    DO 20 J=1,LE2
    DO 10 I=J,N,LE
    IP=I+LE2
    T=X(I)+X(IP)
    X(IP)=(X(I)-X(IP))*U
10  X(I)=T
20  U=U*S
    ND2=N/2
    NM1=N-1
    J=1
    DO 50 I=1,NM1
    IF(I.GE.J) GO TO 30
    T=X(J)
    X(J)=X(I)
    X(I)=T
30  K=ND2
40  IF(K.GE.J) GO TO 50
    J=J-K
    K=K/2
    GO TO 40
50  J=J+K
    RETURN
    END


TABLE 12-4
1000 'THE FAST FOURIER TRANSFORM
1010 'Upon entry, N% contains the number of points in the DFT, REX[ ] and
1020 'IMX[ ] contain the real and imaginary parts of the input.  Upon return,
1030 'REX[ ] & IMX[ ] contain the DFT output. All signals run from 0 to N%-1.
1040 '
1050 PI = 3.14159265                    'Set constants
1060 NM1% = N%-1
1070 ND2% = N%/2
1080 M% = CINT(LOG(N%)/LOG(2))
1090 J% = ND2%
1100 '
1110 FOR I% = 1 TO N%-2                 'Bit reversal sorting
1120   IF I% >= J% THEN GOTO 1190
1130   TR = REX[J%]
1140   TI = IMX[J%]
1150   REX[J%] = REX[I%]
1160   IMX[J%] = IMX[I%]
1170   REX[I%] = TR
1180   IMX[I%] = TI
1190   K% = ND2%
1200   IF K% > J% THEN GOTO 1240
1210   J% = J%-K%
1220   K% = K%/2
1230   GOTO 1200
1240   J% = J%+K%
1250 NEXT I%
1260 '
1270 FOR L% = 1 TO M%                   'Loop for each stage
1280   LE% = CINT(2^L%)
1290   LE2% = LE%/2
1300   UR = 1
1310   UI = 0
1320   SR =  COS(PI/LE2%)               'Calculate sine & cosine values
1330   SI  = -SIN(PI/LE2%)
1340   FOR J% = 1 TO LE2%               'Loop for each sub DFT