share/fftpack5/fftpack5-interface.lisp

   1 ;;; -*- Mode: lisp -*-
   2
   3 (in-package #:maxima)
   4
   5 (in-package :fftpack5)
   6
   7 (defun mlist->complex-cl-vector (object)
   8   "Convert Maxima list to a (simple-array (complex double-float) (*))
   9   with the same values."
  10   (let ((z (make-array (1- (length object)) :element-type '(complex double-float))))
  11     (loop for k of-type fixnum from 0
  12           for x in (cdr object)
  13           while x
  14           do
  15              (let ((fl (maxima::risplit (maxima::$float x))))
  16                (setf (aref z k) (complex (car fl) (cdr fl)))))
  17     z))
  18
  19 (defun complex-cl-vector->mlist (array)
  20   "Convert a CL array of complex values to a Maxima list containing
  21   the same values."
  22   (declare (type (simple-array (complex double-float) (*)) array))
  23   (cons '(maxima::mlist maxima::simp)
  24         (loop for w of-type (complex double-float) across array
  25               collect (maxima::add (realpart w)
  26                                    (maxima::mul (imagpart w) 'maxima::$%i)))))
  27
  28 (defun array->complex-cl-vector (array)
  29   "Convert a CL vector of maxima values to a (simsple-array (complex
  30   double-float) (*)) with the same complex values."
  31   (let* ((n (length array))
  32          (z (make-array n :element-type '(complex double-float))))
  33     (dotimes (k n)
  34       (let ((fl (maxima::risplit (maxima::$float (aref array k)))))
  35         (setf (aref z k) (complex (car fl) (cdr fl)))))
  36     z))
  37
  38 (defun complex-cl-vector->vector (array)
  39   "Convert (simple-array (complex double-float) (*)) to a CL vector of
  40   maxima expressions of the same value."
  41   (declare (type (simple-array (complex double-float) (*)) array))
  42   (let* ((n (length array))
  43          (z (make-array n)))
  44     (dotimes (k n)
  45       (let ((item (aref array k)))
  46         (setf (aref z k) (maxima::add (realpart item)
  47                                       (maxima::mul 'maxima::$%i (imagpart item))))))
  48     z))
  49
  50 (defun find-complex-fft-converters (object &optional (maxima-function-name "fftpack5"))
  51   "Convert a Maxima OBJECT to a specialized CL vector suitable for the
  52   complex FFT routines.  Two values are returned: The specialized
  53   vector and a function that will convert the result of the FFT
  54   routine to a Maxima object of the same type as
  55   OBJECT. MAXIMA-FUNCTION-NAME is a string to use for error messages."
  56   ;;(format t "object = ~S~%" object)
  57   ;;(format t "arrayp = ~S~%" (arrayp object))
  58   ;;(describe 'arrayp)
  59   (cond ((maxima::$listp object)
  60          (values
  61           (mlist->complex-cl-vector object)
  62           #'complex-cl-vector->mlist))
  63         ((arrayp object)
  64          (values
  65           (array->complex-cl-vector object)
  66           #'complex-cl-vector->vector))
  67         ((and (symbolp object) (maxima::symbol-array (maxima::mget object 'maxima::array)))
  68          (values
  69           (let* ((sym-array (maxima::symbol-array (maxima::mget object 'maxima::array))))
  70             (array->complex-cl-vector sym-array))
  71           #'(lambda (array)
  72               (let ((ar (complex-cl-vector->vector array))
  73                     (sym (maxima::meval `((maxima::$array)
  74                                           ,(intern (symbol-name (gensym "$G")) "MAXIMA")
  75                                           maxima::$float
  76                                           ,(1- (length array))))))
  77                 (setf (maxima::symbol-array (maxima::mget sym 'maxima::array))
  78                       ar)
  79                 sym))))
  80         (t
  81          (maxima::merror "~A: input is not a list or an array: ~M." maxima-function-name object))))
  82
  83 (defun mlist->rfft-vector (object)
  84   "Convert Maxima list to a (simple-array (complex double-float) (*))
  85   with the same values suitable for as the input for the RFFT routine.
  86   The maxima list is assumed to consist only of real values."
  87   (let* ((n (length (cdr object)))
  88          (z (make-array n :element-type 'double-float)))
  89     (loop for k from 0
  90           for x in (cdr (maxima::$float object))
  91           do
  92              (setf (aref z k) x))
  93     z))
  94
  95 (defun array->rfft-vector (array)
  96   "Convert CL array of real (Maxima) values to a (simple-array
  97   (complex double-float) (*)) with the same values suitable for as the
  98   input for the RFFT routine.  The CL array is assumed to consist only
  99   of real values."
 100   (let* ((n (length array))
 101          (z (make-array n :element-type 'double-float)))
 102     (loop for k from 0 below n
 103           for m from 0
 104           do
 105              (setf (aref z m) (maxima::$float (aref array k))))
 106     z))
 107
 108 (defun find-rfft-converters (object)
 109   "Convert a Maxima OBJECT to a specialized CL vector suitable for the
 110   RFFT routines.  Two values are returned: The specialized vector and
 111   a function that will convert the result of the RFFT routine to a
 112   Maxima object of the same type as OBJECT. MAXIMA-FUNCTION-NAME is a
 113   string to use for error messages."
 114   (cond ((maxima::$listp object)
 115          (values
 116           (mlist->rfft-vector object)
 117           #'complex-cl-vector->mlist))
 118         ((arrayp object)
 119          (values
 120           (let ((z (make-array  (length object) :element-type '(complex double-float))))
 121             (loop for k of-type fixnum from 0
 122                   for x across object
 123                   do
 124                      (let ((fl (maxima::risplit (maxima::$float x))))
 125                        (setf (aref z k) (complex (car fl) (cdr fl)))))
 126             z)
 127           #'complex-cl-vector->vector))
 128         ((and (symbolp object) (maxima::symbol-array (maxima::mget object 'maxima::array)))
 129          (let ((sym-array (maxima::symbol-array (maxima::mget object 'maxima::array))))
 130            (values
 131             (array->rfft-vector sym-array)
 132             #'(lambda (array)
 133                 (let ((ar (complex-cl-vector->vector array))
 134                       (sym (maxima::meval `((maxima::$array)
 135                                     ,(intern (symbol-name (gensym "$G")))
 136                                     maxima::$float
 137                                     ,(1- (length array))))))
 138                   (setf (maxima::symbol-array (maxima::mget sym 'maxima::array))
 139                         ar))))))
 140         (t
 141          (maxima::merror "real_fft: input is not a list or an array: ~M." object))))
 142
 143 (defun find-irfft-converters (object)
 144   "Convert a Maxima OBJECT to a specialized CL vector suitable for the
 145   inverse RFFT routines.  Two values are returned: The specialized
 146   vector and a function that will convert the result of the inverse
 147   RFFT routine to a Maxima object of the same type as
 148   OBJECT. MAXIMA-FUNCTION-NAME is a string to use for error messages."
 149   (cond ((maxima::$listp object)
 150          (values
 151           (mlist->complex-cl-vector object)
 152           #'(lambda (obj)
 153               (let (result)
 154                 (map nil #'(lambda (z)
 155                              (push z result))
 156                      obj)
 157                 (cons '(maxima::mlist maxima::simp) (nreverse result))))))
 158         ((arrayp object)
 159          (values
 160           (array->complex-cl-vector object)
 161           #'(lambda (obj)
 162               (let ((result (make-array (length obj))))
 163                 (loop for k from 0 by 2
 164                       for z across obj
 165                       do
 166                          (progn
 167                            (setf (aref result k) z)))
 168                 result))))
 169         (t
 170          (maxima::merror "inverse_real_fft: input is not a list or an array: ~M." object))))
 171
 172 (in-package :maxima)
 173
 174 (defmfun $fftpack5_fft (input)
 175   "fftpack5_fft(z)
 176
 177   Computes the FFT of the complex-valued vector <input>.  The <input>
 178   may have any length, but the FFT is most efficient if the length has
 179   the form 2^r*3^s*5^t.
 180
 181   Let x be the input and y be the transformed value.  Then the FFT is
 182   defined by
 183
 184   y[k] = (1/n) sum(x[j] exp(+2 %i %pi j k / n), j, 0, n - 1)"
 185
 186   (multiple-value-bind (z from-lisp)
 187       (fftpack5::find-complex-fft-converters input)
 188     ;; To match the fft package, we want to do the inverse FFT here
 189     ;; and then we need to scale the result by 1/N.
 190     (fftpack5:inverse-cfft z)
 191     (let* ((n (length z)) (n-float (coerce n 'double-float)))
 192       (dotimes (k n)
 193         (setf (aref z k) (/ (aref z k) n-float)))
 194       (funcall from-lisp z))))
 195
 196 (defmfun $fftpack5_inverse_fft (input)
 197   "fftpack5_inverse_fft(z)
 198
 199   Computes the inverse FFT of the complex-valued vector <input>.  The <input>
 200   may have any length, but the FFT is most efficient if the length has
 201   the form 2^r*3^s*5^t.
 202
 203   Let y be the input and x be the inverse-transformed value.  Then the
 204   inverse FFT is defined by
 205
 206   x[j] = sum(y[k] exp(-2 %i %pi j k / n), k, 0, n - 1)"
 207
 208   (multiple-value-bind (z from-lisp)
 209       (fftpack5::find-complex-fft-converters input)
 210     ;; To match the fft package, we want to do the forward FFT here
 211     ;; and scale the result by N.
 212     (fftpack5:cfft z)
 213     (let* ((n (length z)) (n-float (coerce n 'double-float)))
 214       (dotimes (k n)
 215         (setf (aref z k) (* n-float (aref z k))))
 216       (funcall from-lisp z))))
 217
 218 (defmfun $fftpack5_real_fft (input)
 219   "fftpack5_real_fft(input)
 220
 221   Computes the FFT of the real-valued vector <input>.  The <input>
 222   may have any length, but the FFT is most efficient if the length has
 223   the form 2^r*3^s*5^t.
 224
 225   No check is made to ensure that the input contains only real values.
 226
 227   Let x be the input and y be the transformed value.  Then the FFT is
 228   defined by
 229
 230   y[k] = (1/n) sum(x[j] exp(+2 %i %pi j k / n), j, 0, n - 1)
 231
 232   for k = 0, 1, ... floor(n/2).  Because of the symmetries of the FFT,
 233   y[0] is real.  If n is even y[n/2] is real.  If n is odd,
 234   y[floor(n/2)] may be complex.
 235
 236   Use fftpack5_inverse_rfft to compute the inverse FFT"
 237
 238   (multiple-value-bind (z from-lisp)
 239       (fftpack5::find-rfft-converters input)
 240     (funcall from-lisp (fftpack5::rfft z))))
 241
 242 (defmfun $fftpack5_inverse_real_fft (input len)
 243   "fftpack5_inverse_real_fft(input, len)
 244
 245   Computes the inverse FFT of <input>.  The <input> must be in the
 246   form returned by fftpack5_real_fft.  Note that <len> MUST be
 247   provided and should be the length of the inverse transform (which
 248   must be the length of the original vector that was transformed).
 249   This is needed because the length of <input> cannot not uniquely
 250   determine the length of the output.
 251 "
 252
 253   (multiple-value-bind (z from-lisp)
 254       (fftpack5::find-irfft-converters input)
 255     (funcall from-lisp (fftpack5::inverse-rfft z len))))