home *** CD-ROM | disk | FTP | other *** search
/ Usenet 1994 October / usenetsourcesnewsgroupsinfomagicoctober1994disk2.iso / misc / volume1 / 8707 / 49 / fft.cl < prev    next >
Lisp/Scheme  |  1990-07-13  |  3KB  |  97 lines
  1. ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  2. ; File:         fft.cl
  3. ; Description:  FFT benchmark from the Gabriel tests.
  4. ; Author:       Harry Barrow
  5. ; Created:      8-Apr-85
  6. ; Modified:     6-May-85 09:29:22 (Bob Shaw)
  7. ; Language:     Common Lisp
  8. ; Package:      User
  9. ; Status:       Public Domain
  10. ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  11.  
  12. ;;; FFT -- This is an FFT benchmark written by Harry Barrow.
  13. ;;; It tests a variety of floating point operations,
  14. ;;; including array references.
  15.  
  16. (defvar *re* (make-array 1025 :element-type 'single-float
  17.                   :initial-element 0.0))
  18.  
  19. (defvar *im* (make-array 1025 :element-type 'single-float
  20.                   :initial-element 0.0))
  21.  
  22. (defun fft (areal aimag)
  23.   (declare (vector areal aimag))
  24.   (prog (ar ai i j k m n le le1 ip nv2 nm1 ur ui wr wi tr ti)
  25.     (declare (vector ai ai) (fixnum i j k m n le le1 nv2 nm1))
  26.     ;; initialize
  27.     (setq ar areal
  28.           ai aimag
  29.       n (array-dimension ar 0)
  30.       n (1- n)
  31.       nv2 (floor n 2)
  32.       nm1 (1- n)
  33.       m 0                    ;compute m = log(n)
  34.       i 1)
  35.  l1 (cond ((< i n) (setq m (1+ m) i (+ i i)) (go l1)))
  36.     (cond ((not (equal n (expt 2 m)))
  37.        (princ "error ... array size not a power of two.")
  38.        (read)
  39.        (return (terpri))))
  40.     ;; interchange elements in bit-reversed order
  41.     (setq j 1 i 1)
  42.  l3 (cond ((< i j)
  43.        (setq tr (aref ar j) ti (aref ai j))
  44.        (setf (aref ar j) (aref ar i))
  45.        (setf (aref ai j) (aref ai i))
  46.        (setf (aref ar i) tr)
  47.        (setf (aref ai i) ti)))
  48.     (setq k nv2)
  49.  l6 (cond ((< k j) 
  50.        (setq j (- j k) k (/ k 2))
  51.        (go l6)))
  52.     (setq j (+ j k) i (1+ i))
  53.     (cond ((< i n)
  54.        (go l3)))
  55.     (do ((l 1 (1+ l)))            ;loop thru stages (syntax converted
  56.     ((> l m))                       ; from old MACLISP style \bs)
  57.     (declare (fixnum l))
  58.     (setq le (expt 2 l)
  59.           le1 (floor le 2)
  60.           ur 1.0
  61.           ui 0.
  62.           wr (cos (/ pi (float le1)))
  63.           wi (sin (/ pi (float le1))))
  64.     ;; loop thru butterflies
  65.     (do ((j 1 (1+ j)))
  66.         ((> j le1))
  67.         (declare (fixnum j))
  68.         ;; do a butterfly
  69.         (do ((i j (+ i le)))
  70.         ((> i n))
  71.         (declare (fixnum i))
  72.         (setq ip (+ i le1)
  73.               tr (- (* (aref ar ip) ur)
  74.                 (* (aref ai ip) ui))
  75.               ti (+ (* (aref ar ip) ui)
  76.                 (* (aref ai ip) ur)))
  77.         (setf (aref ar ip) (- (aref ar i) tr))
  78.         (setf (aref ai ip) (- (aref ai i) ti))
  79.         (setf (aref ar i) (+ (aref ar i) tr))
  80.         (setf (aref ai i) (+ (aref ai i) ti))))
  81.     (setq tr (- (* ur wr) (* ui wi))
  82.           ti (+ (* ur wi) (* ui wr))
  83.           ur tr
  84.           ui ti))
  85.     (return t)))
  86.  
  87. ;;; the timer which does 10 calls on fft
  88.  
  89. (defmacro fft-bench ()
  90.   `(do ((ntimes 0 (1+ ntimes)))
  91.       ((= ntimes 10))
  92.     (fft *re* *im*)))
  93.  
  94. ;;; call:  (fft-bench)
  95.  
  96. (run-benchmark "FFT" '(fft-bench))
  97.