jurjen@cwi.nl (Jurjen NE Bos) (02/18/91)
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 <name> <value> <name> <value> ... 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