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