transform.lisp

   1 ;;;; -*- Mode: Lisp; Syntax: ANSI-Common-Lisp; Base: 10 -*-
   2 ;;;; *************************************************************************
   3 ;;;; FILE IDENTIFICATION
   4 ;;;;
   5 ;;;; Name:          transform.lisp
   6 ;;;; Purpose:       Transformation functions
   7 ;;;; Programmer:    Kevin M. Rosenberg
   8 ;;;; Date Started:  March 2009
   9 ;;;;
  10 ;;;; This file and CL-FFTW3 are Copyright (c) 2009-2011 by Kevin M. Rosenberg
  11 ;;;;
  12 ;;;; FFTW3 users are granted the rights to distribute and use this software
  13 ;;;; as governed by the terms of the Lisp Lesser GNU Public License
  14 ;;;; (http://opensource.franz.com/preamble.html), also known as the LLGPL.
  15 ;;;; *************************************************************************
  16
  17 (in-package #:fftw3)
  18
  19 (defconstant +size-double+ 8)
  20
  21 (defun fftw-r2hc-1d (in &key (start 0) (count nil) (normalize t))
  22   (declare #.*standard-optimize-settings*)
  23   (multiple-value-bind (start count) (normalize-input-range in start count)
  24     (declare (fixnum start count))
  25
  26     (let ((in-pos start)
  27           (out-c (make-array count :element-type 'double-float))
  28           (in-cf (fftw-malloc (* count +size-double+)))
  29           (out-cf (fftw-malloc (* count +size-double+))))
  30       (declare (fixnum in-pos))
  31       (dotimes (i count)
  32         (declare (fixnum i))
  33         (setf (cffi:mem-aref in-cf :double i) (coerce (aref in in-pos) 'double-float))
  34         (incf in-pos))
  35
  36       (let ((p (fftw-plan-r2r-1d count in-cf out-cf +fftw-r2hc+
  37                                  (logior +fftw-estimate+ +fftw-destroy-input+))))
  38         (fftw-execute p)
  39         (fftw-destroy-plan p))
  40
  41       (if normalize
  42           (let ((factor (coerce (/ 1 count) 'double-float)))
  43             (declare (double-float factor))
  44             (dotimes (i count)
  45               (declare (fixnum i))
  46               (setf (aref out-c i) (* factor (cffi:mem-aref out-cf :double i)))))
  47           (dotimes (i count)
  48             (declare (fixnum i))
  49             (setf (aref out-c i) (cffi:mem-aref out-cf :double i))))
  50
  51       (fftw-free in-cf)
  52       (fftw-free out-cf)
  53
  54       out-c)))
  55
  56 (defun make-fftw-r2r-1d-multi (count type)
  57   (declare #.*standard-optimize-settings*)
  58   (declare (fixnum count type))
  59   (let* ((out-c (make-array count :element-type 'double-float))
  60          (in-cf (fftw-malloc (* count +size-double+)))
  61          (out-cf (fftw-malloc (* count +size-double+)))
  62          (plan (fftw-plan-r2r-1d count in-cf out-cf type
  63                                  (logior +fftw-estimate+ +fftw-destroy-input+))))
  64     (make-instance 'fftw-r2r-1d-multi
  65                    :in-count count :out-c out-c :plan plan
  66                    :in-cf in-cf :out-cf out-cf)))
  67
  68 (defun fftw-r2r-1d-multi (multi in &key (start 0) (count nil) (normalize t))
  69   (declare #.*standard-optimize-settings*)
  70   (multiple-value-bind (start count) (normalize-input-range in start count)
  71     (declare (fixnum start count))
  72     (unless (equal (in-count multi) count)
  73       (error "Different plan and vector lengths."))
  74
  75     (let ((in-cf (in-cf multi))
  76           (out-cf (out-cf multi))
  77           (out-c (out-c multi))
  78           (in-pos start))
  79       (declare (fixnum in-pos))
  80
  81       (dotimes (i count)
  82         (declare (fixnum i))
  83         (setf (cffi:mem-aref in-cf :double i) (coerce (aref in in-pos) 'double-float))
  84         (incf in-pos))
  85
  86       (fftw-execute (plan multi))
  87       (if normalize
  88           (let ((factor (coerce (/ 1 count) 'double-float)))
  89             (declare (double-float factor))
  90             (dotimes (i count)
  91               (declare (fixnum i))
  92               (setf (aref out-c i) (* factor (cffi:mem-aref out-cf :double i)))))
  93           (dotimes (i count)
  94             (declare (fixnum i))
  95             (setf (aref out-c i) (cffi:mem-aref out-cf :double i))))
  96
  97       out-c)))
  98
  99
 100 (defun fftw-hc2r-1d (in)
 101   (declare #.*standard-optimize-settings*)
 102   (let* ((count (length in))
 103          (out-c (make-array count :element-type 'double-float))
 104          (in-cf (fftw-malloc (* count +size-double+)))
 105          (out-cf (fftw-malloc (* count +size-double+))))
 106     (declare (fixnum count))
 107
 108     (dotimes (i count)
 109       (declare (fixnum i))
 110       (setf (cffi:mem-aref in-cf :double i) (aref in i)))
 111
 112     (let ((p (fftw-plan-r2r-1d count in-cf out-cf +fftw-hc2r+
 113                                (logior +fftw-estimate+ +fftw-destroy-input+))))
 114       (fftw-execute p)
 115       (fftw-destroy-plan p))
 116
 117     (dotimes (i count)
 118       (declare (fixnum i))
 119       (setf (aref out-c i) (cffi:mem-aref out-cf :double i)))
 120     (fftw-free in-cf)
 121     (fftw-free out-cf)
 122     out-c))
 123
 124 (defun fftw-r2c-1d (in &key (start 0) (count nil) (normalize t))
 125   (declare #.*standard-optimize-settings*)
 126   (multiple-value-bind (start count) (normalize-input-range in start count)
 127     (declare (fixnum start count))
 128
 129     (let* ((out-n (1+ (floor count 2)))
 130            (in-cf (fftw-malloc (* count +size-double+)))
 131            (out-cf (fftw-malloc (* 2 out-n +size-double+)))
 132            (in-pos start))
 133       (declare (fixnum in-pos out-n))
 134
 135       (dotimes (i count)
 136         (declare (fixnum i))
 137         (setf (cffi:mem-aref in-cf :double i) (coerce (aref in in-pos) 'double-float))
 138         (incf in-pos))
 139
 140       (let ((p (fftw-plan-dft-r2c-1d count in-cf out-cf
 141                                      (logior +fftw-estimate+ +fftw-destroy-input+))))
 142         (fftw-execute p)
 143         (fftw-destroy-plan p))
 144
 145       (let ((out (make-array out-n :element-type '(or rational complex)))
 146             (opos 0))
 147         (declare (fixnum opos))
 148         (if normalize
 149             (let ((factor (coerce (/ 1 count) 'double-float)))
 150               (declare (double-float factor))
 151               (dotimes (i out-n)
 152                 (declare (fixnum i))
 153                 (let ((impos (1+ opos)))
 154                   (declare (fixnum impos))
 155                   (setf (aref out i) (complex (* (cffi:mem-aref out-cf :double opos) factor)
 156                                             (* (cffi:mem-aref out-cf :double impos) factor)))
 157                   (incf opos 2))))
 158             (dotimes (i out-n)
 159               (declare (fixnum i))
 160               (let ((impos (1+ opos)))
 161                 (declare (fixnum impos))
 162                 (setf (aref out i) (complex (cffi:mem-aref out-cf :double opos)
 163                                             (cffi:mem-aref out-cf :double impos)))
 164                 (incf opos 2))))
 165         (fftw-free in-cf)
 166         (fftw-free out-cf)
 167         out))))
 168
 169 (defun make-fftw-r2c-1d-multi (count &key (complex-output t))
 170   (declare #.*standard-optimize-settings*)
 171   (declare (fixnum count))
 172   (let* ((out-n (1+ (floor count 2)))
 173          (in-cf (fftw-malloc (* count +size-double+)))
 174          (out-cf (fftw-malloc (* 2 out-n +size-double+)))
 175          (plan (fftw-plan-dft-r2c-1d count in-cf out-cf
 176                                      (logior +fftw-estimate+ +fftw-destroy-input+))))
 177     (declare (fixnum out-n))
 178
 179     (unless complex-output
 180       (incf out-n out-n))
 181
 182     (make-instance 'fftw-r2c-1d-multi
 183                    :in-count count :out-count out-n :plan plan
 184                    :in-cf in-cf
 185                    :out-cf out-cf
 186                    :complex-output complex-output
 187                    :out
 188                    (if complex-output
 189                      (make-array out-n :element-type '(or double-float complex) :initial-element 0d0)
 190                      (make-array out-n :element-type 'double-float)))))
 191
 192 (defgeneric destroy-fftw-multi (multi))
 193 (defmethod destroy-fftw-multi ((multi fftw-multi))
 194   (declare #.*standard-optimize-settings*)
 195   (when (plan multi)
 196     (fftw-destroy-plan (plan multi))
 197     (setf (plan multi) nil))
 198   (dolist (slot '(in-cf out-cf))
 199     (when (slot-value multi slot)
 200       (fftw-free (slot-value multi slot))
 201       (setf (slot-value multi slot) nil)))
 202   (dolist (slot '(plan-even plan-odd))
 203     (when (and (slot-exists-p multi slot)
 204                (slot-value multi slot))
 205       (fftw-destroy-plan (slot-value multi slot))
 206       (setf (slot-value multi slot) nil))))
 207
 208
 209 (defun fftw-r2c-1d-multi (multi in &key (start 0) (normalize t))
 210   (declare #.*standard-optimize-settings*)
 211   (multiple-value-bind (start count) (normalize-input-range in start (in-count multi))
 212     (declare (fixnum start count))
 213
 214     (unless (equal (length in) (in-count multi))
 215       (error "In count of multi plan doesn't equal length of in vector."))
 216
 217     (let ((in-cf (in-cf multi))
 218           (out-cf (out-cf multi))
 219           (out (out multi))
 220           (in-pos start)
 221           (opos 0)
 222           (in-count (in-count multi))
 223           (out-count (out-count multi)))
 224       (declare (fixnum in-pos opos in-count out-count))
 225
 226       (dotimes (i in-count)
 227         (declare (fixnum i))
 228         (setf (cffi:mem-aref in-cf :double i) (coerce (aref in in-pos) 'double-float))
 229         (incf in-pos))
 230
 231       (fftw-execute (plan multi))
 232
 233       (setf opos 0)
 234       (cond
 235        ((complex-output multi)
 236         (if normalize
 237             (let ((factor (coerce (/ 1 count) 'double-float)))
 238               (declare (double-float factor))
 239               (dotimes (i out-count)
 240                 (declare (fixnum i))
 241                 (let ((impos (1+ opos)))
 242                   (declare (fixnum impos))
 243                   (setf (aref out i) (complex (* (cffi:mem-aref out-cf :double opos) factor)
 244                                               (* (cffi:mem-aref out-cf :double impos) factor)))
 245                   (incf opos 2))))
 246             (dotimes (i out-count)
 247               (declare (fixnum i))
 248               (let ((impos (1+ opos)))
 249                 (declare (fixnum impos))
 250                 (setf (aref out i) (complex (cffi:mem-aref out-cf :double opos)
 251                                             (cffi:mem-aref out-cf :double impos)))
 252                 (incf opos 2)))))
 253        (t
 254         (when normalize
 255           (let ((o2 (* out-count 2))
 256                 (factor (coerce (/ 1 count) 'double-float)))
 257             (declare (fixnum o2)
 258                      (double-float factor))
 259             (dotimes (i o2)
 260               (declare (fixnum i))
 261               (setf (aref out i) (* factor (aref out i))))))))
 262       out)))
 263
 264 (defun fftw-c2r-1d (in)
 265   (declare #.*standard-optimize-settings*)
 266   (let* ((count (length in))
 267          (out-n (if (complex-is-even-p in)
 268                   (* 2 (1- count))
 269                   (1+ (* 2 (1- count)))))
 270          (out-c (make-array out-n :element-type 'double-float))
 271          (in-cf (fftw-malloc (* 2 count +size-double+)))
 272          (out-cf (fftw-malloc (* 2 count +size-double+)))
 273          (pos 0))
 274     (declare (fixnum count out-n pos))
 275
 276     (setq pos 0)
 277     (dotimes (i count)
 278       (declare (fixnum i))
 279       (let ((c (aref in i)))
 280         (declare (complex c))
 281         (setf (cffi:mem-aref in-cf :double pos) (coerce (realpart c) 'double-float))
 282         (incf pos)
 283         (setf (cffi:mem-aref in-cf :double pos) (coerce (imagpart c) 'double-float))
 284         (incf pos)))
 285
 286     (let ((p (fftw-plan-dft-c2r-1d out-n in-cf out-cf
 287                                    (logior +fftw-estimate+ +fftw-destroy-input+))))
 288       (fftw-execute p)
 289       (fftw-destroy-plan p))
 290     (dotimes (i out-n)
 291       (declare (fixnum i))
 292       (setf (aref out-c i) (cffi:mem-aref out-cf :double i)))
 293
 294     (fftw-free in-cf)
 295     (fftw-free out-cf)
 296     out-c))
 297
 298 (defun make-fftw-c2r-1d-multi (count)
 299   (declare #.*standard-optimize-settings*)
 300   (declare (fixnum count))
 301   (let* ((out-count-even (* 2 (1- count)))
 302          (out-count-odd (1+ (* 2 (1- count))))
 303          (in-cf (fftw-malloc (* 2 count +size-double+)))
 304          (out-cf (fftw-malloc (* 2 (max out-count-odd out-count-even) +size-double+))) ;; use larger size
 305          (plan-odd (fftw-plan-dft-c2r-1d out-count-odd in-cf out-cf
 306                                          (logior +fftw-estimate+ +fftw-destroy-input+)))
 307          (plan-even (fftw-plan-dft-c2r-1d out-count-even in-cf out-cf +fftw-measure+)))
 308     (declare (fixnum out-count-even out-count-odd))
 309
 310     (make-instance 'fftw-c2r-1d-multi
 311                    :in-count count
 312                    :plan-even plan-even
 313                    :plan-odd plan-odd
 314                    :in-cf in-cf
 315                    :out-cf out-cf
 316                    :out-count-even out-count-even
 317                    :out-count-odd out-count-odd
 318                    :out-even (make-array out-count-even :element-type 'double-float)
 319                    :out-odd (make-array out-count-odd :element-type 'double-float))))
 320
 321 (defun fftw-c2r-1d-multi (multi in)
 322   (declare #.*standard-optimize-settings*)
 323   (let* ((is-even (complex-is-even-p in))
 324          (out (if is-even (out-even multi) (out-odd multi)))
 325          (out-count (if is-even (out-count-even multi) (out-count-odd multi)))
 326          (plan (if is-even (plan-even multi) (plan-odd multi)))
 327          (count (in-count multi))
 328          (in-cf (in-cf multi))
 329          (out-cf (out-cf multi))
 330          (pos 0))
 331     (declare (fixnum count out-count pos))
 332
 333     (setq pos 0)
 334     (dotimes (i count)
 335       (declare (fixnum i))
 336       (let ((c (aref in i)))
 337         (setf (cffi:mem-aref in-cf :double pos) (coerce (realpart c) 'double-float))
 338         (incf pos)
 339         (setf (cffi:mem-aref in-cf :double pos) (coerce (imagpart c) 'double-float))
 340         (incf pos)))
 341
 342     (fftw-execute plan)
 343
 344     (dotimes (i out-count)
 345       (declare (fixnum i))
 346       (setf (aref out i) (cffi:mem-aref out-cf :double i)))
 347
 348     out))
 349
 350 (defun fftw-c-1d (in &key (direction +fftw-forward+)
 351                 (start 0) (count nil) (normalize t))
 352   (declare #.*standard-optimize-settings*)
 353   (multiple-value-bind (start count) (normalize-input-range in start count)
 354     (declare (fixnum start count))
 355
 356     (let* ((out-n count)
 357            (out (make-array out-n :element-type '(or rational complex)))
 358            (in-cf (fftw-malloc (* 2 count +size-double+)))
 359            (out-cf (fftw-malloc (* 2 out-n +size-double+)))
 360            (opos 0)
 361            (ipos 0))
 362       (declare (fixnum out-n ipos opos))
 363
 364       (setf opos 0)
 365       (setf ipos start)
 366       (dotimes (i count)
 367         (declare (fixnum i))
 368         (let ((c (aref in ipos)))
 369           (incf ipos)
 370           (setf (cffi:mem-aref in-cf :double opos) (coerce (realpart c) 'double-float))
 371           (incf opos)
 372           (setf (cffi:mem-aref in-cf :double opos) (coerce (imagpart c) 'double-float))
 373           (incf opos)))
 374
 375       (let ((p (fftw-plan-dft-1d count in-cf out-cf direction
 376                                  (logior +fftw-estimate+ +fftw-destroy-input+))))
 377         (fftw-execute p)
 378         (fftw-destroy-plan p))
 379
 380       (setf opos 0)
 381       (if normalize
 382           (let ((factor (coerce (/ 1 count) 'double-float)))
 383             (declare (double-float factor))
 384             (dotimes (i out-n)
 385               (declare (fixnum i))
 386               (let ((impos (1+ opos)))
 387                 (declare (fixnum impos))
 388                 (setf (aref out i) (complex (* factor (cffi:mem-aref out-cf :double opos))
 389                                             (* factor (cffi:mem-aref out-cf :double impos)))))
 390               (incf opos 2)))
 391           (dotimes (i out-n)
 392               (declare (fixnum i))
 393               (let ((impos (1+ opos)))
 394                 (declare (fixnum impos))
 395                 (setf (aref out i) (complex (cffi:mem-aref out-cf :double opos)
 396                                             (cffi:mem-aref out-cf :double impos))))
 397               (incf opos 2)))
 398
 399       (fftw-free in-cf)
 400       (fftw-free out-cf)
 401       out)))