fftw-fft.lisp

   1 (require :vol)
   2 (declaim (optimize (speed 2) (debug 3) (safety 3)))
   3 (in-package :vol)
   4
   5 (defconstant +forward+ 1)
   6 (defconstant +backward+ -1)
   7 (defconstant +estimate+ (ash 1 6))
   8
   9 (defmacro def-ffi (&optional (library "libfftw3.so") (prefix "fftw"))
  10   `(progn
  11      (load-shared-object ,(format nil "~a" library))
  12      (define-alien-routine (,(format nil "~a_execute" prefix)
  13                              ,(vol:format-symbol "~a-execute" prefix))
  14          void
  15        (plan (* int)))
  16      ,@(loop for rank in '(1 2 3) collect
  17             `(define-alien-routine
  18                  (,(format nil "~a_plan_dft_~ad" prefix rank)
  19                    ,(format-symbol "~a-plan-dft-~ad" prefix rank))
  20                  (* int)
  21                ,@(loop for i below rank collect
  22                       `(,(format-symbol "n~d" i) int))
  23                (in (* double-float))
  24                (out (* double-float))
  25                (sign int)
  26                (flags unsigned-int)))))
  27
  28 ;; define the foreign functions
  29 ;; fftw{f}-execute, fftw{f}-plan-dft-{1,2,3}d
  30 (def-ffi "libfftw3.so.3" "fftw")
  31 (def-ffi "libfftw3f.so.3" "fftwf")
  32
  33 #+nil
  34 (progn
  35   (load-shared-object "/usr/lib/libfftw3_threads.so")
  36
  37   (define-alien-routine ("fftw_init_threads" init-threads)
  38       int)
  39   (define-alien-routine ("fftw_plan_with_nthreads" plan-with-nthreads)
  40       void
  41     (nthreads int))
  42
  43   (defun init-ft ()
  44     (init-threads)
  45     (plan-with-nthreads 4)))
  46
  47 ;; to clean up completely call void fftw_cleanup_threads(void)
  48
  49
  50
  51 (def-generator (ft (rank type))
  52   (let* ((prefix (ecase type
  53                   (cdf 'fftw)
  54                   (csf 'fftwf)))
  55          (rlong-type (ecase type
  56                        (csf (get-long-type 'sf))
  57                        (cdf (get-long-type 'df))))
  58          (plan (format-symbol "~a-plan-dft-~ad" prefix rank))
  59          (execute (format-symbol "~a-execute" prefix)))
  60     `(defun ,name (in &key (forward t))
  61        (declare ((simple-array ,long-type ,rank) in)
  62                 (boolean forward)
  63                (values (simple-array ,long-type ,rank) &optional))
  64       (let* ((dims (array-dimensions in))
  65              (out (make-array dims :element-type ',long-type))
  66              (in-sap (sb-sys:vector-sap
  67                       (sb-ext:array-storage-vector in)))
  68              (out-sap (sb-sys:vector-sap
  69                        (sb-ext:array-storage-vector out)))
  70              (dir (if forward +forward+ +backward+)))
  71         (sb-sys:with-pinned-objects (in out)
  72           ,(ecase rank
  73                   (1 `(destructuring-bind (x) dims
  74                         (let ((p (,plan x in-sap out-sap
  75                                         dir +estimate+)))
  76                           (,execute p))))
  77                   (2 `(destructuring-bind (y x) dims
  78                         (let ((p (,plan y x in-sap out-sap
  79                                         dir +estimate+)))
  80                           (,execute p))))
  81                   (3 `(destructuring-bind (z y x) dims
  82                         (let ((p (,plan z y x in-sap out-sap
  83                                         dir +estimate+)))
  84                           (,execute p))))))
  85         (when forward ;; normalize if forward
  86           (s* (/ ,(coerce 1 rlong-type) (array-total-size out)) out))
  87         out))))
  88
  89 #+nil
  90 (def-ft-rank-type 1 csf)
  91
  92 (defmacro def-ft-functions (ranks types)
  93   (let* ((specifics nil)
  94          (cases nil)
  95          (name (format-symbol "ft")))
  96     (loop for rank in ranks do
  97          (loop for type in types do
  98               (let ((def-name (format-symbol "def-~a-rank-type" name))
  99                     (specific-name (format-symbol "~a-~a-~a" name rank type)))
 100                 (push `(,def-name ,rank ,type) specifics)
 101                 (push `((simple-array ,(get-long-type type) ,rank)
 102                         (,specific-name a :forward forward))
 103                       cases))))
 104     (store-new-function name)
 105     `(progn ,@specifics
 106             (defun ,name (a &key forward)
 107                (etypecase a
 108                  ,@cases
 109                  (t (error "The given type can't be handled with a generic ~a function." ',name)))))))
 110
 111 (def-ft-functions (1 2 3) (csf cdf))
 112
 113 (defmacro ift (in)
 114   `(ft ,in :forward nil))
 115
 116 #+nil
 117 (progn
 118   (time
 119    (let* ((nx 128)
 120           (ny nx)
 121           (nz ny)
 122           (a (convert-3-ub8/csf-mul
 123               (draw-sphere-ub8 20.0 nz ny nx))))
 124      (write-pgm "/home/martin/tmp/fftw.pgm"
 125                 (normalize-2-cdf/ub8-abs
 126                  (cross-section-xz-cdf (ft a)))))))