C C----------------------------------------------------------------------- C SUBROUTINE: FOUREA C PERFORMS COOLEY-TUKEY FAST FOURIER TRANSFORM C----------------------------------------------------------------------- C SUBROUTINE FOUREA(ID1,ID2,IC,IR,IVN,ISI) C ID1,ID2 = COORDS OF FIRST CELL. IC AND IR ARE 0, OR 1 C ONLY ONE OF IC, IR MAY BE NONZERO. (FLAGS HORIZ/VERTICAL C DATA AREA) C C THE COOLEY-TUKEY FAST FOURIER TRANSFORM IN ANSI FORTRAN C C DATA IS A ONE-DIMENSIONAL COMPLEX ARRAY WHOSE LENGTH, N, IS A C POWER OF TWO. ISI IS +1 FOR AN INVERSE TRANSFORM AND -1 FOR A C FORWARD TRANSFORM. TRANSFORM VALUES ARE RETURNED IN THE INPUT C ARRAY, REPLACING THE INPUT. C TRANSFORM(J)=SUM(DATA(I)*W**((I-1)*(J-1))), WHERE I AND J RUN C FROM 1 TO N AND W = EXP (ISI*2*PI*SQRT(-1)/N). PROGRAM ALSO C COMPUTES INVERSE TRANSFORM, FOR WHICH THE DEFINING EXPRESSION C IS INVTR(J)=(1/N)*SUM(DATA(I)*W**((I-1)*(J-1))). C RUNNING TIME IS PROPORTIONAL TO N*LOG2(N), RATHER THAN TO THE C CLASSICAL N**2. C AFTER PROGRAM BY BRENNER, JUNE 1967. THIS IS A VERY SHORT VERSION C OF THE FFT AND IS INTENDED MAINLY FOR DEMONSTRATION. PROGRAMS C ARE AVAILABLE IN THIS COLLECTION WHICH RUN FASTER AND ARE NOT C RESTRICTED TO POWERS OF 2 OR TO ONE-DIMENSIONAL ARRAYS. C SEE -- IEEE TRANS AUDIO (JUNE 1967), SPECIAL ISSUE ON FFT. C C COMPLEX*16 WC C ASSUMES THAT FIRST N/2 ELEMENTS ARE REAL, SECOND COMPLEX... C STORES DATA THAT WAY ALSO... C C COMPLEX DATA(1) C COMPLEX TEMP, W C MAKE THIS A REAL FFT, NOT COMPLEX... REAL*8 DATA(1),TEMP,W,TEMP2,TEMPI,WI INTEGER*2 ID1,ID2,IC,IR,IRX,IRXX,IVN,N C SET UP STMT FUNCTIONS... ID1F(K)=ID1+IC*(K-1) ID2F(K)=ID2+IR*(K-1) N=IVN C C CHECK FOR POWER OF TWO UP TO 14 C C INITIALLY SAY ALL OK NN = 1 DO 10 I=1,14 M = I NN = NN*2 IF (NN.EQ.N) GO TO 20 IF(NN.GT.N)GOTO 11 10 CONTINUE 11 CONTINUE N=NN/2 C USE NEXT SMALLER POWER OF 2 ARRAY... C RETURN C HERE BEGINNETH ACTUAL WORK. C SET UP DATA COORDS ON THE FLY. NORMALLY I,J RUN IN RANGE 1 TO N C SO WHERE K=(I OR J) (I.E., ONE OF THE TWO) WE USE A RELATION C ID1V=ID1+IC*(K-1) AND ID2V=ID2+IR*(K-1). WE USE STMT FUNCTIONS C ID1F AND ID2F FOR THIS. 20 CONTINUE NOV2=N/2 C C PI = 4.*ATAN(1.) PI=3.14159265358979323846264 FN = NOV2 C C THIS SECTION PUTS DATA IN BIT-REVERSED ORDER C J = 1 DO 80 I=1,NOV2 C C AT THIS POINT, I AND J ARE A BIT REVERSED PAIR (EXCEPT FOR THE C DISPLACEMENT OF +1) C IF(I.GE.J)GOTO 40 C C EXCHANGE DATA(I) WITH DATA(J) IF I.LT.J. C 30 CONTINUE C EXCHANGE DATA(J), DATA(I) CALL XVBLGT(ID1F(J),ID2F(J),TEMP) CALL XVBLGT(ID1F(I),ID2F(I),DATA(1)) CALL XVBLST(ID1F(J),ID2F(J),DATA(1)) CALL XVBLST(ID1F(I),ID2F(I),TEMP) C FLIP BOTH REAL AND COMPLEX PARTS OF DATA CALL XVBLGT(ID1F(J+NOV2),ID2F(J+NOV2),TEMP) CALL XVBLGT(ID1F(I+NOV2),ID2F(I+NOV2),DATA(1)) CALL XVBLST(ID1F(J+NOV2),ID2F(J+NOV2),DATA(1)) CALL XVBLST(ID1F(I+NOV2),ID2F(I+NOV2),TEMP) C 30 TEMP = DATA(J) C DATA(J) = DATA(I) C DATA(I) = TEMP C C IMPLEMENT J=J+1, BIT-REVERSED COUNTER C 40 M = NOV2/2 50 IF (J.LE.M) GOTO 70 60 J = J - M M = (M+1)/2 GO TO 50 70 J = J + M 80 CONTINUE C C NOW COMPUTE THE BUTTERFLIES C MMAX = 1 90 IF (MMAX.GE.NOV2)GOTO 130 100 ISTEP = 2*MMAX DO 120 M=1,MMAX THETA = PI*FLOAT(ISI*(M-1))/FLOAT(MMAX) W = COS(THETA) WI = SIN(THETA) C W = CMPLX(COS(THETA),SIN(THETA)) DO 110 I=M,NOV2,ISTEP J = I + MMAX C GET REAL AND IMAG HALVES OF NUMBER... CALL XVBLGT(ID1F(J),ID2F(J),TEMP) CALL XVBLGT(ID1F(J+NOV2),ID2F(J+NOV2),TEMPI) C DO COMPLEX MULTIPLICATION BY HAND TO AVOID LARGE RUNTIME SYSTEM C ROUTINE INCLUSION. TEMP2=W*TEMP-WI*TEMPI TEMPI=WI*TEMP+W*TEMPI TEMP=TEMP2 C TEMP = W*DATA(J) C DATA(J) = DATA(I) - TEMP C DATA(I) = DATA(I) + TEMP CALL XVBLGT(ID1F(I),ID2F(I),DATA(1)) TEMP2=DATA(1)+TEMP DATA(1)=DATA(1) - TEMP CALL XVBLST(ID1F(J),ID2F(J),DATA(1)) CALL XVBLST(ID1F(I),ID2F(I),TEMP2) C COMPLEX PART CALL XVBLGT(ID1F(I+NOV2),ID2F(I+NOV2),DATA(1)) TEMP2=DATA(1)+TEMPI DATA(1)=DATA(1) - TEMPI CALL XVBLST(ID1F(J+NOV2),ID2F(J+NOV2),DATA(1)) CALL XVBLST(ID1F(I+NOV2),ID2F(I+NOV2),TEMP2) 110 CONTINUE 120 CONTINUE MMAX = ISTEP GO TO 90 130 IF (ISI.LT.0) GOTO 160 C C FOR INV TRANS -- ISI=1 -- MULTIPLY OUTPUT BY 1/N C 140 DO 150 I=1,N C DATA(I) = DATA(I)/FN CALL XVBLGT(ID1F(I),ID2F(I),TEMP) TEMP=TEMP/FN CALL XVBLST(ID1F(I),ID2F(I),TEMP) 150 CONTINUE 160 RETURN END