      SUBROUTINE FOURT (DATA,N,NDIM,ISIGN,IFORM,WORK,NWORK) 
C     * APR 21/83 - R.LAPRISE.
C     COOLEY-TUKEY FAST FOURIER TRANSFORM IN USASI BASIC FORTRAN. 
C     MULTI-DIMENSIONAL TRANSFORM, DIMENSIONS OF ARBITRARY SIZE,
C     COMPLEX OR REAL DATA.  N POINTS CAN BE TRANSFORMED IN TIME
C     PROPORTIONAL TO N*LOG(N), WHEREAS OTHER METHODS TAKE N**2 TIME. 
C     FURTHERMORE, LESS ERROR IS BUILT UP.  WRITTEN BY NORMAN BRENNER 
C     OF MIT LINCOLN LABORATORY, JUNE 1968. 
C 
C     DIMENSION DATA(N(1),N(2),...),TRANSFORM(N(1),N(2),...),N(NDIM)
C     TRANSFORM(K1,K2,...) = SUM(DATA(J1,J2,...)*EXP(ISIGN*2*PI*SQRT(-1)
C     *((J1-1)*(K1-1)/N(1)+(J2-1)*(K2-1)/N(2)+...))), SUMMED FOR ALL
C     J1 AND K1 FROM 1 TO N(1), J2 AND K2 FROM 1 TO N(2), ETC. FOR ALL
C     NDIM SUBSCRIPTS.  NDIM MUST BE POSITIVE AND EACH N(IDIM) MAY BE 
C     ANY INTEGER.  ISIGN IS +1 OR -1.  LET NTOT = N(1)*N(2)... 
C     ...*N(NDIM).  THEN A -1 TRANSFORM FOLLOWED BY A +1 ONE
C     (OR VICE VERSA) RETURNS NTOT TIMES THE ORIGINAL DATA. 
C     IFORM = 1, 0 OR -1, AS DATA IS COMPLEX, REAL OR THE 
C     FIRST HALF OF A COMPLEX ARRAY.  TRANSFORM VALUES ARE
C     RETURNED TO ARRAY DATA.  THEY ARE COMPLEX, REAL OR
C     THE FIRST HALF OF A COMPLEX ARRAY, AS IFORM = 1, -1 OR 0. 
C     THE TRANSFORM OF A REAL ARRAY (IFORM = 0) DIMENSIONED N(1) BY N(2)
C     BY ... WILL BE RETURNED IN THE SAME ARRAY, NOW CONSIDERED TO
C     BE COMPLEX OF DIMENSIONS N(1)/2+1 BY N(2) BY ....  NOTE THAT IF 
C     IFORM = 0 OR -1, N(1) MUST BE EVEN, AND ENOUGH ROOM MUST BE 
C     RESERVED.  THE MISSING VALUES MAY BE OBTAINED BY COMPLEX CONJU- 
C     GATION.  THE REVERSE TRANSFORMATION, OF A HALF COMPLEX ARRAY
C     DIMENSIONED N(1)/2+1 BY N(2) BY ..., IS ACCOMPLISHED SETTING IFORM
C     TO -1.  IN THE N ARRAY, N(1) MUST BE THE TRUE N(1), NOT N(1)/2+1. 
C     THE TRANSFORM WILL BE REAL AND RETURNED TO THE INPUT ARRAY. 
C     WORK IS A ONE-DIMENSIONAL COMPLEX ARRAY USED FOR WORKING STORAGE. 
C     ITS LENGTH, NWORK, NEED NEVER BE LARGER THAN THE LARGEST N(IDIM)
C     AND FREQUENTLY MAY BE MUCH SMALLER.  FOURT COMPUTES THE MINIMUM 
C     LENGTH WORKING STORAGE REQUIRED AND CHECKS THAT NWORK IS AT LEAST 
C     AS LONG.  THIS MINIMUM LENGTH IS COMPUTED AS SHOWN BELOW. 
C 
C     FOR EXAMPLE-- 
C     DIMENSION DATA(1960),WORK(10) 
C     COMPLEX DATA,WORK 
C     CALL FOURT(DATA,1960,1,-1,+1,WORK,10) 
C 
C     THE MULTI-DIMENSIONAL TRANSFORM IS BROKEN DOWN INTO ONE-DIMEN-
C     SIONAL TRANSFORMS OF LENGTH N(IDIM).  THESE ARE FURTHER BROKEN
C     DOWN INTO TRANSFORMS OF LENGTH IFACT(IF), WHERE THESE ARE THE 
C     PRIME FACTORS OF N(IDIM).  FOR EXAMPLE, N(1) = 1960, IFACT(IF) =
C     2, 2, 2, 5, 7 AND 7.  THE RUNNING TIME IS PROPORTIONAL TO NTOT *
C     SUM(IFACT(IF)), THOUGH FACTORS OF TWO AND THREE WILL RUN ESPE-
C     CIALLY FAST.  NAIVE TRANSFORM PROGRAMS WILL RUN IN TIME NTOT**2.
C     ARRAYS WHOSE SIZE NTOT IS PRIME WILL RUN MUCH SLOWER THAN THOSE 
C     WITH COMPOSITE NTOT.  FOR EXAMPLE, NTOT = N(1) = 1951 (A PRIME),
C     RUNNING TIME WILL BE 1951*1951, WHILE FOR NTOT = 1960, IT WILL
C     BE 1960*(2+2+2+5+7+7), A SPEEDUP OF EIGHTY TIMES.  NAIVE CALCUL-
C     ATION WILL RUN BOTH IN THE SLOWER TIME.  IF AN ARRAY IS OF
C     INCONVENIENT LENGTH, SIMPLY ADD ZEROES TO PAD IT OUT.  THE RESULTS
C     WILL BE INTERPOLATED ACCORDING TO THE NEW LENGTH (SEE BELOW). 
C 
C     A FOURIER TRANSFORM OF LENGTH IFACT(IF) REQUIRES A WORK ARRAY 
C     OF THAT LENGTH.  THEREFORE, NWORK MUST BE AS BIG AS THE LARGEST 
C     PRIME FACTOR.  FURTHER, WORK IS NEEDED FOR DIGIT REVERSAL-- 
C     EACH N(IDIM) (BUT N(1)/2 IF IFORM = 0 OR -1) IS FACTORED SYMMETRI-
C     CALLY, AND NWORK MUST BE AS BIG AS THE CENTER FACTOR.  (TO FACTOR 
C     SYMMETRICALLY, SEPARATE PAIRS OF IDENTICAL FACTORS TO THE FLANKS, 
C     COMBINING ALL LEFTOVERS IN THE CENTER.)  FOR EXAMPLE, N(1) = 1960 
C     =2*2*2*5*7*7=2*7*10*7*2, SO NWORK MUST AT LEAST MAX(7,10) = 10. 
C 
C     AN UPPER BOUND FOR THE RMS RELATIVE ERROR IS GIVEN BY GENTLEMAN 
C     AND SANDE (3)-- 3 * 2**(-B) * SUM(F**1.5), WHERE 2**(-B) IS THE 
C     SMALLEST BIT IN THE FLOATING POINT FRACTION AND THE SUM IS OVER 
C     THE PRIME FACTORS OF NTOT.
C 
C     IF THE INPUT DATA ARE A TIME SERIES, WITH INDEX J REPRESENTING
C     A TIME (J-1)*DELTAT, THEN THE CORRESPONDING INDEX K IN THE
C     TRANSFORM REPRESENTS THE FREQUENCY (K-1)*2*PI/(N*DELTAT), WHICH 
C     BY PERIODICITY, IS THE SAME AS FREQUENCY -(N-K+1)*2*PI/(N*DELTAT).
C     THIS IS TRUE FOR N = EACH N(IDIM) INDEPENDENTLY.
C 
C     REFERENCES--
C     1.  COOLEY, J.W. AND TUKEY, J.W., AN ALGORITHM FOR THE MACHINE
C     CALCULATION OF COMPLEX FOURIER SERIES.  MATH. COMP., 19, 90,
C     (APRIL 1967), 297-301.
C     2.  RADER, C., ET AL., WHAT IS THE FAST FOURIER TRANSFORM, IEEE 
C     TRANSACTIONS ON AUDIO AND ELECTROACOUSTICS, AU-15, 2 (JUNE 1967). 
C     (SPECIAL ISSUE ON THE FAST FOURIER TRANSFORM AND ITS APPLICATIONS)
C     3.  GENTLEMAN, W.M. AND SANDE, G., FAST FOURIER TRANSFORMS--
C     FOR FUN AND PROFIT.  1966 FALL JOINT COMP. CONF., SPARTAN BOOKS,
C     WASHINGTON, 1966. 
C     4.  GOERTZEL, G., AN ALGORITHM FOR THE EVALUATION OF FINITE 
C     TRIGONOMETRIC SERIES.  AM. MATH. MO., 65, (1958), 34-35.
C     5.  SINGLETON, R.C., A METHOD FOR COMPUTING THE FAST FOURIER
C     TRANSFORM WITH AUXILIARY MEMORY AND LIMITED HIGH-SPEED STORAGE. 
C     IN (2). 
      IMPLICIT REAL (A-H,O-Z),
     +INTEGER (I-N)
      DIMENSION DATA(1), N(1), WORK(1), IFSYM(32), IFCNT(10), IFACT(32) 
      IF (IFORM) 10,10,40 
 10   IF (N(1)-2*(N(1)/2)) 20,40,20 
 20   WRITE (6,30) IFORM,(N(IDIM),IDIM=1,NDIM)
 30   FORMAT ('0ERROR IN FOURT.  IFORM = ',I2, 
     1 ' (REAL OR HALF-CMPLEX)',
     2 ', BUT N(1) IS NOT EVEN.',/,' DIMENSIONS = ',20I5) 
      RETURN
 40   NTOT=1
      DO 50 IDIM=1,NDIM 
 50   NTOT=NTOT*N(IDIM) 
      NREM=NTOT 
      IF (IFORM) 60,70,70 
 60   NREM=1
      NTOT=(NTOT/N(1))*(N(1)/2+1) 
