extra/math/transforms/fft/fft.factor

   1 ! Copyright (c) 2012 John Benediktsson
   2 ! See http://factorcode.org/license.txt for BSD license.
   3 USING: arrays kernel locals math math.constants math.functions
   4 math.vectors sequences sequences.extras sequences.private ;
   5 IN: math.transforms.fft
   6
   7 <PRIVATE
   8
   9 :: (slow-fft) ( seq inverse? -- seq' )
  10     seq length :> N
  11     inverse? 1 -1 ? 2pi * i* N / :> O
  12     N iota [| k |
  13         0 seq [ O k * * e^ * + ] each-index
  14         inverse? [ N / ] when
  15     ] map ; inline
  16
  17 :: (fft) ( seq inverse? -- seq' )
  18     seq length :> N
  19     N 1 = [ seq ] [
  20         inverse? 1 -1 ? 2pi * i* N / :> O
  21         N 0 <array> :> X
  22         N 2/ :> M
  23         seq even-indices inverse? (fft)
  24         seq odd-indices inverse? (fft)
  25         [ [ [ O * e^ * + inverse? [ 2 / ] when ] [ X set-nth-unsafe ] bi ] 2each-index ]
  26         [ [ [ O * e^ * - inverse? [ 2 / ] when ] [ M + X set-nth-unsafe ] bi ] 2each-index ]
  27         2bi
  28         X
  29     ] if ; inline recursive
  30
  31 PRIVATE>
  32
  33 ERROR: not-enough-data ;
  34
  35 : fft ( seq -- seq' )
  36     [ not-enough-data ] [
  37         f over length even? [ (fft) ] [ (slow-fft) ] if
  38     ] if-empty ;
  39
  40 : ifft ( seq -- seq' )
  41     [ not-enough-data ] [
  42         t over length even? [ (fft) ] [ (slow-fft) ] if
  43     ] if-empty ;
  44
  45 : correlate ( x y -- z )
  46     [ fft ] [ reverse fft ] bi* v* ifft ;