SUBROUTINE FOUR2 (DATA,N,NDIM,ISIGN,IFORM) FF2 1 C COOLEY-TUKEY FAST FOURIER TRANSFORM IN USASI BASIC FORTRAN. FF2 2 C MULTI-DIMENSIONAL TRANSFORM, EACH DIMENSION A POWER OF TWO, FF2 3 C COMPLEX OR REAL DATA. FF2 4 C TRANSFORM(K1,K2,...) = SUM(DATA(J1,J2,...)*EXP(ISIGN*2*PI*SQRT(-1)FF2 5 C *((J1-1)*(K1-1)/N(1)+(J2-1)*(K2-1)/N(2)+...))), SUMMED FOR ALL FF2 6 C J1 AND K1 FROM 1 TO N(1), J2 AND K2 FROM 1 TO N(2), ETC. FOR ALL FF2 7 C NDIM SUBSCRIPTS. NDIM MUST BE POSITIVE AND EACH N(IDIM) MUST BE FF2 8 C A POWER OF TWO. ISIGN IS +1 OR -1. LET NTOT = N(1)*N(2)... FF2 9 C ...*N(NDIM). THEN A -1 TRANSFORM FOLLOWED BY A +1 ONE FF2 10 C (OR VICE VERSA) RETURNS NTOT (NTOT/2 IF IFORM = 0 OR FF2 11 C -1) TIMES THE ORIGINAL DATA. IFORM = 1, 0 OR -1, AS DATA IS FF2 12 C COMPLEX, REAL OR THE FIRST HALF OF A COMPLEX ARRAY. TRANSFORM FF2 13 C VALUES ARE RETURNED TO ARRAY DATA. THEY ARE COMPLEX, REAL OR FF2 14 C THE FIRST HALF OF A COMPLEX ARRAY, AS IFORM = 1, -1 OR 0. FF2 15 C THE TRANSFORM OF A REAL ARRAY (IFORM = 0) DIMENSIONED N(1) BY N(2)FF2 16 C BY ... WILL BE RETURNED IN THE SAME ARRAY, NOW CONSIDERED TO FF2 17 C BE COMPLEX OF DIMENSIONS N(1)/2+1 BY N(2) BY .... NOTE THAT IF FF2 18 C IFORM = 0 OR -1, N(1) MUST BE EVEN, AND ENOUGH ROOM MUST BE FF2 19 C RESERVED. THE MISSING VALUES MAY BE OBTAINED BY COMPLEX CONJU- FF2 20 C GATION. THE REVERSE TRANSFORMATION, OF A HALF COMPLEX ARRAY FF2 21 C DIMENSIONED N(1)/2+1 BY N(2) BY ..., IS ACCOMPLISHED SETTING IFORMFF2 22 C TO -1. IN THE N ARRAY, N(1) MUST BE THE TRUE N(1), NOT N(1)/2+1. FF2 23 C THE TRANSFORM WILL BE REAL AND RETURNED TO THE INPUT ARRAY. FF2 24 C RUNNING TIME IS PROPORTIONAL TO NTOT*LOG2(NTOT), RATHER THAN FF2 25 C THE NAIVE NTOT**2. FF2 26 C WRITTEN BY NORMAN BRENNER OF MIT LINCOLN LABORATORY, JUNE 1968. FF2 27 C SEE-- IEEE AUDIO TRANSACTIONS (JUNE 1967), SPECIAL ISSUE ON FFT. FF2 28 DIMENSION DATA(1), N(1) FF2 29 NTOT=1 FF2 30 DO 10 IDIM=1,NDIM FF2 31 10 NTOT=NTOT*N(IDIM) FF2 32 IF (IFORM) 70,20,20 FF2 33 20 NREM=NTOT FF2 34 DO 60 IDIM=1,NDIM FF2 35 NREM=NREM/N(IDIM) FF2 36 NPREV=NTOT/(N(IDIM)*NREM) FF2 37 NCURR=N(IDIM) FF2 38 IF (IDIM-1+IFORM) 30,30,40 FF2 39 30 NCURR=NCURR/2 FF2 40 40 CALL BITRV (DATA,NPREV,NCURR,NREM) FF2 41 CALL COOL2 (DATA,NPREV,NCURR,NREM,ISIGN) FF2 42 IF (IDIM-1+IFORM) 50,50,60 FF2 43 50 CALL FIXRL (DATA,N(1),NREM,ISIGN,IFORM) FF2 44 NTOT=(NTOT/N(1))*(N(1)/2+1) FF2 45 60 CONTINUE FF2 46 RETURN FF2 47 70 NTOT=(NTOT/N(1))*(N(1)/2+1) FF2 48 NREM=1 FF2 49 DO 100 JDIM=1,NDIM FF2 50 IDIM=NDIM+1-JDIM FF2 51 NCURR=N(IDIM) FF2 52 IF (IDIM-1) 80,80,90 FF2 53 80 NCURR=NCURR/2 FF2 54 CALL FIXRL (DATA,N(1),NREM,ISIGN,IFORM) FF2 55 NTOT=NTOT/(N(1)/2+1)*N(1) FF2 56 90 NPREV=NTOT/(N(IDIM)*NREM) FF2 57 CALL BITRV (DATA,NPREV,NCURR,NREM) FF2 58 CALL COOL2 (DATA,NPREV,NCURR,NREM,ISIGN) FF2 59 100 NREM=NREM*N(IDIM) FF2 60 RETURN FF2 61 END FF2 62- SUBROUTINE FIXRL (DATA,N,NREM,ISIGN,IFORM) FIX 1 C FOR IFORM = 0, CONVERT THE TRANSFORM OF A DOUBLED-UP REAL ARRAY, FIX 2 C CONSIDERED COMPLEX, INTO ITS TRUE TRANSFORM. SUPPLY ONLY THE FIX 3 C FIRST HALF OF THE COMPLEX TRANSFORM, AS THE SECOND HALF HAS FIX 4 C CONJUGATE SYMMETRY. FOR IFORM = -1, CONVERT THE FIRST HALF FIX 5 C OF THE TRUE TRANSFORM INTO THE TRANSFORM OF A DOUBLED-UP REAL FIX 6 C ARRAY. N MUST BE EVEN. FIX 7 C USING COMPLEX NOTATION AND SUBSCRIPTS STARTING AT ZERO, THE FIX 8 C TRANSFORMATION IS-- FIX 9 C DIMENSION DATA(N,NREM) FIX 10 C ZSTP = EXP(ISIGN*2*PI*I/N) FIX 11 C DO 10 I2=0,NREM-1 FIX 12 C DATA(0,I2) = CONJ(DATA(0,I2))*(1+I)/(1-IFORM) FIX 13 C DO 10 I1=1,N/4 FIX 14 C Z = (1+(2*IFORM+1)*I*ZSTP**I1)/2 FIX 15 C I1CNJ = N/2-I1 FIX 16 C DIF = DATA(I1,I2)-CONJ(DATA(I1CNJ,I2)) FIX 17 C TEMP = Z*DIF FIX 18 C DATA(I1,I2) = DATA(I1,I2)-TEMP FIX 19 C 10 DATA(I1CNJ,I2) = DATA(I1CNJ,I2)+CONJ(TEMP) FIX 20 C IF I1=I1CNJ, THE CALCULATION FOR THAT VALUE COLLAPSES INTO FIX 21 C A SIMPLE CONJUGATION OF DATA(I1,I2). FIX 22 DIMENSION DATA(1) FIX 23 TWOPI=6.283185307*FLOAT(ISIGN) FIX 24 IP0=2 FIX 25 IP1=IP0*(N/2) FIX 26 IP2=IP1*NREM FIX 27 IF (IFORM) 10,60,60 FIX 28 10 J1=IP1+1 FIX 29 DATA(2)=DATA(J1) FIX 30 J1=J1+IP0 FIX 31 I2MIN=IP1+1 FIX 32 DO 50 I2=I2MIN,IP2,IP1 FIX 33 DATA(I2)=DATA(J1) FIX 34 J1=J1+IP0 FIX 35 IF (N-2) 40,40,20 FIX 36 20 I1MIN=I2+IP0 FIX 37 I1MAX=I2+IP1-IP0 FIX 38 DO 30 I1=I1MIN,I1MAX,IP0 FIX 39 DATA(I1)=DATA(J1) FIX 40 DATA(I1+1)=DATA(J1+1) FIX 41 30 J1=J1+IP0 FIX 42 40 DATA(I2+1)=DATA(J1) FIX 43 50 J1=J1+IP0 FIX 44 60 DO 80 I2=1,IP2,IP1 FIX 45 TEMPR=DATA(I2) FIX 46 DATA(I2)=DATA(I2)+DATA(I2+1) FIX 47 DATA(I2+1)=TEMPR-DATA(I2+1) FIX 48 IF (IFORM) 70,80,80 FIX 49 70 DATA(I2)=DATA(I2)/2. FIX 50 DATA(I2+1)=DATA(I2+1)/2. FIX 51 80 CONTINUE FIX 52 IF (N-2) 170,170,90 FIX 53 90 THETA=TWOPI/FLOAT(N) FIX 54 SINTH=SIN(THETA/2.) FIX 55 ZSTPR=-2.*SINTH*SINTH FIX 56 ZSTPI=SIN(THETA) FIX 57 ZR=(1.-ZSTPI)/2. FIX 58 ZI=(1.+ZSTPR)/2. FIX 59 IF (IFORM) 100,110,110 FIX 60 100 ZR=1.-ZR FIX 61 ZI=-ZI FIX 62 110 I1MIN=IP0+1 FIX 63 I1MAX=IP0*(N/4)+1 FIX 64 DO 160 I1=I1MIN,I1MAX,IP0 FIX 65 DO 150 I2=I1,IP2,IP1 FIX 66 I2CNJ=N+IP0-2*I1+I2 FIX 67 IF (I2-I2CNJ) 140,120,120 FIX 68 120 IF (ISIGN*(2*IFORM+1)) 130,150,150 FIX 69 130 DATA(I2+1)=-DATA(I2+1) FIX 70 GO TO 150 FIX 71 140 DIFR=DATA(I2)-DATA(I2CNJ) FIX 72 DIFI=DATA(I2+1)+DATA(I2CNJ+1) FIX 73 TEMPR=DIFR*ZR-DIFI*ZI FIX 74 TEMPI=DIFR*ZI+DIFI*ZR FIX 75 DATA(I2)=DATA(I2)-TEMPR FIX 76 DATA(I2+1)=DATA(I2+1)-TEMPI FIX 77 DATA(I2CNJ)=DATA(I2CNJ)+TEMPR FIX 78 DATA(I2CNJ+1)=DATA(I2CNJ+1)-TEMPI FIX 79 150 CONTINUE FIX 80 TEMPR=ZR-.5 FIX 81 ZR=ZSTPR*TEMPR-ZSTPI*ZI+ZR FIX 82 160 ZI=ZSTPR*ZI+ZSTPI*TEMPR+ZI FIX 83 170 IF (IFORM) 240,180,180 FIX 84 180 I2=IP2+1 FIX 85 I1=I2 FIX 86 J1=IP0*(N/2+1)*NREM+1 FIX 87 GO TO 220 FIX 88 190 DATA(J1)=DATA(I1) FIX 89 DATA(J1+1)=DATA(I1+1) FIX 90 I1=I1-IP0 FIX 91 J1=J1-IP0 FIX 92 200 IF (I2-I1) 190,210,210 FIX 93 210 DATA(J1)=DATA(I1) FIX 94 DATA(J1+1)=0. FIX 95 220 I2=I2-IP1 FIX 96 J1=J1-IP0 FIX 97 DATA(J1)=DATA(I2+1) FIX 98 DATA(J1+1)=0. FIX 99 I1=I1-IP0 FIX 100 J1=J1-IP0 FIX 101 IF (I2-1) 230,230,200 FIX 102 230 DATA(2)=0. FIX 103 240 RETURN FIX 104 END FIX 105- SUBROUTINE COOL2 (DATA,NPREV,N,NREM,ISIGN) CO2 1 C FOURIER TRANSFORM OF LENGTH N BY THE COOLEY-TUKEY ALGORITHM. CO2 2 C BIT-REVERSED TO NORMAL ORDER. CO2 3 C DIMENSION DATA(NPREV,N,NREM) CO2 4 C COMPLEX DATA CO2 5 C DATA(I1,J2,I3) = SUM(DATA(I1,I2,I3)*EXP(ISIGN*2*PI*I*((I2-1)* CO2 6 C (J2-1)/N))), SUMMED OVER I2 = 1 TO N FOR ALL I1 FROM 1 TO NPREV, CO2 7 C J2 FROM 1 TO N AND I3 FROM 1 TO NREM. N MUST BE A POWER OF TWO. CO2 8 C FACTORING N BY FOURS SAVES ABOUT TWENTY FIVE PERCENT OVER FACTOR- CO2 9 C ING BY TWOS. CO2 10 C NOTE--IT IS NOT NECESSARY TO REWRITE THIS SUBROUTINE INTO COMPLEX CO2 11 C NOTATION SO LONG AS THE FORTRAN COMPILER USED STORES REAL AND CO2 12 C IMAGINARY PARTS IN ADJACENT STORAGE LOCATIONS. IT MUST ALSO CO2 13 C STORE ARRAYS WITH THE FIRST SUBSCRIPT INCREASING FASTEST. CO2 14 DIMENSION DATA(1) CO2 15 TWOPI=6.2831853072*FLOAT(ISIGN) CO2 16 IP0=2 CO2 17 IP1=IP0*NPREV CO2 18 IP4=IP1*N CO2 19 IP5=IP4*NREM CO2 20 IP2=IP1 CO2 21 NPART=N CO2 22 10 IF (NPART-2) 50,30,20 CO2 23 20 NPART=NPART/4 CO2 24 GO TO 10 CO2 25 C DO A FOURIER TRANSFORM OF LENGTH TWO CO2 26 30 IP3=IP2*2 CO2 27 DO 40 I1=1,IP1,IP0 CO2 28 DO 40 I5=I1,IP5,IP3 CO2 29 J0=I5 CO2 30 J1=J0+IP2 CO2 31 TEMPR=DATA(J1) CO2 32 TEMPI=DATA(J1+1) CO2 33 DATA(J1)=DATA(J0)-TEMPR CO2 34 DATA(J1+1)=DATA(J0+1)-TEMPI CO2 35 DATA(J0)=DATA(J0)+TEMPR CO2 36 40 DATA(J0+1)=DATA(J0+1)+TEMPI CO2 37 GO TO 140 CO2 38 C DO A FOURIER TRANSFORM OF LENGTH FOUR (FROM BIT REVERSED ORDER) CO2 39 50 IP3=IP2*4 CO2 40 THETA=TWOPI/FLOAT(IP3/IP1) CO2 41 SINTH=SIN(THETA/2.) CO2 42 WSTPR=-2.*SINTH*SINTH CO2 43 C COS(THETA)-1, FOR ACCURACY. CO2 44 WSTPI=SIN(THETA) CO2 45 WR=1. CO2 46 WI=0. CO2 47 DO 130 I2=1,IP2,IP1 CO2 48 IF (I2-1) 70,70,60 CO2 49 60 W2R=WR*WR-WI*WI CO2 50 W2I=2.*WR*WI CO2 51 W3R=W2R*WR-W2I*WI CO2 52 W3I=W2R*WI+W2I*WR CO2 53 70 I1MAX=I2+IP1-IP0 CO2 54 DO 120 I1=I2,I1MAX,IP0 CO2 55 DO 120 I5=I1,IP5,IP3 CO2 56 J0=I5 CO2 57 J1=J0+IP2 CO2 58 J2=J1+IP2 CO2 59 J3=J2+IP2 CO2 60 IF (I2-1) 90,90,80 CO2 61 C APPLY THE PHASE SHIFT FACTORS CO2 62 80 TEMPR=DATA(J1) CO2 63 DATA(J1)=W2R*TEMPR-W2I*DATA(J1+1) CO2 64 DATA(J1+1)=W2R*DATA(J1+1)+W2I*TEMPR CO2 65 TEMPR=DATA(J2) CO2 66 DATA(J2)=WR*TEMPR-WI*DATA(J2+1) CO2 67 DATA(J2+1)=WR*DATA(J2+1)+WI*TEMPR CO2 68 TEMPR=DATA(J3) CO2 69 DATA(J3)=W3R*TEMPR-W3I*DATA(J3+1) CO2 70 DATA(J3+1)=W3R*DATA(J3+1)+W3I*TEMPR CO2 71 90 T0R=DATA(J0)+DATA(J1) CO2 72 T0I=DATA(J0+1)+DATA(J1+1) CO2 73 T1R=DATA(J0)-DATA(J1) CO2 74 T1I=DATA(J0+1)-DATA(J1+1) CO2 75 T2R=DATA(J2)+DATA(J3) CO2 76 T2I=DATA(J2+1)+DATA(J3+1) CO2 77 T3R=DATA(J2)-DATA(J3) CO2 78 T3I=DATA(J2+1)-DATA(J3+1) CO2 79 DATA(J0)=T0R+T2R CO2 80 DATA(J0+1)=T0I+T2I CO2 81 DATA(J2)=T0R-T2R CO2 82 DATA(J2+1)=T0I-T2I CO2 83 IF (ISIGN) 100,100,110 CO2 84 100 T3R=-T3R CO2 85 T3I=-T3I CO2 86 110 DATA(J1)=T1R-T3I CO2 87 DATA(J1+1)=T1I+T3R CO2 88 DATA(J3)=T1R+T3I CO2 89 120 DATA(J3+1)=T1I-T3R CO2 90 TEMPR=WR CO2 91 WR=WSTPR*TEMPR-WSTPI*WI+TEMPR CO2 92 130 WI=WSTPR*WI+WSTPI*TEMPR+WI CO2 93 140 IP2=IP3 CO2 94 IF (IP3-IP4) 50,150,150 CO2 95 150 RETURN CO2 96 END CO2 97- SUBROUTINE BITRV (DATA,NPREV,N,NREM) BIT 1 C SHUFFLE THE DATA BY @BIT REVERSAL@. BIT 2 C DIMENSION DATA(NPREV,N,NREM) BIT 3 C DATA(I1,I2REV,I3) = DATA(I1,I2,I3), FOR ALL I1 FROM 1 TO NPREV, BIT 4 C ALL I2 FROM 1 TO N (WHICH MUST BE A POWER OF TWO), AND ALL I3 BIT 5 C FROM 1 TO NREM, WHERE I2REV-1 IS THE BITWISE REVERSAL OF I2-1. BIT 6 C FOR EXAMPLE, N = 32, I2-1 = 10011 AND I2REV-1 = 11001. BIT 7 DIMENSION DATA(1) BIT 8 IP0=2 BIT 9 IP1=IP0*NPREV BIT 10 IP4=IP1*N BIT 11 IP5=IP4*NREM BIT 12 I4REV=1 BIT 13 DO 60 I4=1,IP4,IP1 BIT 14 IF (I4-I4REV) 10,30,30 BIT 15 10 I1MAX=I4+IP1-IP0 BIT 16 DO 20 I1=I4,I1MAX,IP0 BIT 17 DO 20 I5=I1,IP5,IP4 BIT 18 I5REV=I4REV+I5-I4 BIT 19 TEMPR=DATA(I5) BIT 20 TEMPI=DATA(I5+1) BIT 21 DATA(I5)=DATA(I5REV) BIT 22 DATA(I5+1)=DATA(I5REV+1) BIT 23 DATA(I5REV)=TEMPR BIT 24 20 DATA(I5REV+1)=TEMPI BIT 25 30 IP2=IP4/2 BIT 26 40 IF (I4REV-IP2) 60,60,50 BIT 27 50 I4REV=I4REV-IP2 BIT 28 IP2=IP2/2 BIT 29 IF (IP2-IP1) 60,40,40 BIT 30 60 I4REV=I4REV+IP2 BIT 31 RETURN BIT 32 END BIT 33-