Update of /project/sapaclisp/cvsroot/sapaclisp In directory common-lisp.net:/tmp/cvs-serv32056/sapaclisp

Modified Files: dft-and-fft.lisp Log Message: Substituted the recurrence relation used for computing the sines and cosines needed by the FFT. The new one is more stable. (Thanks to B. Lucier)

Date: Mon May 16 11:59:12 2005 Author: mmommer

Index: sapaclisp/dft-and-fft.lisp diff -u sapaclisp/dft-and-fft.lisp:1.2 sapaclisp/dft-and-fft.lisp:1.3 --- sapaclisp/dft-and-fft.lisp:1.2 Mon May 16 11:07:35 2005 +++ sapaclisp/dft-and-fft.lisp Mon May 16 11:59:12 2005 @@ -337,21 +337,30 @@ ;------------------------------------------------------------------------------- ;;; Everything below here consists of internal symbols in the SAPA package ;;; and should be regarded as "dirty laundry" ... -;------------------------------------------------------------------------------- -;------------------------------------------------------------------------------- -(defun pre-fft - (n - &key - (complex-exp-vector nil)) - (if (power-of-2 n) +;----------------------------------------------------------------------------- +;----------------------------------------------------------------------------- +(defun pre-fft (n &key (complex-exp-vector nil)) + "Computes a vector with the cosines and sines needed for the FFT" + (when (power-of-2 n) (let* ((the-vector (if complex-exp-vector complex-exp-vector - (make-array n))) + (make-array n + :element-type '(complex double-float)))) (s (/ (* 2 pi) n)) - (c1 (complex (cos s) (- (sin s)))) - (c2 (complex 1.0d0 0.0d0))) + (2sdh (* 2 (expt (sin (/ s 2)) 2))) + (sd (sin s)) + (sk 0.0d0) + (ck 1.0d0)) + (declare (double-float ck sk sd 2sdh)) + (dotimes (i n the-vector) - (setf (aref the-vector i) c2 - c2 (* c2 c1)))))) + (setf (aref the-vector i) + (complex ck sk) + + (values ck sk) + (values (- ck (* sk sd) (* ck 2sdh)) + (+ sk (* ck sd) (* (- sk) 2sdh)))))))) + +

;------------------------------------------------------------------------------- ;;; saves results up to size (expt 2 31) = 2147483648