ycy@walt.cc.utexas.edu (Joseph Yip) (06/14/90)
Hi, I am looking for C source for Fast Fourier Transform (any dimensions and power of 2). The netlib only has the Fortran version and I do not want to use f2c. If you have a good FFT program (fast!), that is excellent. Thank you. Joseph Yip
ycy@walt.cc.utexas.edu (Joseph Yip) (06/14/90)
Hi all those interested people, Most people recommended me to look at "Numerical Recipe in C". Two people sent me their FFT source. I am going to share with you guys what I got. Here is, From black@seismo.CSS.GOV Wed Jun 13 18:16:08 1990 Received: from seismo.CSS.GOV by walt.cc.utexas.edu (5.61/1.34/CCWF 1.5) id AA00510; Wed, 13 Jun 90 18:15:53 -0500 Received: from beno.CSS.GOV by seismo.CSS.GOV (5.61/1.14) id AA26957; Wed, 13 Jun 90 19:15:46 -0400 Received: by beno.CSS.GOV (4.0/SMI-4.0) id AA06806; Wed, 13 Jun 90 19:15:45 EDT Date: Wed, 13 Jun 90 19:15:45 EDT From: black@seismo.CSS.GOV (Mike Black) Message-Id: <9006132315.AA06806@beno.CSS.GOV> To: ycy@ccwf.cc.utexas.edu Subject: Re: C code for FFT Newsgroups: sci.math,comp.sources.wanted,comp.lang.c In-Reply-To: <31567@ut-emx.UUCP> Organization: Center for Seismic Studies, Arlington, VA Cc: Status: R In article <31567@ut-emx.UUCP> you write: >Hi, > >I am looking for C source for Fast Fourier Transform (any dimensions and >power of 2). The netlib only has the Fortran version and I do not want >to use f2c. Here you go, a uuencoded compressed tar file: ----CUT HERE---- begin 644 fft.tar.Z M'YV.9LS0<3$&@,&#"!,J7,BPH<.'$",B!$'1!@T:( !0C!%CQHV,%"G.@ $R MY,89-&R0!)%R1HT:-F3$H!$C9(T8,&)HE,BSI\^?0(,*'4JTJ-&C2),JA?A" MA0(0*A($'#CFJ5.H":JX28,'A)4R<N:D>>,&A P7&,7D 3&%3AD[9=B&:0-G M#ED6(*#4$<,FS1@01-ZT"9/&#=XI9>"X;2,&+(@8.7#@L&HU 14T:>: @"/G MS1DY<S=W)E-G3!G-84 8D5,F3ITR;L:L#3RX, @RF>&P";/63.<V(.B@B4LE M39NXM F[J4PF#)W4A>'4H;,#1!VQ;LX$'ZXZS!PZJM_4D9/&,170;N:8>2.G MC0O*6"_'E5(F#)O;SJ%K/@,;K/,R9("@UG9Q%=8"'&AX%Q<*2:0 PA@)NN%& M&6S@%88; 0I71F5)#'9&86'(L59SSX&0F6AOD&9:@ -J"()K89 !&AWCE7%@ M@G,L&(6#$%XXX7UOF%%9:H'!H9MC4L08(@BLF98&7'*X0)%\(+21!AEDL!&7 M&6FPX98<)FJF!UAO@("$'E*", 0;;V!W1F7/R<$?':BQ=EMNN^4!8'!ENDC> M&6C082&&(#P!!VR%:1?GG)K1T2=W6@KT'E25.?&&6SJ <$4:PH'081@?NA'B MB/F!D"-XCH(P9F<$@B >'=*!=V)J5LK1&9C&@5K&I%=%5<4<NNH@E4"FSG6D M9M%-Y^ITL5:60 Q2 DN7ELBB!H(=(:81!E]QP?'&'8X%&=P=;_":0 )G>>I& MK&&"@(*TQT)E:AIC!HD"&>)QFP*/"8(VQI=SF)O #%(^P:RR)[YK++7RBE6O M&??FJ^6^#_8;QK]@!>PL#5*V=6%S<@0(UJWM#L<&'*Z"24898M1Q1JC:'3<' ML/QI3"E6B,5EGUV9LA;C<2[@.P9>TK&1(T&#!HC=KF/P^H(""C2E[AALU+&R M9BH\K< (A5%M=5P\?(=;N6CXL'4:9JQL!@A!:$%%$4-\,<39L56],@AATT%> M=@&7O36%.=+M]=T\#,8&FTW[/0)LN G)==U?XSV8<"[X';4*@)7!Y818:SV" MVH658=D53WP!11()H' 629@3/%,,->0@@PTUN(1##C?D,(,,-*20@-3V::O9 M#"K=5H9V66\-^H0)#(%&'6ZLD?H4"Y_6,+UEV(OO7A.G /7GFH>>P-(HX.%[ M E.$OD7Y*72A?/C,C]%FZ.RCKSX*[(.P@KO43VN]#SXP"\5*(!?_S<%]EP/! M$=@DAIWAQ0CLN4.( K0RJH7(.6-)#U2T]AT,_J4P=$!?]>: %RA\"RP[4$ ' MZ>"7VTA,="J@CWWPHH)/G0$O30@#'@Q3)5!MA48K0T$*4KC"%MKA#5<:WQCL M4P8AXF4J3NPA'KY@K"CB9@ZZX4T41_8%UGA+#G008@I5^!P6?M -(3Q4>Z;3 MQ"&2T8,NY!X,TS<A(I;1B$@D0P):UB4R?$%O:3B#&*%6Q+\8(0E,*$("5& $ M.*2!AHU\PQ@32"7.> 8TP$F> FKC!A2$Z Q# \$G[> ]$";@DV-(88_DL$@5 MC#*%>U! M#F+E2"( 0]8(GWSL5%+[(GC#!PHRS[%Z]<;FL.GG14&E"@@A6L M8)3[2F$"3 @N,.6236? "?ZVIR]B4LM!+P !-K6Y.F'^#G-!.-P;ENB6*I6A M#>SI#7M:A2WR;$M+HK05;SH'M5FN#7\RO$\N(R9'J*3 <(A#@?.@MP:*]2"7 M3J@"$YBPRP3T,C&_1$$,S$E+_-D0! /E)CY5<%#[)'2AT7,H1"5*45E:U%9= MQ"@85<?1?Z* CG$)Z0L-BM!U*O1Y*7700T$0T8E6]*)?#.,,S"FU-QRJ+"XB M46JXQ#!-^M-=*&AD&D *@O4\E9G.'"5>1" '$:B4J"P]*DQ].=/>C?&J^(LD M5[T*&[ ^4TYV&.L=S"K4E1K5I4C-: V8BCF?!<@V42V5;P33JF2!AZIQL:IO MZD,&% 04+P[+'L1$6H84'&:$#W2D&Z5B6 9Y"+/8T]X+/5O 8X4V#:,U ]7: MU$2MFI./;/ C( 4YVCDL44N#'%881WO%+.9AD+*4VAW(TTX73=8UL)&-"SF) MGQ*ERD4=Q- $EP4K94EVN9QJXF7G]3""ZNNS!@SM&V(+7K>8%E2H+2]G6>O- MTZ@WMK/-45:]1=BNLB8N=4#9<> I(E$22DMT.,%^R@">!HYA#7Q:%IC,$(8C M,C<N<S#9?23[7_RMTJ !C6V'4?!ADMI0Q&5H(HG[95"<FK,,>. 4"H*9PCY M[8UFO%8>01#8F;JAE&A,@!M@*<LYW(%3$'+7CRD2RP0L,4<@@($.7&H&SH 0 M8F(;V5A_%:RN$@M>#"O,%R#+W3%WJ0Q<<(-9I2D5*Z,1RW18F:W&>H7AV G, MUIO59DX()G'1@5QI7K-+Q>"S-4CSR7&)P93/5>6]T0'.<I;#6*>2*8.!1UP# MCB<(YL3"[#"I/O?!XL70K.;1[K'0AU:061;=9D=#6LL@$ &E"Z4L3+]3TYQ. ME(D\9"HXC#K0IIZLBDN<@A"SF=#U,;0L$0V"&;"ZT5=^5YQA+6N!5+K6:\MT M@7/MZ:7U^M>E9K.P/<QBDAJ;RB,F]HF/C>IEJYH&SW;SHZ4=:194FPZ9TLH] MXY(JI\*FNM!9E[(@"VQQI[O<Q09UL ^^)!-[:.$I)G?#4^!B=B<[U5"N0;Q= M36\MWSO?HN)6A/U=%JEREUT$#S>5\UM;T1H\XBN>^+D9S7 P.1Q4$!\VPBL^ MZ'8GP,:RA+&,-UKC&Q=2QU<RU6_;J( FL^9#WW&,*=GDU+>N!TPHH#K*<@F# MZF@=;ZVEEM??X%1G.@@$34Y 0+>@]2Z X 6YK.\<V&Q#MI,=#FZ'>]A/(TT; M QW'>$PZ%+WG].-EYDLF"K+6\:+U&##^[C*0)O3$<H8)Z3$-,=C!N<ZY:21Z M.E5_"I2I,$,LX9PH>;,,\BQE@)<TS*#U-,!+'J3)V7.%P?&B9+V <"\&W>L! M]WJ(?#\QSU6Y@P" CY%F&F3 U<SW4VIJO#IP7-0R.GS)#&S(@X(WV,^KN^OK M7!^[4\%.312*\^YF9[)+7<_5KJ\?([G$_%NE,L^LW[TFX3^_4VO" Q-%7O]P M$ /IAW8NE0!KD4MJU 9LA +L%TX,B$8I@'FFE@"_QU4/,C\3@@)Y,('!QU4M M0#X;.'_TAW5:QWSQ-P/B!P?,UW]I0 ,I* ,#&$L%>$KX!P)K5X)=P&:V9X*> MXB%VYU0RD(/]M'EB4!-<]7LQ!&H_J((K@'ENUP*J(@,UY(,EV(0Q((2;MT?, M=X12>(.09X5NMS](6'=5Z(0BJ'9*B(-<=7OZLWLZF !D"'EN9TS,MS^]=X9> M"(1@N(8U 85%^(9QJ(=.N(;,YX?"MWE MWGLMP(#17P\T']$5X MJ#^-R'^0 M:&J).$LU 4"YY'SGLGQXPP.=V'<* '1-41E$@">\T2K/]1JQ,2*"H1R\4AE4 M0E4 DT\Z@Q[\H32/XAR ,03M@EC<8258@D^V"!;F<AG7 0(M4"1'@BN:01:M M$BF"HC_.J"5]!E7<$7H#814<=$=_<41)5UQY(D1-)TM/=WA2IWAWAQ=!8@:G MHGRJQT-X\%:#,455-%IE$7<C='P!Y'RRY'WV-WZY-&0 "'9REX(QZ%)=T8@0 MB (>\D-?,Y!PX""8$SLN0!+AE$,[9&KR!C$B4 )DP 5TP <B\'CC!X4_QF;O M>"KMQV9W@!GXA +.U)(,AC>Y9#X%&"NKA (GT ,G8$X)P)/]XI-I%I2D&)#U M!WY1EH)@9Y!:MY#GTI#N D(I )$^Q"D3J746"0(8J9$@P)$KZ5(?B0(A.9(E M>9('N3]C>2XV"1[AYU(Q>6;N4I,!X9*B" (ZN7E$&2(^"91"V9=R8)1N@)2R MY'=0<XI8802%$2#W:!QU !P+HY4KXXV$!(Y(YY@Z1$5T88Z%!W6(-W5W1WLO M="X>\E9BV7X9*4T"R91==Y#]EY *)5PA97:D4L129E-Q)4.$D!B65$)D)JX M"2JD:(J]TCSV,09UL!OM5'Y@DD.5IYN6"7A_P5D]%)U U$0_I@"F-&3G*&1< ME8 +V):L02-R4!9HD ?>$D9KYP9=@!=UYYYN9)R5(02<\FE0 F4NX@:1V1AR M !]1X3P7PA^/P1$P<* '*@81AJ 'RA$QH*#S9 =^82/Y&0;3>70@M!E@H8#N MU9@P!F0AY*'U^)WI&'6XPHY5)TN3%TB6)V2>:$H)P!ISL)S4T4]N4(,)B8ZG M0:,OV7U+>7<O"9MYP6<*N0)GEW8WZH^CN'[_)*)/N5$%*#]H5!BO,7]2HQ5< M 0(X(DK!(0?0PTX9A"+< AQ'UBG><@="Q'TN):,\RHA5^9!GJCH9B1?F-3%< MV889.5ALYJ0M4)">F(CE.1YEP:9>4G274QE0P!HM(#]TP4:MXFT?<X'>]AS< M,@<7BIGBR")UT$=_]">>J:.@N8XAI'6D*4>VEQU:TGJQP1K'@49O-:#XQ'6K M*4M=PZJP 9<@0 6D8SI)\';N,E\YJI0D"*3Y]W6Q.4)%>J0NA5-+.(?STDFP MVEEL%JV4F'AC8*NN>IBE"#5+T:W>^JW@&JX1431'4Q#B>JY $1(6@1$:4: > M41(B01+M:A+-)@,'2A$M\1(QP1$801$U0 ,RL1/H.K $6[ &>[ (F[!$H9A1 M0:Y, Z!9L15=\15A$:9GD19KT19O$1?$9!<\I!=\T4+)41B'D1B+X9^/$1F3 M<3-1<1DG8DF?$1HG<AU[DBK\,2$SPF_&$1?X0ET."W CZ&54X33?TS5V S9B M,Q:58S9<DS::PS9N S=R(SA'BS??L3=GT#=,"S@;\CB# S8]E3A;BR%H\SW+ M(SHH)3TWU8^8DUFJ)4?[<IEP5'OU54)W@!?%<1PT-%Y32)S4F7@AE*-'EZDO M)0<QE52#='2'E$B+U$B3Q$F>)">A!$W<&62HI$K]TDJO1**&9Z* FP"D.GS_ M9$NXU&QJ9;AL!4SF9'S&) ;(Y!Q(9%?05%-8%5!<5:>155+J- 8_Q5!G550M MQ4MK)5-A%(FT:E,?I5,%15)AV[M!!5)^%;R%>[@9)0/F5)8BD 0"!QY4HA=& M$Q>;0@:=<E-E(*5D, <.\@.QUEM+Y 8@60)L8 9J:0)0< >CA;W:RRY4XDVZ MEK<+@ACFB[X@H+Z"ID07\K[Q.[_U94[^JYKN%[04&:3&NG=S-YM&JGZ?^$\- MW'_U"YQYB'?--ZN\]+T%^,'."@.SZE*!",*R^L"6L;/56E])69N--%<D)[MX MI5=\!;UH]5?"B[K$BP)+977M)5Z@%E^:A;OT!5JJ49&L6<3O=4/DE<3SA5ZN MU<3XQ2;ZU4C7.QT]>0)':4[C%G,VEW#V@6(Z-W'KUD]"I[J&^K>9RF/#B[C; MV9U$-CYEFF0HL&0$Z&2J)F545I99-F<BP&7\D2D_BV?1>#!T8&9:4G"!S'&# M+&DB4&=@@6']J&=G&BYK\V=O ,GG@FQAH&Q^#&6*%LG1-LGVYK#7=FG9=FO; MQF"=IAT^$VJ^9AJ@?&H7YVY0)@,;E\K31LBL3&NN[$X$MA;<IAVYHAVBALLJ M1W,P1VPSI\NCC'%QX6RH_&8=)\S+F2,@MV\1=JWUT4Y%W%5G!LJ)V,8:-9^) MB3E%\ 1&P'T& ;RG2[TS9;UOA;_;FZO<X;U0%K[C"\!D<;[IN[[2Y%L';);P M*[]X0;_V*TWZK+_<P;^>UL#D&\ %7< ([;X*G< -O<!OU< MS)H_2I!-*:2R M&947W,=PM<%Y8;\EG(9WYZS0XL)<&]/VT:P./'\K?,(B_,+'$</5,\,=M5\V M_%7-=%=GD%>QME>_FU: -<<9M30$H;!6/:_K6A(<\:[SVFSR2J\B ;#]FJ\P M(1.S8Q.\(P,">]5LW=9N_=9P[=8,.SZA0Q 0BZ43FS$6BQ8"DK%N 1>MY;$E MQ#TB&XLDRQ8F^TXH"QF2 ;$NJQDPBTGM0K,94B8WZQ_-!<,]JQS/6@82!-A2 M51G>-Q4$D29)L#9Y(!X@($%H9&"N\E2ZMB@WJ9+'\Q\$ B>G\5BL 5VO2+14 M&SEY,S9+>S9.NS9M\S9Q,S=>6[5Y@[5:^S<DO-S ';;#O3ADZSA&.]W.@0;# M_3Z;(SJZ6CJGDSJK(R^N0Q.Q,SNU4P.WDSN[TSO>'3IJ E1J:WQMFUJ;M5K> M\[=T6ST/M-MXN[-[>\1084-V!$>F)+B8FD?3F[J)BYF+JTB,! >/JQR1"TH6 M@E<@>DJ2B[DAHKEX=<<E&IHH2N&B6TN2>TNYM%11#<2(2V/]Q+JBY+K)%+M) M/;MO5=2VJ[PBQ[PFY5-I^]0^W.!!;+RUF;R_NE,^OKO.VU!]U</2VV-AA,_] M%-'*0B54D-NJL=NN*%TH<"89S;X)'9(>#0(FL!JM<;\<E[W[O+^3Z6EHSMM> M#N8#;-!%UKX(S-!F#M+#5#TA#)8RW$\B'64BW)K$>M(3G-+HM])I9\(6B!TH M$-Z\*B]Q+B_^:VH][8'R@TR2?CJ4OMN6OK.FUL!N&N@_=^)QM75==<,WGL-- MO<-#1<\M;L]*)<90/%YNF]]P:\74$EKL=6%1C,1O>UX4[.NLR7+[U<5T\,5A M;'4U!V(*QYK/?G-G\&(QYL;;*K<Y%L=2OL<;[IU-9F1(A@9*IJRE'!> S&B" M',R4;,AED"G>IL@G-QV.3&H%W&K 7&^5;&>8;$#MDAJ;W&>=#&C/3,VDS&RG MK.Z2S.[VMC2M["JO?,R;)LNZ5LO?YLSW+LH'KVJ^G,WSILHBX/#$#/'&C&L4 M[VG+?/'VGG,25\;3K/'6W&R_K,T@+_+Z)G*I(LZV7<XI5\#I?.WK;*A2\\[Q MC'I%C5.W:YU+?E+T+>11+M4S-<3/ASDDURHF!UG\=+Q8I55'75>MOM1C559. M7\\.[E:HOE]EDDMTU4E?S]0BX-1/+NL_3.LHH*=3_VDQDGA6KUAR$-=^__> M'_B"/_B$7_B&?_B(G_B*O_B,3[!4*(>$V(9W.(3G\L$P.(C&U(=N.(,K?/E7 M&/F&J(.9.$LS4*T,:(F/,8&3Z*:G'XJISV:92'R<F'SKMX)Y^:?9/M>IB$5Y MPHI<'EVPR$FS&!_<<8QA@8OY-*![8A?;X8M$ (PG(HQQ08Q9LB5=\B7)^#R: MT8Q.]8SM(HTN0HUXL0+7&"[:&!?<Z-N#NV/DJ$6$!ZKJ>**CVHX0#X\,)H\A M)&1X,:*R=(^<"0?!91]1,']$^X35]SET4 E('2L#DJPP6 *@2@^(#EREW)2= M*%)7^DJ^ZC>Q&>PEDDB224))*$,EF9JWU*/.Q5R:27;)_H&'O+27SH5@^DN& MB2QYL:($Q@K3:$E$ALZDO:8)E@#+#J-C2+?+*MDFB11$[M1%R@$9*0/JD+8T ME-9<!TQ+(+ -,4$2&)=,H$Q:$"D0+^4DX/0"?U(,=($ST"_50,-TZNC38FI, M4@0R2:9I(9W4U/I+.OXO'W&N4"7_0-=HDB6U)P&<IGXBG C= \N!JFX'*D * MU@!96E'C-<,).VVENT,Q?-,2!$Y]< ^6P796&89 <EI.MLTYA:6LE)TNU=S: M*11P(M6Q( /N9(D %$_NA3PQ&$&EI=33I; L2L@]P2<?))^$WG&R3^"!->0G M?L,=^!-C L02T!EAT1CH!B4@DH5#"HG< 0(!28DE&EH 17*$^:8#(4*FXB( MVG B2L1U+A(W_U)4 EA1E0= N"C\=RX(58TRA3@*6>FH&>4E>I0!C&#%:A Z MIT*(I#:1#UA2&>Q-K8RNT']N%'"24IVF2@VA*R6QM)1WX%)ZXTMA$&ED2<;4 M:N,4Y"Y.622M40YW5#MD?4(P3JT..C5?[M3^R%.FAD_YJ1D6J,[3IV&'-<H, M1H5$92,85:S@A1A&OD6J32??*!6UD(:!1U-QJMWRJ6+4-A15=M ;YL%HI:K$ M6:LZAZ<$5>44/RB/C.*MXBJ=KE<Y(&"E#M^A:W)*"_!8T,-EI3[:#E>!=-'* MU% K-U6KWLFM*D[<JO&IQ<3GL*K:6@P*ZNHB:+6.\!&ZVDB 5R9A=]PKEF # M7$)9XU<V 6#IA(SP%@NC83R,@&^NM<6JP+(B5I:B6&)!&EVLOL86_AK',A:" M+2\0MK\PLG@(8E ,BLTQ,+:5=9P>&XJ(6<!A9N6(RC;Q<)9M8R%!;;/9AI\E M58(6:6L:TTFZ$8ZD13:8%MI0&U +N4TMW8BT G?J&XWC3A*CA\GMO[&=3-; M\ -MT;?IP;:FV+#K'MJNA?0W_V&W IS>@@I\R\#]K02G#M^@'NEV;D1Q(1() MY[AN#.1")1ENJ6VXRZ4 5@F(6VK:D Y^KM"E]2[<7RA=+&[N.3@8U^?\7>MZ M7<H$A\E'V@500$W24W*ZB^GYKK@'U01DD6.04>Q!+J\(">2:'H4<<MV.RLD2 M*\>]^EDW U]7(J"5KX$FP B8F.-H9$[/.30U%VW8G$2[C&RPHL&PB\8BP]Q! MPW,=34;R.5DRZ%#8'RQI@1 K$D(+9NY:&@SC8#!M\SBZ3O33;IJ4E&E.Q:>Y B'Q7V^++D3G,II(X?^0^B9E-JF-IC=6'%U;T]6!>]R%X0 _33 end
ycy@walt.cc.utexas.edu (Joseph Yip) (06/14/90)
Here is the second program I got. From tpf@jdyx.atlanta.ga.us Wed Jun 13 21:44:53 1990 Received: from emory.mathcs.emory.edu by walt.cc.utexas.edu (5.61/1.34/CCWF 1.5) id AA01107; Wed, 13 Jun 90 21:44:40 -0500 Received: from jdyx.UUCP by emory.mathcs.emory.edu (5.59/1.17.EUCC) via UUCP id AA00433 ; Wed, 13 Jun 90 22:44:37 EDT Return-Path: tpf@jdyx.atlanta.ga.us Received: by jdyx.UUCP (5.51/smail2.2/06-30-87) id AA12786; Wed, 13 Jun 90 21:47:09 EDT Date: Wed, 13 Jun 90 21:47:09 EDT From: tpf@jdyx.atlanta.ga.us (Tom Friedel) Message-Id: <9006140147.AA12786@jdyx.UUCP> To: ycy@ccwf.cc.utexas.edu Subject: Re: C code for FFT Newsgroups: sci.math,comp.sources.wanted,comp.lang.c References: <31567@ut-emx.UUCP> Status: R In comp.sources.wanted you write: >Hi, >I am looking for C source for Fast Fourier Transform (any dimensions and >power of 2). The netlib only has the Fortran version and I do not want >to use f2c. >If you have a good FFT program (fast!), that is excellent. >Thank you. >Joseph Yip Maybe this will help. tom Tom Friedel JDyx Enterprises (404) 320-7624 tpf@jdyx.UUCP Unix BBS: (404) 325-1719 <= 2400 ; (404) 321-5020 >= 2400 "Live simply, so that others may simply live." Well, according to the amount of E-mail I received, and as I cannot answer to all of them, I am posting the 2D FFT program I have written in C. How to use it : first, it can be used only on matrix of the type complexe, as defined at the begining of the program, with a size in 2^maxlevel (2 power maxlevel). One just have to call the program as following : InvFFT(matrix,maxlevel); As you have noticed, it is in fact an Inverse Fourier Transform (sampling matrix in the frequency space -> matrix of heights) which is performed, but just an initialisation line at the begining of the procedure InvFFT has to be change to obtain the FFT : Omega = 2 * M_PI / dimension; must be change into Omega = - 2 * M_PI / dimension; I am clearly aware that this program has not been optimised as it could be, mainly in order to keep the the clarity of the algorithm, but I am open to any suggestion that you could make in order to improve the process. I have tested it in a fractal landscape project and it seems to work correctly (when there were bugs in it such as forgetting the multiplication of omega by two at the recursive calls, I had funny landscapes looking like a micro-chip seen through an electronic microscop). Well, I think that is all, for the moment. Here comes the program. Note : you will have probably to cut the signature at the end of the file as well. Jean-Marc ----cut here -----cut here -----cut here -----cut here -----cut here ----- /* It is a pity I have to add this notice, but it must be done. */ /* NOTICE * * Copyright (c) 1990, Jean-Marc de Lauwereyns * * for the procedures mirror, RecFFT, InvFFT. * * These procedures can be used and given freely to anybody but can not * * be sold or used in any commercial program or in any book without the * * permission of the original author. * * These procedures can be modified at the condition that the nature of * * modification is mentioned in comments at the place of the modification * * with the original lines accompanied by the name of the original author.* * This notice itself can not be removed or modified except by the * * original author himself. */ /* you must have included the file /usr/include/math.h at the begining of * * your program */ typedef struct { double a,b; } complexe; /****************************************************/ /* function mirror just transform the binary */ /* representation of an integer into its mirror */ /* binary representation on maxlevel bits */ /* eg : with maxlevel = 6, 000101 becomes 101000 */ /****************************************************/ int mirror(index,maxlevel) int index; int maxlevel; { register int i; int res; res = 0; for (i=0;i<maxlevel;i++) { res = (res << 1) | (index & (int)1); index = (index >> 1); } return(res); } /*****************************************************/ /* recursive procedure to calculate the FFT on 2^N */ /*****************************************************/ void RecFFT(A,k,l,p,q,Omega) complexe **A; int k,l; /* indexes for the begin and the end on the lines taken */ int p,q; /* indexes for the begin and the end of the columns taken */ double Omega; { int i,j,m; complexe a,b,c,d; double argument[3]; if (k != l) { m = (l-k)/2; for (i = k;i < k + m;i++) { argument[0] = Omega * (double)(i-k); for (j = p;j < p + m;j++) { argument[1] = Omega * (double)(j-p); argument[2] = Omega * (double)(i-k + j-p); a.a = A[i][j].a;a.b = A[i][j].b; b.a = A[i+m][j].a;b.b = A[i+m][j].b; c.a = A[i][j+m].a;c.b = A[i][j+m].b; d.a = A[i+m][j+m].a;d.b = A[i+m][j+m].b; A[i][j].a = a.a + b.a + c.a + d.a; A[i][j].b = a.b + b.b + c.b + d.b; A[i+m][j].a = (a.a - b.a + c.a - d.a)*cos(argument[0]) - (a.b - b.b + c.b - d.b)*sin(argument[0]); A[i+m][j].b = (a.a - b.a + c.a - d.a)*sin(argument[0]) + (a.b - b.b + c.b - d.b)*cos(argument[0]); A[i][j+m].a = (a.a + b.a - c.a - d.a)*cos(argument[1]) - (a.b + b.b - c.b - d.b)*sin(argument[1]); A[i][j+m].b = (a.a + b.a - c.a - d.a)*sin(argument[1]) + (a.b + b.b - c.b - d.b)*cos(argument[1]); A[i+m][j+m].a = (a.a - b.a - c.a + d.a)*cos(argument[2]) - (a.b - b.b - c.b + d.b)*sin(argument[2]); B[i+m][j+m].b = (a.a - b.a - c.a + d.a)*sin(argument[2]) + (a.b - b.b - c.b + d.b)*cos(argument[2]); } } if (m != 0) { RecFFT(k,k+m,p,p+m,Omega*2); RecFFT(k+m,l,p,p+m,Omega*2); RecFFT(k,k+m,p+m,q,Omega*2); RecFFT(k+m,l,p+m,q,Omega*2); } } } /********************************************************/ /* InvFFT is a procedure which produce the inverse of */ /* Fourier Transformation. */ /********************************************************/ void InvFFT(A,interlevel) complexe **A; int interlevel; { int dimension; int i,j,i1,j1; double coeff,Omega; dimension = (1 << interlevel); Omega = 2 * M_PI / dimension; RecFFT(A,0,dimension,0,dimension,Omega); for (i=0;i<=(dimension - 1);i++) /* swap every rows i with its mirror image */ { i1 = mirror(i,interlevel); if (i1 > i) /* if i1 <= i, the swaaping have been already made */ for (j=0;j<dimension;j++) { coeff = A[i][j].a; A[i][j].a = A[i1][j].a; A[i1][j].a = coeff; /* this may be remove if you are sure that the result is a * * matrix of real which means that the original matrix must * * verify the following conditions : * * ______ ______ * * A(N-i,N-j) = A(i,j) , A(0,N-j) = A(0,j) * * ______ * * A(N-i,0) = A(i,0) for i,j > 0 and A(0,0) is a pure * * real */ coeff = A[i][j].b; A[i][j].a = A[i1][j].a; A[i1][j].a = coeff; } } for (j=0;j<=(dimension - 1);j++) /* idem for the columns */ { j1 = mirror(j,interlevel); if (j1 >j) for (i=0;i<dimension;i++) { coeff = A[i][j].a; A[i][j].a = A[i][j1].a; A[i][j1].a = coeff; /* part which may be removed if (see above) */ coeff = A[i][j].b; A[i][j].b = A[i][j1].b; A[i][j1].b = coeff; } } } ----cut here -----cut here -----cut here -----cut here -----cut here ---- Jean-Marc de Lauwereyns | ____ | e-mail addresses : Heriot-Watt University | |\ /| | \ | JANET: lauwer@uk.ac.hw.cs Computer Science Department | | \/ | |___/ | ARPA.: lauwer@cs.hw.ac.uk Edinburgh | \___/ |___| \ |