Commit 4f8b8372 authored by Raymond Toy's avatar Raymond Toy
Browse files

Factor out common routine for getting the wsave entry

parent 54c40a35
...@@ -6,11 +6,13 @@ ...@@ -6,11 +6,13 @@
(defvar *wsave-cache* (defvar *wsave-cache*
(make-hash-table) (make-hash-table)
"Cache for different wsave tables. The key is the FFT size; the "Cache for different wsave tables. The key is the FFT size; the
value is the wsave table needed by the FFT routines.") value is a the wsave table needed by the FFT routines.")
(defun convert-rfft (x) (defun convert-rfft (x)
(declare (type (simple-array single-float (*)) x)) "Convert the output of FFTPACK RFFT1F (forward real FFT) into a more
(let* ((n (length x)) user-friendly format with complex values"
(declare (type
(simple-array single-float (*)) x)) (let* ((n (length x))
(nhalf (floor (/ n 2))) (nhalf (floor (/ n 2)))
(out (make-array (+ 1 nhalf) :element-type '(complex single-float)))) (out (make-array (+ 1 nhalf) :element-type '(complex single-float))))
(setf (aref out 0) (complex (aref x 0) 0.0)) (setf (aref out 0) (complex (aref x 0) 0.0))
...@@ -24,23 +26,34 @@ ...@@ -24,23 +26,34 @@
(complex (aref x (1- n)) 0.0))) (complex (aref x (1- n)) 0.0)))
out)) out))
(defun rfft (x) (defun get-wsave-entry (n)
(declare (type (simple-array single-float (*)) x)) "Get the wsave array and it's length that is needed for computing
;; Initialize the wave table if needed the (forward and inverse) FFTs. The value is cached, so if it's the
(let* ((n (length x)) cache, return it. Otherwise compute a new value and save it in the
(lensav (+ n (floor (log n 2)) 4)) cache."
(work (make-array n :element-type 'single-float))) (let ((entry (gethash n *wsave-cache*)))
(unless (gethash n *wsave-cache*) (if entry
(let* ((wsave (make-array lensav :element-type 'single-float))) entry
(let* ((lensav (+ n 4 (floor (log n 2))))
(wsave (make-array lensav :element-type 'single-float)))
(multiple-value-bind (ignore-0 ignore-1 ignore-2 ier) (multiple-value-bind (ignore-0 ignore-1 ignore-2 ier)
(rfft1i n wsave lensav 0) (rfft1i n wsave lensav 0)
(declare (ignore ignore-0 ignore-1 ignore-2)) (declare (ignore ignore-0 ignore-1 ignore-2))
(unless (zerop ier) (unless (zerop ier)
(error "lensav is not big enough")) (error "lensav is not big enough"))
(setf (gethash n *wsave-cache*) wsave)))) (setf (gethash n *wsave-cache*) wsave)
wsave)))))
(defun rfft (x)
(declare (type (simple-array single-float (*)) x))
;; Initialize the wave table if needed
(let* ((n (length x))
(work (make-array n :element-type 'single-float))
(wsave (get-wsave-entry n))
(lensav (length wsave)))
(let ((ier (let ((ier
(nth-value 8 (nth-value 8
(rfft1f n 1 x n (gethash n *wsave-cache*) lensav work n 0)))) (rfft1f n 1 x n wsave lensav work n 0))))
(unless (zerop ier) (unless (zerop ier)
(error "rfft1f failed with code ~A" ier)) (error "rfft1f failed with code ~A" ier))
;; If X is the transformed value, then the output from rfftf is: ;; If X is the transformed value, then the output from rfftf is:
...@@ -72,21 +85,14 @@ ...@@ -72,21 +85,14 @@
(declare (type (simple-array (complex single-float) (*)) x)) (declare (type (simple-array (complex single-float) (*)) x))
;; Initialize the wave table if needed ;; Initialize the wave table if needed
(let* ((inv (convert-inverse-rfft x n)) (let* ((inv (convert-inverse-rfft x n))
(lensav (+ n (floor (log n 2)) 4)) (work (make-array n :element-type 'single-float))
(work (make-array n :element-type 'single-float))) (wsave (get-wsave-entry n))
(unless (gethash n *wsave-cache*) (lensav (length wsave)))
(let* ((wsave (make-array lensav :element-type 'single-float)))
(multiple-value-bind (ignore-0 ignore-1 ignore-2 ier)
(rfft1i n wsave lensav 0)
(declare (ignore ignore-0 ignore-1 ignore-2))
(unless (zerop ier)
(error "lensav is not big enough"))
(setf (gethash n *wsave-cache*) wsave))))
(let* ((ier (let* ((ier
(progn (progn
(format t "converted ~A~%" inv) (format t "converted ~A~%" inv)
(nth-value 8 (nth-value 8
(rfft1b n 1 inv n (gethash n *wsave-cache*) lensav work n 0))))) (rfft1b n 1 inv n wsave lensav work n 0)))))
(unless (zerop ier) (unless (zerop ier)
(error "rfft1b failed with code ~A" ier)) (error "rfft1b failed with code ~A" ier))
inv))) inv)))
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment