tarticle-tgtimes-Relics-of-Fast-Fourrier-Transform.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-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.