extra/math/transforms/fft/fft.factor

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