Mixed radix, begining rader

src/fft/rader.rs

@@ -0,0 +1,160 @@
+// Implementation of raders's fft for prime sized ffts
+
+use std::f32::consts::PI;
+
+use crate::{complex::Complex32, fft::{dft::NaiveDFT, mixed_radix::is_prime, DFT}};
+use super::mixed_radix;
+
+pub struct RaderFFT
+{
+    input_buffer: Box<[Complex32]>, 
+    output_buffer: Box<[Complex32]>, 
+    
+    size: usize,
+    g: usize,
+
+    // Fourrier transform of the exponential terms
+    convolution_operand: Box<[Complex32]>,
+    convolution_fft: NaiveDFT, // TODO: Use fft
+}
+
+impl DFT for RaderFFT
+{
+    fn create(size: usize) -> Self
+        where Self: Sized {
+            let g = compute_prime_primitive_root(size);
+            RaderFFT
+            {
+                input_buffer: vec![Complex32::zero(); size].into_boxed_slice(), 
+                output_buffer: vec![Complex32::zero(); size].into_boxed_slice(), 
+
+                size,
+                g,
+                
+                convolution_operand: compute_convolution_operand(size, g),
+                convolution_fft: NaiveDFT::create(size - 1),
+            }
+    }
+
+    fn execute(&mut self) {
+
+        // Copy to convolution fft with permutation
+        for i in 0..(self.size - 1)
+        {
+            self.convolution_fft.get_input()[i] = self.input_buffer[exp_mod(self.g, self.size - 1 - i - 1, self.size)];
+        }
+
+        self.convolution_fft.execute();
+
+        // Convolve (use output buffer as staging buffer)
+        self.output_buffer
+            .iter_mut()
+            .skip(1)
+            .zip(self.convolution_operand.iter())
+            .zip(self.convolution_fft.get_output().iter())
+            .for_each(|((dest, a), b)| *dest = *a * *b);
+
+        // Add first sample as DC-offset
+        //self.output_buffer[1] = self.output_buffer[1] + self.input_buffer[0];
+
+        // Copy to input
+        self
+            .convolution_fft
+            .get_input()
+            .iter_mut()
+            .zip(self.output_buffer.iter().skip(1))
+            .for_each(|(dest, x)| *dest = - *x);
+
+        self.convolution_fft.execute(); // Inverse fft
+
+        // Copy to output buffer : n - 1 terms to copy
+        for i in 1..(self.size)
+        {
+            self.output_buffer[i] = 
+                self.convolution_fft.get_output()[exp_mod(self.g, i, self.size) - 1] / (self.size as f32 - 1.) + self.input_buffer[0];
+        }
+        self.output_buffer[0] = self.input_buffer.iter().copied().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, g: usize) -> Box<[Complex32]>
+{
+    println!("TODO: Change to better fft");
+    let mut fft = NaiveDFT::create(n - 1); //TODO: Use fft
+
+    fft.get_input().iter_mut().enumerate()
+        .for_each(|(i, x)| *x = Complex32::cexp(- 2. * PI * (exp_mod(g, i + 1, n) as f32) / (n as f32)));
+    fft.execute();
+    fft.get_output().iter().map(|x| *x).collect::<Vec<_>>().into_boxed_slice()
+}
+
+
+
+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(n: usize, exp: usize, m: usize) -> usize
+{
+    let mut num = n % m;
+    let mut acc = 1;
+    let mut exp = exp;
+    if exp == 0
+    {
+        return 1;
+    }
+
+    while exp != 1
+    {
+        if exp % 2 == 0
+        {
+            num = (num * num) % m;
+            exp /= 2;
+        }
+        else
+        {
+            acc = (acc * n) % m;
+            exp -= 1;
+            num = num * num % m;
+            exp /= 2;
+        }
+    }
+
+    (num * acc) % m
+}