common.lisp

   1 ;;;; -*- Mode: Lisp; Syntax: ANSI-Common-Lisp; Base: 10 -*-
   2 ;;;; *************************************************************************
   3 ;;;; FILE IDENTIFICATION
   4 ;;;;
   5 ;;;; Name:          common.lisp
   6 ;;;; Purpose:       Common functions for FFTW3 package
   7 ;;;; Programmer:    Kevin M. Rosenberg
   8 ;;;; Date Started:  Jan 2009
   9 ;;;;
  10 ;;;; $Id$
  11 ;;;;
  12 ;;;; This file, part of FFTW3, is Copyright (c) 2009 by Kevin M. Rosenberg
  13 ;;;;
  14 ;;;; FFTW3 users are granted the rights to distribute and use this software
  15 ;;;; as governed by the terms of the Lisp Lesser GNU Public License
  16 ;;;; (http://opensource.franz.com/preamble.html), also known as the LLGPL.
  17 ;;;; *************************************************************************
  18
  19 (in-package #:fftw3)
  20
  21 (defclass fftw-multi ()
  22   ((plan :initarg :plan :initform nil :accessor plan)))
  23
  24 (defclass fftw-r2r-1d-multi (fftw-multi)
  25   ((in-count :initarg :in-count :initform 0 :type fixnum :accessor in-count)
  26    (out-c :initarg :out-c :initform nil :accessor out-c)
  27    (in-cf :initarg :in-cf :initform nil :accessor in-cf)
  28    (out-cf :initarg :out-cf :initform nil :accessor out-cf)))
  29
  30 (defclass fftw-r2c-1d-multi (fftw-multi)
  31   ((in-count :initarg :in-count :initform 0 :type fixnum :accessor in-count)
  32    (out-count :initarg :out-count :initform 0 :type fixnum :accessor out-count)
  33    (complex-output :initarg :complex-output :initform t :type boolean :accessor complex-output)
  34    (in-cf :initarg :in-cf :initform nil :accessor in-cf)
  35    (out-cf :initarg :out-cf :initform nil :accessor out-cf)
  36    (out :initarg :out :initform nil :accessor out)))
  37
  38 (defclass fftw-1d-multi (fftw-multi)
  39   ((in-count :initarg :in-count :initform 0 :type fixnum :accessor in-count)
  40    (out-count :initarg :out-count :initform 0 :type fixnum :accessor out-count)
  41    (in-cf :initarg :in-cf :initform nil :accessor in-cf)
  42    (out-cf :initarg :out-cf :initform nil :accessor out-cf)
  43    (out :initarg :out :initform nil :accessor out)))
  44
  45 (defclass fftw-c2r-1d-multi (fftw-multi)
  46   ((plan-even :initarg :plan-even :initform nil :accessor plan-even)
  47    (plan-odd :initarg :plan-odd :initform nil :accessor plan-odd)
  48    (in-count :initarg :in-count :initform 0 :type fixnum :accessor in-count)
  49    (out-count-odd :initarg :out-count-odd :initform 0 :type fixnum :accessor out-count-odd)
  50    (out-count-even :initarg :out-count-even :initform 0 :type fixnum :accessor out-count-even)
  51    (in-cf :initarg :in-cf :initform nil :accessor in-cf)
  52    (out-cf :initarg :out-cf :initform nil :accessor out-cf)
  53    (out-odd :initarg :out-odd :initform nil :accessor out-odd)
  54    (out-even :initarg :out-even :initform nil :accessor out-even)))
  55
  56 (defun normalize-input-range (in start count)
  57   "Returns the start and count of an input range with range checking."
  58   (declare #.*standard-optimize-settings*)
  59   (unless (vectorp in)
  60     (error "Input needs to be a vector."))
  61   (when (or (not (integerp start)) (minusp start))
  62     (error "Start must be an integer, got ~A." start))
  63   (let ((len (length in)))
  64     (declare (fixnum len))
  65     (when (>= start len)
  66       (error "Start (~D) greater than last element (~D)." start (1- len)))
  67
  68     (cond
  69      ((null count)
  70       (setq count (- len start)))
  71      ((or (not (integerp count)) (minusp count))
  72       (error "Count must be positive integer or nil, got ~A." count))
  73      ((> (+ start count) len)
  74       (error "Requesting element (~D) past end of vector (~D)." (+ start count) len)))
  75     (values start count)))
  76
  77 (defun complex-is-even-p (v)
  78   "Returns T if first and last members have zero imaginary component."
  79   ;; Do not use *standard-optimize-settings-here* -- they fail on lispworks 5.1.2
  80   (declare (optimize (speed 3) (space 0)))
  81   (and (zerop (imagpart (aref v 0)))
  82        (zerop (imagpart (aref v (1- (length v)))))))
  83
  84 (defun complex-to-mag-phase (cmplx &optional (elements-complex t))
  85   (declare #.*standard-optimize-settings*)
  86   (let* ((len (length cmplx))
  87          (out-len (if elements-complex len (/ len 2)))
  88          (mag (make-array out-len :element-type 'double-float))
  89          (phase (make-array out-len :element-type 'double-float)))
  90     (declare (fixnum len out-len))
  91     (cond
  92       (elements-complex
  93        (dotimes (i out-len)
  94          (declare (fixnum i))
  95          (let* ((e (aref cmplx i)))
  96            (setf (aref mag i) (abs e))
  97            (setf (aref phase i) (phase e)))))
  98       (t
  99        (dotimes (i out-len)
 100          (declare (fixnum i))
 101          (let* ((base (+ i i))
 102                 (r (aref cmplx base))
 103                 (i (aref cmplx (1+ base))))
 104            (setf (aref mag i) (sqrt (+ (* r r) (* i i))))
 105            (setf (aref phase i) (atan i r))))))
 106     (values mag phase)))
 107
 108 (defun hc-to-mag-phase (hc)
 109   "Turns half-complex into magniture and phase vectors."
 110   (declare #.*standard-optimize-settings*)
 111   (let* ((n (length hc))
 112          (out-n (1+ (floor n 2)))
 113          (mag (make-array out-n :element-type 'double-float))
 114          (phase (make-array out-n :element-type 'double-float)))
 115     (declare (fixnum n out-n))
 116     (dotimes (i out-n)
 117       (cond
 118        ((or (zerop i)
 119             (and (evenp n) (eql i (1- out-n))))
 120         (setf (aref mag i) (aref hc (/ i 2)))
 121         (setf (aref phase i) 0d0))
 122        (t
 123         (let ((re (aref hc i))
 124               (im (aref hc (- n i))))
 125           (setf (aref mag i) (sqrt (+ (* re re) (* im im))))
 126           (setf (aref phase i) (atan im re))))))
 127     (values mag phase)))
 128