gistfile1.rkt

#lang racket

(provide (all-defined-out))

(require ffi/unsafe
         ffi/unsafe/define)

;; on systems machines use this:
;;(define-ffi-definer define-fftw (ffi-lib "libfftw3.so.3"))
(define-ffi-definer define-fftw (ffi-lib "libfftw3"))

(define _fftw_complex (_array _double 2))

;; creates a c array from a racket vector
(define (complex-vector->fftw_complex* vs)
  (let* ((N (vector-length vs))
         (block (fftw_malloc (* N (ctype-sizeof _fftw_complex)))))
    (for ((v (in-vector vs))
          (n (in-naturals)))
      (ptr-set! block _double (+ (* n 2) 0) (exact->inexact (real-part v)))
      (ptr-set! block _double (+ (* n 2) 1) (exact->inexact (imag-part v))))
    block))

;; creates a racket vector from a c array given some length
(define (fftw_complex*->complex-vector ptr N)
  (for/vector ((n (in-range 0 N)))
    (make-rectangular (array-ref (ptr-ref ptr _fftw_complex n) 0)
                      (array-ref (ptr-ref ptr _fftw_complex n) 1))))

;; some built-in FFTW malloc? not sure why I can't use malloc
(define-fftw fftw_malloc (_fun (N : _int) -> _pointer))
;; produces FFTW plans
(define-fftw fftw_plan_dft_1d
  (_fun (N : _int)
        (in : _pointer) ;; fftw_complex*
        (out : _pointer) ;; fftw_complex*
        (sign : _int)
        (flags : _uint)
        ->
        (plan : _pointer)))
;; runs a plan
(define-fftw fftw_execute
  (_fun _pointer -> _void))

(define (fft v)
  (let* ((N (vector-length v))
         (in (complex-vector->fftw_complex* v))
         (out (fftw_malloc (* N (ctype-sizeof _fftw_complex)))))

    (displayln (fftw_complex*->complex-vector in N))
    (fftw_execute (fftw_plan_dft_1d N
                                    in
                                    out
                                    -1
                                    #b100000))
    (displayln (fftw_complex*->complex-vector in N))
    (fftw_complex*->complex-vector out N)))