|
|
tarticle-tgtimes-Relics-of-Fast-Fourrier-Transform.mw - tgtimes - The Gopher Times |
|
|
 |
git clone git://bitreich.org/tgtimes git://enlrupgkhuxnvlhsf6lc3fziv5h2hhfrinws65d7roiv6bfj7d652fid.onion/tgtimes (git://bitreich.org) |
|
|
 |
Log |
|
|
|
Files |
|
|
|
Refs |
|
|
|
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. |
|