C     LOOP OVER ALL DIMENSIONS. 
 70   DO 230 JDIM=1,NDIM
      IF (IFORM) 80,90,90 
 80   IDIM=NDIM+1-JDIM
      GO TO 100 
 90   IDIM=JDIM 
      NREM=NREM/N(IDIM) 
 100  NCURR=N(IDIM) 
      IF (IDIM-1) 110,110,140 
 110  IF (IFORM) 120,130,140
 120  CALL FIXRL (DATA,N(1),NREM,ISIGN,IFORM) 
      NTOT=(NTOT/(N(1)/2+1))*N(1) 
 130  NCURR=NCURR/2 
 140  IF (NCURR-1) 190,190,150
C     FACTOR N(IDIM), THE LENGTH OF THIS DIMENSION. 
 150  CALL FACTR (NCURR,IFACT,NFACT)
      IFMAX=IFACT(NFACT)
C     ARRANGE THE FACTORS SYMMETRICALLY FOR SIMPLER DIGIT REVERSAL. 
      CALL SMFAC (IFACT,NFACT,ISYM,IFSYM,NFSYM,ICENT,IFCNT,NFCNT) 
      IFMAX=MAX(IFMAX,ICENT) 
      IF (IFMAX-NWORK) 180,180,160
 160  WRITE (6,170) NWORK,IDIM,NCURR,ICENT,(IFACT(IF),IF=1,NFACT) 
 170  FORMAT ('0ERROR IN FOURT.  NWORK = ',I4,' IS TOO SMALL FOR N(', 
     XI1,') = ',I5,', WHOSE CENTER = ',I4,', AND WHOSE PRIME FACTORS 
     XARE--'/(1X,20I5))
      RETURN
 180  NPREV=NTOT/(N(IDIM)*NREM) 
C     DIGIT REVERSE ON SYMMETRIC FACTORS, FOR EXAMPLE 2*7*6*7*2.
      CALL SYMRV (DATA,NPREV,NCURR,NREM,IFSYM,NFSYM)
C     DIGIT REVERSE THE ASYMMETRIC CENTER, FOR EXAMPLE, ON 6 = 2*3. 
      CALL ASMRV (DATA,NPREV*ISYM,ICENT,ISYM*NREM,IFCNT,NFCNT,WORK) 
C     FOURIER TRANSFORM ON EACH FACTOR, FOR EXAMPLE, ON 2,7,2,3,7 AND 2.
      CALL COOL (DATA,NPREV,NCURR,NREM,ISIGN,IFACT,WORK)
 190  IF (IFORM) 200,210,230
 200  NREM=NREM*N(IDIM) 
      GO TO 230 
 210  IF (IDIM-1) 220,220,230 
 220  CALL FIXRL (DATA,N(1),NREM,ISIGN,IFORM) 
      NTOT=NTOT/N(1)*(N(1)/2+1) 
 230  CONTINUE
      RETURN
      END 
