Simple unoptimized FFT implementation

fft.rs

use num::complex::Complex;
use std::f64;

fn ditfft2(x: &[Complex<f64>], n: usize, s: usize) -> Vec<Complex<f64>> {
    if n == 1 {
        return vec![x[0]];
    }

    let mut out: Vec<_> = ditfft2(x, n / 2, s * 2)
        .into_iter()
        .chain(ditfft2(&x[s..], n / 2, s * 2))
        .collect();

    for k in 0..n / 2 {
        let t = out[k];
        let e = Complex::new(0.0, -2.0 * f64::consts::PI * k as f64 / n as f64).exp();
        out[k] = t + e * out[k + n / 2];
        out[k + n / 2] = t - e * out[k + n / 2];
    }

    out
}

fn main() {
    let x = vec![
        Complex::new(1.0, 0.0),
        Complex::new(2.0, 0.0),
        Complex::new(3.0, 0.0),
        Complex::new(4.0, 0.0),
        Complex::new(3.0, 0.0),
        Complex::new(2.0, 0.0),
        Complex::new(0.0, 0.0),
        Complex::new(0.0, 0.0),
    ];
    
    let mut y = ditfft2(&x, 8, 1);
    
    println!("{:?}", y);
    
    // three different approaches to compute the inverse (https://en.wikipedia.org/wiki/Discrete_Fourier_transform#Expressing_the_inverse_DFT_in_terms_of_the_DFT)
    // 1. swap imaginary and real parts in input and output
    // 2. conjugate input and output
    // 3. reverse all but the first elements
    
    /*for v in &mut y {
        // *v = Complex::new(v.im, v.re);
        *v = v.conj();
    }*/
    
    y[1..].reverse();
    
    let mut z = ditfft2(&y, 8, 1);
    
    /*for v in &mut z {
        // *v = Complex::new(v.im, v.re);
        *v = v.conj();
        
        *v /= x.len() as f64;
    }*/
    
    for v in &mut z {
        *v /= x.len() as f64;
    }
    
    println!("{:?}", z);
}