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 fry kernel locals math
   5 math.functions math.vectors sequences sequences.private
   6 specialized-arrays ;
   7 SPECIALIZED-ARRAY: double
   8 SPECIALIZED-ARRAY: fftw_complex
   9
  10 IN: fftw
  11
  12 <PRIVATE
  13
  14 : <fftw-array> ( length -- array )
  15     [ fftw_complex heap-size * fftw_malloc &fftw_free ] keep
  16     fftw_complex-array boa ;
  17
  18 : >fftw-array ( seq -- array )
  19     [ length <fftw-array> ] keep over '[
  20         [ >rect 0 1 ] [ _ nth ] bi*
  21         [ set-nth-unsafe ] curry bi-curry@ bi*
  22     ] each-index ;
  23
  24 : fftw-array> ( array -- seq )
  25     [ first2 rect> ] { } map-as ;
  26
  27 :: (fft1d) ( seq sign -- seq' )
  28     seq length :> n
  29     [
  30         n
  31         seq >fftw-array
  32         n <fftw-array> [
  33             sign FFTW_ESTIMATE fftw_plan_dft_1d
  34             [ fftw_execute ] [ fftw_destroy_plan ] bi
  35         ] keep fftw-array>
  36     ] with-destructors ;
  37
  38 PRIVATE>
  39
  40 : fft1d ( seq -- seq' ) FFTW_FORWARD (fft1d) ;
  41
  42 : ifft1d ( seq -- seq' )
  43     [ FFTW_BACKWARD (fft1d) ] [ length v/n ] bi ;
  44
  45 : correlate1d ( x y -- z )
  46     [ fft1d ] [ reverse fft1d ] bi* v* ifft1d ;