programs.txt

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

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[ ] & B[ ]
690 END


TABLE 20-5
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 24-1
100 CONVENTIONAL IMAGE CONVOLUTION
110 '
120 DIM X[99,99]         'holds the input image,   100 100 pixels
130 DIM H[28,28]         'holds the filter kernel,   29 29 pixels 
140 DIM Y[127,127]       'holds the output image, 128 128 pixels
150 '
160 FOR R% = 0 TO 127    'loop through each row and column in the output
170 FOR C% = 0 TO 127    'image calculating the pixel value via Eq. 24-3
180 '
190   Y[R%,C%] = 0       'zero the pixel so it can be used as an accumulator
200 '
210   FOR J% = 0 TO 28   'multiply each pixel in the kernel by the 
220   FOR K% = 0 TO 28   'corresponding pixel in the input image, and add to
225                      'the accumulator
230      Y[R%,C%] = Y[R%,C%] + H[J%,K%] * X[R%-J%,C%-J%]
240   NEXT K%
250   NEXT J%
260 '
270 NEXT C%
280 NEXT R%
290 '   
300 END


TABLE 25-1
100 'SKELETONIZATION PROGRAM
110 'Object pixels have a value of 0 (displayed as black)
120 'Background pixels have a value of 255 (displayed as white)
130 '
140 DIM X%[149,149]        'X%[  ,  ] holds the image being processed
150 '
160 GOSUB XXXX             'Mythical subroutine to load X%[ , ]
170 '
180 FOR ITER% = 0 TO 5     'Run through six iteration loops
190 '
200   FOR R% = 1 TO 148    'Loop through each pixel in the image.
210   FOR C% = 1 TO 148    'Subroutine 5000 (Table 25-2) indicates which 
220     GOSUB 5000         'pixels can be changed from black to white,
230   NEXT C%              'by marking the pixels with a value of 1.
240   NEXT R%
250 '
260   FOR R% = 0 TO 149    'Loop through each pixel in the image changing
270   FOR C% = 0 TO 149    'the marked pixels from black to white.
280     IF X%(R%,C%) = 1 THEN X%(R%,C%) = 255
290   NEXT C%
300   NEXT R%
310 '
320 NEXT ITER%
330 '
340 END


TABLE 25-2
5000 'Subroutine to determine if the pixel at X%[R%,C%] can be removed.
5010 'If all four of the rules are satisfied, then X%(R%,C%], is set to a 
5020 'value of 1, indicating it should be removed at the end of the iteration.
5030 '
5040 'RULE #1:  Do nothing if the pixel already white
5050 IF X%(R%,C%) = 255 THEN RETURN
5060 '
5070 '
5080 'RULE #2:  Do nothing if all of the close neighbors are black
5090 IF  X%[R% -1,C%  ] <> 255 AND X%[R%  ,C%+1] <> 255 AND 
     X%[R%+1,C%  ] <> 255 AND X%[R%  ,C% -1] <> 255  THEN RETURN
5100 '
5110 '
5120 'RULE #3:  Do nothing if only a single neighbor pixel is black
5130 COUNT% = 0
5140 IF X%[R%-1,C%-1]  = 0 THEN COUNT% = COUNT% + 1
5150 IF X%[R%-1,C%  ]  = 0 THEN COUNT% = COUNT% + 1
5160 IF X%[R%-1,C%+1]  = 0 THEN COUNT% = COUNT% + 1
5170 IF X%[R%  ,C%+1]  = 0 THEN COUNT% = COUNT% + 1
5180 IF X%[R%+1,C%+1]  = 0 THEN COUNT% = COUNT% + 1
5190 IF X%[R%+1,C%  ]  = 0 THEN COUNT% = COUNT% + 1
5200 IF X%[R%+1,C%-1]  = 0 THEN COUNT% = COUNT% + 1
5210 IF X%[R%  ,C%-1]  = 0 THEN COUNT% = COUNT% + 1
5220 IF COUNT% = 1 THEN RETURN
5230 '
5240 '
5250 'RULE 4:  Do nothing if the neighbors are unconnected. 
5260 'Determine this by counting the black-to-white transitions
5270 'while moving clockwise through the 8 neighboring pixels.
5280 COUNT% = 0
5290 IF X%[R%-1,C%-1]  = 0 AND X%[R%-1,C%  ]  > 0  THEN COUNT% = COUNT% +
1
5300 IF X%[R%-1,C%  ]  = 0 AND X%[R%-1,C%+1]  > 0  AND X%[R%  ,C%+1] > 0  
     THEN COUNT% = COUNT% + 1
