From 4c0778836b4800966653f4b2e8aeb761fcb5cd93 Mon Sep 17 00:00:00 2001 From: Kevin Rosenberg Date: Sat, 2 May 2009 10:20:50 -0600 Subject: [PATCH] Initial import --- c-tests/Makefile | 8 + c-tests/sizes.c | 7 + cl-fftw3-tests.asd | 26 +++ cl-fftw3.asd | 39 +++++ common.lisp | 128 ++++++++++++++ ffi.lisp | 253 ++++++++++++++++++++++++++++ package.lisp | 59 +++++++ specials.lisp | 47 ++++++ tests.lisp | 231 ++++++++++++++++++++++++++ transform.lisp | 403 +++++++++++++++++++++++++++++++++++++++++++++ wisdom.lisp | 35 ++++ 11 files changed, 1236 insertions(+) create mode 100644 c-tests/Makefile create mode 100644 c-tests/sizes.c create mode 100644 cl-fftw3-tests.asd create mode 100644 cl-fftw3.asd create mode 100644 common.lisp create mode 100644 ffi.lisp create mode 100644 package.lisp create mode 100644 specials.lisp create mode 100644 tests.lisp create mode 100644 transform.lisp create mode 100644 wisdom.lisp diff --git a/c-tests/Makefile b/c-tests/Makefile new file mode 100644 index 0000000..38aac38 --- /dev/null +++ b/c-tests/Makefile @@ -0,0 +1,8 @@ +all: sizes + +sizes: sizes.c + gcc -I/opt/local/include sizes.c -L/opt/local/lib -lfftw3 -lm -o sizes + +clean: + rm -f *.o *~ sizes + diff --git a/c-tests/sizes.c b/c-tests/sizes.c new file mode 100644 index 0000000..134304a --- /dev/null +++ b/c-tests/sizes.c @@ -0,0 +1,7 @@ +#include +#include + +main () { + printf("Size of fftw_complex: %d\n", sizeof (fftw_complex)); + printf("Size of fftw_plan: %d\n", sizeof (fftw_plan)); +} diff --git a/cl-fftw3-tests.asd b/cl-fftw3-tests.asd new file mode 100644 index 0000000..1f2f59f --- /dev/null +++ b/cl-fftw3-tests.asd @@ -0,0 +1,26 @@ +;;;; -*- Mode: Lisp; Syntax: ANSI-Common-Lisp; Base: 10 -*- +;;;; ************************************************************************* +;;;; FILE IDENTIFICATION +;;;; +;;;; Name: fftw3-tests.asd +;;;; Purpose: ASDF system definitionf for fftw3 testing package +;;;; Author: Kevin M. Rosenberg +;;;; Date Started: Jan 2009 +;;;; +;;;; $Id$ +;;;; ************************************************************************* + +(defpackage #:fftw3-tests-system + (:use #:asdf #:cl)) +(in-package #:fftw3-tests-system) + +(defsystem cl-fftw3-tests + :depends-on (rt cl-fftw3) + :components + ((:file "tests"))) + +(defmethod perform ((o test-op) (c (eql (find-system 'cl-fftw3-tests)))) + (or (funcall (intern (symbol-name '#:do-tests) + (find-package '#:regression-test))) + (error "test-op failed"))) + diff --git a/cl-fftw3.asd b/cl-fftw3.asd new file mode 100644 index 0000000..846d2bd --- /dev/null +++ b/cl-fftw3.asd @@ -0,0 +1,39 @@ +;;;; -*- Mode: Lisp; Syntax: ANSI-Common-Lisp; Base: 10 -*- +;;;; ************************************************************************* +;;;; FILE IDENTIFICATION +;;;; +;;;; Name: cl-fftw3.asd +;;;; Purpose: ASDF system definition for FFTW3 package +;;;; Programmer: Kevin M. Rosenberg +;;;; Date Started: Jan 2009 +;;;; +;;;; $Id$ +;;;; +;;;; This file, part of FFTW3, is Copyright (c) 2009 by Kevin M. Rosenberg +;;;; +;;;; FFTW3 users are granted the rights to distribute and use this software +;;;; as governed by the terms of the Lisp Lesser GNU Public License +;;;; (http://opensource.franz.com/preamble.html), also known as the LLGPL. +;;;; ************************************************************************* + +(in-package #:cl-user) +(defpackage #:fftw3-system (:use #:asdf #:cl)) +(in-package #:fftw3-system) + +(defsystem cl-fftw3 + :name "fftw3" + :author "Kevin M. Rosenberg " + :maintainer "Kevin M. Rosenberg " + :licence "LLGPL" + :depends-on (:kmrcl :cffi) + :components + ((:file "package") + (:file "specials" :depends-on ("package")) + (:file "common" :depends-on ("specials")) + (:file "ffi" :depends-on ("common")) + (:file "transform" :depends-on ("ffi")) + (:file "wisdom" :depends-on ("ffi")))) + +(defmethod perform ((o test-op) (c (eql (find-system 'cl-fftw3)))) + (operate 'load-op 'cl-fftw3-tests) + (operate 'test-op 'cl-fftw3-tests :force t)) diff --git a/common.lisp b/common.lisp new file mode 100644 index 0000000..a82b9de --- /dev/null +++ b/common.lisp @@ -0,0 +1,128 @@ +;;;; -*- Mode: Lisp; Syntax: ANSI-Common-Lisp; Base: 10 -*- +;;;; ************************************************************************* +;;;; FILE IDENTIFICATION +;;;; +;;;; Name: common.lisp +;;;; Purpose: Common functions for FFTW3 package +;;;; Programmer: Kevin M. Rosenberg +;;;; Date Started: Jan 2009 +;;;; +;;;; $Id$ +;;;; +;;;; This file, part of FFTW3, is Copyright (c) 2009 by Kevin M. Rosenberg +;;;; +;;;; FFTW3 users are granted the rights to distribute and use this software +;;;; as governed by the terms of the Lisp Lesser GNU Public License +;;;; (http://opensource.franz.com/preamble.html), also known as the LLGPL. +;;;; ************************************************************************* + +(in-package #:fftw3) + +(defclass fftw-multi () + ((plan :initarg :plan :initform nil :accessor plan))) + +(defclass fftw-r2r-1d-multi (fftw-multi) + ((in-count :initarg :in-count :initform 0 :type fixnum :accessor in-count) + (out-c :initarg :out-c :initform nil :accessor out-c) + (in-cf :initarg :in-cf :initform nil :accessor in-cf) + (out-cf :initarg :out-cf :initform nil :accessor out-cf))) + +(defclass fftw-r2c-1d-multi (fftw-multi) + ((in-count :initarg :in-count :initform 0 :type fixnum :accessor in-count) + (out-count :initarg :out-count :initform 0 :type fixnum :accessor out-count) + (complex-output :initarg :complex-output :initform t :type boolean :accessor complex-output) + (in-cf :initarg :in-cf :initform nil :accessor in-cf) + (out-cf :initarg :out-cf :initform nil :accessor out-cf) + (out :initarg :out :initform nil :accessor out))) + +(defclass fftw-1d-multi (fftw-multi) + ((in-count :initarg :in-count :initform 0 :type fixnum :accessor in-count) + (out-count :initarg :out-count :initform 0 :type fixnum :accessor out-count) + (in-cf :initarg :in-cf :initform nil :accessor in-cf) + (out-cf :initarg :out-cf :initform nil :accessor out-cf) + (out :initarg :out :initform nil :accessor out))) + +(defclass fftw-c2r-1d-multi (fftw-multi) + ((plan-even :initarg :plan-even :initform nil :accessor plan-even) + (plan-odd :initarg :plan-odd :initform nil :accessor plan-odd) + (in-count :initarg :in-count :initform 0 :type fixnum :accessor in-count) + (out-count-odd :initarg :out-count-odd :initform 0 :type fixnum :accessor out-count-odd) + (out-count-even :initarg :out-count-even :initform 0 :type fixnum :accessor out-count-even) + (in-cf :initarg :in-cf :initform nil :accessor in-cf) + (out-cf :initarg :out-cf :initform nil :accessor out-cf) + (out-odd :initarg :out-odd :initform nil :accessor out-odd) + (out-even :initarg :out-even :initform nil :accessor out-even))) + +(defun normalize-input-range (in start count) + "Returns the start and count of an input range with range checking." + (declare #.*standard-optimize-settings*) + (unless (vectorp in) + (error "Input needs to be a vector.")) + (when (or (not (integerp start)) (minusp start)) + (error "Start must be an integer, got ~A." start)) + (let ((len (length in))) + (declare (fixnum len)) + (when (>= start len) + (error "Start (~D) greater than last element (~D)." start (1- len))) + + (cond + ((null count) + (setq count (- len start))) + ((or (not (integerp count)) (minusp count)) + (error "Count must be positive integer or nil, got ~A." count)) + ((> (+ start count) len) + (error "Requesting element (~D) past end of vector (~D)." (+ start count) len))) + (values start count))) + +(defun complex-is-even-p (v) + "Returns T if first and last members have zero imaginary component." + ;; Do not use *standard-optimize-settings-here* -- they fail on lispworks 5.1.2 + (declare (optimize (speed 3) (space 0))) + (and (zerop (imagpart (aref v 0))) + (zerop (imagpart (aref v (1- (length v))))))) + +(defun complex-to-mag-phase (cmplx &optional (elements-complex t)) + (declare #.*standard-optimize-settings*) + (let* ((len (length cmplx)) + (out-len (if elements-complex len (/ len 2))) + (mag (make-array out-len :element-type 'double-float)) + (phase (make-array out-len :element-type 'double-float))) + (declare (fixnum len out-len)) + (cond + (elements-complex + (dotimes (i out-len) + (declare (fixnum i)) + (let* ((e (aref cmplx i))) + (setf (aref mag i) (abs e)) + (setf (aref phase i) (phase e))))) + (t + (dotimes (i out-len) + (declare (fixnum i)) + (let* ((base (+ i i)) + (r (aref cmplx base)) + (i (aref cmplx (1+ base)))) + (setf (aref mag i) (sqrt (+ (* r r) (* i i)))) + (setf (aref phase i) (atan i r)))))) + (values mag phase))) + +(defun hc-to-mag-phase (hc) + "Turns half-complex into magniture and phase vectors." + (declare #.*standard-optimize-settings*) + (let* ((n (length hc)) + (out-n (1+ (floor n 2))) + (mag (make-array out-n :element-type 'double-float)) + (phase (make-array out-n :element-type 'double-float))) + (declare (fixnum n out-n)) + (dotimes (i out-n) + (cond + ((or (zerop i) + (and (evenp n) (eql i (1- out-n)))) + (setf (aref mag i) (aref hc (/ i 2))) + (setf (aref phase i) 0d0)) + (t + (let ((re (aref hc i)) + (im (aref hc (- n i)))) + (setf (aref mag i) (sqrt (+ (* re re) (* im im)))) + (setf (aref phase i) (atan im re)))))) + (values mag phase))) + diff --git a/ffi.lisp b/ffi.lisp new file mode 100644 index 0000000..8d16045 --- /dev/null +++ b/ffi.lisp @@ -0,0 +1,253 @@ +;;;; -*- Mode: Lisp; Syntax: ANSI-Common-Lisp; Base: 10 -*- +;;;; ************************************************************************* +;;;; FILE IDENTIFICATION +;;;; +;;;; Name: ffi.lisp +;;;; Purpose: CFFI interface for FFTW3 package +;;;; Programmer: Kevin M. Rosenberg +;;;; Date Started: Mar 2009 +;;;; +;;;; $Id$ +;;;; +;;;; This file, part of FFTW3, is Copyright (c) 2009 by Kevin M. Rosenberg +;;;; +;;;; FFTW3 users are granted the rights to distribute and use this software +;;;; as governed by the terms of the Lisp Lesser GNU Public License +;;;; (http://opensource.franz.com/preamble.html), also known as the LLGPL. +;;;; ************************************************************************* + +(in-package #:fftw3) + +(cffi:define-foreign-library fftw3 + ((or :darwin :macosx) (:or #P"/opt/local/lib/libfftw3.dylib" + #P"/usr/lib/libfftw3.dylib")) + (:linux (:or #P"/usr/lib/libfftw3.so" + #P"/usr/local/lib/libfftw3.so")) + (t (:default "libfftw3"))) + +(cffi:use-foreign-library fftw3) + +(cffi:defcstruct fftw-complex-struct + "Complex number structure." + (re :double) + (im :double)) + +(cffi:defctype fftw-plan :pointer) + +(declaim (inline fftw-plan-dft-1d)) +(cffi:defcfun ("fftw_plan_dft_1d" fftw-plan-dft-1d) fftw-plan + (n :int) + (in (:pointer fftw-complex-struct)) + (out (:pointer fftw-complex-struct)) + (sign :int) + (flags :uint)) + +(declaim (inline fftw-plan-dft-2d)) +(cffi:defcfun ("fftw_plan_dft_2d" fftw-plan-dft-2d) fftw-plan + (n0 :int) + (n1 :int) + (in (:pointer fftw-complex-struct)) + (out (:pointer fftw-complex-struct)) + (sign :int) + (flags :uint)) + +(declaim (inline fftw-plan-dft-3d)) +(cffi:defcfun ("fftw_plan_dft_3d" fftw-plan-dft-3d) fftw-plan + (n0 :int) + (n1 :int) + (n2 :int) + (in (:pointer fftw-complex-struct)) + (out (:pointer fftw-complex-struct)) + (sign :int) + (flags :uint)) + +(declaim (inline fftw-plan-dft)) +(cffi:defcfun ("fftw_plan_dft" fftw-plan-dft) fftw-plan + (rank :int) + (n (:pointer :int)) + (in (:pointer fftw-complex-struct)) + (out (:pointer fftw-complex-struct)) + (sign :int) + (flags :uint)) + +(declaim (inline fftw-plan-r2r-1d)) +(cffi:defcfun ("fftw_plan_r2r_1d" fftw-plan-r2r-1d) fftw-plan + (n :int) + (in :pointer) + (out :pointer) + (kind :int) + (flags :uint)) + +(declaim (inline fftw-plan-r2r-2d)) +(cffi:defcfun ("fftw_plan_r2r_2d" fftw-plan-r2r-2d) fftw-plan + (n0 :int) + (n1 :int) + (in :pointer) + (out :pointer) + (kind :int) + (flags :uint)) + +(declaim (inline fftw-plan-r2r-3d)) +(cffi:defcfun ("fftw_plan_r2r_3d" fftw-plan-r2r-3d) fftw-plan + (n0 :int) + (n1 :int) + (n2 :int) + (in :pointer) + (out :pointer) + (kind :int) + (flags :uint)) + +(declaim (inline fftw-plan-dft-r2c-1d)) +(cffi:defcfun ("fftw_plan_dft_r2c_1d" fftw-plan-dft-r2c-1d) fftw-plan + (n :int) + (in :pointer) + (out :pointer) + (flags :uint)) + +(declaim (inline fftw-plan-dft-r2c-2d)) +(cffi:defcfun ("fftw_plan_dft_r2c_2d" fftw-plan-dft-r2c-2d) fftw-plan + (n0 :int) + (n1 :int) + (in :pointer) + (out :pointer) + (flags :uint)) + +(declaim (inline fftw-plan-dft-r2c-3d)) +(cffi:defcfun ("fftw_plan_dft_r2c_3d" fftw-plan-dft-r2c-3d) fftw-plan + (n0 :int) + (n1 :int) + (n2 :int) + (in :pointer) + (out :pointer) + (flags :uint)) + +(declaim (inline fftw-plan-dft-r2c)) +(cffi:defcfun ("fftw_plan_dft_r2c" fftw-plan-dft-r2c) fftw-plan + (rank :int) + (n (:pointer :int)) + (in :pointer) + (out :pointer) + (flags :uint)) + +(declaim (inline fftw-plan-dft-c2r-1d)) +(cffi:defcfun ("fftw_plan_dft_c2r_1d" fftw-plan-dft-c2r-1d) fftw-plan + (n :int) + (in :pointer) + (out :pointer) + (flags :uint)) + +(declaim (inline fftw-plan-dft-c2r-2d)) +(cffi:defcfun ("fftw_plan_dft_c2r_2d" fftw-plan-dft-c2r-2d) fftw-plan + (n0 :int) + (n1 :int) + (in :pointer) + (out :pointer) + (flags :uint)) + +(declaim (inline fftw-plan-dft-c2r-3d)) +(cffi:defcfun ("fftw_plan_dft_c2r_3d" fftw-plan-dft-c2r-3d) fftw-plan + (n0 :int) + (n1 :int) + (n2 :int) + (in :pointer) + (out :pointer) + (flags :uint)) + +(declaim (inline fftw-plan-dft-c2r)) +(cffi:defcfun ("fftw_plan_dft_c2r" fftw-plan-dft-c2r) fftw-plan + (rank :int) + (n (:pointer :int)) + (in :pointer) + (out :pointer) + (flags :uint)) + +(declaim (inline fftw-plan-many-dft)) +(cffi:defcfun ("fftw_plan_many_dft" fftw-plan-many-dft) fftw-plan + (rank :int) + (n (:pointer :int)) + (howmany :int) + (in (:pointer fftw-complex-struct)) + (inembed (:pointer :int)) + (istride :int) + (idist :int) + (out (:pointer fftw-complex-struct)) + (onembed (:pointer :int)) + (ostride :int) + (odist :int) + (sign :int) + (flags :uint)) + +(declaim (inline fftw-plan-many-dft_r2c)) +(cffi:defcfun ("fftw_plan_many_dft_r2c" fftw-plan-many-dft-r2c) fftw-plan + (rank :int) + (n (:pointer :int)) + (howmany :int) + (in (:pointer :double)) + (inembed (:pointer :int)) + (istride :int) + (idist :int) + (out (:pointer fftw-complex-struct)) + (onembed (:pointer :int)) + (ostride :int) + (odist :int) + (flags :uint)) + +(declaim (inline fftw-plan-many-dft_c2r)) +(cffi:defcfun ("fftw_plan_many_dft_c2r" fftw-plan-many-dft-c2r) fftw-plan + (rank :int) + (n (:pointer :int)) + (howmany :int) + (in (:pointer fftw-complex-struct)) + (inembed (:pointer :int)) + (istride :int) + (idist :int) + (out (:pointer :double)) + (onembed (:pointer :int)) + (ostride :int) + (odist :int) + (flags :uint)) + +(declaim (inline fftw-plan-many-dft_r2r)) +(cffi:defcfun ("fftw_plan_many_dft_r2r" fftw-plan-many-dft-r2r) fftw-plan + (rank :int) + (n (:pointer :int)) + (howmany :int) + (in (:pointer :double)) + (inembed (:pointer :int)) + (istride :int) + (idist :int) + (out (:pointer :double)) + (onembed (:pointer :int)) + (ostride :int) + (odist :int) + (kind (:pointer :int)) + (flags :uint)) + + +(declaim (inline fftw-malloc)) +(cffi:defcfun ("fftw_malloc" fftw-malloc) (:pointer fftw-complex-struct) + (n :int)) + +(declaim (inline fftw-execute)) +(cffi:defcfun ("fftw_execute" fftw-execute) :void + (p fftw-plan)) + +(declaim (inline fftw-destroy-plan)) +(cffi:defcfun ("fftw_destroy_plan" fftw-destroy-plan) :void + (p fftw-plan)) + +(declaim (inline fftw-free)) +(cffi:defcfun ("fftw_free" fftw-free) :void + (p :pointer)) + + +;;; Wisdom functions + +(cffi:defcfun ("fftw_import_wisdom_from_string" fftw-import-wisdom-from-string) :int + (input-string :string)) + +(cffi:defcfun ("fftw_export_wisdom_to_string" fftw-export-wisdom-to-string) :string+ptr + ) + +(cffi:defcfun free :void + (ptr :pointer)) diff --git a/package.lisp b/package.lisp new file mode 100644 index 0000000..6ed7ede --- /dev/null +++ b/package.lisp @@ -0,0 +1,59 @@ +;;;; -*- Mode: Lisp; Syntax: ANSI-Common-Lisp; Base: 10 -*- +;;;; ************************************************************************* +;;;; FILE IDENTIFICATION +;;;; +;;;; Name: package.lisp +;;;; Purpose: Package definition for fftw3 package +;;;; Programmer: Kevin M. Rosenberg +;;;; Date Started: Jan 2009 +;;;; +;;;; $Id$ +;;;; +;;;; This file, part of FFTW3, is Copyright (c) 2009 by Kevin M. Rosenberg +;;;; +;;;; FFTW3 users are granted the rights to distribute and use this software +;;;; as governed by the terms of the Lisp Lesser GNU Public License +;;;; (http://opensource.franz.com/preamble.html), also known as the LLGPL. +;;;; ************************************************************************* + +(in-package #:cl-user) + +(defpackage #:cl-fftw3 + (:nicknames #:fftw3) + (:use #:common-lisp) + (:export + ;; transform.lisp + #:fftw-r2hc-1d + #:make-fftw-r2r-1d-multi + #:fftw-r2r-1d-multi + #:hc-to-mag-phase + #:fftw-hc2r-1d + #:fftw-r2c-1d + #:make-fftw-r2c-1d-multi + #:make-fftw-c2r-1d-multi + #:destroy-fftw-multi + #:fftw-r2c-1d-multi + #:fftw-c2r-1d-multi + #:complex-to-mag-phase + #:fftw-c2r-1d + #:fftw-c-1d + + ;; wisdom.lisp + #:import-user-wisdom + #:export-user-wisdom + + ;; specials + #:+fftw-forward+ + #:+fftw-backward+ + #:+fftw-hc2r+ + #:+fftw-r2hc+ + + #:+fftw-measure+ + #:+fftw-destroy-input+ + #:+fftw-unaligned+ + #:+fftw-conservative-memory+ + #:+fftw-exhaustive+ + #:+fftw-preserve-input+ + #:+fftw-patiento+ + #:+fftw-estimate+ + )) diff --git a/specials.lisp b/specials.lisp new file mode 100644 index 0000000..53ae7d0 --- /dev/null +++ b/specials.lisp @@ -0,0 +1,47 @@ +;;;; -*- Mode: Lisp; Syntax: ANSI-Common-Lisp; Base: 10 -*- +;;;; ************************************************************************* +;;;; FILE IDENTIFICATION +;;;; +;;;; Name: specials.lisp +;;;; Purpose: Special data declarations for FFTW3 package +;;;; Programmer: Kevin M. Rosenberg +;;;; Date Started: March 2009 +;;;; +;;;; $Id$ +;;;; +;;;; This file, part of FFTW3, is Copyright (c) 2009 by Kevin M. Rosenberg +;;;; +;;;; FFTW3 users are granted the rights to distribute and use this software +;;;; as governed by the terms of the Lisp Lesser GNU Public License +;;;; (http://opensource.franz.com/preamble.html), also known as the LLGPL. +;;;; ************************************************************************* + +(in-package #:fftw3) + +(defconstant +fftw-forward+ -1) +(defconstant +fftw-backward+ +1) + +(defconstant +fftw-r2hc+ 0) +(defconstant +fftw-hc2r+ 1) + +(defconstant +fftw-measure+ 0) +(defconstant +fftw-destroy-input+ (ash 1 0)) +(defconstant +fftw-unaligned+ (ash 1 1)) +(defconstant +fftw-conservative-memory+ (ash 1 2)) +(defconstant +fftw-exhaustive+ (ash 1 3)) +(defconstant +fftw-preserve-input+ (ash 1 4)) +(defconstant +fftw-patiento+ (ash 1 5)) +(defconstant +fftw-estimate+ (ash 1 6)) + +(defparameter *user-wisdom-file* + (concatenate 'string (namestring (user-homedir-pathname)) ".fftw3-wisdom")) + +(defvar *standard-optimize-settings* + '(optimize + speed + (safety 0) + (space 0) + (debug 1) + (compilation-speed 0) + #+:lispworks (hcl:fixnum-safety 0)) + "The standard optimize settings used by most declaration expressions.") diff --git a/tests.lisp b/tests.lisp new file mode 100644 index 0000000..0803415 --- /dev/null +++ b/tests.lisp @@ -0,0 +1,231 @@ +(defpackage #:cl-fftw3-tests + (:nicknames #:fftw3-tests) + (:use #:cl #:rtest #:cl-fftw3)) + +(in-package #:cl-fftw3-tests) + +(defun double-float-close-p (d1 d2 &key (epsilon 1d-7)) + (when (and (zerop d1) (zerop d2)) + (return-from double-float-close-p t)) + + (let* ((delta (abs (- d1 d2))) + (eps (/ delta (max (abs d1) (abs d2))))) + (when (or (and (zerop (min (abs d1) (abs d2))) (< (max (abs d1) (abs d2)) epsilon)) + (< eps epsilon)) + t))) + +(defun vectors-close-p (v1 v2) + (unless (eql (length v1) (length v2)) + (return-from vectors-close-p nil)) + + (dotimes (i (length v1)) + (unless (and (double-float-close-p (realpart (aref v1 i)) (realpart (aref v2 i))) + (double-float-close-p (imagpart (aref v1 i)) (imagpart (aref v2 i)))) + (return-from vectors-close-p nil))) + t) + +(rem-all-tests) + +(defparameter *a7* #(0 1 2 3 2 1 0)) +(defparameter *a8* #(0 1 2 3 3 2 1 0)) +(defparameter *ac7* #(#C(0 0) #C(1 1) #C(2 2) #C(3 3) #C(2 2) #C(1 1) #C(0 0))) +(defparameter *ac8* #(#C(0 0) #C(1 1) #C(2 2) #C(3 3) #C(3 3) #C(2 2) #C(1 1) #C(0 0))) +(defparameter *cc7* #(#C(0 1) #C(1 1) #C(2 2) #C(3 3) #C(2 2) #C(1 1) #C(0 0))) +(defparameter *cc8* #(#C(0 1) #C(1 1) #C(2 2) #C(3 3) #C(3 3) #C(2 2) #C(1 1) #C(0 0))) + +(defparameter *b7* #(3 2 1 0 1 2 3)) +(defparameter *b8* #(3 2 1 0 0 1 2 3)) +(defparameter *bc7* #(#C(3 3) #C(2 2) #C(1 1) #C(0 0) #C(1 1) #C(2 2) #C(3 3))) +(defparameter *bc8* #(#C(3 3) #C(2 2) #C(1 1) #C(0 0) #C(0 0) #C(1 1) #C(2 2) #C(3 3))) + +(defparameter *a7hc* (fftw-r2hc-1d *a7*)) +(defparameter *a8hc* (fftw-r2hc-1d *a8*)) +(defparameter *a7c* (fftw-r2c-1d *a7*)) +(defparameter *a8c* (fftw-r2c-1d *a8*)) +(defparameter *ac7c* (fftw-c-1d *a7*)) +(defparameter *ac8c* (fftw-c-1d *a8*)) + +(defparameter *a7c-octave* #(#C(1.28571428571428580945d0 0d0) + #C(-0.64984533421747225912d0 -0.31294901910963018876d0) + #C(0.02743163880429938875d0 0.03439818705768063478d0) + #C(-0.02044344744397007946d0 -0.08956859554733599682d0))) + +(defparameter *a8c-octave* #(#C(1.50000000000000000000d0 0.00000000000000000000d0) + #C(-0.72855339059327373086d0 -0.30177669529663686543d0) + #C(0.00000000000000000000d0 0.00000000000000000000d0) + #C(-0.02144660940672621363d0 -0.05177669529663689318d0) + #C(0.00000000000000000000d0 0.00000000000000000000d0))) + +(defparameter *ac7c-octave* #(#C(1.28571428571428580945d0 1.28571428571428580945d0) + #C(-0.33689631510784201485d0 -0.96279435332710250339d0) + #C(-0.00696654825338124516d0 0.06182982586198002700d0) + #C(0.06912514810336592430d0 -0.11001204299130606934d0) + #C(-0.11001204299130606934d0 0.06912514810336592430d0) + #C(0.06182982586198002700d0 -0.00696654825338124516d0) + #C(-0.96279435332710250339d0 -0.33689631510784201485d0))) + +(defparameter *cc7c-octave* #(#C(1.28571428571428580945d0 1.42857142857142860315d0) + #C(-0.33689631510784201485d0 -0.81993721046995970969d0) + #C(-0.00696654825338124516d0 0.20468696871912286928d0) + #C(0.06912514810336592430d0 0.03284509986583677293d0) + #C(-0.11001204299130606934d0 0.21198229096050877351d0) + #C(0.06182982586198002700d0 0.13589059460376159971d0) + #C(-0.96279435332710250339d0 -0.19403917225069916563d0))) + +(defparameter *ac8c-octave* #(#C(1.50000000000000000000d0 1.50000000000000000000d0) + #C(-0.42677669529663686543d0 -1.03033008588991070731d0) + #C(0.00000000000000000000d0 0.00000000000000000000d0) + #C(0.03033008588991070731d0 -0.07322330470336310682d0) + #C(0.00000000000000000000d0 0.00000000000000000000d0) + #C(-0.07322330470336310682d0 0.03033008588991070731d0) + #C(0.00000000000000000000d0 0.00000000000000000000d0) + #C(-1.03033008588991070731d0 -0.42677669529663686543d0))) + +(defparameter *cc8c-octave* #(#C(1.50000000000000000000d0 1.62500000000000000000d0) + #C(-0.42677669529663686543d0 -0.90533008588991070731) + #C(0.00000000000000000000d0 0.12500000000000000000d0) + #C(0.03033008588991070731d0 0.05177669529663689318d0) + #C(0.00000000000000000000d0 0.12500000000000000000d0) + #C(-0.07322330470336310682d0 0.15533008588991070731d0) + #C(0.00000000000000000000d0 0.12500000000000000000d0) + #C(-1.03033008588991070731d0 -0.30177669529663686543d0))) + +(deftest :len.1 (length *a7*) 7) +(deftest :len.2 (length *a8*) 8) +(deftest :len.3 (length *a7hc*) 7) +(deftest :len.4 (length *a8hc*) 8) +(deftest :len.5 (length *a7c*) 4) +(deftest :len.6 (length *a8c*) 5) +(deftest :len.7 (length *b7*) 7) +(deftest :len.8 (length *b8*) 8) + +(deftest :eql.1 (vectors-close-p *a7* (fftw-hc2r-1d *a7hc*)) t) +(deftest :eql.2 (vectors-close-p *a8* (fftw-hc2r-1d *a8hc*)) t) +(deftest :eql.3 (vectors-close-p *a7* (fftw-c2r-1d *a7c*)) t) +(deftest :eql.4 (vectors-close-p *a8* (fftw-c2r-1d *a8c*)) t) +(deftest :eql.5 (vectors-close-p *a7* (fftw-c-1d *ac7c* :normalize nil :direction +fftw-backward+)) t) +(deftest :eql.6 (vectors-close-p *a8* (fftw-c-1d *ac8c* :normalize nil :direction +fftw-backward+)) t) + +(deftest :eqlc.1 (vectors-close-p + *ac7* + (fftw-c-1d (fftw-c-1d *ac7* :normalize t :direction +fftw-forward+) + :normalize nil :direction +fftw-backward+)) + t) + +(deftest :multi.1 + (let* ((multi-r2hc (make-fftw-r2r-1d-multi (length *a7*) +fftw-r2hc+)) + (multi-hc2r (make-fftw-r2r-1d-multi (length *a7*) +fftw-hc2r+)) + (a7hc (copy-seq (fftw-r2r-1d-multi multi-r2hc *a7*))) + (b7hc (copy-seq (fftw-r2r-1d-multi multi-r2hc *b7*))) + (a7r (copy-seq (fftw-r2r-1d-multi multi-hc2r a7hc :normalize nil))) + (b7r (copy-seq (fftw-r2r-1d-multi multi-hc2r b7hc :normalize nil))) + (eq (and (vectors-close-p *a7* a7r) (vectors-close-p *b7* b7r) t))) + (destroy-fftw-multi multi-r2hc) + (destroy-fftw-multi multi-hc2r) + eq) + t) + +(deftest :multi.2 + (let* ((multi-r2hc (make-fftw-r2r-1d-multi (length *a8*) +fftw-r2hc+)) + (multi-hc2r (make-fftw-r2r-1d-multi (length *a8*) +fftw-hc2r+)) + (a8hc (copy-seq (fftw-r2r-1d-multi multi-r2hc *a8*))) + (b8hc (copy-seq (fftw-r2r-1d-multi multi-r2hc *b8*))) + (a8r (copy-seq (fftw-r2r-1d-multi multi-hc2r a8hc :normalize nil))) + (b8r (copy-seq (fftw-r2r-1d-multi multi-hc2r b8hc :normalize nil))) + (eq (and (vectors-close-p *a8* a8r) (vectors-close-p *b8* b8r) t))) + (destroy-fftw-multi multi-r2hc) + (destroy-fftw-multi multi-hc2r) + eq) + t) + +(deftest :multi.3 + (let* ((multi-r2c (make-fftw-r2c-1d-multi (length *a7*))) + (a7c (copy-seq (fftw-r2c-1d-multi multi-r2c *a7*))) + (b7c (copy-seq (fftw-r2c-1d-multi multi-r2c *b7*))) + (a7 (fftw-c2r-1d a7c)) + (b7 (fftw-c2r-1d b7c)) + (eq (and (vectors-close-p *a7* a7) + (vectors-close-p *b7* b7)))) + (destroy-fftw-multi multi-r2c) + eq) + t) + +(deftest :multi.4 + (let* ((multi-r2c (make-fftw-r2c-1d-multi (length *a8*))) + (a8c (copy-seq (fftw-r2c-1d-multi multi-r2c *a8*))) + (b8c (copy-seq (fftw-r2c-1d-multi multi-r2c *b8*))) + (a8 (fftw-c2r-1d a8c)) + (b8 (fftw-c2r-1d b8c)) + (eq (and (vectors-close-p *a8* a8) + (vectors-close-p *b8* b8)))) + (destroy-fftw-multi multi-r2c) + eq) + t) + +(deftest :multi.5 + (let* ((multi-r2c (make-fftw-r2c-1d-multi (length *a7*))) + (a7c (copy-seq (fftw-r2c-1d-multi multi-r2c *a7*))) + (b7c (copy-seq (fftw-r2c-1d-multi multi-r2c *b7*))) + (multi-c2r (make-fftw-c2r-1d-multi (length a7c))) + (a7r (copy-seq (fftw-c2r-1d-multi multi-c2r a7c))) + (b7r (copy-seq (fftw-c2r-1d-multi multi-c2r b7c))) + (eq (and (vectors-close-p *a7* a7r) (vectors-close-p *b7* b7r) t))) + (destroy-fftw-multi multi-r2c) + (destroy-fftw-multi multi-c2r) + eq) + t) + +(deftest :multi.6 + (let* ((multi-r2c (make-fftw-r2c-1d-multi (length *a8*))) + (a8c (copy-seq (fftw-r2c-1d-multi multi-r2c *a8*))) + (b8c (copy-seq (fftw-r2c-1d-multi multi-r2c *b8*))) + (multi-c2r (make-fftw-c2r-1d-multi (length a8c))) + (a8r (copy-seq (fftw-c2r-1d-multi multi-c2r a8c))) + (b8r (copy-seq (fftw-c2r-1d-multi multi-c2r b8c))) + (eq (and (vectors-close-p *a8* a8r) (vectors-close-p *b8* b8r) t))) + (destroy-fftw-multi multi-r2c) + (destroy-fftw-multi multi-c2r) + eq) + t) + +(deftest :c2c.1 + (let* ((ac7c (fftw-c-1d *ac7*)) + (bc7c (fftw-c-1d *bc7*)) + (cc7c (fftw-c-1d *cc7*)) + (ac7i (fftw-c-1d ac7c :direction +fftw-backward+ :normalize nil)) + (bc7i (fftw-c-1d bc7c :direction +fftw-backward+ :normalize nil)) + (cc7i (fftw-c-1d cc7c :direction +fftw-backward+ :normalize nil)) + (eq (and (vectors-close-p *ac7* ac7i) (vectors-close-p *bc7* bc7i) + (vectors-close-p *cc7* cc7i)))) + eq) + t) + +(deftest :c2c.2 + (let* ((ac8c (fftw-c-1d *ac8*)) + (bc8c (fftw-c-1d *bc8*)) + (cc8c (fftw-c-1d *cc8*)) + (ac8i (fftw-c-1d ac8c :direction +fftw-backward+ :normalize nil)) + (bc8i (fftw-c-1d bc8c :direction +fftw-backward+ :normalize nil)) + (cc8i (fftw-c-1d cc8c :direction +fftw-backward+ :normalize nil)) + (eq (and (vectors-close-p *ac8* ac8i) (vectors-close-p *bc8* bc8i) + (vectors-close-p *cc8* cc8i)))) + eq) + t) + +(deftest :r2c-octave.1 + (vectors-close-p *a7c* *a7c-octave*) t) + +(deftest :r2c-octave.2 + (vectors-close-p *a8c* *a8c-octave*) t) + +(deftest :c2c-octave.1 + (vectors-close-p (fftw-c-1d *ac7*) *ac7c-octave*) t) + +(deftest :c2c-octave.2 + (vectors-close-p (fftw-c-1d *cc7*) *cc7c-octave*) t) + +(deftest :c2c-octave.3 + (vectors-close-p (fftw-c-1d *ac8*) *ac8c-octave*) t) + +(deftest :c2c-octave.4 + (vectors-close-p (fftw-c-1d *cc8*) *cc8c-octave*) t) diff --git a/transform.lisp b/transform.lisp new file mode 100644 index 0000000..adb5b31 --- /dev/null +++ b/transform.lisp @@ -0,0 +1,403 @@ +;;;; -*- Mode: Lisp; Syntax: ANSI-Common-Lisp; Base: 10 -*- +;;;; ************************************************************************* +;;;; FILE IDENTIFICATION +;;;; +;;;; Name: transform.lisp +;;;; Purpose: Transformation functions +;;;; Programmer: Kevin M. Rosenberg +;;;; Date Started: March 2009 +;;;; +;;;; $Id$ +;;;; +;;;; This file, part of FFTW3, is Copyright (c) 2009 by Kevin M. Rosenberg +;;;; +;;;; FFTW3 users are granted the rights to distribute and use this software +;;;; as governed by the terms of the Lisp Lesser GNU Public License +;;;; (http://opensource.franz.com/preamble.html), also known as the LLGPL. +;;;; ************************************************************************* + +(in-package #:fftw3) + +(defconstant +size-double+ 8) + +(defun fftw-r2hc-1d (in &key (start 0) (count nil) (normalize t)) + (declare #.*standard-optimize-settings*) + (multiple-value-bind (start count) (normalize-input-range in start count) + (declare (fixnum start count)) + + (let ((in-pos start) + (out-c (make-array count :element-type 'double-float)) + (in-cf (fftw-malloc (* count +size-double+))) + (out-cf (fftw-malloc (* count +size-double+)))) + (declare (fixnum in-pos)) + (dotimes (i count) + (declare (fixnum i)) + (setf (cffi:mem-aref in-cf :double i) (coerce (aref in in-pos) 'double-float)) + (incf in-pos)) + + (let ((p (fftw-plan-r2r-1d count in-cf out-cf +fftw-r2hc+ + (logior +fftw-estimate+ +fftw-destroy-input+)))) + (fftw-execute p) + (fftw-destroy-plan p)) + + (if normalize + (let ((factor (coerce (/ 1 count) 'double-float))) + (declare (double-float factor)) + (dotimes (i count) + (declare (fixnum i)) + (setf (aref out-c i) (* factor (cffi:mem-aref out-cf :double i))))) + (dotimes (i count) + (declare (fixnum i)) + (setf (aref out-c i) (cffi:mem-aref out-cf :double i)))) + + (fftw-free in-cf) + (fftw-free out-cf) + + out-c))) + +(defun make-fftw-r2r-1d-multi (count type) + (declare #.*standard-optimize-settings*) + (declare (fixnum count type)) + (let* ((out-c (make-array count :element-type 'double-float)) + (in-cf (fftw-malloc (* count +size-double+))) + (out-cf (fftw-malloc (* count +size-double+))) + (plan (fftw-plan-r2r-1d count in-cf out-cf type + (logior +fftw-estimate+ +fftw-destroy-input+)))) + (make-instance 'fftw-r2r-1d-multi + :in-count count :out-c out-c :plan plan + :in-cf in-cf :out-cf out-cf))) + +(defun fftw-r2r-1d-multi (multi in &key (start 0) (count nil) (normalize t)) + (declare #.*standard-optimize-settings*) + (multiple-value-bind (start count) (normalize-input-range in start count) + (declare (fixnum start count)) + (unless (equal (in-count multi) count) + (error "Different plan and vector lengths.")) + + (let ((in-cf (in-cf multi)) + (out-cf (out-cf multi)) + (out-c (out-c multi)) + (in-pos start)) + (declare (fixnum in-pos)) + + (dotimes (i count) + (declare (fixnum i)) + (setf (cffi:mem-aref in-cf :double i) (coerce (aref in in-pos) 'double-float)) + (incf in-pos)) + + (fftw-execute (plan multi)) + (if normalize + (let ((factor (coerce (/ 1 count) 'double-float))) + (declare (double-float factor)) + (dotimes (i count) + (declare (fixnum i)) + (setf (aref out-c i) (* factor (cffi:mem-aref out-cf :double i))))) + (dotimes (i count) + (declare (fixnum i)) + (setf (aref out-c i) (cffi:mem-aref out-cf :double i)))) + + out-c))) + + +(defun fftw-hc2r-1d (in) + (declare #.*standard-optimize-settings*) + (let* ((count (length in)) + (out-c (make-array count :element-type 'double-float)) + (in-cf (fftw-malloc (* count +size-double+))) + (out-cf (fftw-malloc (* count +size-double+)))) + (declare (fixnum count)) + + (dotimes (i count) + (declare (fixnum i)) + (setf (cffi:mem-aref in-cf :double i) (aref in i))) + + (let ((p (fftw-plan-r2r-1d count in-cf out-cf +fftw-hc2r+ + (logior +fftw-estimate+ +fftw-destroy-input+)))) + (fftw-execute p) + (fftw-destroy-plan p)) + + (dotimes (i count) + (declare (fixnum i)) + (setf (aref out-c i) (cffi:mem-aref out-cf :double i))) + (fftw-free in-cf) + (fftw-free out-cf) + out-c)) + +(defun fftw-r2c-1d (in &key (start 0) (count nil) (normalize t)) + (declare #.*standard-optimize-settings*) + (multiple-value-bind (start count) (normalize-input-range in start count) + (declare (fixnum start count)) + + (let* ((out-n (1+ (floor count 2))) + (in-cf (fftw-malloc (* count +size-double+))) + (out-cf (fftw-malloc (* 2 out-n +size-double+))) + (in-pos start)) + (declare (fixnum in-pos out-n)) + + (dotimes (i count) + (declare (fixnum i)) + (setf (cffi:mem-aref in-cf :double i) (coerce (aref in in-pos) 'double-float)) + (incf in-pos)) + + (let ((p (fftw-plan-dft-r2c-1d count in-cf out-cf + (logior +fftw-estimate+ +fftw-destroy-input+)))) + (fftw-execute p) + (fftw-destroy-plan p)) + + (let ((out (make-array out-n :element-type '(or rational complex))) + (opos 0)) + (declare (fixnum opos)) + (if normalize + (let ((factor (coerce (/ 1 count) 'double-float))) + (declare (double-float factor)) + (dotimes (i out-n) + (declare (fixnum i)) + (let ((impos (1+ opos))) + (declare (fixnum impos)) + (setf (aref out i) (complex (* (cffi:mem-aref out-cf :double opos) factor) + (* (cffi:mem-aref out-cf :double impos) factor))) + (incf opos 2)))) + (dotimes (i out-n) + (declare (fixnum i)) + (let ((impos (1+ opos))) + (declare (fixnum impos)) + (setf (aref out i) (complex (cffi:mem-aref out-cf :double opos) + (cffi:mem-aref out-cf :double impos))) + (incf opos 2)))) + (fftw-free in-cf) + (fftw-free out-cf) + out)))) + +(defun make-fftw-r2c-1d-multi (count &key (complex-output t)) + (declare #.*standard-optimize-settings*) + (declare (fixnum count)) + (let* ((out-n (1+ (floor count 2))) + (in-cf (fftw-malloc (* count +size-double+))) + (out-cf (fftw-malloc (* 2 out-n +size-double+))) + (plan (fftw-plan-dft-r2c-1d count in-cf out-cf + (logior +fftw-estimate+ +fftw-destroy-input+)))) + (declare (fixnum out-n)) + + (unless complex-output + (incf out-n out-n)) + + (make-instance 'fftw-r2c-1d-multi + :in-count count :out-count out-n :plan plan + :in-cf in-cf + :out-cf out-cf + :complex-output complex-output + :out + (if complex-output + (make-array out-n :element-type '(or double-float complex) :initial-element 0d0) + (make-array out-n :element-type 'double-float))))) + +(defgeneric destroy-fftw-multi (multi)) +(defmethod destroy-fftw-multi ((multi fftw-multi)) + (declare #.*standard-optimize-settings*) + (when (plan multi) + (fftw-destroy-plan (plan multi)) + (setf (plan multi) nil)) + (dolist (slot '(in-cf out-cf)) + (when (slot-value multi slot) + (fftw-free (slot-value multi slot)) + (setf (slot-value multi slot) nil))) + (dolist (slot '(plan-even plan-odd)) + (when (and (slot-exists-p multi slot) + (slot-value multi slot)) + (fftw-destroy-plan (slot-value multi slot)) + (setf (slot-value multi slot) nil)))) + + +(defun fftw-r2c-1d-multi (multi in &key (start 0) (normalize t)) + (declare #.*standard-optimize-settings*) + (multiple-value-bind (start count) (normalize-input-range in start (in-count multi)) + (declare (fixnum start count)) + + (unless (equal (length in) (in-count multi)) + (error "In count of multi plan doesn't equal length of in vector.")) + + (let ((in-cf (in-cf multi)) + (out-cf (out-cf multi)) + (out (out multi)) + (in-pos start) + (opos 0) + (in-count (in-count multi)) + (out-count (out-count multi))) + (declare (fixnum in-pos opos in-count out-count)) + + (dotimes (i in-count) + (declare (fixnum i)) + (setf (cffi:mem-aref in-cf :double i) (coerce (aref in in-pos) 'double-float)) + (incf in-pos)) + + (fftw-execute (plan multi)) + + (setf opos 0) + (cond + ((complex-output multi) + (if normalize + (let ((factor (coerce (/ 1 count) 'double-float))) + (declare (double-float factor)) + (dotimes (i out-count) + (declare (fixnum i)) + (let ((impos (1+ opos))) + (declare (fixnum impos)) + (setf (aref out i) (complex (* (cffi:mem-aref out-cf :double opos) factor) + (* (cffi:mem-aref out-cf :double impos) factor))) + (incf opos 2)))) + (dotimes (i out-count) + (declare (fixnum i)) + (let ((impos (1+ opos))) + (declare (fixnum impos)) + (setf (aref out i) (complex (cffi:mem-aref out-cf :double opos) + (cffi:mem-aref out-cf :double impos))) + (incf opos 2))))) + (t + (when normalize + (let ((o2 (* out-count 2)) + (factor (coerce (/ 1 count) 'double-float))) + (declare (fixnum o2) + (double-float factor)) + (dotimes (i o2) + (declare (fixnum i)) + (setf (aref out i) (* factor (aref out i)))))))) + out))) + +(defun fftw-c2r-1d (in) + (declare #.*standard-optimize-settings*) + (let* ((count (length in)) + (out-n (if (complex-is-even-p in) + (* 2 (1- count)) + (1+ (* 2 (1- count))))) + (out-c (make-array out-n :element-type 'double-float)) + (in-cf (fftw-malloc (* 2 count +size-double+))) + (out-cf (fftw-malloc (* 2 count +size-double+))) + (pos 0)) + (declare (fixnum count out-n pos)) + + (setq pos 0) + (dotimes (i count) + (declare (fixnum i)) + (let ((c (aref in i))) + (declare (complex c)) + (setf (cffi:mem-aref in-cf :double pos) (coerce (realpart c) 'double-float)) + (incf pos) + (setf (cffi:mem-aref in-cf :double pos) (coerce (imagpart c) 'double-float)) + (incf pos))) + + (let ((p (fftw-plan-dft-c2r-1d out-n in-cf out-cf + (logior +fftw-estimate+ +fftw-destroy-input+)))) + (fftw-execute p) + (fftw-destroy-plan p)) + (dotimes (i out-n) + (declare (fixnum i)) + (setf (aref out-c i) (cffi:mem-aref out-cf :double i))) + + (fftw-free in-cf) + (fftw-free out-cf) + out-c)) + +(defun make-fftw-c2r-1d-multi (count) + (declare #.*standard-optimize-settings*) + (declare (fixnum count)) + (let* ((out-count-even (* 2 (1- count))) + (out-count-odd (1+ (* 2 (1- count)))) + (in-cf (fftw-malloc (* 2 count +size-double+))) + (out-cf (fftw-malloc (* 2 (max out-count-odd out-count-even) +size-double+))) ;; use larger size + (plan-odd (fftw-plan-dft-c2r-1d out-count-odd in-cf out-cf + (logior +fftw-estimate+ +fftw-destroy-input+))) + (plan-even (fftw-plan-dft-c2r-1d out-count-even in-cf out-cf +fftw-measure+))) + (declare (fixnum out-count-even out-count-odd)) + + (make-instance 'fftw-c2r-1d-multi + :in-count count + :plan-even plan-even + :plan-odd plan-odd + :in-cf in-cf + :out-cf out-cf + :out-count-even out-count-even + :out-count-odd out-count-odd + :out-even (make-array out-count-even :element-type 'double-float) + :out-odd (make-array out-count-odd :element-type 'double-float)))) + +(defun fftw-c2r-1d-multi (multi in) + (declare #.*standard-optimize-settings*) + (let* ((is-even (complex-is-even-p in)) + (out (if is-even (out-even multi) (out-odd multi))) + (out-count (if is-even (out-count-even multi) (out-count-odd multi))) + (plan (if is-even (plan-even multi) (plan-odd multi))) + (count (in-count multi)) + (in-cf (in-cf multi)) + (out-cf (out-cf multi)) + (pos 0)) + (declare (fixnum count out-count pos)) + + (setq pos 0) + (dotimes (i count) + (declare (fixnum i)) + (let ((c (aref in i))) + (setf (cffi:mem-aref in-cf :double pos) (coerce (realpart c) 'double-float)) + (incf pos) + (setf (cffi:mem-aref in-cf :double pos) (coerce (imagpart c) 'double-float)) + (incf pos))) + + (fftw-execute plan) + + (dotimes (i out-count) + (declare (fixnum i)) + (setf (aref out i) (cffi:mem-aref out-cf :double i))) + + out)) + +(defun fftw-c-1d (in &key (direction +fftw-forward+) + (start 0) (count nil) (normalize t)) + (declare #.*standard-optimize-settings*) + (multiple-value-bind (start count) (normalize-input-range in start count) + (declare (fixnum start count)) + + (let* ((out-n count) + (out (make-array out-n :element-type '(or rational complex))) + (in-cf (fftw-malloc (* 2 count +size-double+))) + (out-cf (fftw-malloc (* 2 out-n +size-double+))) + (opos 0) + (ipos 0)) + (declare (fixnum out-n ipos opos)) + + (setf opos 0) + (setf ipos start) + (dotimes (i count) + (declare (fixnum i)) + (let ((c (aref in ipos))) + (incf ipos) + (setf (cffi:mem-aref in-cf :double opos) (coerce (realpart c) 'double-float)) + (incf opos) + (setf (cffi:mem-aref in-cf :double opos) (coerce (imagpart c) 'double-float)) + (incf opos))) + + (let ((p (fftw-plan-dft-1d count in-cf out-cf direction + (logior +fftw-estimate+ +fftw-destroy-input+)))) + (fftw-execute p) + (fftw-destroy-plan p)) + + (setf opos 0) + (if normalize + (let ((factor (coerce (/ 1 count) 'double-float))) + (declare (double-float factor)) + (dotimes (i out-n) + (declare (fixnum i)) + (let ((impos (1+ opos))) + (declare (fixnum impos)) + (setf (aref out i) (complex (* factor (cffi:mem-aref out-cf :double opos)) + (* factor (cffi:mem-aref out-cf :double impos))))) + (incf opos 2))) + (dotimes (i out-n) + (declare (fixnum i)) + (let ((impos (1+ opos))) + (declare (fixnum impos)) + (setf (aref out i) (complex (cffi:mem-aref out-cf :double opos) + (cffi:mem-aref out-cf :double impos)))) + (incf opos 2))) + + (fftw-free in-cf) + (fftw-free out-cf) + out))) diff --git a/wisdom.lisp b/wisdom.lisp new file mode 100644 index 0000000..0b79318 --- /dev/null +++ b/wisdom.lisp @@ -0,0 +1,35 @@ +;;;; -*- Mode: Lisp; Syntax: ANSI-Common-Lisp; Base: 10 -*- +;;;; ************************************************************************* +;;;; FILE IDENTIFICATION +;;;; +;;;; Name: wisdowm.lisp +;;;; Purpose: Functions for handling FFTW wisdom +;;;; Programmer: Kevin M. Rosenberg +;;;; Date Started: March 2009 +;;;; +;;;; $Id$ +;;;; +;;;; This file, part of FFTW3, is Copyright (c) 2009 by Kevin M. Rosenberg +;;;; +;;;; FFTW3 users are granted the rights to distribute and use this software +;;;; as governed by the terms of the Lisp Lesser GNU Public License +;;;; (http://opensource.franz.com/preamble.html), also known as the LLGPL. +;;;; ************************************************************************* + +(in-package #:fftw3) + +(defun import-user-wisdom () + (let ((file-contents (kmrcl:read-file-to-string *user-wisdom-file*))) + (when (and (stringp file-contents) (plusp (length file-contents))) + (fftw-import-wisdom-from-string file-contents)))) + +(defun export-user-wisdom () + (let ((str+ptr (fftw-export-wisdom-to-string))) + (when (probe-file *user-wisdom-file*) + (delete-file *user-wisdom-file*)) + (with-open-file (out *user-wisdom-file* :direction :output :if-exists :overwrite + :if-does-not-exist :create) + (format out "~A" (first str+ptr))) + (free (second str+ptr)))) + + -- 2.34.1