8-AUG-1997 20:47:53 DEC Fortran V6.3-711 Page 1 6-AUG-1997 12:22:28 DEPT1:[EDMUNDS.SPARE_6]FORWARD_FOURIER.FOR;1 1 C 2 C This subroutine thakes the data from the integer input array 3 C that is 1024 locations long and puts this data into the real 4 C array that will be transformed. These 1024 points are forward 5 C transformed into 1024 complex pairs. 6 C 7 C Next these 1024 complex pairs are used to calculate 1024 magnitude 8 C values. These magnitude values are then normalized to 1.0 9 C This final result is stored at the odd index locations in the 10 C real array that is 2048 locations long. 11 C 12 C 13 C Created 4-AUG-1997 Daniel Edmunds 14 C 15 16 SUBROUTINE FORWARD_FOURIER 17 18 C ***************************************** 19 C Define all Variables and Arrays 20 C ***************************************** 21 22 IMPLICIT NONE 23 24 REAL*8 CDB(1:2048) 25 INTEGER*2 ICDB(1:1024) 26 COMMON ICDB, CDB 27 28 INTEGER*2 I, INX, INXI, INXR, IStep, J, M, MMax, N 29 30 REAL*8 CDBMax, TempR, TempI, Theta, SinTh, WI, WR, WStpR, WStpI 31 32 CHARACTER FuncName*32 33 34 35 C Put the fixed point discrete time sample ADC data array 36 C into the floating point array that will be transformed. 37 38 DO 20 INX=1,1024,1 39 INXR=(INX*2)-1 40 INXI=INX*2 41 CDB(INXR) = DFloat(ICDB(INX)) 42 CDB(INXI) = 0.0 43 20 CONTINUE 44 45 46 C Now we want to center this input data about zero before 47 C multiplying it by a windowing function. Note that we are 48 C only working on the entries in the array into which input 49 C data was just loaded, i.e. the odd entries. 50 51 TempR = 0.0 52 53 DO 25 INX=1,2048,2 54 TempR = TempR + CDB(INX) 55 25 CONTINUE 56 57 TempR = TempR / 1024.0 FORWARD_FOURIER 8-AUG-1997 20:47:53 DEC Fortran V6.3-711 Page 2 6-AUG-1997 12:22:28 DEPT1:[EDMUNDS.SPARE_6]FORWARD_FOURIER.FOR;1 58 59 DO 26 INX=1,2048,2 60 CDB(INX) = CDB(INX) - TempR 61 26 CONTINUE 62 63 64 65 66 67 C 68 C We now have a chance to apply a "Windowing Function" to the 69 C input data before we transform it to frequency space. 70 C We can select from the following Window Functions: 71 C Rectangular i.e. no window function 72 C Hamming 0.08 + 0.46( 1 - cos 2pi n / N ) 73 C Blackman-Harris 0.423 - 0.497 cos 2pi n/N + 0.079 cos 4pi n/N 74 C Rosenfeld 0.724 - 0.950 cos 2pi n/N + 0.226 cos 4pi n/N 75 C 76 500 WRITE ( 6, 505 ) 77 505 FORMAT ( // 78 & ' A Windowing Function may be applied to the input ', / 79 & ' data before the transform to frequency space is ', / 80 & ' performed. Pick from the following window functions: ', / 81 & ' Rectangular = rec Hamming = ham ', / 82 & ' Rosenfeld = rose Blackman-Harris = black ', / 83 & ' enter: rec, ham, rose, black : ', $ ) 84 85 READ ( 5, 510 ) FuncName 86 510 FORMAT ( A ) 87 88 89 IF ( FuncName .EQ. 'rec ' ) THEN 90 WRITE ( 6, 550 ) 91 550 FORMAT ( // ' You selected a Rectangular Window ', // ) 92 GOTO 700 93 94 95 ELSE IF ( FuncName .EQ. 'ham' ) THEN 96 WRITE ( 6, 560 ) 97 560 FORMAT ( // ' You selected a Hamming Window ', // ) 98 99 C Apply the Hamming window to the input data. 100 C This will be a 0.46 * cosine on an 8% pedestal. 101 102 TempR = 2.0 * 3.1415926 / 1024.0 103 DO 600 INX=1,1024,1 104 INXR=(INX*2)-1 105 TempI = 0.46 * ( 1.0 - COS( TempR * DFloat(INX) )) 106 CDB(INXR) = CDB(INXR) * ( 0.08 + TempI ) 107 600 CONTINUE 108 109 GOTO 700 110 111 112 ELSE IF ( FuncName .EQ. 'black' ) THEN 113 WRITE ( 6, 570 ) 114 570 FORMAT ( // ' You selected a Blackman-Harris Window ', // ) FORWARD_FOURIER 8-AUG-1997 20:47:53 DEC Fortran V6.3-711 Page 3 6-AUG-1997 12:22:28 DEPT1:[EDMUNDS.SPARE_6]FORWARD_FOURIER.FOR;1 115 116 C Apply a Blackman-Harris window to the input data. 117 C This will be 3 terms with coefficients +0.42323, -0.49755, 118 C +0.07922 119 120 TempR = 2.0 * 3.1415926 / 1024.0 121 122 DO 610 INX=1,1024,1 123 INXR=(INX*2)-1 124 TempI = 0.42323 125 TempI = TempI - (0.49755 * COS( 1.0 * TempR * DFloat(INX) )) 126 TempI = TempI + (0.07922 * COS( 2.0 * TempR * DFloat(INX) )) 127 CDB(INXR) = CDB(INXR) * TempI 128 610 CONTINUE 129 130 GOTO 700 131 132 133 ELSE IF ( FuncName .EQ. 'rose' ) THEN 134 WRITE ( 6, 580 ) 135 580 FORMAT ( // ' You selected a Rosenfeld Window ', // ) 136 137 C Apply a Rosenfeld window to the input data. 138 C This will be 3 terms with coefficients +0.724, -0.950, 139 C +0.226 140 141 TempR = 2.0 * 3.1415926 / 1024.0 142 143 DO 620 INX=1,1024,1 144 INXR=(INX*2)-1 145 TempI = 0.724 146 TempI = TempI - (0.950 * COS( 1.0 * TempR * DFloat(INX) )) 147 TempI = TempI + (0.226 * COS( 2.0 * TempR * DFloat(INX) )) 148 CDB(INXR) = CDB(INXR) * TempI 149 620 CONTINUE 150 151 GOTO 700 152 153 154 ELSE 155 WRITE ( 6, 590 ) 156 590 FORMAT ( // ' You must select from: ', / 157 & ' rec, ham, black, rose ', // ) 158 GOTO 500 159 160 END IF 161 162 163 164 165 166 167 168 169 700 CONTINUE 170 171 C Now do the actual Fourier Transform. This is a 1024 point transform. FORWARD_FOURIER 8-AUG-1997 20:47:53 DEC Fortran V6.3-711 Page 4 6-AUG-1997 12:22:28 DEPT1:[EDMUNDS.SPARE_6]FORWARD_FOURIER.FOR;1 172 C i.e. N = 2 * the number of points in the transform 173 174 N = 2048 175 J = 1 176 177 DO 5 I=1,N,2 178 179 IF (I-J) 1,2,2 180 181 1 TempR = CDB(J) 182 TempI = CDB(J+1) 183 184 CDB(J) = CDB(I) 185 CDB(J+1) = CDB(I+1) 186 187 CDB(I) = TempR 188 CDB(I+1) = TempI 189 190 2 M = N/2 191 192 3 IF (J-M) 5,5,4 193 194 4 J = J - M 195 M = M / 2 196 197 IF (M-2) 5,3,3 198 199 5 J = J + M 200 201 202 MMax = 2 203 6 IF (MMax - N) 7,10,10 204 7 IStep = 2 * MMax 205 206 Theta = 6.2831853 / DFloat(MMax) 207 SinTh = SIN(Theta/2.0) 208 209 WStpR = -2.0 * SinTh * SinTh 210 WStpI = SIN(Theta) 211 212 WR = 1.0 213 WI = 0.0 214 215 DO 9 M=1,MMax,2 216 217 DO 8 I=M,N,IStep 218 219 J = I + MMax 220 221 TempR = WR * CDB(J) - WI * CDB(J+1) 222 TempI = WR * CDB(J+1) + WI * CDB(J) 223 224 CDB(J) = CDB(I) - TempR 225 CDB(J+1) = CDB(I+1) - TempI 226 CDB(I) = CDB(I) + TempR 227 8 CDB(I+1) = CDB(I+1) + TempI 228 FORWARD_FOURIER 8-AUG-1997 20:47:53 DEC Fortran V6.3-711 Page 5 6-AUG-1997 12:22:28 DEPT1:[EDMUNDS.SPARE_6]FORWARD_FOURIER.FOR;1 229 TempR = WR 230 WR = WR + WR * WStpR - WI * WstpI 231 9 WI = WI + WI * WStpR + TempR * WStpI 232 233 MMax = IStep 234 GOTO 6 235 236 10 CONTINUE 237 238 C 239 C Now calculate the magnitude of each point in the transform. 240 C Put these magnitude values into every other location in the 241 C array, i.e. the index odd locations. From now on we will 242 C work only with these magnitude values. 243 244 DO 40 INX=1,2048,2 245 TempR = ( CDB(INX) * CDB(INX) ) + ( CDB(INX+1) * CDB(INX+1) ) 246 CDB(INX) = SQRT( TempR ) 247 40 CONTINUE 248 249 250 C 251 C Now zero the first point in the array and find the maximum 252 C magnitude value. 253 C 254 CDB(1) = 0.0 255 CDBMax = 0.0 256 DO 60 INX=1,2048,2 257 IF (CDB(INX)-CDBMax) 70,70,80 258 80 CDBMax = CDB(INX) 259 70 CONTINUE 260 60 CONTINUE 261 262 263 C 264 C Now normalize the array to a value of 1.0 Once again it 265 C is only the magnitude values that we are working with. 266 C If after normalization any value is less than 0.000001, then 267 C set it to 0.000001 268 269 DO 90 INX=1,2048,2 270 271 CDB(INX) = CDB(INX) / CDBMax 272 273 IF ( CDB(INX) .LE. 0.000001 ) THEN 274 CDB(INX) = 0.000001 275 END IF 276 277 90 CONTINUE 278 279 RETURN 280 281 END FORWARD_FOURIER 8-AUG-1997 20:47:53 DEC Fortran V6.3-711 Page 6 Symbol Table 6-AUG-1997 12:22:28 DEPT1:[EDMUNDS.SPARE_6]FORWARD_FOURIER.FOR;1 PROGRAM SECTIONS Name Bytes Attributes 1 $BSS$ 176 NOPIC CON REL LCL NOSHR NOEXE RD WRT NOSEC OCTA 2 $IODATA$ 760 NOPIC CON REL LCL NOSHR NOEXE RD WRT NOSEC OCTA 3 $CODE$ 5192 PIC CON REL LCL SHR EXE NORD NOWRT NOSEC OCTA 4 $LINK$ 336 NOPIC CON REL LCL NOSHR NOEXE RD NOWRT NOSEC OCTA 5 $BLANK 18432 NOPIC OVR REL GBL NOSHR NOEXE RD WRT NOSEC OCTA Total Space Allocated 24896 ENTRY POINTS Address Type Name 3-00000000 FORWARD_FOURIER VARIABLES Address Type Name Address Type Name Address Type Name Address Type Name 1-00000048 R*8 CDBMAX 1-00000018 I*2 INXR 1-00000040 I*2 N 1-00000070 R*8 WI 1-00000090 CHAR FUNCNAME 1-00000020 I*2 ISTEP 1-00000068 R*8 SINTH 1-00000078 R*8 WR 1-00000000 I*2 I 1-00000028 I*2 J 1-00000058 R*8 TEMPI 1-00000088 R*8 WSTPI 1-00000008 I*2 INX 1-00000030 I*2 M 1-00000050 R*8 TEMPR 1-00000080 R*8 WSTPR 1-00000010 I*2 INXI 1-00000038 I*2 MMAX 1-00000060 R*8 THETA ARRAYS Address Type Name Bytes Dimensions 5-00000800 R*8 CDB 16384 (2048) 5-00000000 I*2 ICDB 2048 (1024) LABELS Address Label Address Label Address Label Address Label Address Label 3-0000091C 1 3-00000BEC 7 3-00000000 26 3-000002B4 500 ** 580 3-00000AAC 2 3-00000F70 8 3-00000000 40 ** 505 ** 590 3-00000AE4 3 3-00001054 9 3-000012B8 60 ** 510 3-00000000 600 3-00000AFC 4 3-000010DC 10 3-000012B8 70 ** 550 3-00000000 610 3-00000B6C 5 3-00000000 20 3-00001284 80 ** 560 3-00000000 620 3-00000BD4 6 3-00000000 25 3-000013C8 90 ** 570 3-000008C8 700 FORWARD_FOURIER 8-AUG-1997 20:47:53 DEC Fortran V6.3-711 Page 7 Symbol Table 6-AUG-1997 12:22:28 DEPT1:[EDMUNDS.SPARE_6]FORWARD_FOURIER.FOR;1 282 +---------------------------------------------------+ | KEY TO ADDRESS CODE FORMATS | | ppp-oooooooo - In Psect ppp, Offset oooooooo | | ***-******** - External | | # - Suffix: Also In Registers | | REG-rrrrrrrr - In Register rrrrrrrr | | REG-######## - In Various Registers | | ** - Not Used; Not Allocated | +---------------------------------------------------+ COMMAND QUALIFIERS /ALIGNMENT=(COMMONS=(PACKED,NOMULTILANGUAGE),RECORDS=NATURAL) /ASSUME=(ACCURACY_SENSITIVE,BACKSLASH,NODUMMY_ALIASES,NOUNDERSCORE,NOSOURCE_INCLUDE) /CHECK=(ASSERTIONS,BOUNDS,FORMAT,FP_EXCEPTIONS,OUTPUT_CONVERSION,OVERFLOW,UNDERFLOW) /DEBUG=(PARAMETERS=NONE,SYMBOLS,TRACEBACK) /DESIGN=(NOCOMMENTS) /SHOW=(NODICTIONARY,NOINCLUDE,MAP,NOPREPROCESSOR,SINGLE) /STANDARD=(NOMIA,NOSEMANTIC,NOSOURCE_FORM,NOSYNTAX) /WARNINGS=(ALIGNMENTS,NOARGUMENT_CHECKING,NODECLARATIONS,GENERAL,INFORMATIONAL,NOTRUNCATED_SOURCE, UNCALLED,UNINITIALIZED,UNREACHABLE,NOUNUSED,USAGE) /NOAUTOMATIC /BLAS=NOMAPPED /CONVERT=NATIVE /NOCROSS_REFERENCE /NOD_LINES /DOUBLE_SIZE=64 /ERROR_LIMIT=30 /NOEXTEND_SOURCE /F77 /FLOAT=G_FLOAT /GRANULARITY=QUADWORD /IEEE_MODE=FAST /INSTRUCTION_SET=FLOATING_POINT /INTEGER_SIZE=32 /NOMACHINE_CODE /MATH_LIBRARY=ACCURATE /NAMES=UPPERCASE /OPTIMIZE=(INLINE=NONE,LEVEL=0,SPECULATE=NONE,TUNE=GENERIC,UNROLL=0) /NOPAD_SOURCE /REAL_SIZE=32 /NORECURSIVE /ROUNDING_MODE=NEAREST /NOSEPARATE_COMPILATION /NOSYNCHRONOUS_EXCEPTIONS /NOSYNTAX_ONLY /TERMINAL=NOSTATISTICS /NOTIE /VMS /NOANALYSIS_DATA /NODIAGNOSTICS /INCLUDE=(.FOR,.f,FORT$INCLUDE:.FOR,FORT$INCLUDE:.f) /LIST=DEPT1:[EDMUNDS.SPARE_6]FORWARD_FOURIER.LIS;1 /OBJECT=DEPT1:[EDMUNDS.SPARE_6]FORWARD_FOURIER.OBJ;1 /NOLIBRARY sys$lib=SYS$COMMON:[SYSLIB]FORSYSDEF.TLB;1 COMPILER: DEC Fortran V6.3-711-293V COMPILATION STATISTICS CPU time: 0.47 seconds Elapsed time: 0.88 seconds Pagefaults: 488 I/O Count: 14