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