SUBROUTINE FOUR1(DATA,NN,ISIGN) C THE COOLEY-TUKEY FAST ROURIER TRANSFORM IN USASI BASIC FORTRAN C TRANSFORM(J) = SUM(DATA(I)+W**((I-1)*(J-1)). WHERE I AND J RUN C FROM 1 TO NN AND W = EXP(ISIGN*2*PI+SQRT(-1)/NN). DATA IS ONE- C DIMENSIONAL COMPLEX ARRAY (I.E.: THE REAL AND IMAGINARY PARTS OF C THE DATA ARE LOCATE IMMEDIATELY ADJACENT IN STORAGE, SUCH AS C FORTRAN IV PLACES THEM) WHOSE LENGTH NN IS A POWER OF TWO. ISIGN C IS +1 OR -1, GIVING THE SIGN OF THE TRANSFORM, TRANSFORM VALUES C ARE RETURNED IN ARRAY DATA, REPLACING THE INPUT DATA. THE TIME IS C PROPORTIONAL TO N*LOG2(N), RATHER THAN THE USUAL N**2. WRITTEN BY C NORMAN BRENNER, JUNE 1967, THIS IS THE SHOURTEST VERSION C OF FFT KNOWN THE THE AUTHOR, AND IS INTENDED MAINLY FOR C DEMONSTRATION. PROGRAMS FOUR2 AND FOURT ARE AVAILABLE THAT RUN C TWICE AS FAST AND OPERATE ON MULTIDIMENSIONAL ARRAYS WHOSE C DIMENSIONS ARE NOT RESTRICTED TO POWERS OR TWO. (LOOKING UP SINES C AND COSINES IN A TABLE WILL CUT RUNNING TIME OF FOUR1 BY A THIRD.) C SEE-- IEEE AUDIO TRANSACTIONS (JUNE 1967), SPECIAL ISSUE ON FFT. DIMENSION DATA(1) N=2*NN J=1 DO 5 I=1,N,2 IF(I-J)1,2,2 1 TEMPR=DATA(J) TEMPI=DATA(J+1) DATA(J)=DATA(I) DATA(J+1)=DATA(I+1) DATA(I)=TEMPR DATA(I+1)=TEMPI 2 M=N/2 3 IF(J-M)5,5,4 4 J=J-M M=M/2 IF(M-2)5,3,3 5 J=J+M C FFT ;) MMAX=2 6 IF(MMAX-N)7,9,9 7 ISTEP=2*MMAX DO 8 M=1,MMAX,2 THETA=3.1415926535*FLOAT(ISIGN*(M-1))/FLOAT(MMAX) WR=COS(THETA) WI=SIN(THETA) DO 8 I=M,N,ISTEP J=I+MMAX TEMPR=WR*DATA(J)-WI*DATA(J+1) TEMPI=WR*DATA(J+1)+WI*DATA(J) DATA(J)=DATA(I)-TEMPR DATA(J+1)=DATA(I+1)-TEMPI DATA(I)=DATA(I)+TEMPR 8 DATA(I+1)=DATA(I+1)+TEMPI MMAX=ISTEP GO TO 6 9 RETURN END C ------------------------------------------------------------------- DIMENSION DATA(1) N = 2 * NN J = 1 DO 5 I=1,N,2 IF(I-J)1,2,2 1 TEMPR = DATA(J) TEMPI = DATA(J+1) DATA(J) = DATA(I) DATA(J+1) = DATA(I+1) DATA(I) = TEMPR DATA(I+1) = TEMPI 2 M = N / 2 3 IF(J-M)5,5,4 4 J = J - M M = M / 2 IF(M-2)5,3,3 5 J = J + M MMAX=2 6 IF(MMAX-N) 7, 9, 9 7 ISTEP = 2 * MMAX DO 8 M=1, MMAX, 2 // will equal 2, 4, 8, 16.... THETA = 3.1415926535*FLOAT(ISIGN*(M-1))/FLOAT(MMAX) // M/MMAX will only ever take this up to 0.5 WR = COS(THETA) WI = SIN(THETA) DO 8 I=M, N, ISTEP J = I + MMAX TEMPR = WR * DATA(J) - WI * DATA(J+1) TEMPI = WR * DATA(J+1) + WI * DATA(J) DATA(J) = DATA(I) - TEMPR DATA(J+1) = DATA(I+1) - TEMPI DATA(I) = DATA(I) + TEMPR 8 DATA(I+1) = DATA(I+1) + TEMPI MMAX = ISTEP GO TO 6 9 RETURN c -------------------------------------------------------------------- n = nn*2 ; j = 1; for ( mmax = 2; while(mmax < n) { istep = mmax << 1; for( m = 1; m <= mmax; m+=2) { theta = M_PI*((1-m)/mmax); wr = cos(theta); wi = sin(theta); for (i = m; i <= n; i+=istep) { j = i + mmax; tempr = wr * data[j] - wi * data[j+1]; tempi = wr * data[j+1] + wi * data[j]; data[j] = data[i] - tempr; data[j+1] = data[i+1] - tempi; data[i] += tempr; data[i+1] += tempi; } } mmax = istep; }