source: project/release/3/crunch/trunk/fft.scm @ 9918

Last change on this file since 9918 was 9918, checked in by Kon Lovett, 12 years ago

Using canonical directory structure.

File size: 3.4 KB
Line 
1;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2;;; File:         fft.sc
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;;;               11-Aug-87 (Will Clinger)
8;;;               16-Nov-94 (Qobi)
9;;;               31-Mar-98 (Qobi)
10;;;               26-Mar-00 (flw)
11;;; Language:     Scheme
12;;; Status:       Public Domain
13;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
14
15(define pi (atan2 0 -1))
16
17(define count 0)                        ;(felix)
18
19;;; FFT -- This is an FFT benchmark written by Harry Barrow.
20;;; It tests a variety of floating point operations,
21;;; including array references.
22
23(define *re* (make-f64vector 1025 0.0))
24
25(define *im* (make-f64vector 1025 0.0))
26
27(define (fft areal aimag)
28 (let ((ar areal)                       ;Qobi
29       (ai aimag)                       ;Qobi
30       (i 0)
31       (j 0)
32       (k 0)
33       (m 0)
34       (n 0)
35       (le 0)
36       (le1 0)
37       (ip 0)
38       (nv2 0)
39       (nm1 0)
40       (ur 0.0)                         ;Qobi
41       (ui 0.0)                         ;Qobi
42       (wr 0.0)                         ;Qobi
43       (wi 0.0)                         ;Qobi
44       (tr 0.0)                         ;Qobi
45       (ti 0.0))                        ;Qobi
46  ;; initialize
47  (set! ar areal)
48  (set! ai aimag)
49  (set! n 1025)
50  (set! n (- n 1))
51  (set! nv2 (quotient n 2))             ;(felix) was "quotient"
52  (set! nm1 (- n 1))
53  (set! m 0)                            ;compute m = log(n)
54  (set! i 1)
55  (let loop ()
56   (if (< i n)
57       (begin (set! m (+ m 1))
58              (set! i (+ i i))
59              (loop))))
60  (cond ((not (= n (let ((px 0))
61                     (let loop ((i m) (p 1)) ;Qobi
62                       (if (zero? i) 
63                           (set! px p)  ;(felix) omit loop result
64                           (loop (- i 1) (* 2 p))))
65                     px) ) )
66         (display "array size not a power of two.")
67         (newline)))
68  ;; interchange elements in bit-reversed order
69  (set! j 1)
70  (set! i 1)
71  (let l3 ()
72   (cond ((< i j)
73          (set! tr (f64vector-ref ar j))
74          (set! ti (f64vector-ref ai j))
75          (f64vector-set! ar j (f64vector-ref ar i))
76          (f64vector-set! ai j (f64vector-ref ai i))
77          (f64vector-set! ar i tr)
78          (f64vector-set! ai i ti)))
79   (set! k nv2)
80   (let l6 ()
81    (cond ((< k j)
82           (set! j (- j k))
83           (set! k (quotient k 2))      ;Qobi: was / but this violates R4RS
84           (l6))))
85   (set! j (+ j k))
86   (set! i (+ i 1))
87   (cond ((< i n) (l3))))
88  ;; loop thru stages (syntax converted from old MACLISP style \bs)
89  (do ((l 1 (+ l 1))) ((> l m))
90   (let loop ((i l) (p 1))      ;Qobi
91     (if (zero? i) 
92         (set! le p)                    ;(felix) omit loop result
93         (loop (- i 1) (* 2 p))))
94   (set! le1 (quotient le 2))
95   (set! ur 1.0)
96   (set! ui 0.0)
97   (set! wr (cos (/ pi le1)))
98   (set! wi (sin (/ pi le1)))
99   ;; loop thru butterflies
100   (do ((j 1 (+ j 1))) ((> j le1))
101    ;; do a butterfly
102    (do ((i j (+ i le))) ((> i n))
103      (set! count (+ count 1))          ;(felix)
104     (set! ip (+ i le1))
105     (set! tr (- (* (f64vector-ref ar ip) ur) (* (f64vector-ref ai ip) ui)))
106     (set! ti (+ (* (f64vector-ref ar ip) ui) (* (f64vector-ref ai ip) ur)))
107     (f64vector-set! ar ip (- (f64vector-ref ar i) tr))
108     (f64vector-set! ai ip (- (f64vector-ref ai i) ti))
109     (f64vector-set! ar i (+ (f64vector-ref ar i) tr))
110     (f64vector-set! ai i (+ (f64vector-ref ai i) ti))))
111   (set! tr (- (* ur wr) (* ui wi)))
112   (set! ti (+ (* ur wi) (* ui wr)))
113   (set! ur tr)
114   (set! ui ti))
115  #t))
116
117;;; the timer which does 10 calls on fft
118
119(define (fft-bench)
120 (do ((ntimes 0 (+ ntimes 1))) ((= ntimes 10))
121  (fft *re* *im*)))
122
123(fft-bench)
124(display count)                         ;(felix)
125(newline)
Note: See TracBrowser for help on using the repository browser.