math.transforms.fft: simple, fairly slow, fft.

extra/math/transforms/fft/authors.txt

@@ -0,0 +1 @@
+John Benediktsson

extra/math/transforms/fft/fft-docs.factor

@@ -0,0 +1,10 @@
+USING: help.markup help.syntax sequences ;
+IN: math.transforms.fft
+
+HELP: fft
+{ $values { "seq" sequence } { "seq'" sequence } }
+{ $description "Fast Fourier transform function." } ;
+
+HELP: ifft
+{ $values { "seq" sequence } { "seq'" sequence } }
+{ $description "Inverse fast Fourier transform function." } ;

extra/math/transforms/fft/fft-tests.factor

@@ -0,0 +1,26 @@
+USING: math.vectors tools.test ;
+IN: math.transforms.fft
+
+! even lengths
+
+{ t } [
+    { C{ 10 0 } C{ -2 2 } C{ -2 0 } C{ -2 -2 } }
+    { 1 2 3 4 } fft 1e-12 v~
+] unit-test
+
+{ t } [
+    { C{ 2+1/2 0 } C{ -1/2 -1/2 } C{ -1/2 0 } C{ -1/2 1/2 } }
+    { 1 2 3 4 } ifft 1e-12 v~
+] unit-test
+
+! odd lengths
+
+{ t } [
+    { C{ 5 0 } C{ -1 0 } C{ -1 0 } }
+    { 1 2 2 } fft 1e-12 v~
+] unit-test
+
+{ t } [
+    { C{ 1+2/3 0 } C{ -1/3 0 } C{ -1/3 0 } }
+    { 1 2 2 } ifft 1e-12 v~
+] unit-test

extra/math/transforms/fft/fft.factor

@@ -0,0 +1,43 @@
+! Copyright (c) 2012 John Benediktsson
+! See http://factorcode.org/license.txt for BSD license.
+USING: arrays kernel locals math math.constants math.functions
+math.vectors sequences sequences.extras sequences.private ;
+IN: math.transforms.fft
+
+<PRIVATE
+
+:: (slow-fft) ( seq inverse? -- seq' )
+    seq length :> N
+    inverse? 1 -1 ? 2pi * i* N / :> O
+    N iota [| k |
+        0 seq [ O k * * e^ * + ] each-index
+        inverse? [ N / ] when
+    ] map ; inline
+
+:: (fft) ( seq inverse? -- seq' )
+    seq length :> N
+    N 1 = [ seq ] [
+        inverse? 1 -1 ? 2pi * i* N / :> O
+        N 0 <array> :> X
+        N 2/ :> M
+        seq even-indices inverse? (fft)
+        seq odd-indices inverse? (fft)
+        [ [ [ O * e^ * + inverse? [ 2 / ] when ] [ X set-nth-unsafe ] bi ] 2each-index ]
+        [ [ [ O * e^ * - inverse? [ 2 / ] when ] [ M + X set-nth-unsafe ] bi ] 2each-index ]
+        2bi
+        X
+    ] if ; inline recursive
+
+PRIVATE>
+
+ERROR: not-enough-data ;
+
+: fft ( seq -- seq' )
+    [ not-enough-data ] [
+        f over length even? [ (fft) ] [ (slow-fft) ] if
+    ] if-empty ;
+
+: ifft ( seq -- seq' )
+    [ not-enough-data ] [
+        t over length even? [ (fft) ] [ (slow-fft) ] if
+    ] if-empty ;

extra/math/transforms/fft/summary.txt

@@ -0,0 +1 @@
+Fast fourier transform