C C This subroutine thakes the data from the integer input array C that is 1024 locations long and puts this data into the real C array that will be transformed. These 1024 points are forward C transformed into 1024 complex pairs. C C Next these 1024 complex pairs are used to calculate 1024 magnitude C values. These magnitude values are then normalized to 1.0 C This final result is stored at the odd index locations in the C real array that is 2048 locations long. C C C Created 4-AUG-1997 Daniel Edmunds C SUBROUTINE FORWARD_FOURIER C ***************************************** C Define all Variables and Arrays C ***************************************** IMPLICIT NONE REAL*8 CDB(1:2048) INTEGER*2 ICDB(1:1024) COMMON ICDB, CDB INTEGER*2 I, INX, INXI, INXR, IStep, J, M, MMax, N REAL*8 CDBMax, TempR, TempI, Theta, SinTh, WI, WR, WStpR, WStpI CHARACTER FuncName*32 C Put the fixed point discrete time sample ADC data array C into the floating point array that will be transformed. DO 20 INX=1,1024,1 INXR=(INX*2)-1 INXI=INX*2 CDB(INXR) = DFloat(ICDB(INX)) CDB(INXI) = 0.0 20 CONTINUE C Now we want to center this input data about zero before C multiplying it by a windowing function. Note that we are C only working on the entries in the array into which input C data was just loaded, i.e. the odd entries. TempR = 0.0 DO 25 INX=1,2048,2 TempR = TempR + CDB(INX) 25 CONTINUE TempR = TempR / 1024.0 DO 26 INX=1,2048,2 CDB(INX) = CDB(INX) - TempR 26 CONTINUE C C We now have a chance to apply a "Windowing Function" to the C input data before we transform it to frequency space. C We can select from the following Window Functions: C Rectangular i.e. no window function C Hamming 0.08 + 0.46( 1 - cos 2pi n / N ) C Blackman-Harris 0.423 - 0.497 cos 2pi n/N + 0.079 cos 4pi n/N C Rosenfeld 0.724 - 0.950 cos 2pi n/N + 0.226 cos 4pi n/N C 500 WRITE ( 6, 505 ) 505 FORMAT ( // & ' A Windowing Function may be applied to the input ', / & ' data before the transform to frequency space is ', / & ' performed. Pick from the following window functions: ', / & ' Rectangular = rec Hamming = ham ', / & ' Rosenfeld = rose Blackman-Harris = black ', / & ' enter: rec, ham, rose, black : ', $ ) READ ( 5, 510 ) FuncName 510 FORMAT ( A ) IF ( FuncName .EQ. 'rec ' ) THEN WRITE ( 6, 550 ) 550 FORMAT ( // ' You selected a Rectangular Window ', // ) GOTO 700 ELSE IF ( FuncName .EQ. 'ham' ) THEN WRITE ( 6, 560 ) 560 FORMAT ( // ' You selected a Hamming Window ', // ) C Apply the Hamming window to the input data. C This will be a 0.46 * cosine on an 8% pedestal. TempR = 2.0 * 3.1415926 / 1024.0 DO 600 INX=1,1024,1 INXR=(INX*2)-1 TempI = 0.46 * ( 1.0 - COS( TempR * DFloat(INX) )) CDB(INXR) = CDB(INXR) * ( 0.08 + TempI ) 600 CONTINUE GOTO 700 ELSE IF ( FuncName .EQ. 'black' ) THEN WRITE ( 6, 570 ) 570 FORMAT ( // ' You selected a Blackman-Harris Window ', // ) C Apply a Blackman-Harris window to the input data. C This will be 3 terms with coefficients +0.42323, -0.49755, C +0.07922 TempR = 2.0 * 3.1415926 / 1024.0 DO 610 INX=1,1024,1 INXR=(INX*2)-1 TempI = 0.42323 TempI = TempI - (0.49755 * COS( 1.0 * TempR * DFloat(INX) )) TempI = TempI + (0.07922 * COS( 2.0 * TempR * DFloat(INX) )) CDB(INXR) = CDB(INXR) * TempI 610 CONTINUE GOTO 700 ELSE IF ( FuncName .EQ. 'rose' ) THEN WRITE ( 6, 580 ) 580 FORMAT ( // ' You selected a Rosenfeld Window ', // ) C Apply a Rosenfeld window to the input data. C This will be 3 terms with coefficients +0.724, -0.950, C +0.226 TempR = 2.0 * 3.1415926 / 1024.0 DO 620 INX=1,1024,1 INXR=(INX*2)-1 TempI = 0.724 TempI = TempI - (0.950 * COS( 1.0 * TempR * DFloat(INX) )) TempI = TempI + (0.226 * COS( 2.0 * TempR * DFloat(INX) )) CDB(INXR) = CDB(INXR) * TempI 620 CONTINUE GOTO 700 ELSE WRITE ( 6, 590 ) 590 FORMAT ( // ' You must select from: ', / & ' rec, ham, black, rose ', // ) GOTO 500 END IF 700 CONTINUE C Now do the actual Fourier Transform. This is a 1024 point transform. C i.e. N = 2 * the number of points in the transform N = 2048 J = 1 DO 5 I=1,N,2 IF (I-J) 1,2,2 1 TempR = CDB(J) TempI = CDB(J+1) CDB(J) = CDB(I) CDB(J+1) = CDB(I+1) CDB(I) = TempR CDB(I+1) = TempI 2 M = N/2 3 IF (J-M) 5,5,4 4 J = J - M M = M / 2 IF (M-2) 5,3,3 5 J = J + M MMax = 2 6 IF (MMax - N) 7,10,10 7 IStep = 2 * MMax Theta = 6.2831853 / DFloat(MMax) SinTh = SIN(Theta/2.0) WStpR = -2.0 * SinTh * SinTh WStpI = SIN(Theta) WR = 1.0 WI = 0.0 DO 9 M=1,MMax,2 DO 8 I=M,N,IStep J = I + MMax TempR = WR * CDB(J) - WI * CDB(J+1) TempI = WR * CDB(J+1) + WI * CDB(J) CDB(J) = CDB(I) - TempR CDB(J+1) = CDB(I+1) - TempI CDB(I) = CDB(I) + TempR 8 CDB(I+1) = CDB(I+1) + TempI TempR = WR WR = WR + WR * WStpR - WI * WstpI 9 WI = WI + WI * WStpR + TempR * WStpI MMax = IStep GOTO 6 10 CONTINUE C C Now calculate the magnitude of each point in the transform. C Put these magnitude values into every other location in the C array, i.e. the index odd locations. From now on we will C work only with these magnitude values. DO 40 INX=1,2048,2 TempR = ( CDB(INX) * CDB(INX) ) + ( CDB(INX+1) * CDB(INX+1) ) CDB(INX) = SQRT( TempR ) 40 CONTINUE C C Now zero the first point in the array and find the maximum C magnitude value. C CDB(1) = 0.0 CDBMax = 0.0 DO 60 INX=1,2048,2 IF (CDB(INX)-CDBMax) 70,70,80 80 CDBMax = CDB(INX) 70 CONTINUE 60 CONTINUE C C Now normalize the array to a value of 1.0 Once again it C is only the magnitude values that we are working with. C If after normalization any value is less than 0.000001, then C set it to 0.000001 DO 90 INX=1,2048,2 CDB(INX) = CDB(INX) / CDBMax IF ( CDB(INX) .LE. 0.000001 ) THEN CDB(INX) = 0.000001 END IF 90 CONTINUE RETURN END