share/fftpack5/lisp/fftpack5.lisp

   1 (defpackage :fftpack5
   2   (:use :common-lisp)
   3   (:export "RFFT"
   4            "INVERSE-RFFT"))
   5
   6 (in-package :fftpack5)
   7
   8 (defvar *wsave-cache*
   9   (make-hash-table)
  10   "Cache for different wsave tables.  The key is the FFT size; the
  11   value is a the wsave table needed by the FFT routines.")
  12
  13 (defun get-wsave-entry (n)
  14   "Get the wsave array and it's length that is needed for computing
  15   the (forward and inverse) FFTs.  The value is cached, so if it's the
  16   cache, return it.  Otherwise compute a new value and save it in the
  17   cache."
  18   (let ((entry (gethash n *wsave-cache*)))
  19     (if entry
  20         entry
  21         (let* ((lensav (+ n 4 (floor (log n 2))))
  22                (wsave (make-array lensav :element-type 'single-float)))
  23           (multiple-value-bind (ignore-0 ignore-1 ignore-2 ier)
  24               (rfft1i n wsave lensav 0)
  25             (declare (ignore ignore-0 ignore-1 ignore-2))
  26             (unless (zerop ier)
  27               ;; This shouldn't really ever happen.
  28               (error "lensav is not big enough"))
  29             (setf (gethash n *wsave-cache*) wsave)
  30             wsave)))))
  31
  32 (defun convert-rfft (x)
  33   "Convert the output of FFTPACK RFFT1F (forward real FFT) into a more
  34   user-friendly format with complex values"
  35   (declare (type (simple-array single-float (*)) x))
  36   (let* ((n (length x))
  37          (nhalf (floor (/ n 2)))
  38          (out (make-array (+ 1 nhalf) :element-type '(complex single-float))))
  39     ;; If X is the transformed value, then the output from rfftf is:
  40     ;; 0: X(0) / N
  41     ;; 1: 2*realpart(X(1))/N
  42     ;; 2: 2*imagpart(X(1))/N
  43     ;; 3: 2*realpart(X(2))/N
  44     ;; 4: 2*realpart(X(2))/N
  45     ;; ...
  46     ;; N-1: X(N-1)/N
  47     ;; The last term exists only if N is even.
  48     (setf (aref out 0) (complex (aref x 0) 0.0))
  49     (loop for j from 1 to (if (evenp n) (1- nhalf) nhalf)
  50           for k from 1 by 2
  51           do
  52              ;; Need to remove the factor of 2 that rfftf added.
  53              (setf (aref out j) (complex (* 0.5 (aref x k))
  54                                          (* 0.5 (aref x (1+ k))))))
  55     (when (evenp n)
  56       (setf (aref out nhalf)
  57             (complex (aref x (1- n)) 0.0)))
  58     out))
  59
  60 (defun convert-inverse-rfft (x n)
  61   "Convert the complex-valued input, X, (that was produced by rfft)
  62   into the form needed by rfft1b.  The length of the transform, N, is
  63   needed because it cannot be uniquely determined from the length of
  64   X."
  65   (declare (type (simple-array (complex single-float) (*)) x))
  66   (let ((res (make-array n :element-type 'single-float)))
  67     (setf (aref res 0) (realpart (aref x 0)))
  68     (loop for j from 1 below (if (evenp n) (1- (length x)) (length x))
  69           for k from 1 by 2
  70           do
  71              (let ((z (aref x j)))
  72                ;; Put back the factor of 2 that we removed in rfft but
  73                ;; is needed by rfft1b.
  74                (setf (aref res k) (* 2 (realpart z)))
  75                (setf (aref res (1+ k)) (* 2 (imagpart z)))))
  76     (when (evenp n)
  77       (setf (aref res (1- n)) (realpart (aref x (1- (length x))))))
  78     res))
  79
  80 (defun rfft (x)
  81   "Compute the real FFT of X.
  82
  83    Let N be the length of X.  The FFT is:
  84
  85     Y[n] = 1/N*sum(x[k] * exp(2*%i*%pi*k*n/N), k, 0, N-1)
  86
  87    for n = 0, 1,...,floor(N/2)+1
  88
  89    NOTE: This differs from the typical engineering definition of the
  90    forward FFT, where the transform is often not normalized by the
  91    length and the argument to exp has a negative sign.  This is
  92    generally considered the inverse FFT in engineering"
  93
  94   (declare (type (simple-array single-float (*)) x))
  95   (let* ((n (length x))
  96          (work (make-array n :element-type 'single-float))
  97          (wsave (get-wsave-entry n))
  98          (lensav (length wsave)))
  99     (let ((ier
 100             (nth-value 8
 101                        (rfft1f n 1 x n wsave lensav work n 0))))
 102       (unless (zerop ier)
 103         ;; This should never happen because we always allocate enough
 104         ;; space for x, wsave, and the work array.
 105         (error "rfft1f failed with code ~A" ier))
 106       (convert-rfft x))))
 107
 108 (defun inverse-rfft (x n)
 109   "Compute the inverse real FFT of X. N is the length of the transform
 110   (and not the length of X!).
 111
 112    See RFFT for the definition of the FFT used in these routines."
 113   (declare (type (simple-array (complex single-float) (*)) x))
 114
 115   ;; Check to see if N is consistent with the length of X.  However,
 116   ;; this isn't foolproof.  For a real FFT of length 5, the FFT array
 117   ;; has length 3.  But a real FFT of length 4 also has an array
 118   ;; length of 3.
 119   (unless (= (length x) (1+ (floor n 2)))
 120     (error "Length of X (~A) is not compatible with N (~D)"
 121            (length x) n))
 122
 123   (let* ((inv (convert-inverse-rfft x n))
 124          (work (make-array n :element-type 'single-float))
 125          (wsave (get-wsave-entry n))
 126          (lensav (length wsave)))
 127     (let* ((ier
 128              (progn
 129                (nth-value 8
 130                           (rfft1b n 1 inv n wsave lensav work n 0)))))
 131       (unless (zerop ier)
 132         ;; This should never happen because we should have always
 133         ;; allocated the correct size of the arrays.
 134         (error "rfft1b failed with code ~A" ier))
 135       inv)))
 136
 137 (defun test-rfft (n &key verbose)
 138   "Test the rfft and inverse-rfft routines using a simple ramp of
 139   length n"
 140   ;; For testing, use a simple ramp: 1, 2, 3, ..., N.
 141   ;; With some care, we can derive the FFT for this is
 142   ;;
 143   ;; X[0] = (N + 1) / 2
 144   ;; X[n] = -1/2 - cot(%pi*n/N)/2, n = 1, N
 145
 146   (let* ((x (make-array n :element-type 'single-float)))
 147     (loop for k from 0 below n
 148           do
 149              (setf (aref x k) (float (+ k 1) 1f0)))
 150     (let* ((xfrm (rfft x))
 151            (xfrm-len (length xfrm))
 152            (expected (make-array xfrm-len :element-type '(complex single-float)))
 153            (noise-pwr 0.0)
 154            (signal-pwr 0.0))
 155       (declare (single-float noise-pwr signal-pwr))
 156
 157       (setf (aref expected 0) (complex (* 0.5 (+ n 1)) 0.0))
 158       (loop for k from 1 below xfrm-len
 159             with omega = (coerce (/ pi n) 'single-float)
 160             do
 161                (setf (aref expected k)
 162                      (complex -0.5 (/ -0.5 (tan (* omega k))))))
 163       (when verbose
 164         (format t "Forward transform; actual vs expected~%")
 165         (loop for k from 0 below xfrm-len
 166               do
 167                  (format t "~4d: ~A ~A~%" k (aref xfrm k) (aref expected k))))
 168       (incf noise-pwr (expt (abs (- (aref xfrm 0) (* 0.5 (+ n 1)))) 2))
 169       (incf signal-pwr (expt (* 0.5 (+ n 1)) 2))
 170       (loop for k from 1 below (length xfrm)
 171             do
 172                (incf noise-pwr (expt (abs (- (aref xfrm k)
 173                                              (aref expected k)))
 174                                      2))
 175                (incf signal-pwr (expt (abs (aref expected k)) 2)))
 176       (let ((inv (inverse-rfft xfrm n))
 177             (inv-noise-pwr 0.0)
 178             (inv-signal-pwr 0.0))
 179         (when verbose
 180           (format t "Inverse transform; actual vs expected~%")
 181           (loop for k from 0 below n
 182                 do
 183                    (format t "~4d: ~A ~A~%" k (aref inv k) (+ k 1))))
 184         (loop for k from 0 below n
 185               do
 186                  (incf inv-noise-pwr (expt (- (aref inv k) (+ k 1)) 2))
 187                  (incf inv-signal-pwr (expt (coerce (+ k 1) 'single-float) 2)))
 188         (flet ((db (s n)
 189                  (if (zerop n)
 190                      1000.0
 191                      (* 10 (log (/ s n) 10)))))
 192           (values (db signal-pwr noise-pwr) (db inv-signal-pwr inv-noise-pwr) noise-pwr signal-pwr inv-noise-pwr inv-signal-pwr))))))
 193