From d46617e323c649e127f28546b404ddbcf10d3e1a Mon Sep 17 00:00:00 2001
From: John Benediktsson <mrjbq7@gmail.com>
Date: Wed, 8 Aug 2012 20:28:14 -0700
Subject: [PATCH] math.transforms.fft: simple, fairly slow, fft.

---
 extra/math/transforms/fft/authors.txt      |  1 +
 extra/math/transforms/fft/fft-docs.factor  | 10 +++++
 extra/math/transforms/fft/fft-tests.factor | 26 +++++++++++++
 extra/math/transforms/fft/fft.factor       | 43 ++++++++++++++++++++++
 extra/math/transforms/fft/summary.txt      |  1 +
 5 files changed, 81 insertions(+)
 create mode 100644 extra/math/transforms/fft/authors.txt
 create mode 100644 extra/math/transforms/fft/fft-docs.factor
 create mode 100644 extra/math/transforms/fft/fft-tests.factor
 create mode 100644 extra/math/transforms/fft/fft.factor
 create mode 100644 extra/math/transforms/fft/summary.txt

diff --git a/extra/math/transforms/fft/authors.txt b/extra/math/transforms/fft/authors.txt
new file mode 100644
index 0000000000..e091bb8164
--- /dev/null
+++ b/extra/math/transforms/fft/authors.txt
@@ -0,0 +1 @@
+John Benediktsson
diff --git a/extra/math/transforms/fft/fft-docs.factor b/extra/math/transforms/fft/fft-docs.factor
new file mode 100644
index 0000000000..f1cabf7b60
--- /dev/null
+++ b/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." } ;
diff --git a/extra/math/transforms/fft/fft-tests.factor b/extra/math/transforms/fft/fft-tests.factor
new file mode 100644
index 0000000000..fdfeab5f7a
--- /dev/null
+++ b/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
diff --git a/extra/math/transforms/fft/fft.factor b/extra/math/transforms/fft/fft.factor
new file mode 100644
index 0000000000..16c202bd3b
--- /dev/null
+++ b/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 ;
diff --git a/extra/math/transforms/fft/summary.txt b/extra/math/transforms/fft/summary.txt
new file mode 100644
index 0000000000..3d71dfa49a
--- /dev/null
+++ b/extra/math/transforms/fft/summary.txt
@@ -0,0 +1 @@
+Fast fourier transform