home *** CD-ROM | disk | FTP | other *** search
/ HAKERIS 11 / HAKERIS 11.ISO / linux / system / LinuxConsole 0.4 / linuxconsole0.4install-en.iso / guile0.4.lcm / share / guile / slib / fft.scm < prev    next >
Encoding:
Text File  |  2004-01-06  |  2.3 KB  |  71 lines
  1. ;;;"fft.scm" Fast Fourier Transform
  2. ;Copyright (C) 1999 Aubrey Jaffer
  3. ;
  4. ;Permission to copy this software, to redistribute it, and to use it
  5. ;for any purpose is granted, subject to the following restrictions and
  6. ;understandings.
  7. ;
  8. ;1.  Any copy made of this software must include this copyright notice
  9. ;in full.
  10. ;
  11. ;2.  I have made no warrantee or representation that the operation of
  12. ;this software will be error-free, and I am under no obligation to
  13. ;provide any services, by way of maintenance, update, or otherwise.
  14. ;
  15. ;3.  In conjunction with products arising from the use of this
  16. ;material, there shall be no use of my name in any advertising,
  17. ;promotional, or sales literature without prior written consent in
  18. ;each case.
  19.  
  20. ;;;; See:
  21. ;;; Introduction to Algorithms (MIT Electrical
  22. ;;;    Engineering and Computer Science Series)
  23. ;;; by Thomas H. Cormen, Charles E. Leiserson (Contributor),
  24. ;;;    Ronald L. Rivest (Contributor)
  25. ;;; MIT Press; ISBN: 0-262-03141-8 (July 1990)
  26.  
  27. ;;; http://www.astro.virginia.edu/~eww6n/math/DiscreteFourierTransform.html
  28. ;;; differs in the direction of rotation of the complex unit vectors.
  29.  
  30. (require 'array)
  31.  
  32. (define (fft:shuffled&scaled ara n scale)
  33.   (define lgn (integer-length (+ -1 n)))
  34.   (define new (apply make-array 0 (array-dimensions ara)))
  35.   (define bit-reverse (lambda (width in)
  36.             (if (zero? width) 0
  37.                 (+ (bit-reverse (+ -1 width) (quotient in 2))
  38.                    (ash (modulo in 2) (+ -1 width))))))
  39.   (if (not (eqv? n (expt 2 lgn)))
  40.       (slib:error 'fft "array length not power of 2" n))
  41.   (do ((k 0 (+ 1 k)))
  42.       ((>= k n) new)
  43.     (array-set! new (* (array-ref ara k) scale) (bit-reverse lgn k))))
  44.  
  45. (define (dft! ara n dir)
  46.   (define lgn (integer-length (+ -1 n)))
  47.   (define pi2i (* 0+8i (atan 1)))
  48.   (do ((s 1 (+ 1 s)))
  49.       ((> s lgn) ara)
  50.     (let* ((m (expt 2 s))
  51.        (w_m (exp (* dir (/ pi2i m))))
  52.        (m/2-1 (+ (quotient m 2) -1)))
  53.       (do ((j 0 (+ 1 j))
  54.        (w 1 (* w w_m)))
  55.       ((> j m/2-1))
  56.     (do ((k j (+ m k)))
  57.         ((>= k n))
  58.       (let* ((k+m/2 (+ k m/2-1 1))
  59.          (t (* w (array-ref ara k+m/2)))
  60.          (u (array-ref ara k)))
  61.         (array-set! ara (+ u t) k)
  62.         (array-set! ara (- u t) k+m/2)))))))
  63.  
  64. (define (fft ara)
  65.   (define n (car (array-dimensions ara)))
  66.   (dft! (fft:shuffled&scaled ara n 1) n 1))
  67.  
  68. (define (fft-1 ara)
  69.   (define n (car (array-dimensions ara)))
  70.   (dft! (fft:shuffled&scaled ara n (/ n)) n -1))
  71.