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