tarticle-tgtimes-Relics-of-Fast-Fourrier-Transform.mw - tgtimes - The Gopher Times
 (HTM) git clone git://bitreich.org/tgtimes git://enlrupgkhuxnvlhsf6lc3fziv5h2hhfrinws65d7roiv6bfj7d652fid.onion/tgtimes
 (DIR) Log
 (DIR) Files
 (DIR) Refs
 (DIR) Tags
       ---
       tarticle-tgtimes-Relics-of-Fast-Fourrier-Transform.mw (3217B)
       ---
            1 .SH tgtimes
            2 Relics of Fast Fourrier Transform
            3 .2C 3v
            4 .
            5 .PP
            6 In 1967, the Kooley-Tukey FFT algorythm (the one we all use now) was written in Fortran.
            7 What the hell were they running it on, and what damned data were they feeding into it?!
            8 .
            9 .1C
           10 .
           11 .DS
           12       SUBROUTINE FOUR1(DATA,NN,ISIGN)
           13 C     THE COOLEY-TUKEY FAST ROURIER TRANSFORM IN USASI BASIC FORTRAN
           14 C     TRANSFORM(J) = SUM(DATA(I)+W**((I-1)*(J-1)). WHERE I AND J RUN
           15 C     FROM 1 TO NN AND W = EXP(ISIGN*2*PI+SQRT(-1)/NN). DATA IS ONE-
           16 C     DIMENSIONAL COMPLEX ARRAY (I.E.: THE REAL AND IMAGINARY PARTS OF
           17 C     THE DATA ARE LOCATE IMMEDIATELY ADJACENT IN STORAGE, SUCH AS
           18 C     FORTRAN IV PLACES THEM) WHOSE LENGTH NN IS A POWER OF TWO. ISIGN
           19 C     IS +1 OR -1, GIVING THE SIGN OF THE TRANSFORM, TRANSFORM VALUES
           20 C     ARE RETURNED IN ARRAY DATA, REPLACING THE INPUT DATA. THE TIME IS
           21 C     PROPORTIONAL TO N*LOG2(N), RATHER THAN THE USUAL N**2. WRITTEN BY
           22 C     NORMAN BRENNER, JUNE 1967, THIS IS THE SHOURTEST VERSION
           23 C     OF FFT KNOWN THE THE AUTHOR, AND IS INTENDED MAINLY FOR
           24 C     DEMONSTRATION. PROGRAMS FOUR2 AND FOURT ARE AVAILABLE THAT RUN
           25 C     TWICE AS FAST AND OPERATE ON MULTIDIMENSIONAL ARRAYS WHOSE
           26 C     DIMENSIONS ARE NOT RESTRICTED TO POWERS OR TWO. (LOOKING UP SINES
           27 C     AND COSINES IN A TABLE WILL CUT RUNNING TIME OF FOUR1 BY A THIRD.)
           28 C     SEE-- IEEE AUDIO TRANSACTIONS (JUNE 1967), SPECIAL ISSUE ON FFT.
           29       DIMENSION DATA(1)
           30       N=2*NN
           31       J=1
           32       DO 5 I=1,N,2
           33       IF(I-J)1,2,2
           34 1     TEMPR=DATA(J)
           35       TEMPI=DATA(J+1)
           36       DATA(J)=DATA(I)
           37       DATA(J+1)=DATA(I+1)
           38       DATA(I)=TEMPR
           39       DATA(I+1)=TEMPI
           40 2     M=N/2
           41 3     IF(J-M)5,5,4
           42 4     J=J-M
           43       M=M/2
           44       IF(M-2)5,3,3
           45 5     J=J+M     
           46       MMAX=2
           47 6     IF(MMAX-N)7,9,9
           48 7     ISTEP=2*MMAX
           49       DO 8 M=1,MMAX,2
           50       THETA=3.1415926535*FLOAT(ISIGN*(M-1))/FLOAT(MMAX)
           51       WR=COS(THETA)
           52       WI=SIN(THETA)
           53       DO 8 I=M,N,ISTEP
           54       J=I+MMAX
           55       TEMPR=WR*DATA(J)-WI*DATA(J+1)
           56       TEMPI=WR*DATA(J+1)+WI*DATA(J)
           57       DATA(J)=DATA(I)-TEMPR
           58       DATA(J+1)=DATA(I+1)-TEMPI
           59       DATA(I)=DATA(I)+TEMPR
           60 8     DATA(I+1)=DATA(I+1)+TEMPI
           61       MMAX=ISTEP
           62       GO TO 6
           63 9     RETURN
           64       END
           65 .DE
           66 .
           67 .1C
           68 .
           69 .PP
           70 And no, you \fBcannot\fR get the IEEE document because IEEE broke it up into pages and sells each page individually.
           71 .
           72 .1C
           73 .
           74 .DS
           75 "PROGRAMS FOUR2 AND FOURT ARE AVAILABLE THAT RUN
           76 C     TWICE AS FAST AND OPERATE ON MULTIDIMENSIONAL ARRAYS WHOSE
           77 C     DIMENSIONS ARE NOT RESTRICTED TO POWERS OR TWO."
           78 .DE
           79 .
           80 .2C 15v
           81 .
           82 .PP
           83 But, this code was easy to port because it was small, so, to this day, we use it.
           84 It was ported from Fortran to BASIC, then to C, then to C++ and everything else.
           85 .
           86 .PP
           87 Nobody ever actually understood it, so they didn't fix anything.
           88 You see, Fortran has no bitwise operateors, so alot of the acrobatics
           89 in that code are just doing bitwise operations in regular math.
           90 Its absolutely amazing when you tear it apart.
           91 .
           92 .PP
           93 I got the code from a bad scan of a document off a military ftp site.
           94 What I love, and find halarious, is that this code has been ported and hacked a million times since it was written.
           95 .
           96 .PP
           97 But, from the comments, it, itself, is a hack.
           98 It is a mash up of cooley and tukeys code.
           99 It is a hack, from 1967.