(* ;;-*-LISP-*- KEEP EMACS HAPPY ******************************** * * FFT -- FAST FOURIER TRANSFORM * *(1) Stolen from <LISPCORE>BENCHMARKS>BARROW>FFTI.LSP ****************************************************************) (DEFPROP FASTELT DMACRO (OPENLAMBDA (A N) (\GETBASEPTR (ARRAYP.BASE A) (LLSH (1- N) 1)))) (DEFPROP FASTSETA DMACRO (OPENLAMBDA (A N V) (\PUTBASEPTR (ARRAYP.BASE A) (LLSH (1- N) 1) V))) (DEFEXPR (FFT AREAL AIMAG) (* Fast Fourier Transform AREAL = real part%, AIMAG = imaginary part) (PROG (AR AI PI I J K M N LE LE1 IP NV2 NM1 UR UI WR WI TR TI) (* Initialize) (SETQ AR AREAL) (SETQ AI AIMAG) (SETQ PI 3.141593) (SETQ N (ARRAYSIZE AR)) (SETQ NV2 (/ N 2)) (SETQ NM1 (1- N)) (* Compute M = log (N)) (SETQ M 0) (SETQ I 1) L1 (COND ((< I N) (SETQ M (1+ M)) (SETQ I (+ I I)) (GO L1))) (COND ((NOT (= N (EXPT 2.0 M))) (PRIN1 "Error ... array size not a power of two.") (HELP) (RETURN (TERPRI)))) (* Interchange elements) (SETQ J 1) (* in bit-reversed order) (SETQ I 1) L3 (COND ((< I J) (SETQ TR (FASTELT AR J)) (SETQ TI (FASTELT AI J)) (FASTSETA AR J (FASTELT AR I)) (FASTSETA AI J (FASTELT AI I)) (FASTSETA AR I TR) (FASTSETA AI I TI))) (SETQ K NV2) L6 (COND ((< K J) (SETQ J (- J K)) (SETQ K (/ K 2)) (GO L6))) (SETQ J (+ J K)) (SETQ I (1+ I)) (COND ((< I N) (GO L3))) (FOR L FROM 1 TO M DO (* Loop thru stages) (SETQ LE (EXPT 2.0 L)) (SETQ LE1 (/ LE 2)) (SETQ UR 1.0) (SETQ UI 0.0) (SETQ WR (COS (/$ PI (FLOAT LE1)))) (SETQ WI (SIN (/$ PI (FLOAT LE1)))) (FOR J FROM 1 TO LE1 DO (* Loop thru butterflies) (FOR I FROM J BY LE TO N DO (* Do a butterfly) (SETQ IP (+ I LE1)) (SETQ TR (-$ (x$ (FASTELT AR IP) UR) (x$ (FASTELT AI IP) UI))) (SETQ TI (+$ (x$ (FASTELT AR IP) UI) (x$ (FASTELT AI IP) UR))) (FASTSETA AR IP (-$ (FASTELT AR I) TR)) (FASTSETA AI IP (-$ (FASTELT AI I) TI)) (FASTSETA AR I (+$ (FASTELT AR I) TR)) (FASTSETA AI I (+$ (FASTELT AI I) TI))) (SETQ TR (-$ (x$ UR WR) (x$ UI WI))) (SETQ TI (+$ (x$ UR WI) (x$ UI WR))) (SETQ UR TR) (SETQ UI TI))) (RETURN T) )) (DEFEXPR (FFT.TEST (OPTIONAL SIZE 1024)) (PROG (RE IM) (SETQ RE (ARRAY SIZE)) (SETQ IM (ARRAY SIZE)) (FOR I FROM 1 TO SIZE DO (SETF (ELT RE I) 0.0) (SETF (ELT IM I) 0.0)) (RETURN (TIMEALL (FFT RE IM))) )) STOP