Commit eaf7f881 authored by Raymond Toy's avatar Raymond Toy
Browse files

Clean up, adding comments and docstrings

parent 4f8b8372
......@@ -8,24 +8,6 @@
"Cache for different wsave tables. The key is the FFT size; the
value is a the wsave table needed by the FFT routines.")
(defun convert-rfft (x)
"Convert the output of FFTPACK RFFT1F (forward real FFT) into a more
user-friendly format with complex values"
(declare (type
(simple-array single-float (*)) x)) (let* ((n (length x))
(nhalf (floor (/ n 2)))
(out (make-array (+ 1 nhalf) :element-type '(complex single-float))))
(setf (aref out 0) (complex (aref x 0) 0.0))
(loop for j from 1 to (if (evenp n) (1- nhalf) nhalf)
for k from 1 by 2
do
(setf (aref out j) (complex (* 0.5 (aref x k))
(* 0.5 (aref x (1+ k))))))
(when (evenp n)
(setf (aref out nhalf)
(complex (aref x (1- n)) 0.0)))
out))
(defun get-wsave-entry (n)
"Get the wsave array and it's length that is needed for computing
the (forward and inverse) FFTs. The value is cached, so if it's the
......@@ -40,13 +22,54 @@
(rfft1i n wsave lensav 0)
(declare (ignore ignore-0 ignore-1 ignore-2))
(unless (zerop ier)
;; This shouldn't really ever happen.
(error "lensav is not big enough"))
(setf (gethash n *wsave-cache*) wsave)
wsave)))))
(defun convert-rfft (x)
"Convert the output of FFTPACK RFFT1F (forward real FFT) into a more
user-friendly format with complex values"
(declare (type (simple-array single-float (*)) x))
(let* ((n (length x))
(nhalf (floor (/ n 2)))
(out (make-array (+ 1 nhalf) :element-type '(complex single-float))))
;; If X is the transformed value, then the output from rfftf is:
;; 0: X(0) / N
;; 1: 2*realpart(X(1))/N
;; 2: 2*imagpart(X(1))/N
;; 3: 2*realpart(X(2))/N
;; 4: 2*realpart(X(2))/N
;; ...
;; N-1: X(N-1)/N
;; The last term exists only if N is even.
(setf (aref out 0) (complex (aref x 0) 0.0))
(loop for j from 1 to (if (evenp n) (1- nhalf) nhalf)
for k from 1 by 2
do
;; Need to remove the factor of 2 that rfftf added.
(setf (aref out j) (complex (* 0.5 (aref x k))
(* 0.5 (aref x (1+ k))))))
(when (evenp n)
(setf (aref out nhalf)
(complex (aref x (1- n)) 0.0)))
out))
(defun rfft (x)
"Compute the real FFT of X.
Let N be the length of X. The FFT is:
Y[n] = 1/N*sum(x[k] * exp(2*%i*%pi*k*n/N), n, 0, N-1)
for n = 0, 1,...,floor(N/2)+1
NOTE: This differs from the typical engineering definition of the
forward FFT, where the transfrom is often not normalized by the
length and the argument to exp has a negative sign. This is
generally considered the inverse FFT in engineering"
(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))
......@@ -55,19 +78,16 @@
(nth-value 8
(rfft1f n 1 x n wsave lensav work n 0))))
(unless (zerop ier)
;; This should never happen because we always allocate enough
;; space for x, wsave, and the work array.
(error "rfft1f failed with code ~A" ier))
;; If X is the transformed value, then the output from rfftf is:
;; 0: X(0) / N
;; 1: 2*realpart(X(1))/N
;; 2: 2*imagpart(X(1))/N
;; 3: 2*realpart(X(2))/N
;; 4: 2*realpart(X(2))/N
;; ...
;; N-1: X(N-1)/N
;; The last term exists only if N is even.
(convert-rfft x))))
(defun convert-inverse-rfft (x n)
"Convert the complex-valued input, X, (that was produced by rfft)
into the form needed by rfft1b. The length of the transform, N, is
needed because it cannot be uniquely determined from the length of
X."
(declare (type (simple-array (complex single-float) (*)) x))
(let ((res (make-array n :element-type 'single-float)))
(setf (aref res 0) (realpart (aref x 0)))
......@@ -75,6 +95,8 @@
for k from 1 by 2
do
(let ((z (aref x j)))
;; Put back the factor of 2 that we removed in rfft but
;; is needed by rfft1b.
(setf (aref res k) (* 2 (realpart z)))
(setf (aref res (1+ k)) (* 2 (imagpart z)))))
(when (evenp n)
......@@ -82,17 +104,21 @@
res))
(defun inverse-rfft (x n)
"Compute the inverse real FFT of X. N is the length of the transform
(and not the length of X!).
See RFFT for the definition of the FFT used in these routines."
(declare (type (simple-array (complex single-float) (*)) x))
;; Initialize the wave table if needed
(let* ((inv (convert-inverse-rfft x n))
(work (make-array n :element-type 'single-float))
(wsave (get-wsave-entry n))
(lensav (length wsave)))
(let* ((ier
(progn
(format t "converted ~A~%" inv)
(nth-value 8
(rfft1b n 1 inv n wsave lensav work n 0)))))
(unless (zerop ier)
;; This should never happen because we should have always
;; allocated the correct size of the arrays.
(error "rfft1b failed with code ~A" ier))
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