5310 IF X%[R%-1,C%+1]  = 0 AND X%[R%  ,C%+1]  > 0  THEN COUNT% = COUNT%
+ 1
5320 IF X%[R%  ,C%+1]  = 0 AND X%[R%+1,C%+1]  > 0  AND X%[R%+1,C%  ] > 0  
     THEN COUNT% = COUNT% + 1
5330 IF X%[R%+1,C%+1]  = 0 AND X%[R%+1,C%  ]  > 0  THEN COUNT% =
COUNT% + 1
5340 IF X%[R%+1,C%  ]  = 0 AND X%[R%+1,C%-1]  > 0  AND X%[R%  ,C%-1] > 0  
     THEN COUNT% = COUNT% + 1
5350 IF X%[R%+1,C%-1]  = 0 AND X%[R%  ,C%-1]  > 0  THEN COUNT% = COUNT%
+ 1
5360 IF X%[R%  ,C%-1]  = 0 AND X%[R%-1,C%-1]  > 0  AND X%[R%-1,C%  ] > 0  
     THEN COUNT% = COUNT% + 1
5370 IF COUNT% > 1 THEN RETURN
5380 '
5390 '
5400 'If all rules are satisfied, mark the pixel to be set to white at the 
5410 X%(R%,C%) = 1              'end of the iteration
5420 '
5430 RETURN


TABLE 26-1
100 'NEURAL NETWORK (FOR THE FLOW DIAGRAM IN FIG. 26-5)
110 '
120 DIM X1[15]      'holds the input values
130 DIM X2[4]       'holds the values exiting the hidden layer
140 DIM X3[2]       'holds the values exiting the output layer
150 DIM WH[4,15]    'holds the hidden layer weights
160 DIM WO[2,4]     'holds the output layer weights
170 '
180 GOSUB XXXX      'mythical sub to load X1[ ] with the input data
190 GOSUB XXXX      'mythical sub to load the weights, WH[ , ] & W0[ , ]
200 '
210 '                                     'FIND THE HIDDEN NODE VALUES, X2[  ]
220 FOR J% = 1 TO 4                       'loop for each hidden layer node
230   ACC = 0                             'clear accumulator variable, ACC
240   FOR I% = 1 TO 15                    'weight and sum each input node  
250     ACC = ACC + X1[I%] * WH[J%,I%]
260   NEXT I%
270   X2[J%] = 1 / (1 + EXP(-ACC) )       'pass value through the sigmoid
280 NEXT J%
290 '
300 '                                     'FIND THE OUTPUT NODE VALUES, X3[]
310 FOR J% = 1 TO 2                       'loop for each output layer node
320   ACC = 0                             'clear accumulator variable, ACC
330   FOR I% = 1 TO 4                     'weight and sum each hidden node  
340     ACC = ACC + X2[I%] * WO[J%,I%]
350   NEXT I%
360   X3[J%] = 1 / (1 + EXP(-ACC) )       'pass value through the sigmoid
370 NEXT J%
380 '
390 END


TABLE 26-2
100 'NEURAL NETWORK TRAINING (Determination of weights)
110 '
120                                 'INITIALIZE
130 MU = .000005                    'iteration step size
140 DIM X1[101]                     'holds the input layer signal + bias term
150 DIM X2[10]                      'holds the hidden layer signal
160 DIM WH[10,101]                  'holds hidden layer weights
170 DIM WO[10]                      'holds output layer weights
180 '
190 FOR H% = 1 TO 10                'SET WEIGHTS TO RANDOM VALUES
200   WO[H%] = (RND-0.5)            'output layer weights: -0.5 to 0.5
210   FOR I% = 1 TO 101             'hidden layer weights: -0.0005 to 0.0005
220     WH[H%,I%] = (RND-0.5)/1000
230   NEXT I%
240 NEXT H%
250 '
260 '                               'ITERATION LOOP
270 FOR ITER% = 1 TO 800            'loop for 800 iterations
280 '
290   ESUM = 0                      'clear the error accumulator, ESUM
300   '
310   FOR LETTER% = 1 TO 260        'loop for each letter in the training set
320     GOSUB 1000                  'load X1[ ] with training set
330     GOSUB 2000                  'find the error for this letter, ELET
340     ESUM = ESUM + ELET^2        'accumulate error for this iteration
350     GOSUB 3000                  'find the new weights
360   NEXT LETTER%
370   '
380   PRINT  ITER%  ESUM            'print the progress to the video screen
390 '
400 NEXT ITER%
410 '
420 GOSUB XXXX                      'mythical subroutine to save the weights
430 END


