[sapaclisp-cvs] CVS update: sapaclisp/dft-and-fft.lisp

[Mario Mommer]

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