1 files changed, 59 insertions, 1 deletions
diff --git a/tools/tfft.scm b/tools/tfft.scm
index 9c7bbaa..8da1272 100644
--- a/tools/tfft.scm
+++ b/tools/tfft.scm
@@ -84,12 +84,70 @@
       (set! ui ti))
     #t))
  
+(if (not (defined? 'complex))
+  (define complex make-rectangular))
+(if (not (defined? 'when))
+  (define-macro (when test . body) `(if ,test (begin ,@body))))
+
+(define* (cfft data n (dir 1)) ; complex data
+  (unless n (set! n (length data)))
+  (do ((i 0 (+ i 1))
+       (j 0))
+      ((= i n))
+    (if (> j i)
+	(let ((temp (data j)))
+	  (set! (data j) (data i))
+	  (set! (data i) temp)))
+    (do ((m (/ n 2) (/ m 2)))
+        ((or (< m 2) 
+             (< j m))
+         (set! j (+ j m)))
+     (set! j (- j m))))
+  (do ((ipow (floor (log n 2)))
+       (prev 1)
+       (lg 0 (+ lg 1))
+       (mmax 2 (* mmax 2))
+       (pow (/ n 2) (/ pow 2))
+       (theta (complex 0.0 (* pi dir)) (* theta 0.5)))
+      ((= lg ipow))
+    (do ((wpc (exp theta))
+         (wc 1.0)
+         (ii 0 (+ ii 1)))
+	((= ii prev)
+	 (set! prev mmax))
+      (do ((jj 0 (+ jj 1))
+           (i ii (+ i mmax))
+           (j (+ ii prev) (+ j mmax)))
+          ((>= jj pow))
+        (let ((tc (* wc (data j))))
+          (set! (data j) (- (data i) tc))
+          (set! (data i) (+ (data i) tc))))
+      (set! wc (* wc wpc))))
+  data)
+
+(when (defined? 'equivalent?)
+  (unless (equivalent? (cfft (vector 0.0 1+i 0.0 0.0)) #(1+1i -1+1i -1-1i 1-1i))
+    (format *stderr* "cfft 1: ~S~%" (cfft (vector 0.0 1+i 0.0 0.0))))
+  (let-temporarily (((*s7* 'equivalent-float-epsilon) 1e-14))
+    (unless (equivalent? (cfft (vector 0 0 1+i 0 0 0 1-i 0)) #(2 -2 -2 2 2 -2 -2 2))
+      (format *stderr* "cfft 2: ~S~%" (cfft (vector 0 0 1+i 0 0 0 1-i 0))))))
+
 (define (fft-bench)
   (let ((*re* (make-float-vector (+ size 1) 0.10))
 	(*im* (make-float-vector (+ size 1) 0.10)))
     (do ((ntimes 0 (+ ntimes 1)))
 	((= ntimes times))
-      (b_fft *re* *im*))))
+      (b_fft *re* *im*)))
+
+  (let* ((n 256)
+	 (cdata (make-vector n 0.0)))
+    (do ((i 0 (+ i 1)))
+	((= i times))
+      (fill! cdata 0.0)
+      (vector-set! cdata 2 1+i)
+      (vector-set! cdata (- n 1) 1-i)
+      (cfft cdata)))
+  )
  
 (fft-bench)