SUBROUTINE FOURG (DATA,N,ISIGN,WORK) FFG 1
C COOLEY-TUKEY FAST FOURIER TRANSFORM IN USASI BASIC FORTRAN. FFG 2
C ONE-DIMENSIONAL TRANSFORM OF COMPLEX DATA, ARBITRARY NUMBER OF FFG 3
C POINTS. N POINTS CAN BE TRANSFORMED IN TIME PROPORTIONAL TO FFG 4
C N*LOG(N) (FOR N NON-PRIME), WHEREAS OTHER METHODS TAKE N**2 TIME. FFG 5
C FURTHERMORE, BECAUSE FEWER ARITHMETIC OPERATIONS ARE PERFORMED, FFG 6
C LESS ERROR IS BUILT UP. THE TRANSFORM DONE IS-- FFG 7
C DIMENSION DATA(N),TRANSFORM(N),WORK(N) FFG 8
C COMPLEX DATA,TRANSFORM,WORK FFG 9
C TRANSFORM(K) = SUM(DATA(J)*EXP(ISIGN*2*PI*I*(J-1)*(K-1)/N)), FFG 10
C SUMMED FROM J = 1 TO N FOR ALL K FROM 1 TO N. THE TRANSFORM FFG 11
C VALUES ARE RETURNED TO DATA, REPLACING THE INPUT. N MAY BE ANY FFG 12
C POSITIVE NUMBER, BUT IT SHOULD BE NON-PRIME FOR SPEED. ISIGN = FFG 13
C +1 OR -1. A -1 TRANSFORM FOLLOWED BY A +1 ONE (OR VICE VERSA) FFG 14
C RETURNS N TIMES THE ORIGINAL DATA. WORK IS A ONE-DIMENSIONAL FFG 15
C COMPLEX ARRAY OF LENGTH N USED FOR WORKING STORAGE. FFG 16
C RUNNING TIME IS PROPORTIONAL TO N * (SUM OF THE PRIME FACTORS OF FFG 17
C N). FOR EXAMPLE, N = 1960, TIME IS T0 * 1960 * (2+2+2+5+7+7). FFG 18
C NAIVE METHODS DIRECTLY IMPLEMENTING THE SUMMATION RUN IN TIME FFG 19
C PROPORTIONAL TO N**2. AN UPPER BOUND FOR THE RMS RELATIVE ERROR FFG 20
C IS 3 * 2**(-B) * SUM(F**1.5), WHERE B IS THE NUMBER OF BITS IN FFG 21
C THE FLOATING POINT FRACTION AND THE SUM IS OVER THE PRIME FFG 22
C FACTORS OF N. WRITTEN BY NORMAN BRENNER, MIT LINCOLN LABORATORY, FFG 23
C AUGUST 1968. SEE--IEEE TRANSACTIONS ON AUDIO AND ELECTROACOUSTICSFFG 24
C (JUNE 1967), SPECIAL ISSUE ON THE FAST FOURIER TRANSFORM. FFG 25
DIMENSION DATA(1), WORK(1), IFACT(32) FFG 26
TWOPI=6.283185307*FLOAT(ISIGN) FFG 27
C FACTOR N INTO ITS PRIME FACTORS, NFACT IN NUMBER. FOR EXAMPLE, FFG 28
C FOR N = 1960, NFACT = 6 AND IFACT(IF) = 2, 2, 2, 5, 7 AND 7. FFG 29
IF=0 FFG 30
NPART=N FFG 31
DO 50 ID=1,N,2 FFG 32
IDIV=ID FFG 33
IF (ID-1) 10,10,20 FFG 34
10 IDIV=2 FFG 35
20 IQUOT=NPART/IDIV FFG 36
IF (NPART-IDIV*IQUOT) 40,30,40 FFG 37
30 IF=IF+1 FFG 38
IFACT(IF)=IDIV FFG 39
NPART=IQUOT FFG 40
GO TO 20 FFG 41
40 IF (IQUOT-IDIV) 60,60,50 FFG 42
50 CONTINUE FFG 43
60 IF (NPART-1) 80,80,70 FFG 44
70 IF=IF+1 FFG 45
IFACT(IF)=NPART FFG 46
80 NFACT=IF FFG 47
C SHUFFLE THE DATA ARRAY BY REVERSING THE DIGITS OF THE INDEX. FFG 48
C REPLACE DATA(I) BY DATA(IREV) FOR ALL I FROM 1 TO N. IREV-1 IS FFG 49
C THE INTEGER WHOSE DIGIT REPRESENTATION IN THE MULTI-RADIX FFG 50
C NOTATION OF FACTORS IFACT(IF) IS THE REVERSE OF THE REPRESENTATIONFFG 51
C OF I-1. FOR EXAMPLE, IF ALL IFACT(IF) = 2, THEN FOR I-1 = 11001, FFG 52
C IREV-1 = 10011. A WORK ARRAY OF LENGTH N IS NEEDED. FFG 53
IP0=2 FFG 54
IP3=IP0*N FFG 55
IWORK=1 FFG 56
I3REV=1 FFG 57
DO 110 I3=1,IP3,IP0 FFG 58
WORK(IWORK)=DATA(I3REV) FFG 59
WORK(IWORK+1)=DATA(I3REV+1) FFG 60
IP2=IP3 FFG 61
DO 100 IF=1,NFACT FFG 62
IP1=IP2/IFACT(IF) FFG 63
I3REV=I3REV+IP1 FFG 64
IF (I3REV-IP2) 110,110,90 FFG 65
90 I3REV=I3REV-IP2 FFG 66
100 IP2=IP1 FFG 67
110 IWORK=IWORK+IP0 FFG 68
IWORK=1 FFG 69
DO 120 I3=1,IP3,IP0 FFG 70
DATA(I3)=WORK(IWORK) FFG 71
DATA(I3+1)=WORK(IWORK+1) FFG 72
120 IWORK=IWORK+IP0 FFG 73
C PHASE-SHIFTED FOURIER TRANSFORM OF LENGTH IFACT(IF). FFG 74
C IPROD=IP1/IP0 FFG 75
C IREM=N/(IFACT(IF)*IPROD) FFG 76
C DIMENSION DATA(IPROD,IFACT(IF),IREM),WORK(IFACT(IF)) FFG 77
C COMPLEX DATA,WORK FFG 78
C DATA(I1,J2,I3) = SUM(DATA(I1,I2,I3) * W**(I2-1)), SUMMED OVER FFG 79
C I2 = 1 TO IFACT(IF) FOR ALL I1 FROM 1 TO IPROD, J2 FROM 1 TO FFG 80
C IFACT(IF) AND I3 FROM 1 TO IREM. FFG 81
C W = EXP(ISIGN*2*PI*I*(I1-1+IPROD*(J2-1))/(IPROD*IFACT(IF))). FFG 82
IF=0 FFG 83
IP1=IP0 FFG 84
130 IF (IP1-IP3) 140,240,240 FFG 85
140 IF=IF+1 FFG 86
IFCUR=IFACT(IF) FFG 87
IP2=IP1*IFCUR FFG 88
THETA=TWOPI/FLOAT(IFCUR) FFG 89
SINTH=SIN(THETA/2.) FFG 90
ROOTR=-2.*SINTH*SINTH FFG 91
C COS(THETA)-1, FOR ACCURACY FFG 92
ROOTI=SIN(THETA) FFG 93
THETA=TWOPI/FLOAT(IP2/IP0) FFG 94
SINTH=SIN(THETA/2.) FFG 95
WSTPR=-2.*SINTH*SINTH FFG 96
WSTPI=SIN(THETA) FFG 97
WMINR=1. FFG 98
WMINI=0. FFG 99
DO 230 I1=1,IP1,IP0 FFG 100
IF (IFCUR-2) 150,150,170 FFG 101
150 DO 160 I3=I1,IP3,IP2 FFG 102
J0=I3 FFG 103
J1=I3+IP1 FFG 104
TEMPR=WMINR*DATA(J1)-WMINI*DATA(J1+1) FFG 105
TEMPI=WMINR*DATA(J1+1)+WMINI*DATA(J1) FFG 106
DATA(J1)=DATA(J0)-TEMPR FFG 107
DATA(J1+1)=DATA(J0+1)-TEMPI FFG 108
DATA(J0)=DATA(J0)+TEMPR FFG 109
160 DATA(J0+1)=DATA(J0+1)+TEMPI FFG 110
GO TO 220 FFG 111
170 IWMAX=IP0*IFCUR FFG 112
DO 210 I3=I1,IP3,IP2 FFG 113
I2MAX=I3+IP2-IP1 FFG 114
WR=WMINR FFG 115
WI=WMINI FFG 116
DO 200 IWORK=1,IWMAX,IP0 FFG 117
I2=I2MAX FFG 118
SUMR=DATA(I2) FFG 119
SUMI=DATA(I2+1) FFG 120
180 I2=I2-IP1 FFG 121
TEMPR=SUMR FFG 122
SUMR=WR*SUMR-WI*SUMI+DATA(I2) FFG 123
SUMI=WR*SUMI+WI*TEMPR+DATA(I2+1) FFG 124
IF (I2-I3) 190,190,180 FFG 125
190 WORK(IWORK)=SUMR FFG 126
WORK(IWORK+1)=SUMI FFG 127
TEMPR=WR FFG 128
WR=WR*ROOTR-WI*ROOTI+WR FFG 129
200 WI=TEMPR*ROOTI+WI*ROOTR+WI FFG 130
IWORK=1 FFG 131
DO 210 I2=I3,I2MAX,IP1 FFG 132
DATA(I2)=WORK(IWORK) FFG 133
DATA(I2+1)=WORK(IWORK+1) FFG 134
210 IWORK=IWORK+IP0 FFG 135
220 TEMPR=WMINR FFG 136
WMINR=WMINR*WSTPR-WMINI*WSTPI+WMINR FFG 137
230 WMINI=TEMPR*WSTPI+WMINI*WSTPR+WMINI FFG 138
IP1=IP2 FFG 139
GO TO 130 FFG 140
240 RETURN FFG 141
END FFG 142-