diff options
Diffstat (limited to 'fft.scm')
-rw-r--r-- | fft.scm | 94 |
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 |