From c47196e3e91382faa62200487f07c4014daa7065 Mon Sep 17 00:00:00 2001
From: Ivan Ukhov <ivan.ukhov@gmail.com>
Date: Sun, 21 Jun 2015 06:46:14 -0400
Subject: [PATCH] =?UTF-8?q?Add=20an=20implementation=20based=20on=20Brenne?=
 =?UTF-8?q?r=E2=80=99s=20(1976)?=
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

---
 src/lib.rs | 68 ++++++++++++++++++++++++++++++++++++++++++++++++++----
 1 file changed, 63 insertions(+), 5 deletions(-)

diff --git a/src/lib.rs b/src/lib.rs
index 208ba55..e3f8b2b 100644
--- a/src/lib.rs
+++ b/src/lib.rs
@@ -1,11 +1,69 @@
+#![feature(step_by)]
+
+/// A direction of the Fourier transform.
 pub enum Direction {
+    /// From the time domain to the frequency domain.
     Forward,
-    Backward,
+    /// From the frequency domain to the time domain.
+    Inverse,
 }
 
-pub fn transform(data: &mut [f64], _: Direction) {
-    let n = data.len();
-    if n < 2 || n & (n - 1) != 0 {
-        panic!("the data size should be a power of two");
+/// Perform the Fourier transform.
+pub fn transform(data: &mut [f64], direction: Direction) {
+    use std::f64::consts::PI;
+
+    let n = data.len() / 2;
+    let nn = n << 1;
+
+    let isign = match direction {
+        Direction::Forward => 1.0,
+        Direction::Inverse => -1.0,
+    };
+
+    let mut j = 1;
+    for i in (1..nn).step_by(2) {
+        if j > i {
+            data.swap(j - 1, i - 1);
+            data.swap(j, i);
+        }
+        let mut m = n;
+        while m >= 2 && j > m {
+            j -= m;
+            m >>= 1;
+        }
+        j += m;
+    }
+
+    let mut mmax = 2;
+    while nn > mmax {
+        let istep = mmax << 1;
+        let theta = isign * (2.0 * PI / mmax as f64);
+        let wtemp = (0.5 * theta).sin();
+        let wpr = -2.0 * wtemp * wtemp;
+        let wpi = theta.sin();
+        let mut wr = 1.0;
+        let mut wi = 0.0;
+        for m in (1..mmax).step_by(2) {
+            for i in (m..(nn + 1)).step_by(istep) {
+                let j = i + mmax;
+                let tempr = wr * data[j - 1] - wi * data[j];
+                let tempi = wr * data[j] + wi * data[j - 1];
+                data[j - 1] = data[i - 1] - tempr;
+                data[j] = data[i] - tempi;
+                data[i - 1] += tempr;
+                data[i] += tempi;
+            }
+            let wtemp = wr;
+            wr = wr * wpr - wi * wpi + wr;
+            wi = wi * wpr + wtemp * wpi + wi;
+        }
+        mmax = istep;
+    }
+
+    if let Direction::Inverse = direction {
+        let scale = 1.0 / n as f64;
+        for i in 0..(2 * n) {
+            data[i] *= scale;
+        }
     }
 }