1(define times 1000) 2(define size 1024) 3 4(define (b_fft areal aimag) 5 (let ((ar 0) 6 (ai 0) 7 (i 1) 8 (j 1) 9 (k 0) 10 (m 0) 11 (n 0) 12 (le 1) 13 (le1 0) (le2 0) 14 (ip 0) 15 (nv2 0) 16 (nm1 0) 17 (ur 0.0) 18 (ui 0.0) 19 (wr 0.0) 20 (wi 0.0) 21 (tr 0.0) 22 (ti 0.0)) 23 ;; initialize 24 (set! ar areal) 25 (set! ai aimag) 26 (set! n (length ar)) 27 (set! n (- n 1)) 28 (set! nv2 (quotient n 2)) 29 (set! nm1 (- n 1)) 30 31 (do () 32 ((>= i n)) 33 (set! m (+ m 1)) 34 (set! i (+ i i))) 35 (set! i 1) 36 37 (do () 38 ((>= i n)) 39 (when (< i j) 40 (set! tr (float-vector-ref ar j)) 41 (set! ti (float-vector-ref ai j)) 42 (float-vector-set! ar j (float-vector-ref ar i)) 43 (float-vector-set! ai j (float-vector-ref ai i)) 44 (float-vector-set! ar i tr) 45 (float-vector-set! ai i ti)) 46 (set! k nv2) 47 (do () 48 ((>= k j)) 49 (set! j (- j k)) 50 (set! k (quotient k 2))) 51 52 (set! j (+ j k)) 53 (set! i (+ i 1))) 54 55 (do ((l 1 (+ l 1))) 56 ((> l m)) 57 (set! le1 le) 58 (set! le2 (+ le 1)) 59 (set! le (* le 2)) 60 (set! ur 1.0) 61 (set! ui 0.) 62 (set! wr (cos (/ pi le1))) 63 (set! wi (sin (/ pi le1))) 64 (do ((j1 1 (+ j1 1))) 65 ((= j1 le2)) 66 (do ((i1 j1 (+ i1 le))) 67 ((> i1 n)) 68 (set! ip (+ i1 le1)) 69 (set! tr (- (* (float-vector-ref ar ip) ur) 70 (* (float-vector-ref ai ip) ui))) 71 (set! ti (+ (* (float-vector-ref ar ip) ui) 72 (* (float-vector-ref ai ip) ur))) 73 (float-vector-set! ar ip (- (float-vector-ref ar i1) tr)) 74 (float-vector-set! ai ip (- (float-vector-ref ai i1) ti)) 75 (float-vector-set! ar i1 (+ (float-vector-ref ar i1) tr)) 76 (float-vector-set! ai i1 (+ (float-vector-ref ai i1) ti)))) 77 (set! tr (- (* ur wr) (* ui wi))) 78 (set! ti (+ (* ur wi) (* ui wr))) 79 (set! ur tr) 80 (set! ui ti)) 81 #t)) 82 83(if (not (defined? 'complex)) 84 (define complex make-rectangular)) 85(if (not (defined? 'when)) 86 (define-macro (when test . body) `(if ,test (begin ,@body)))) 87 88(define* (cfft data n (dir 1)) ; complex data 89 (unless n (set! n (length data))) 90 (do ((i 0 (+ i 1)) 91 (j 0)) 92 ((= i n)) 93 (if (> j i) 94 (let ((temp (data j))) 95 (set! (data j) (data i)) 96 (set! (data i) temp))) 97 (do ((m (/ n 2) (/ m 2))) 98 ((or (< m 2) 99 (< j m)) 100 (set! j (+ j m))) 101 (set! j (- j m)))) 102 (do ((ipow (floor (log n 2))) 103 (prev 1) 104 (lg 0 (+ lg 1)) 105 (mmax 2 (* mmax 2)) 106 (pow (/ n 2) (/ pow 2)) 107 (theta (complex 0.0 (* pi dir)) (* theta 0.5))) 108 ((= lg ipow)) 109 (do ((wpc (exp theta)) 110 (wc 1.0) 111 (ii 0 (+ ii 1))) 112 ((= ii prev) 113 (set! prev mmax)) 114 (do ((jj 0 (+ jj 1)) 115 (i ii (+ i mmax)) 116 (j (+ ii prev) (+ j mmax))) 117 ((>= jj pow)) 118 (let ((tc (* wc (data j)))) 119 (set! (data j) (- (data i) tc)) 120 (set! (data i) (+ (data i) tc)))) 121 (set! wc (* wc wpc)))) 122 data) 123 124(when (defined? 'equivalent?) 125 (unless (equivalent? (cfft (vector 0.0 1+i 0.0 0.0)) #(1+1i -1+1i -1-1i 1-1i)) 126 (format *stderr* "cfft 1: ~S~%" (cfft (vector 0.0 1+i 0.0 0.0)))) 127 (let-temporarily (((*s7* 'equivalent-float-epsilon) 1e-14)) 128 (unless (equivalent? (cfft (vector 0 0 1+i 0 0 0 1-i 0)) #(2 -2 -2 2 2 -2 -2 2)) 129 (format *stderr* "cfft 2: ~S~%" (cfft (vector 0 0 1+i 0 0 0 1-i 0)))))) 130 131(define (fft-bench) 132 (let ((*re* (make-float-vector (+ size 1) 0.10)) 133 (*im* (make-float-vector (+ size 1) 0.10))) 134 (do ((ntimes 0 (+ ntimes 1))) 135 ((= ntimes times)) 136 (b_fft *re* *im*))) 137 138 (let* ((n 256) 139 (cdata (make-vector n 0.0))) 140 (do ((i 0 (+ i 1))) 141 ((= i times)) 142 (fill! cdata 0.0) 143 (vector-set! cdata 2 1+i) 144 (vector-set! cdata (- n 1) 1-i) 145 (cfft cdata)))) 146 147(fft-bench) 148 149(when (> (*s7* 'profile) 0) 150 (show-profile 200)) 151(exit) 152