Google
 

Trailing-Edge - PDP-10 Archives - decus_20tap2_198111 - decus/20-0051/ctfft.for
There is 1 other file named ctfft.for in the archive. Click here to see a list.
      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-