ikarus/benchmarks/r6rs-benchmarks/todo-src/fft.scm

80 lines
2.8 KiB
Scheme

;;; FFT - Fast Fourier Transform, translated from "Numerical Recipes in C"
(define (four1 data)
(let ((n (FLOATvector-length data))
(pi*2 6.28318530717959)) ; to compute the inverse, negate this value
; bit-reversal section
(let loop1 ((i 0) (j 0))
(if (< i n)
(begin
(if (< i j)
(begin
(let ((temp (FLOATvector-ref data i)))
(FLOATvector-set! data i (FLOATvector-ref data j))
(FLOATvector-set! data j temp))
(let ((temp (FLOATvector-ref data (+ i 1))))
(FLOATvector-set! data (+ i 1) (FLOATvector-ref data (+ j 1)))
(FLOATvector-set! data (+ j 1) temp))))
(let loop2 ((m (quotient n 2)) (j j))
(if (and (>= m 2) (>= j m))
(loop2 (quotient m 2) (- j m))
(loop1 (+ i 2) (+ j m)))))))
; Danielson-Lanczos section
(let loop3 ((mmax 2))
(if (< mmax n)
(let* ((theta
(FLOAT/ pi*2 (exact->inexact mmax)))
(wpr
(let ((x (FLOATsin (FLOAT* 0.5 theta))))
(FLOAT* -2.0 (FLOAT* x x))))
(wpi
(FLOATsin theta)))
(let loop4 ((wr 1.0) (wi 0.0) (m 0))
(if (< m mmax)
(begin
(let loop5 ((i m))
(if (< i n)
(let* ((j
(+ i mmax))
(tempr
(FLOAT-
(FLOAT* wr (FLOATvector-ref data j))
(FLOAT* wi (FLOATvector-ref data (+ j 1)))))
(tempi
(FLOAT+
(FLOAT* wr (FLOATvector-ref data (+ j 1)))
(FLOAT* wi (FLOATvector-ref data j)))))
(FLOATvector-set! data j
(FLOAT- (FLOATvector-ref data i) tempr))
(FLOATvector-set! data (+ j 1)
(FLOAT- (FLOATvector-ref data (+ i 1)) tempi))
(FLOATvector-set! data i
(FLOAT+ (FLOATvector-ref data i) tempr))
(FLOATvector-set! data (+ i 1)
(FLOAT+ (FLOATvector-ref data (+ i 1)) tempi))
(loop5 (+ j mmax)));***))
(loop4 (FLOAT+ (FLOAT- (FLOAT* wr wpr) (FLOAT* wi wpi)) wr)
(FLOAT+ (FLOAT+ (FLOAT* wi wpr) (FLOAT* wr wpi)) wi)
(+ m 2)))))
));******
(loop3 (* mmax 2)))))))
(define data
(FLOATmake-vector 1024 0.0))
(define (run data)
(four1 data)
(FLOATvector-ref data 0))
(define (main . args)
(run-benchmark
"fft"
fft-iters
(lambda (result) (equal? result 0.0))
(lambda (data) (lambda () (run data)))
data))