tarticle-tgtimes-fft-hack.mw - tgtimes - The Gopher Times
HTML 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.