summaryrefslogtreecommitdiffstats
path: root/fft.scm
diff options
context:
space:
mode:
Diffstat (limited to 'fft.scm')
-rw-r--r--fft.scm94
1 files changed, 0 insertions, 94 deletions
diff --git a/fft.scm b/fft.scm
deleted file mode 100644
index feefec8..0000000
--- a/fft.scm
+++ /dev/null
@@ -1,94 +0,0 @@
-;;;"fft.scm" Fast Fourier Transform
-;Copyright (C) 1999, 2003 Aubrey Jaffer
-;
-;Permission to copy this software, to modify it, to redistribute it,
-;to distribute modified versions, and to use it for any purpose is
-;granted, subject to the following restrictions and understandings.
-;
-;1. Any copy made of this software must include this copyright notice
-;in full.
-;
-;2. I have made no warranty or representation that the operation of
-;this software will be error-free, and I am under no obligation to
-;provide any services, by way of maintenance, update, or otherwise.
-;
-;3. In conjunction with products arising from the use of this
-;material, there shall be no use of my name in any advertising,
-;promotional, or sales literature without prior written consent in
-;each case.
-
-;;;; See:
-;;; Introduction to Algorithms (MIT Electrical
-;;; Engineering and Computer Science Series)
-;;; by Thomas H. Cormen, Charles E. Leiserson (Contributor),
-;;; Ronald L. Rivest (Contributor)
-;;; MIT Press; ISBN: 0-262-03141-8 (July 1990)
-
-;;; http://www.astro.virginia.edu/~eww6n/math/DiscreteFourierTransform.html
-;;; differs in the direction of rotation of the complex unit vectors.
-
-(require 'array)
-(require 'logical)
-
-;;@code{(require 'fft)}
-;;@ftindex fft
-
-(define (fft:shuffle&scale new ara n scale)
- (define lgn (integer-length (+ -1 n)))
- (if (not (eqv? n (expt 2 lgn)))
- (slib:error 'fft "array length not power of 2" n))
- (do ((k 0 (+ 1 k)))
- ((>= k n) new)
- (array-set! new (* (array-ref ara k) scale) (reverse-bit-field k 0 lgn))))
-
-(define (dft! ara n dir)
- (define lgn (integer-length (+ -1 n)))
- (define pi2i (* 0+8i (atan 1)))
- (do ((s 1 (+ 1 s)))
- ((> s lgn) ara)
- (let* ((m (expt 2 s))
- (w_m (exp (* dir (/ pi2i m))))
- (m/2-1 (+ (quotient m 2) -1)))
- (do ((j 0 (+ 1 j))
- (w 1 (* w w_m)))
- ((> j m/2-1))
- (do ((k j (+ m k)))
- ((>= k n))
- (let* ((k+m/2 (+ k m/2-1 1))
- (t (* w (array-ref ara k+m/2)))
- (u (array-ref ara k)))
- (array-set! ara (+ u t) k)
- (array-set! ara (- u t) k+m/2)))))))
-
-;;@args array
-;;@var{array} is an array of @code{(expt 2 n)} numbers. @code{fft}
-;;returns an array of complex numbers comprising the
-;;@dfn{Discrete Fourier Transform} of @var{array}.
-(define (fft ara)
- (define n (car (array-dimensions ara)))
- (define new (apply make-array ara (array-dimensions ara)))
- (dft! (fft:shuffle&scale new ara n 1) n 1))
-
-;;@args array
-;;@code{fft-1} returns an array of complex numbers comprising the
-;;inverse Discrete Fourier Transform of @var{array}.
-(define (fft-1 ara)
- (define n (car (array-dimensions ara)))
- (define new (apply make-array ara (array-dimensions ara)))
- (dft! (fft:shuffle&scale new ara n (/ n)) n -1))
-
-;;@noindent
-;;@code{(fft-1 (fft @var{array}))} will return an array of values close to
-;;@var{array}.
-;;
-;;@example
-;;(fft '#(1 0+i -1 0-i 1 0+i -1 0-i)) @result{}
-;;
-;;#(0.0 0.0 0.0+628.0783185208527e-18i 0.0
-;; 0.0 0.0 8.0-628.0783185208527e-18i 0.0)
-;;
-;;(fft-1 '#(0 0 0 0 0 0 8 0)) @result{}
-;;
-;;#(1.0 -61.23031769111886e-18+1.0i -1.0 61.23031769111886e-18-1.0i
-;; 1.0 -61.23031769111886e-18+1.0i -1.0 61.23031769111886e-18-1.0i)
-;;@end example