rdsp-experiments/src/fft/rader2.rs

// Implementation of raders's fft for prime sized ffts

/*
use std::{f32::consts::PI, ops::Deref};

use super::mixed_radix;
use crate::{
    complex::Complex32,
    fft::{DFT, FFTDirection, create_fft, dft::NaiveDFT, is_prime, windows},
};

pub struct Rader2FFT {
    input_buffer: Box<[Complex32]>,
    output_buffer: Box<[Complex32]>,

    size: usize,
    sub_size: usize,

    // Fourrier transform of the exponential terms
    convolution_operand: Box<[Complex32]>,
    convolution_fft: Box<dyn DFT>, // TODO: Use fft
    permutation: Box<[usize]>,
}

impl DFT for Rader2FFT {
    fn create(size: usize, direction: FFTDirection) -> Self
    where
        Self: Sized,
    {
        let g = compute_prime_primitive_root(size);
        let permutation: Box<[usize]> = (0..(size - 1)).map(|i| exp_mod(g, i + 1, size)).collect();
        let sub_size = next_pow2((2 * size - 4) - 1);
        Rader2FFT {
            input_buffer: vec![Complex32::zero(); size].into_boxed_slice(),
            output_buffer: vec![Complex32::zero(); sub_size].into_boxed_slice(),

            size,
            sub_size,

            convolution_operand: compute_convolution_operand(size, sub_size, &permutation),
            //convolution_fft: create_fft(next_pow2((2 * size - 4) - 1)),
            convolution_fft: Box::new(NaiveDFT::create(sub_size)),
            permutation,
        }
    }

    fn execute(&mut self, window: fn(f32) -> f32) {
        self.convolution_fft.get_input()[0] = self.input_buffer[self.permutation[self.size - 2]];

        for i in 0..(self.sub_size - self.size + 1) {
            self.convolution_fft.get_input()[i + 1] = Complex32::zero();
        }
        for i in 1..(self.size - 1) {
            let k = self.permutation[self.size - 1 - i - 1];
            self.convolution_fft.get_input()[i + self.sub_size - self.size + 1] =
                self.input_buffer[k] * window(k as f32 / self.size as f32)
        }

        self.convolution_fft.execute(windows::rectanguar);

        // Use output buffer as staging buffer
        for i in 0..(self.sub_size) {
            self.output_buffer[i] =
                self.convolution_fft.get_output()[i] * self.convolution_operand[i];
        }

        for i in 0..(self.sub_size) {
            self.convolution_fft.get_input()[i] = self.output_buffer[i];
        }
        /*
        self.convolution_fft.get_input()[0] =
            self.convolution_fft.get_input()[0] + self.input_buffer[0] * window(0.);
        */

        // Compute ifft to obtain convolution
        self.convolution_fft.execute(window);

        for i in 0..(self.size - 1) {
            self.output_buffer[self.permutation[i]] =
                self.convolution_fft.get_output()[i] / self.sub_size as f32;
        }

        self.output_buffer[0] = self
            .input_buffer
            .iter()
            .copied()
            .enumerate()
            .map(|(i, x)| x * window(i as f32 / self.size as f32))
            .sum();
    }

    fn get_input(&mut self) -> &mut [Complex32] {
        &mut self.input_buffer
    }

    fn get_output(&self) -> &[Complex32] {
        &self.output_buffer
    }
}

pub fn compute_convolution_operand(
    n: usize,
    sub_size: usize,
    permutation: &[usize],
) -> Box<[Complex32]> {
    //let mut fft = create_fft(sub_size);
    let mut fft = NaiveDFT::create(sub_size);

    fft.get_input().iter_mut().enumerate().for_each(|(i, x)| {
        *x = Complex32::cexp(-2. * PI * (permutation[i % (n - 1)] as f32) / (n as f32))
    });
    fft.execute(windows::rectanguar);
    fft.get_output().iter().copied().collect()
}

pub fn compute_prime_primitive_root(n: usize) -> usize {
    assert!(is_prime(n));

    let phi = n - 1; // Euler's totient for n prime

    // Test all candidates
    for i in 1..(n + 1) {
        // Find multiplicative order of i
        let mut val = i;
        let mut order = 1;
        for j in 0..n {
            if val == 1 {
                break;
            }
            val = (val * i) % n;
            order += 1;
        }

        if order == phi {
            return i;
        }
    }

    unreachable!("Prime must have primitive root");
}

pub fn exp_mod(mut n: usize, mut exp: usize, m: usize) -> usize {
    if m == 1 {
        return 0;
    }

    n %= m;
    let mut r = 1;
    while exp > 0 {
        if exp % 2 == 1 {
            r = (r * n) % m;
        }

        n = (n * n) % m;
        exp >>= 1;
    }

    r
}

pub fn next_pow2(mut n: usize) -> usize {
    let mut pow = 0;
    while n > 0 {
        n >>= 1;
        pow += 1;
    }
    1 << pow
}
*/