palmfft 1.1.2

Palm-sized Faster Fourier Transform
Documentation 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69

use crate::{
    algorithms::cooleytukey::CooleyTukey,
    complex::Complex,
    math::{good_size, sincos_2pibyn}
};

pub struct Bluestein {
    n: usize,
    n2: usize,
    plan: CooleyTukey,
    bk: Vec<Complex>,
    bkf: Vec<Complex>,
}

impl Bluestein {
    pub fn new(length: usize) -> Bluestein {
        let n = length;
        let n2 = if length == 0 { 0 } else { good_size(n * 2 - 1) };
        let mut plan = Bluestein {
            n, n2, plan: CooleyTukey::new(n2),
            bk: vec![Complex::new(0.0, 0.0); n],
            bkf: vec![Complex::new(0.0, 0.0); n2]
        };

        if plan.n < 2 { return plan; }

        let mut tmp = vec![Complex::new(0.0, 0.0); n * 2];
        sincos_2pibyn(n * 2, &mut tmp);
        plan.bk[0] = Complex::new(1.0, 0.0);

        let mut coeff = 0;
        for m in 1..n {
            coeff += 2 * m - 1;
            if coeff >= 2 * n { coeff -= 2 * n; }
            plan.bk[m] = tmp[coeff];
        }

        let xn2 = 1.0 / (n2 as f64);
        plan.bkf[0] = plan.bk[0] * xn2;

        for m in 1..n {
            let norm = plan.bk[m] * xn2;
            (plan.bkf[m], plan.bkf[n2 - m]) = (norm, norm);
        }

        for m in n..=(n2 - n) { plan.bkf[m] = Complex::new(0.0, 0.0); }
        plan.plan.forward(&mut plan.bkf, 1.0);
        return plan;
    }

    pub fn forward(&self, data: &mut [Complex], fct: f64) { self.fft(data, fct, -1); }
    pub fn backward(&self, data: &mut [Complex], fct: f64) { self.fft(data, fct, 1); }

    pub fn fft(&self, data: &mut [Complex], fct: f64, sign: i32) {
        if self.n != data.len() { panic!("Provided buffer({}) not compatible with Bluestein({}).", data.len(), self.n); }
        if self.n < 2 { return; }
        let mut akf = vec![Complex::new(0.0, 0.0); self.n2];

        for m in 0..self.n { akf[m] = data[m] * if sign > 0 { self.bk[m] } else { self.bk[m].conj() }; }
        for m in self.n..self.n2 { akf[m] = Complex::new(0.0, 0.0); }

        self.plan.forward(&mut akf, fct);
        for m in 0..self.n2 { akf[m] *= if sign > 0 { self.bkf[m].conj() } else { self.bkf[m] }; }
        self.plan.backward(&mut akf, fct);
        for m in 0..self.n { data[m] = akf[m] * if sign > 0 { self.bk[m] } else { self.bk[m].conj() }; }
    }

    pub fn len(&self) -> usize { self.n }
}