[plt-scheme] fft

Corey Sweeney wrote:
> I'm about to have a need soon to do fourier transformations.  My web 
> searches havn't been finding much.  Does anyone know of anything that 
> can do fft from scheme? (preferably plt)

If you want to do it in Scheme, then take a look at

     <http://www.scheme.com/tspl3/examples.html#./examples:h9>

or at the source at the bottom (taken from the benchmarks of Gambit).

If you need something fast, then take a look at FFTW and the
foreign function interface.

     <http://www.fftw.org/>


Source taken from the Gambit benchmarks:

;;; 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)))
(vector-length (vector));************ remove this!!!!!
           (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)
   (four1 data)
   (FLOATvector-ref data 0))

(define (main . args)
   (run-benchmark
     "fft"
     fft-iters
     (lambda () (run))
     (lambda (result) (equal? result 0.0))))

-- 
Jens Axel Søgaard