Directory of image this file is from
This file as a plain text file
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-