TABLE 26-3
1000 'SUBROUTINE TO LOAD X1[ ] WITH IMAGES FROM THE DATABASE
1010 'Variables entering routine: LETTER%
1020 'Variables exiting routine: X1[1] to X1[100],  X1[101] = 1, CORRECT
1030 '
1040 'The variable, LETTER%, between 1 and 260, indicates which image in the
1045 'database is to be 
1050 'returned in X1[1] to X1[100].  The bias node, X1[101], always has a
1055 'value of one.  The variable, 
1060 'CORRECT, has a value of one if the image being returned is a vowel, and
1065 'zero otherwise.
1070 '(The details of this subroutine are unimportant, and not listed here).
1900 RETURN

2000 'SUBROUTINE TO CALCULATE THE ERROR WITH THE CURRENT WEIGHTS
2010 'Variables entering routine:  X1[ ], X2[ ], WI[ , ], WH[ ], CORRECT
2020 'Variables exiting routine:   ELET
2030 '
2040 '                                     'FIND THE HIDDEN NODE VALUES, X2[ ]
2050 FOR H% = 1 TO 10                      'loop for each hidden nodes
2060   ACC = 0                             'clear the accumulator
2070   FOR I% = 1 TO 101                   'weight and sum each input node  
2080     ACC = ACC + X1[I%] * WH[H%,I%]
2090   NEXT I%
2100   X2[H%] = 1 / (1 + EXP(-ACC) )       'pass summed value through sigmoid
2110 NEXT H%
2120 '
2130 '                                     'FIND THE OUTPUT VALUE: X3
2140 ACC = 0                               'clear the accumulator
2150 FOR H% = 1 TO 10                      'weight and sum each hidden node
2160   ACC = ACC + X2[H%] * WO[H%]
2170 NEXT H%
2180 X3 = 1 / (1 + EXP(-ACC) )             'pass summed value through sigmoid
2190 '
2200 '                                     'FIND ERROR FOR THIS LETTER, ELET
2210 ELET = (CORRECT - X3)                 'find the error
2220 IF CORRECT = 1 THEN ELET = ELET*5     'give extra weight to targets
2230 '
2240 RETURN

3000 'SUBROUTINE TO FIND NEW WEIGHTS                                    
3010 'Variables entering routine:  X1[ ], X2[ ], X3, WI[ , ], WH[ ], ELET, MU
3020 'Variables exiting routine:   WI[ , ], WH[ ]
3030 '
3040 '                                      'FIND NEW WEIGHTS FOR INPUT NODES
3050 FOR H% = 1 TO 10
3060   FOR I% = 1 TO 101
3070     SLOPEO = X3 * (1 - X3)
3080     SLOPEH = X2(H%) * (1 - X2[H%])
3090     DX3DW = X1[I%] * SLOPEH * WO[H%] * SLOPEO
3100     WH[H%,I%] = WH[H%,I%] + DX3DW * ELET * MU
3110   NEXT I%
3120 NEXT H%
3130 '
3140 '                                     'FIND NEW WEIGHTS FOR HIDDEN NODES
3150 FOR H% = 1 TO 10
3160   SLOPEO = X3 * (1 - X3)
3170   DX3DW = X2[H%] * SLOPEO
3180   WO[H%] = WO[H%] + DX3DW * ELET * MU
3190 NEXT H%
3200 '
3210 RETURN


