diff --git a/fast-fourier-transforms/fft-half-complex.lisp b/fast-fourier-transforms/fft-half-complex.lisp index 66a12ba4e97322eb9077c8dbc8857da531ad8ff3..94436bc04be97d73268f0467d38cf6f8c2215530 100644 --- a/fast-fourier-transforms/fft-half-complex.lisp +++ b/fast-fourier-transforms/fft-half-complex.lisp @@ -1,6 +1,6 @@ ;; Functions for fast fourier transforms on real data. ;; Sumant Oemrawsingh, Sat Oct 31 2009 - 20:12 -;; Time-stamp: <2009-11-01 23:04:46EST fft-half-complex.lisp> +;; Time-stamp: <2009-11-03 23:02:17EST fft-half-complex.lisp> ;; /usr/include/gsl/gsl_fft_halfcomplex.h ;; /usr/include/gsl/gsl_fft_halfcomplex_float.h @@ -33,22 +33,6 @@ :documentation "Inverse FFT for a half-complex radix-2 vector") -(defmfun fft-half-complex-radix2-unpack - ((h-complex vector) - &key (stride 1) (n (expt 2 (floor (log (size h-complex) 2)))) - (complex-output (eltcase single-float (make-marray '(complex single-float) :dimensions n) - t (make-marray '(complex double-float) :dimensions n)))) - ("gsl_fft_halfcomplex" :type "_radix2_unpack") - (((c-pointer h-complex) :pointer) ((c-pointer complex-output) :pointer) (stride sizet) (n sizet)) - :definition :generic - :element-types :float - :inputs (h-complex complex-output) - :outputs (complex-output) - :return (complex-output) - :documentation - "Convert an array of half-complex coefficients as returned by - real-fft-radix2-transform, into an ordinary complex array.") - ;; Mixed Radix general-N functions (defmobject fft-half-complex-wavetable-double-float @@ -131,20 +115,3 @@ :documentation "Forward FFT for a half-complex vector") -(defmfun fft-half-complex-unpack - ((h-complex vector) - &key (stride 1) (n (size h-complex)) - (complex-output (eltcase single-float (make-marray '(complex single-float) :dimensions n) - t (make-marray '(complex double-float) :dimensions n)))) - ("gsl_fft_halfcomplex" :type "_unpack") - (((c-pointer h-complex) :pointer) ((c-pointer complex-output) :pointer) (stride sizet) (n sizet)) - :definition :generic - :element-types :float - :inputs (h-complex complex-output) - :outputs (complex-output) - :return (complex-output) - :documentation - "This function converts an array of half-complex - coefficients as returned by fft-real-transform, into an ordinary complex - array. It fills in the complex array using the symmetry z_k = z_{n-k}^* to - reconstruct the redundant elements.") diff --git a/fast-fourier-transforms/fft-real.lisp b/fast-fourier-transforms/fft-real.lisp index 33875e150ad0b51e725dd0fb2a2de78d3daf27fa..c84f286007d0391fe86a65dab9ce93a04bed1f0e 100644 --- a/fast-fourier-transforms/fft-real.lisp +++ b/fast-fourier-transforms/fft-real.lisp @@ -1,6 +1,6 @@ ;; Functions for fast fourier transforms on real data. ;; Sumant Oemrawsingh, Sun Oct 25 2009 - 16:35 -;; Time-stamp: <2009-11-01 22:41:14EST fft-real.lisp> +;; Time-stamp: <2009-11-03 23:02:17EST fft-real.lisp> ;; /usr/include/gsl/gsl_fft_real.h ;; /usr/include/gsl/gsl_fft_real_float.h @@ -65,21 +65,3 @@ :documentation "This function allocates a workspace for a real float transform of length n.") - -(defmfun fft-real-unpack - ((vector vector) - &key (stride 1) (n (size vector)) - (complex-output - (eltcase single-float (make-marray '(complex single-float) :dimensions n) - t (make-marray '(complex double-float) :dimensions n)))) - ("gsl_fft_real" :type "_unpack") - (((c-pointer vector) :pointer) ((c-pointer complex-output) :pointer) (stride sizet) (n sizet)) - :definition :generic - :element-types :float - :inputs (vector complex-output) - :outputs (complex-output) - :return (complex-output) - :documentation - "This function converts a single real array into an equivalent complex - array (with imaginary part set to zero), suitable for fft-complex - routines.") diff --git a/fast-fourier-transforms/unpack.lisp b/fast-fourier-transforms/unpack.lisp new file mode 100644 index 0000000000000000000000000000000000000000..7b58662de83f554fbc9184d78bcd797fc685cf84 --- /dev/null +++ b/fast-fourier-transforms/unpack.lisp @@ -0,0 +1,93 @@ +;; Unpack functions for FFT vectors. +;; Sumant Oemrawsingh, Sun Oct 25 2009 - 16:35 +;; Time-stamp: <2009-11-03 23:20:35EST unpack.lisp> + +(in-package :gsl) + +;; /usr/include/gsl/gsl_fft_real.h +;; /usr/include/gsl/gsl_fft_real_float.h +;; /usr/include/gsl/gsl_fft_complex.h +;; /usr/include/gsl/gsl_fft_complex_float.h + +(defmfun fft-real-unpack + ((vector vector) + &key (stride 1) + (output + (eltcase single-float + (make-marray '(complex single-float) :dimensions (size vector)) + t + (make-marray '(complex double-float) :dimensions (size vector))))) + ("gsl_fft_real" :type "_unpack") + (((c-pointer vector) :pointer) + ((c-pointer output) :pointer) + (stride sizet) ((size vector) sizet)) + :definition :generic + :element-types :float + :inputs (vector output) + :outputs (output) + :return (output) + :export nil + :index unpack + :documentation + "This function converts a single real array into an equivalent complex + array (with imaginary part set to zero), suitable for fft-complex + routines.") + +(defmfun fft-half-complex-radix2-unpack + ((vector vector) + &key (stride 1) + (output + (eltcase single-float + (make-marray '(complex single-float) :dimensions (size vector)) + t + (make-marray '(complex double-float) :dimensions (size vector))))) + ("gsl_fft_halfcomplex" :type "_radix2_unpack") + (((c-pointer vector) :pointer) + ((c-pointer output) :pointer) + (stride sizet) ((size vector) sizet)) + :definition :generic + :element-types :float + :inputs (vector output) + :outputs (output) + :return (output) + :export nil + :index unpack + :documentation + "Convert an array of half-complex coefficients as returned by + real-fft-radix2-transform, into an ordinary complex array.") + +(defmfun fft-half-complex-unpack + ((vector vector) + &key (stride 1) + (output + (eltcase single-float + (make-marray '(complex single-float) :dimensions (size vector)) + t + (make-marray '(complex double-float) :dimensions (size vector))))) + ("gsl_fft_halfcomplex" :type "_unpack") + (((c-pointer vector) :pointer) + ((c-pointer output) :pointer) + (stride sizet) ((size vector) sizet)) + :definition :generic + :element-types :float + :inputs (vector output) + :outputs (output) + :return (output) + :export nil + :index unpack + :documentation + "This function converts an array of half-complex coefficients as + returned by fft-real-transform, into an ordinary complex array. It + fills in the complex array using the symmetry z_k = z_{n-k}^* to + reconstruct the redundant elements.") + +(export 'unpack) +(defun unpack + (vector &rest args &key (stride 1) (unpack-type 'real) &allow-other-keys) + (let ((pass-on-args (copy-list args))) + (remf pass-on-args :unpack-type) + (if (eq unpack-type 'complex) + (if (power-of-2-p (floor (size vector) stride)) + (apply 'fft-half-complex-radix2-unpack pass-on-args) + (apply 'fft-half-complex-unpack pass-on-args)) + (apply 'fft-real-unpack pass-on-args)))) diff --git a/gsll.asd b/gsll.asd index 30a7b8db408f68bfd0c5670372b464d5697ca53d..8c87182549b44cbaf9b0fc932a6e3344bc605e34 100644 --- a/gsll.asd +++ b/gsll.asd @@ -1,6 +1,6 @@ ;; Definition of GSLL system ;; Liam Healy -;; Time-stamp: <2009-11-01 23:01:16EST gsll.asd> +;; Time-stamp: <2009-11-03 23:23:15EST gsll.asd> (when (asdf:find-system :fsbv nil) (pushnew :fsbv *features*)) @@ -140,6 +140,7 @@ (:file "fft-real") (:file "fft-half-complex" :depends-on ("fft-real")) (:file "forward") + (:file "unpack") (:file "fft-interface" :depends-on ("fft-complex" "fft-half-complex" "fft-real")) (:file "fft-interface-example"))) (:module random