article-rue_mohr-Relics-of-Fast-Fourrier-Transform.mw - tgtimes - The Gopher Ti… | |
git clone git://bitreich.org/tgtimes git://enlrupgkhuxnvlhsf6lc3fziv5h2hhfrinws… | |
Log | |
Files | |
Refs | |
Tags | |
README | |
--- | |
article-rue_mohr-Relics-of-Fast-Fourrier-Transform.mw (3218B) | |
--- | |
1 .SH rue_mohr | |
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 wri… | |
7 What the hell were they running it on, and what damned data were they fe… | |
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 … | |
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, w… | |
84 It was ported from Fortran to BASIC, then to C, then to C++ and everythi… | |
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 h… | |
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. |