TABLE 26-4
100 'ITERATIVE DESIGN OF RECURSIVE FILTER
110 '
120                                'INITIALIZE
130 N%  = 256                      'number of points in FFT
140 NP% = 8                        'number of poles in filter
150 DELTA = .00001                 'perturbation increment
160 MU = .2                        'iteration step size
170 DIM REX[255]                   'real part of signal during FFT
180 DIM IMX[255]                   'imaginary part of signal during FFT
190 DIM T[128]                     'desired frequency response (mag only)
200 DIM A[8]                       'the "a" recursion coefficients
210 DIM B[8]                       'the "b" recursion coefficients
220 DIM SA[8]                      'slope for "a" coefficients
230 DIM SB[8]                      'slope for "b" coefficients
240 '
250 GOSUB XXXX                     'mythical subroutine to load T[ ]
260 '
270 FOR P% = 0 TO NP%              'initialize coefficients to identity system
280   A[P%] = 0  
290   B[P%] = 0
300 NEXT P%
310 A[0] = 1
320 '
330 '                                  'ITERATION LOOP
340 FOR ITER% = 1 TO 100               'loop for desired number of iterations
350   GOSUB 2000                       'calculate new coefficients
360   PRINT  ITER%  ENEW  MU           'print current status to video screen
370   IF ENEW > EOLD THEN MU = MU/2    'adjust the value of MU
380 NEXT ITER%
390 '                                                            
400 ' 
410 FOR P% = 0 TO NP%    'PRINT OUT THE COEFFICIENTS
420   PRINT  A[P%]    B[P%]
430 NEXT P%
440 '
450 END


TABLE 26-5
2000 'SUBROUTINE TO CALCULATE THE NEW RECURSION COEFFICIENTS
2010 'Variables entering routine:  A[ ], B[ ], DELTA, MU
2020 'Variables exiting routine:   A[ ], B[ ], EOLD, ENEW
2030 '
2040 GOSUB 3000                     'FIND THE CURRENT ERROR
2050 EOLD = ER                      'store current error in variable, EOLD
2060 '
2070                                'FIND THE ERROR SLOPES 
2080 FOR P% = 0 TO NP%              'loop through each "a" coefficient
2090   A[P%] = A[P%] + DELTA        'add a small increment to the coefficient
2100   GOSUB 3000                   'find the error with the change
2110   SA[P%] = (ER-EOLD)/DELTA     'calculate the error slope, store in SA[ ]
2120   A[P%] = A[P%] - DELTA        'return coefficient to original value
2130 NEXT P%
2140 '
2150 FOR P% = 1 TO NP%              'repeat process for each "b" coefficient
2160   B[P%] = B[P%] + DELTA                 
2170   GOSUB 3000
2180   SB[P%] = (ER-EOLD)/DELTA     'calculate the error slope, store in SB[ ]
2190   B[P%] = B[P%] - DELTA
2200 NEXT P%
2210 '                                 'CALCULATE NEW COEFFICIENTS
2220 FOR P% = 0 TO NP%                 'loop through each coefficient
2230   A[P%] = A[P%] - SA[P%] * MU     'adjust coefficients to move "downhill"
2240   B[P%] = B[P%] - SB[P%] * MU
2250 NEXT P%
2260 '
2270 GOSUB 3000                        'FIND THE NEW ERROR
2280 ENEW = ER                         'store new error in variable, ENEW
2290 '
2300 RETURN

3000 'SUBROUTINE TO CALCULATE THE FREQUENCY DOMAIN ERROR
3010 'Variables entering routine:  A[ ], B[ ], T[ ]
3020 'Variables exiting routine:   ER   
3030 '
3040 FOR I% = 0 TO N%-1                'LOAD SHIFTED  IMPULSE INTO IMX[ ]
3050  REX[I%] = 0
3060  IMX[I%] = 0
3070 NEXT I%
3080 IMX[12] = 1
3090 '                                 'CALCULATE IMPULSE RESPONSE
3100 FOR I% = 12 TO N%-1
3110  FOR J% = 0 TO NP%
3120    REX[I%] = REX[I%]  +  A[J%] * IMX[I%-J%]  +  B[J%] * REX[I%-J%]
3130  NEXT J%
3140 NEXT I%
3150 IMX[12] = 0
3160 '                                 'CALCULATE THE FFT
3170 GOSUB 1000                        'Table 12-4, uses REX[ ], IMX[ ], N%
3180 '
3190                                     'FIND FREQUENCY DOMAIN ERROR
3200 ER = 0                              'zero ER, to use as an accumulator
3210 FOR I% = 0 TO N%/2                  'loop through each positive frequency 
3220   MAG = SQR(REX[I%]^2 + IMX[I%]^2)  'rectangular --> polar conversion
3230   ER = ER + ( MAG - T[I%] )^2       'calculate & accumulate squared error
3240 NEXT I%
3250 ER = SQR( ER/(N%/2+1) )             'finish calculation of error, ER 
3260 '
3270 RETURN