Commit 03f166cf authored by Raymond Toy's avatar Raymond Toy
Browse files

Add support for double-precision verrsion of FFTPACK5

This adds the asd file and the simple test program, based on the
single-precision version, but updated to compile and test the
double-precision version.
parent 289774ae
;;; -*- Mode: lisp -*-
;;;
;; Need f2cl to be loaded before we can even read this file.
(load-system "f2cl")
(in-package :f2cl-asdf)
;;; The files marked out with #+bad-array-type calls functions with
;;; REAL array arguments and the called function was expecting REAL*8
;;; arguments, or vice versa. This isn't supported by f2cl.
;;; See https://www2.cisl.ucar.edu/resources/legacy/fft5/documentation
;;; for documentation of this package.
(defsystem fftpack5-double
:description "F2CL conversion of FFTPACK5."
:class f2cl-system
:default-component-class :fortran-file
:f2cl-options (:package "FFTPACK5" :array-slicing t :array-type :array :common-as-array t :keep-lisp-file t :relaxed-array-decls t :promote-to-double t)
:components
((:module "package"
:pathname "fftpack5"
:components
((:cl-source-file "fftpack5-double")))
(:module "fftpack5"
:depends-on ("package")
:components
((:file "c1f2kb")
(:file "c1f2kf")
(:file "c1f3kb")
(:file "c1f3kf")
(:file "c1f4kb")
(:file "c1f4kf")
(:file "c1f5kb")
(:file "c1f5kf")
(:file "c1fgkb")
(:file "c1fgkf")
(:file "c1fm1b"
:f2cl-options (:copy-array-parameter t))
(:file "c1fm1f"
:f2cl-options (:copy-array-parameter t))
(:file "cfft1b"
:depends-on ("xerfft"))
(:file "cfft1f"
:depends-on ("xerfft" "mcfti1"))
(:file "cfft1i")
(:file "cfft2b"
:depends-on ("xerfft" "cfftmb"))
(:file "cfft2f"
:depends-on ("xerfft" "cfftmf"))
(:file "cfft2i"
:depends-on ("xerfft" "cfftmi"))
(:file "cfftmb"
:depends-on ("xerfft" "xercon" "cmfm1b"))
(:file "cfftmf"
:depends-on ("xerfft" "xercon" "cmfm1f"))
(:file "cfftmi"
:depends-on ("xerfft" "mcfti1"))
(:file "cmf2kb")
(:file "cmf2kf")
(:file "cmf3kb")
(:file "cmf3kf")
(:file "cmf4kb")
(:file "cmf4kf")
(:file "cmf5kb")
(:file "cmf5kf")
(:file "cmfgkb")
(:file "cmfgkf")
(:file "cmfm1b"
:f2cl-options (:copy-array-parameter t)
:depends-on ("cmf2kb" "cmf3kb" "cmf4kb" "cmf5kb" "cmfgkb"))
(:file "cmfm1f"
:f2cl-options (:copy-array-parameter t)
:depends-on ("cmf2kf" "cmf3kf" "cmf4kf" "cmf5kf" "cmfgkf"))
(:file "cosq1b"
:depends-on ("xerfft" "cosqb1"))
(:file "cosq1f"
:depends-on ("xerfft" "cosqf1"))
(:file "cosq1i"
:depends-on ("xerfft" "rfft1i"))
(:file "cosqb1"
:depends-on ("rfft1b" "xerfft"))
(:file "cosqf1"
:depends-on ("rfft1f" "xerfft"))
(:file "cosqmb"
:depends-on ("xerfft" "xercon" "mcsqb1"))
(:file "cosqmf"
:depends-on ("xerfft" "xercon" "mcsqf1"))
(:file "cosqmi"
:depends-on ("xerfft" "rfftmi"))
(:file "cost1b"
:depends-on ("xerfft" "costb1"))
(:file "cost1f"
:depends-on ("xerfft" "costf1"))
(:file "cost1i"
:depends-on ("xerfft" "rfft1i"))
(:file "costb1"
:depends-on ("rfft1f"))
(:file "costf1"
:depends-on ("rfft1f" "xerfft"))
#+bad-array-type
(:file "costmb"
:depends-on ("xerfft" "xercon" "mcstb1"))
#+bad-array-type
(:file "costmf"
:depends-on ("xerfft" "xercon" "mcstf1"))
(:file "costmi"
:depends-on ("xerfft" "rfftmi"))
(:file "factor")
(:file "mcfti1"
:depends-on ("tables" "factor"))
(:file "mcsqb1"
:depends-on ("rfftmb" "xerfft"))
(:file "mcsqf1"
:depends-on ("rfftmf" "xerfft"))
#+bad-array-type
(:file "mcstb1"
:depends-on ("rfftmf" "xerfft"))
#+bad-array-type
(:file "mcstf1"
:depends-on ("rfftmf" "xerfft"))
(:file "mradb2")
(:file "mradb3")
(:file "mradb4")
(:file "mradb5")
(:file "mradbg")
(:file "mradf2")
(:file "mradf3")
(:file "mradf4")
(:file "mradf5")
(:file "mradfg")
(:file "mrftb1"
:depends-on ("mradb4" "mradb2" "mradb3" "mradb5" "mradbg"))
(:file "mrftf1"
:depends-on ("mradf4" "mradf2" "mradf3" "mradf5" "mradfg"))
(:file "mrfti1")
(:file "msntb1"
:depends-on ("rfftmf" "xerfft"))
(:file "msntf1")
(:file "r1f2kb")
(:file "r1f2kf")
(:file "r1f3kb")
(:file "r1f3kf")
(:file "r1f4kb")
(:file "r1f4kf")
(:file "r1f5kb")
(:file "r1f5kf")
(:file "r1fgkb")
(:file "r1fgkf")
(:file "rfft1b"
:depends-on ("xerfft" "rfftb1"))
(:file "rfft1f"
:depends-on ("xerfft" "rfftf1"))
(:file "rfft1i"
:depends-on ("xerfft" "rffti1"))
(:file "rfft2b"
:f2cl-options (:copy-array-parameter t)
:depends-on ("xerfft" "rfftmb"))
(:file "rfft2f"
:f2cl-options (:copy-array-parameter t)
:depends-on ("xerfft" "rfftmf"))
(:file "rfft2i"
:depends-on ("xerfft" "rfftmi"))
(:file "rfftb1"
:depends-on ("r1f4kb" "r1f2kb" "r1f3kb" "r1fgkb" "r1f5kb"))
(:file "rfftf1"
:depends-on ("r1f4kf" "r1f2kf" "r1f3kf" "r1fgkf" "r1f5kf"))
(:file "rffti1")
(:file "rfftmb"
:depends-on ("xerfft" "xercon" "mrftb1"))
(:file "rfftmf"
:depends-on ("xerfft" "xercon" "mrftf1"))
(:file "rfftmi"
:depends-on ("xerfft" "mrfti1"))
(:file "sinq1b"
:depends-on ("xerfft"))
(:file "sinq1f"
:depends-on ("xerfft"))
(:file "sinq1i"
:depends-on ("xerfft"))
(:file "sinqmb"
:depends-on ("xerfft"))
(:file "sinqmf"
:depends-on ("xerfft" "xercon"))
(:file "sinqmi"
:depends-on ("xerfft"))
(:file "sint1b"
:depends-on ("xerfft" "sintb1"))
(:file "sint1f"
:depends-on ("xerfft" "sintf1"))
(:file "sint1i"
:depends-on ("xerfft"))
(:file "sintb1"
:depends-on ("xerfft"))
(:file "sintf1"
:depends-on ("xerfft"))
#+bad-array-type
(:file "sintmb"
:depends-on ("xerfft" "xercon"))
#+bad-array-type
(:file "sintmf"
:depends-on ("xerfft" "xercon"))
(:file "sintmi"
:depends-on ("xerfft"))
(:file "tables")
(:file "xercon")
(:file "xerfft")
))))
\ No newline at end of file
(defpackage :fftpack5
(:use :common-lisp)
(:export "RFFT"
"INVERSE-RFFT"))
(in-package :fftpack5)
(defvar *wsave-cache*
(make-hash-table)
"Cache for different wsave tables. The key is the FFT size; the
value is a the wsave table needed by the FFT routines.")
(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
cache, return it. Otherwise compute a new value and save it in the
cache."
(let ((entry (gethash n *wsave-cache*)))
(if entry
entry
(let* ((lensav (+ n 4 (floor (log n 2))))
(wsave (make-array lensav :element-type 'double-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)
;; 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 double-float (*)) x))
(let* ((n (length x))
(nhalf (floor (/ n 2)))
(out (make-array (+ 1 nhalf) :element-type '(complex double-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 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 double-float) (*)) x))
(let ((res (make-array n :element-type 'double-float)))
(setf (aref res 0) (realpart (aref x 0)))
(loop for j from 1 below (if (evenp n) (1- (length x)) (length x))
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)
(setf (aref res (1- n)) (realpart (aref x (1- (length x))))))
res))
(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), k, 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 double-float (*)) x))
(let* ((n (length x))
(work (make-array n :element-type 'double-float))
(wsave (get-wsave-entry n))
(lensav (length wsave)))
(let ((ier
(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))
(convert-rfft x))))
(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 double-float) (*)) x))
;; Check to see if N is consistent with the length of X. However,
;; this isn't foolproof. For a real FFT of length 5, the FFT array
;; has length 3. But a real FFT of length 4 also has an array
;; length of 3.
(unless (= (length x) (1+ (floor n 2)))
(error "Length of X (~A) is not compatible with N (~D)"
(length x) n))
(let* ((inv (convert-inverse-rfft x n))
(work (make-array n :element-type 'double-float))
(wsave (get-wsave-entry n))
(lensav (length wsave)))
(let* ((ier
(progn
(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)))
(defun test-rfft (n &key verbose)
"Test the rfft and inverse-rfft routines using a simple ramp of
length n"
;; For testing, use a simple ramp: 1, 2, 3, ..., N.
;; With some care, we can derive the FFT for this is
;;
;; X[0] = (N + 1) / 2
;; X[n] = -1/2 - cot(%pi*n/N)/2, n = 1, N
(let* ((x (make-array n :element-type 'double-float)))
(loop for k from 0 below n
do
(setf (aref x k) (float (+ k 1) 1d0)))
(let* ((xfrm (rfft x))
(xfrm-len (length xfrm))
(expected (make-array xfrm-len :element-type '(complex double-float)))
(noise-pwr 0d0)
(signal-pwr 0d0))
(declare (double-float noise-pwr signal-pwr))
(setf (aref expected 0) (complex (* 0.5 (+ n 1)) 0.0))
(loop for k from 1 below xfrm-len
with omega = (coerce (/ pi n) 'double-float)
do
(setf (aref expected k)
(complex -0.5 (/ -0.5 (tan (* omega k))))))
(when verbose
(format t "Forward transform; actual vs expected~%")
(loop for k from 0 below xfrm-len
do
(format t "~4d: ~A ~A~%" k (aref xfrm k) (aref expected k))))
(incf noise-pwr (expt (abs (- (aref xfrm 0) (* 0.5 (+ n 1)))) 2))
(incf signal-pwr (expt (* 0.5 (+ n 1)) 2))
(loop for k from 1 below (length xfrm)
do
(incf noise-pwr (expt (abs (- (aref xfrm k)
(aref expected k)))
2))
(incf signal-pwr (expt (abs (aref expected k)) 2)))
(let ((inv (inverse-rfft xfrm n))
(inv-noise-pwr 0.0)
(inv-signal-pwr 0.0))
(when verbose
(format t "Inverse transform; actual vs expected~%")
(loop for k from 0 below n
do
(format t "~4d: ~A ~A~%" k (aref inv k) (+ k 1))))
(loop for k from 0 below n
do
(incf inv-noise-pwr (expt (- (aref inv k) (+ k 1)) 2))
(incf inv-signal-pwr (expt (coerce (+ k 1) 'double-float) 2)))
(flet ((db (s n)
(if (zerop n)
1000d0
(* 10 (log (/ s n) 10)))))
(values (db signal-pwr noise-pwr) (db inv-signal-pwr inv-noise-pwr) noise-pwr signal-pwr inv-noise-pwr inv-signal-pwr))))))
\ 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