home *** CD-ROM | disk | FTP | other *** search
/ Eagles Nest BBS 8 / Eagles_Nest_Mac_Collection_Disc_8.TOAST / Developer Tools⁄Additions / MacScheme20 / Mathlib / fft.scm < prev    next >
Encoding:
Text File  |  1989-04-27  |  2.2 KB  |  62 lines  |  [TEXT/????]
  1. ;;; $Header: fft.scm,v 1.4 88/01/26 07:43:39 GMT gjs Exp $
  2.  
  3. (if-mit
  4.  (declare (usual-integrations = + - * /
  5.                  zero? 1+ -1+
  6.                  ;; truncate round floor ceiling
  7.                  sqrt exp log sin cos)))
  8.  
  9. ;;;;    Discrete Fourier Transform -- GJS -- 8 April 1987
  10. ;;; Implemented using FFT for records whose length is a power of 2.
  11. ;;; Data records are represented as lists of (complex) numbers.
  12.  
  13. (define (power-of-2? n)
  14.   (and (integer? n) 
  15.        (positive? n)
  16.        (let lp ((n n)) (or (= n 1) (and (even? n) (lp (quotient n 2)))))))
  17.  
  18. (define (make-transform-pair period)
  19.   (define (roots-of-unity n)
  20.     (let ((n/2 (/ n 2)) (2pi/n (/ 2pi n)))
  21.       (let loop ((k 0))
  22.         (if (= k n/2)
  23.             '()
  24.             (cons (make-rectangular (cos (* 2pi/n k)) (sin (* 2pi/n k)))
  25.                   (loop (1+ k)))))))
  26.   (define (evens list)
  27.     (if (null? list)
  28.         '()
  29.         (cons (car list) (evens (cddr list)))))
  30.   (define (odds list)
  31.     (if (null? list)
  32.         '()
  33.         (cons (cadr list) (odds (cddr list)))))
  34.   ;; Below here all arithmetic is assumed to be complex generic.
  35.   (define (fft-internal data ws)
  36.     (if (null? (cdr data))
  37.         data
  38.         (let ((ews (if (null? (cdr ws)) ws (evens ws))))
  39.           (let ((even-fft (fft-internal (evens data) ews))
  40.                 (odd-fft (map * ws (fft-internal (odds data) ews))))
  41.             (append (map + even-fft odd-fft)
  42.                     (map - even-fft odd-fft))))))
  43.   (if (power-of-2? period)
  44.       (let* ((w-list (roots-of-unity period))
  45.          (w*-list (map conjugate w-list)))
  46.         (define (ft data)
  47.           (if (= period (length data))
  48.               (map (lambda (z) (/ z period))
  49.                    (fft-internal data w*-list))
  50.               (error "Wrong length data record -- FT" period)))
  51.         (define (ift data)
  52.           (if (= period (length data))
  53.               (fft-internal data w-list)
  54.               (error "Wrong length data record -- IFT" period)))
  55.     (list period ft ift))
  56.       (error "Period is not a power of 2 -- MAKE-TRANSFORM-PAIR")))
  57.  
  58. ;;; For example, we may:
  59. ;;; (define fts (make-transform-pair 64))
  60. ;;; (define ft (cadr fts)) ; This gets the transform.
  61. ;;; (define ift (caddr fts)) ; This gets the inverse transform.
  62.