      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;
 }
 
 
 
 
