|
|
tarticle-tgtimes-fft-hack.mw - tgtimes - The Gopher Times |
|
|
 |
git clone git://bitreich.org/tgtimes git://enlrupgkhuxnvlhsf6lc3fziv5h2hhfrinws65d7roiv6bfj7d652fid.onion/tgtimes (git://bitreich.org) |
|
|
 |
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 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. |
|