tarticle-tgtimes-fft-hack.mw - tgtimes - The Gopher Times | |
git clone git://bitreich.org/tgtimes git://enlrupgkhuxnvlhsf6lc3fziv5h2hhfrinws… | |
Log | |
Files | |
Refs | |
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 wri… | |
6 What the hell were they running it on, and what damned data were they fe… | |
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 … | |
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, w… | |
75 It was ported from Fortran to BASIC, then to C, then to C++ and everythi… | |
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 h… | |
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. |