extra/fftw/fftw.factor

   1 ! Copyright (c) 2014 John Benediktsson
   2 ! See http://factorcode.org/license.txt for BSD license.
   3
   4 USING: alien.c-types destructors fftw.ffi kernel math
   5 math.vectors sequences sequences.private specialized-arrays ;
   6 SPECIALIZED-ARRAY: double
   7 SPECIALIZED-ARRAY: fftw_complex
   8
   9 IN: fftw
  10
  11 <PRIVATE
  12
  13 : <fftw-array> ( length -- array )
  14     [ fftw_complex heap-size * fftw_malloc &fftw_free ] keep
  15     fftw_complex-array boa ;
  16
  17 : >fftw-array ( seq -- array )
  18     [ length <fftw-array> ] keep over '[
  19         [ >rect 0 1 ] [ _ nth ] bi*
  20         [ set-nth-unsafe ] curry bi-curry@ bi*
  21     ] each-index ;
  22
  23 : fftw-array> ( array -- seq )
  24     [ first2 rect> ] { } map-as ;
  25
  26 :: (fft1d) ( seq sign -- seq' )
  27     seq length :> n
  28     [
  29         n
  30         seq >fftw-array
  31         n <fftw-array> [
  32             sign FFTW_ESTIMATE fftw_plan_dft_1d
  33             [ fftw_execute ] [ fftw_destroy_plan ] bi
  34         ] keep fftw-array>
  35     ] with-destructors ;
  36
  37 PRIVATE>
  38
  39 : fft1d ( seq -- seq' ) FFTW_FORWARD (fft1d) ;
  40
  41 : ifft1d ( seq -- seq' )
  42     [ FFTW_BACKWARD (fft1d) ] [ length v/n ] bi ;
  43
  44 : correlate1d ( x y -- z )
  45     [ fft1d ] [ reverse fft1d ] bi* v* ifft1d ;