share/numeric/fft-interface.lisp

   1 ;; -*- Lisp -*-
   2
   3 ;; Maxima interface to the FFT routines and some various related
   4 ;; utilities.
   5
   6 (in-package :maxima)
   7
   8 ;; These functions perhaps need revision if they are
   9 ;; needed at all. Output ignores the type of input.
  10 (defun $recttopolar (rary iary)
  11   (let ((fast-rary (maxima-fft:fft-arg-check '$recttopolar rary))
  12         (fast-iary (maxima-fft:fft-arg-check '$recttopolar iary)))
  13     (maxima-fft:complex-to-polar fast-rary fast-iary)
  14     (list '(mlist) rary iary)))
  15
  16 (defun $polartorect (rary iary)
  17   (let ((fast-rary (maxima-fft:fft-arg-check '$polartorect rary))
  18         (fast-iary (maxima-fft:fft-arg-check '$polartorect iary)))
  19     (maxima-fft:polar-to-complex fast-rary fast-iary)
  20     (list '(mlist) rary iary)))
  21
  22 ;; Maxima's forward FFT routine.
  23 (defun $fft (input)
  24   (multiple-value-bind (z from-lisp)
  25       (maxima-fft:find-complex-fft-converters input)
  26     (let* ((size (length z))
  27            (order (maxima-fft:log-base2 size)))
  28       (unless (= size (ash 1 order))
  29         (merror "fft: size of input must be a power of 2, not ~M" size))
  30       (when (> (length z) 1)
  31         (setf z (maxima-fft:fft-r2-nn z)))
  32       (funcall from-lisp z))))
  33
  34 ;; Maxima's inverse FFT routine.
  35 (defun $inverse_fft (input)
  36   (multiple-value-bind (z from-lisp)
  37       (maxima-fft:find-complex-fft-converters input "inverse_fft")
  38     (let* ((size (length z))
  39            (order (maxima-fft:log-base2 size)))
  40       (unless (= size (ash 1 order))
  41         (merror "inverse_fft: size of input must be a power of 2, not ~M" size))
  42       (when (> (length z) 1)
  43         (setf z (maxima-fft:fft-r2-nn z :inverse-fft-p t)))
  44       (funcall from-lisp z))))
  45
  46 ;;; Maxima's RFFT function.
  47 ;;;
  48 ;;; The input is assumed to consist of objects that can be converted
  49 ;;; to a real-valued float.  If the input has length N, then the
  50 ;;; output is a complex reuslt of length N/2 + 1 because of the
  51 ;;; symmetry of the FFT for real inputs.
  52 (defun $real_fft (input)
  53   (multiple-value-bind (z from-lisp input-length)
  54       (maxima-fft:find-rfft-converters input)
  55     (declare (type (simple-array (complex double-float) (*)) z))
  56     (unless (= input-length (ash 1 (maxima-fft:log-base2 input-length)))
  57       (merror "real_fft: size of input must be a power of 2, not ~M" input-length))
  58     (let* ((n (ash (length z) 1))
  59            (result (make-array (1+ (length z)) :element-type '(complex double-float))))
  60       ;; This declaration causes CCL to generate incorrect code for at
  61       ;; least version 1.11 DarwinX8664 on OSX 10.11.3.  Disable it.
  62       #-ccl
  63       (declare (type (simple-array (complex double-float) (*)) result))
  64
  65       ;; We don't handle input of length 4 or less.  Use the complex
  66       ;; FFT routine to produce the desired (correct) output.
  67       (when (< n 3)
  68         (return-from $real_fft ($fft input)))
  69       (locally
  70           ;; Setting safety to 0 causes an internal compiler error with CCL 1.11.
  71           (declare (optimize (speed 3) (safety #-ccl 0 #+ccl 1)))
  72         ;; Compute FFT of shorter complex vector.
  73         (setf z (maxima-fft:fft-r2-nn z))
  74
  75         ;; Reconstruct the FFT of the original from the parts
  76         (setf (aref result 0)
  77               (complex (* 0.5
  78                           (+ (realpart (aref z 0))
  79                              (imagpart (aref z 0))))))
  80
  81         (let ((sincos (maxima-fft:sincos-table (maxima-fft:log-base2 n)))
  82               (n/2 (length z)))
  83           (declare (type (simple-array (complex double-float) (*)) sincos))
  84
  85           ;;(format t "n/2 = ~A~%" n/2)
  86           (loop for k of-type fixnum from 1 below (length z) do
  87             (setf (aref result k)
  88                   (* 0.25 (+ (+ (aref z k)
  89                                 (conjugate (aref z (- n/2 k))))
  90                              (* #c(0 -1d0)
  91                                 (aref sincos k)
  92                                 (- (aref z k)
  93                                    (conjugate (aref z (- n/2 k)))))))))
  94           (setf (aref result (length z))
  95                 (complex
  96                  (* 0.5
  97                     (- (realpart (aref z 0))
  98                        (imagpart (aref z 0))))))))
  99         (funcall from-lisp result))))
 100
 101 ;;; Maxima's inverse RFFT routine.
 102 ;;;
 103 ;;; The input is assumed to consist of objects that can be converted
 104 ;;; to a complex-valued floats.  The input MUST have a length that is
 105 ;;; one more than a power of two, as is returned by RFFT.
 106 (defun $inverse_real_fft (input)
 107   (multiple-value-bind (ft from-lisp)
 108       (maxima-fft:find-irfft-converters input)
 109     ;; declarations + (SPEED 3) tickles bug: https://bugs.launchpad.net/sbcl/+bug/1776091
 110     #-sbcl (declare (type (simple-array (complex double-float) (*)) ft))
 111     (let* ((n (1- (length ft))))
 112       (when (< n 2)
 113         ;; Just use the regular inverse fft to compute these values
 114         ;; because inverse_real_fft below doesn't work for these cases.
 115         (return-from $inverse_real_fft ($inverse_fft input)))
 116       (let* ((order (maxima-fft:log-base2 n))
 117              (sincos (maxima-fft:sincos-table (1+ order)))
 118              (z (make-array n :element-type '(complex double-float))))
 119         #-sbcl (declare (type (simple-array (complex double-float) (*)) sincos))
 120
 121         (unless (= n (ash 1 order))
 122           (merror "inverse_real_fft: input length must be one more than a power of two, not ~M" (1+ n)))
 123
 124         (locally
 125             #+sbcl nil #-sbcl (declare (optimize (speed 3)))
 126           (loop for k from 0 below n
 127                 do
 128                    (let ((evenpart (+ (aref ft k)
 129                                       (conjugate (aref ft (- n k)))))
 130                          (oddpart (* (- (aref ft k)
 131                                         (conjugate (aref ft (- n k))))
 132                                      (conjugate (aref sincos k)))))
 133                      (setf (aref z k) (+ evenpart
 134                                          (* #c(0 1d0) oddpart)))))
 135
 136           (setf z (maxima-fft:fft-r2-nn z :inverse-fft-p t)))
 137         (funcall from-lisp z)))))
 138
 139 ;; Bigfloat forward and inverse FFTs.  Length of the input must be a
 140 ;; power of 2.
 141 (defun $bf_fft (input)
 142   (multiple-value-bind (z from-lisp)
 143       (maxima-fft:find-bf-fft-converter input)
 144     (let ((n (length z)))
 145       (unless (= n (ash 1 (maxima-fft:log-base2 n)))
 146         (merror "bf_fft: size of input must be a power of 2, not ~M" n))
 147       (if (> n 1)
 148           (funcall from-lisp (bigfloat::fft-r2-nn z))
 149           (funcall from-lisp z)))))
 150
 151 (defun $bf_inverse_fft (input)
 152   (multiple-value-bind (x unconvert)
 153       (maxima-fft:find-bf-fft-converter input "bf_inverse_fft")
 154     (let ((n (length x)))
 155       (unless (= n (ash 1 (maxima-fft:log-base2 n)))
 156         (merror "bf_inverse_fft: size of input must be a power of 2, not ~M" n))
 157       (if (> n 1)
 158           (funcall unconvert (bigfloat::fft-r2-nn x :inverse-fft-p t))
 159           (funcall unconvert x)))))
 160
 161 (defun $bf_real_fft (input)
 162   (multiple-value-bind (z from-lisp)
 163       (maxima-fft:find-bf-rfft-converters input)
 164     (let* ((n (ash (length z) 1))
 165            (result (make-array (1+ (length z)))))
 166
 167       (when (< n 3)
 168         (return-from $bf_real_fft ($bf_fft input)))
 169
 170       ;; Compute FFT of shorter complex vector.  NOTE: the result
 171       ;; returned by bigfloat:fft has scaled the output by the length of
 172       ;; z.  That is, divided by n/2.  For our final result, we want to
 173       ;;     divide by n, so in the following bits of code, we have an
 174       ;;     extra factor of 2 to divide by.
 175       (setf z (bigfloat::fft-r2-nn z))
 176
 177       ;;(format t "z = ~A~%" z)
 178       ;; Reconstruct the FFT of the original from the parts
 179       (setf (aref result 0)
 180             (bigfloat:* 0.5
 181                         (bigfloat:+ (bigfloat:realpart (aref z 0))
 182                                     (bigfloat:imagpart (aref z 0)))))
 183
 184       (let ((sincos (bigfloat::sincos-table (maxima-fft:log-base2 n)))
 185             (n/2 (length z)))
 186         ;;(format t "n/2 = ~A~%" n/2)
 187         (loop for k from 1 below (length z) do
 188           (setf (aref result k)
 189                 (bigfloat:* 0.25 (bigfloat:+ (bigfloat:+ (aref z k)
 190                                                          (bigfloat:conjugate (aref z (- n/2 k))))
 191                                              (bigfloat:* #c(0 -1)
 192                                                          (aref sincos k)
 193                                                          (bigfloat:- (aref z k)
 194                                                                      (bigfloat:conjugate (aref z (- n/2 k)))))))))
 195         (setf (aref result (length z))
 196               (bigfloat:* 0.5
 197                           (bigfloat:- (bigfloat:realpart (aref z 0))
 198                                       (bigfloat:imagpart (aref z 0)))))
 199         (funcall from-lisp result)))))
 200
 201 (defun $bf_inverse_real_fft (input)
 202   (multiple-value-bind (ft from-lisp)
 203       (maxima-fft:find-bf-irfft-converters input)
 204   (let* ((n (1- (length ft))))
 205     (when (< n 2)
 206       ;; Just use the regular inverse fft to compute these values
 207       ;; because inverse_real_fft below doesn't work for these cases.
 208       (return-from $bf_inverse_real_fft ($bf_inverse_fft input)))
 209
 210     (let* ((z (make-array n))
 211            (order (maxima-fft:log-base2 n))
 212            (sincos (bigfloat::sincos-table (1+ order))))
 213
 214       (unless (= n (ash 1 order))
 215         (merror "bf_inverse_real_fft: input length must be one more than a power of two, not ~M" (1+ n)))
 216
 217       (loop for k from 0 below n
 218             do
 219                (let ((evenpart (bigfloat:+ (aref ft k)
 220                                            (bigfloat:conjugate (aref ft (- n k)))))
 221                      (oddpart (bigfloat:* (bigfloat:- (aref ft k)
 222                                                       (bigfloat:conjugate (aref ft (- n k))))
 223                                           (bigfloat:conjugate (aref sincos k)))))
 224                  (setf (aref z k)
 225                        (bigfloat:+ evenpart
 226                                    (bigfloat:* #c(0 1) oddpart)))))
 227
 228       ;;(format t "z = ~A~%" z)
 229       (let ((inverse (bigfloat::fft-r2-nn z :inverse-fft-p t)))
 230         ;;(format t "inverse = ~A~%" inverse)
 231         (funcall from-lisp inverse))))))
 232
 233