Commit 284a108b authored by Raymond Toy's avatar Raymond Toy
Browse files

Add simpler interface and define package

Put everything in the FFTPACK51 package.

Define a simpler interface for computing the real FFT.  And provide a
function for converting the output of rfft1f to a more obvious complex
result.  This also removes some of the inconsistent, non-obvious,
scaling done in rfft1f.
parent 8198c7f7
......@@ -16,16 +16,15 @@
:description "F2CL conversion of FFTPACK5."
:class f2cl-system
:default-component-class :fortran-file
:f2cl-options (:array-slicing t :array-type :array :common-as-array t :keep-lisp-file t :relaxed-array-decls t)
:f2cl-options (:package "FFTPACK5" :array-slicing t :array-type :array :common-as-array t :keep-lisp-file t :relaxed-array-decls t)
:components
(#+nil
(:module "package"
:pathname "fftpack5"
((:module "package"
:pathname "fftpack51"
:components
((:cl-source-file "fftpack5")))
((:cl-source-file "fftpack51")))
(:module "fftpack51"
:pathname "fftpack51/src"
;;:depends-on ("package")
:depends-on ("package")
:components
(
;;(:file "c1f2kb")
......
(defpackage :fftpack51
(:use :common-lisp))
(in-package :fftpack51)
(defvar *wsave-cache*
(make-hash-table) "Hash table containing cached *wsave* tables used
by the FFT routines. The key is the order of the FFT.")
(defun rfft (x)
(declare (type (simple-array single-float (*)) x))
;; Initialize the wave table if needed
(let* ((n (length x))
(lensav (+ n (floor (log n 2)) 4))
(work (make-array n :element-type 'single-float)))
(unless (gethash n *wsave-cache*)
(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
(nth-value 8
(rfft1f n 1 x n (gethash n *wsave-cache*) lensav work n 0))))
(unless (zerop ier)
(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.
x)))
(defun inverse-rfft (x)
(declare (type (simple-array single-float (*)) x))
;; Initialize the wave table if needed
(let* ((n (length x))
(lensav (+ n (floor (log n 2)) 4))
(work (make-array n :element-type 'single-float)))
(unless (gethash n *wsave-cache*)
(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
(nth-value 8
(rfft1b n 1 x n (gethash n *wsave-cache*) lensav work n 0))))
(unless (zerop ier)
(error "rfft1b failed with code ~A" ier))
x)))
(defun convert-rfft-to-cfft (x)
(declare (type (simple-array single-float (*)) x))
(let* ((n (length x))
(nhalf (floor (/ n 2)))
(out (make-array n :element-type '(complex single-float))))
(setf (aref out 0) (complex (* n (aref x 0)) 0.0))
(loop for j from 1 to (1- nhalf)
do
(setf (aref out j) (complex (* n (aref x (1- (* 2 j))) 0.5)
(* n (aref x (* 2 j)) 0.5)))
(setf (aref out (- n j)) (conjugate (aref out j))))
(when (evenp n)
(setf (aref out nhalf)
(complex (* n (aref x (1- n))) 0.0)))
out))
(defun convert-rfft (x)
(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 (1- nhalf)
do
(setf (aref out j) (complex (* 0.5 (aref x (1- (* 2 j))))
(* 0.5 (aref x (* 2 j))))))
(when (evenp n)
(setf (aref out nhalf)
(complex (aref x (1- n)) 0.0)))
out))
\ No newline at end of file
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