  FFT - Fast Fourier Transform , translated from " Numerical Recipes in C "

(##define-macro (def-macro form . body)
  `(##define-macro ,form (let () ,@body)))

(def-macro (FLOATvector-const . lst) `(f64vector ,@lst))
(def-macro (FLOATvector-ref v i)     `(f64vector-ref ,v ,i))
(def-macro (FLOATvector-set! v i n)  `(f64vector-set! ,v ,i ,n))
(def-macro (FLOATmake-vector l i)    `(make-f64vector ,l ,i))
(def-macro (FLOATvector-length v)    `(f64vector-length ,v))
(def-macro (FLOATsin v)    `(sin ,v))
(def-macro (FLOAT= . lst)    `(= ,@lst))
(def-macro (FLOAT> . lst)    `(> ,@lst))
(def-macro (FLOAT>= . lst)    `(>= ,@lst))
(def-macro (FLOAT< . lst)    `(< ,@lst))
(def-macro (FLOAT<= . lst)    `(<= ,@lst))
(def-macro (FLOAT+ . lst)    `(+ ,@lst))
(def-macro (FLOAT- . lst)    `(- ,@lst))
(def-macro (FLOAT* . lst)    `(* ,@lst))
(def-macro (FLOAT/ . lst)    `(/ ,@lst))

(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)))))))

     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 (run data)
  (four1 data)
  (FLOATvector-ref data 0))

(pp (run (FLOATmake-vector 1024 0.0)))
(pp (run (FLOATmake-vector 1024 42.0)))

;0.
21504 .

 to compute the inverse, negate this value
 bit-reversal section
***))
******
0.
  FFT - Fast Fourier Transform , translated from " Numerical Recipes in C "

(##define-macro (def-macro form . body)
  `(##define-macro ,form (let () ,@body)))

(def-macro (FLOATvector-const . lst) `(f64vector ,@lst))
(def-macro (FLOATvector-ref v i)     `(f64vector-ref ,v ,i))
(def-macro (FLOATvector-set! v i n)  `(f64vector-set! ,v ,i ,n))
(def-macro (FLOATmake-vector l i)    `(make-f64vector ,l ,i))
(def-macro (FLOATvector-length v)    `(f64vector-length ,v))
(def-macro (FLOATsin v)    `(sin ,v))
(def-macro (FLOAT= . lst)    `(= ,@lst))
(def-macro (FLOAT> . lst)    `(> ,@lst))
(def-macro (FLOAT>= . lst)    `(>= ,@lst))
(def-macro (FLOAT< . lst)    `(< ,@lst))
(def-macro (FLOAT<= . lst)    `(<= ,@lst))
(def-macro (FLOAT+ . lst)    `(+ ,@lst))
(def-macro (FLOAT- . lst)    `(- ,@lst))
(def-macro (FLOAT* . lst)    `(* ,@lst))
(def-macro (FLOAT/ . lst)    `(/ ,@lst))

(define (four1 data)
  (let ((n (FLOATvector-length data))


    (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)))))))

     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))
                (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 (run data)
  (four1 data)
  (FLOATvector-ref data 0))

(pp (run (FLOATmake-vector 1024 0.0)))
(pp (run (FLOATmake-vector 1024 42.0)))

21504 .
