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
11 ! Discrete Fourier Transform
12 :: (slow-fft) ( seq inverse? -- seq' )
14 inverse? 1 -1 ? 2pi * N / N iota n*v :> omega
16 0 seq omega [ k * cis * + ] 2each
20 ! Cooley–Tukey Algorithm
21 :: (fast-fft) ( seq inverse? -- 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 ]
31 ] if ; inline recursive
33 : (fft) ( seq inverse? -- seq' )
34 over length power-of-2?
35 [ (fast-fft) ] [ (slow-fft) ] if ; inline
39 ERROR: not-enough-data ;
42 [ not-enough-data ] [ f (fft) ] if-empty ;
44 : ifft ( seq -- seq' )
45 [ not-enough-data ] [ t (fft) ] if-empty ;
47 : correlate ( x y -- z )
48 [ fft ] [ reverse fft ] bi* v* ifft ;