Add an implementation based on Brenner’s (1976)

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;
+        }
     }
 }