(* ;;-*-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)))

(MOVD 'ELT 'FASTELT)
(MOVD 'SETA 'FASTSETA)

(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