From en.ecn.purdue.edu!noose.ecn.purdue.edu!samsung!sdd.hp.com!spool.mu.edu!uunet!mcsun!hp4nl!charon!jurjen 18 Feb 91 10:28:30 GMT Path: en.ecn.purdue.edu!noose.ecn.purdue.edu!samsung!sdd.hp.com!spool.mu.edu!uunet!mcsun!hp4nl!charon!jurjen From: jurjen@cwi.nl (Jurjen NE Bos) Newsgroups: comp.sys.handhelds Subject: HP28/48: Very FFT Keywords: FFT Message-ID: <2956@charon.cwi.nl> Date: 18 Feb 91 10:28:30 GMT Sender: news@cwi.nl Organization: STORC, Veldhoven Lines: 83 Originator: jurjen@lijster.cwi.nl A few friends asked me to write an FFT (fast fourier transform; converts a vector to its "frequency diagram") program. Looking through my archives, I could not find any, although I vaguely remember another program was posted. This program is written with speed in mind. Most computations are done with vectors instead of numbers. At the end of the posting you see a HP48 directory. HP28 users must know the following: The format is DIR ... END. It is sent in translate code 3; that means that \.. codes are to be translated: \<< open program; \>> close program; \|v down arrow (make it a 'd') \Gr greek rho (make it an 'R') \GS greek Sigma (don't type in the routine DFT) @ signifies comment. Don't type in from here to end of line A very nice feature of this program is that it works for ANY vector length. If the vector length is odd, it is only slightly faster than the regular DFT program in the end of the directory. If the length is even, it is slightly faster, and the more factor of two, the faster. If the length is a power of two, the regular FFT algorithm is applied. All this is in one algorithm, without discrimination of all cases! %%HP: T(3)F(,); DIR CST { { "FFT" \<< 1 FFT \>> } { "-FFT" \<< -1 FFT \>> } } FFT @ FFT routine. @ Input: 2: Vector 1: 1 for FFT, -1 for inverse @ Ouput: 1: transformed vector \<< OVER SIZE 1 GET IF OVER 0 < @ Divide by length if inverse FFT THEN ROT OVER / ROT ROT END (0;6,28318530718) ROT * OVER / SWAP DUP WHILE DUP 2 MOD NOT @ Divide out factors of 2 REPEAT 2 / END SWAP OVER / \-> \Gr r p @ p: twopower, r: odd rest @ rho: root of unity \<< { r p } RDM IF r 1 \=/ @ Compute regular r-point FFT p times simultaneously THEN A\->L \Gr p * EXP 1 \-> m \Grp \Grn \<< 1 r START m 1 GET 1 2 r FOR k \Grn * m k GET OVER * ROT + SWAP NEXT DROP OBJ\-> DROP \Grp '\Grn' STO* NEXT \>> { r p } \->ARRY END TRN CONJ @ TRN does transpose & conjugate, so conjugate back WHILE p 1 \=/ @ combine r-point FFTs to 2r-point FFTs REPEAT OBJ\-> DROP 'p' 2 STO/ { p r } \->ARRY TRN A\->L @ from here on, we work with conjugated numbers! \Gr NEG p * EXP 1 \-> m \Grp \Grn \<< { p r } \->ARRY TRN 1 r FOR k m k GET \Grn * OBJ\-> DROP \Grp '\Grn' STO* NEXT \>> { r p } \->ARRY + LASTARG - \-> m \<< OBJ\-> DROP m \>> OBJ\-> DROP 'r' 2 STO* { r p } \->ARRY TRN @ numbers are normal again END { r } RDM \>> \>> A\->L @ convert matrix to list of rows. This doesn't use sigma+, because @ this trick doesn't work with complex numbers :-( \<< OBJ\-> OBJ\-> DROP { } \-> r c u \<< 1 r START c \->ARRY 'u' STO+ NEXT u \>> \>> DFT @ The regular Fourier Transform. Extra short version. \<< SWAP DUP SIZE 1 GET (0;6,28318530718) OVER / 4 PICK * \-> d x s q \<< 0 s 1 - FOR k '\GS(n=0;s-1;x(n+1)*EXP(q*k*n))' \->NUM NEXT s \->ARRY IF d -1 == THEN s / END \>> \>> END