octagonal/rdsp-experiments
1
1
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Finishes fft interface + algorithms

Browse Source
This commit is contained in:
Albin Chaboissier
2025-09-24 21:30:45 +02:00
parent 3cc4144747
commit f62ef05cb8
8 changed files with 169 additions and 88 deletions
Download Patch File Download Diff File Expand all files Collapse all files

21
src/fft/dft.rs
View File

@ -1,31 +1,36 @@
use crate::complex::Complex32;
use crate::fft::{DFT, FFTDirection};
use crate::fft::{DFTAlgorithm, FFTDirection};
use std::f32::consts::PI;
pub struct NaiveDFT {
direction: FFTDirection,
size: usize,
output: Box<[Complex32]>,
}
impl DFT for NaiveDFT {
impl DFTAlgorithm for NaiveDFT {
fn create(size: usize, direction: FFTDirection) -> Self
where
Self: Sized,
{
NaiveDFT { direction, size }
NaiveDFT { direction, size, output: vec![Complex32::zero(); size].into() }
}
fn execute(&mut self, input: &[Complex32], output: &mut [Complex32], window: fn(f32) -> f32) {
for (freq, out) in output.iter_mut().enumerate() {
fn execute(&mut self, input: &[Complex32]) {
for (freq, out) in self.output.iter_mut().enumerate() {
*out = Complex32::zero();
for (i, inp) in input.iter().enumerate() {
*out = *out
+ ((*inp
+ (*inp
* Complex32::cexp(
-2. * self.direction.sign() * PI * (i * freq) as f32 / self.size as f32,
))
* window(i as f32 / self.size as f32));
));
}
}
}
fn get_output(&self) -> &[Complex32]
{
&self.output
}
}

40
src/fft/mixed_radix.rs
View File

@ -4,55 +4,55 @@ use std::f32::consts::PI;
use crate::{
complex::Complex32,
fft::{DFT, FFTDirection, create_fft, dft::NaiveDFT, prime_factors, windows},
fft::{DFTAlgorithm, FFTDirection, create_fft, dft::NaiveDFT, prime_factors, windows},
};
pub struct MixedRadixFFT {
size: usize,
//size: usize, size is implicitely stored in p and q
p: usize,
q: usize,
twiddle_factors: Box<[Complex32]>,
qfft: Box<dyn DFT>,
pfft: Box<dyn DFT>,
qfft: Box<dyn DFTAlgorithm>,
pfft: Box<dyn DFTAlgorithm>,
staging_buffer: Box<[Complex32]>,
pfft_input: Box<[Complex32]>,
output: Box<[Complex32]>
}
impl DFT for MixedRadixFFT {
impl DFTAlgorithm for MixedRadixFFT {
fn create(size: usize, direction: FFTDirection) -> Self {
let q = decide_radix_factor(size);
let p = size / q;
let qfft = create_fft(q, direction);
let pfft = create_fft(p, direction);
//let qfft = Box::new(NaiveDFT::create(q, direction));
//let pfft = Box::new(NaiveDFT::create(p, direction));
MixedRadixFFT {
size,
twiddle_factors: compute_twiddle_factors(size, direction),
qfft,
pfft,
staging_buffer: vec![Complex32::zero(); size].into_boxed_slice(),
pfft_input: vec![Complex32::zero(); p].into_boxed_slice(),
output: vec![Complex32::zero(); size].into_boxed_slice(),
p,
q,
}
}
fn execute(&mut self, input: &[Complex32], output: &mut [Complex32], window: fn(f32) -> f32) {
fn execute(&mut self, input: &[Complex32]) {
// Perform p ffts of size q
for k0 in 0..self.p {
// Copy samples into input buffer
for k1 in 0..self.q {
let k = k1 * self.p + k0;
self.qfft.get_input()[k1] =
self.input_buffer[k] * window(k as f32 / self.size as f32);
// Use output as staging buffer
self.output[k1] = input[k];
}
self.qfft.execute(windows::rectanguar);
self.qfft.execute(&self.output);
for j0 in 0..self.q {
// "Store j0'th of k0'th fft into staging buffer"
@ -65,23 +65,21 @@ impl DFT for MixedRadixFFT {
for j0 in 0..self.q {
// Copy input
for k0 in 0..self.p {
self.pfft.get_input()[k0] = self.staging_buffer[j0 * self.p + k0];
// Use output as staging buffer
self.pfft_input[k0] = self.staging_buffer[j0 * self.p + k0];
}
self.pfft.execute(windows::rectanguar);
self.pfft.execute(&self.pfft_input);
// Actually compute final output
for j1 in 0..self.p {
self.output_buffer[j1 * self.q + j0] = self.pfft.get_output()[j1];
self.output[j1 * self.q + j0] = self.pfft.get_output()[j1];
}
}
}
fn get_input(&mut self) -> &mut [Complex32] {
&mut self.input_buffer
}
fn get_output(&self) -> &[Complex32] {
&self.output_buffer
&self.output
}
}

77
src/fft/rader.rs
View File

@ -1,82 +1,85 @@
// Implementation of raders's fft for prime sized ffts
use std::{f32::consts::PI, ops::Deref};
use std::f32::consts::PI;
use super::mixed_radix;
use crate::{
complex::Complex32,
fft::{DFT, FFTDirection, create_fft, dft::NaiveDFT, is_prime, windows},
fft::{create_fft, is_prime , DFTAlgorithm, FFTDirection},
};
pub struct RaderFFT {
permutations: Box<[usize]>,
convolution_op: Box<[Complex32]>,
staging_buffer: Box<[Complex32]>,
inv_fft: Box<dyn DFT>,
conv_fft: Box<dyn DFT>,
convolution_operand: Box<[Complex32]>,
convolution_ifft: Box<dyn DFTAlgorithm>,
convolution_fft: Box<dyn DFTAlgorithm>,
output: Box<[Complex32]>,
size: usize,
}
impl DFT for RaderFFT {
impl DFTAlgorithm for RaderFFT {
fn create(size: usize, direction: FFTDirection) -> Self
where
Self: Sized,
{
assert!(is_prime(size));
// Primitive root and its powers
let g = compute_prime_primitive_root(size);
let permutations: Box<[usize]> = (0..(size - 1)).map(|i| exp_mod(g, i + 1, size)).collect();
let mut conv_fft = create_fft(size - 1, FFTDirection::Forward);
//let mut conv_fft = create_fft(size - 1);
let mut convolution_op = vec![Complex32::zero(); size - 1];
let conv_fft_input: Vec<Complex32> = (0..(size - 1))
.map(|i| {
Complex32::cexp(
-2. * direction.sign() * PI * (permutations[i] as f32) / (size as f32),
)
})
.collect();
conv_fft.execute(&conv_fft_input, &mut convolution_op, windows::rectangular);
// Compute fourrier transform of twiddle factors
let mut convolution_fft = create_fft(size - 1, FFTDirection::Forward);
let mut convolution_operand = (0..(size - 1))
.map(|i| {Complex32::cexp(-2. * direction.sign() * PI * (permutations[i] as f32) / (size as f32))})
.collect::<Vec<Complex32>>();
convolution_fft.execute(&convolution_operand);
convolution_operand = Vec::from(convolution_fft.get_output());
RaderFFT {
permutations,
convolution_op: convolution_op.into(),
staging_buffer: vec![Complex32::zero(); size - 1].into(),
inv_fft: create_fft(size - 1, FFTDirection::Inverse),
conv_fft,
convolution_operand: convolution_operand.into(),
convolution_ifft: create_fft(size - 1, FFTDirection::Inverse),
convolution_fft,
output: vec![Complex32::zero(); size].into(),
size,
}
}
fn execute(&mut self, input: &[Complex32], output: &mut [Complex32], window: fn(f32) -> f32) {
fn execute(&mut self, input: &[Complex32]) {
// Compute fft of input signal
for i in 0..(self.size - 1) {
let k = self.permutations[self.size - 1 - i - 1];
self.staging_buffer[i] = input[k] * window(k as f32 / (self.size as f32));
// Using output as staging buffer
self.output[i] = input[k];
}
self.conv_fft
.execute(&self.staging_buffer, output, windows::rectangular);
self.convolution_fft.execute(&self.output);
// Compute convolution by multiplying in freq domain
for i in 0..(self.size - 1) {
self.staging_buffer[i] = output[i] * self.convolution_op[i];
// Using output as staging buffer
self.output[i] = self.convolution_fft.get_output()[i] * self.convolution_operand[i];
}
self.inv_fft
.execute(&self.staging_buffer, output, windows::rectangular);
self.convolution_ifft.execute(&self.output);
self.output[0] = input[0];
for i in 0..(self.size - 1) {
// Actually compute the output
let k = self.permutations[i];
self.staging_buffer[k - 1] = output[i];
self.output[k] = (self.convolution_ifft.get_output()[i] / (self.size - 1) as f32) + input[0];
self.output[0] = self.output[0] + input[i + 1];
}
}
output[0] = input[0] * window(0.0);
for i in 0..(self.size - 1) {
output[i + 1] = (self.staging_buffer[i] / (self.size - 1) as f32) + input[0];
output[0] = output[0] + (input[i + 1] * window((i + 1) as f32 / self.size as f32));
}
fn get_output(&self) -> &[Complex32] {
&self.output
}
}
@ -90,7 +93,7 @@ pub fn compute_prime_primitive_root(n: usize) -> usize {
// Find multiplicative order of i
let mut val = i;
let mut order = 1;
for j in 0..n {
for _ in 0..n {
if val == 1 {
break;
}

24
src/fft/radix2.rs
View File

@ -1,16 +1,17 @@
// Cooley-Tukey algorithm
use crate::complex::Complex32;
use crate::fft::{DFT, FFTDirection};
use crate::fft::{DFTAlgorithm, FFTDirection};
use std::f32::consts::PI;
pub struct Radix2FFT {
direction: FFTDirection,
size: usize,
length: usize,
output: Box<[Complex32]>
}
impl DFT for Radix2FFT {
impl DFTAlgorithm for Radix2FFT {
// Size as power of two
fn create(size: usize, direction: FFTDirection) -> Self {
if !is_power_of_two(size) {
@ -21,14 +22,15 @@ impl DFT for Radix2FFT {
size: size.ilog2() as usize,
direction,
length: size,
output: vec![Complex32::zero(); size].into()
}
}
fn execute(&mut self, input: &[Complex32], output: &mut [Complex32], window: fn(f32) -> f32) {
fn execute(&mut self, input: &[Complex32]) {
// Reorder samples
for (i, x) in output.iter_mut().enumerate() {
for (i, x) in self.output.iter_mut().enumerate() {
let k = reverse_bits(i, self.size as u32);
*x = input[k] * window(k as f32 / self.size as f32);
*x = input[k];
}
for step in 1..(self.size + 1) {
@ -37,16 +39,20 @@ impl DFT for Radix2FFT {
for s in (0..(self.length / pol_length)).map(|i| i * pol_length) {
for i in 0..mid_point {
// Compute current polynomial at each unit root
let a = output[s + i];
let b = output[s + i + mid_point];
let a = self.output[s + i];
let b = self.output[s + i + mid_point];
let angle = -2. * self.direction.sign() * PI * (i as f32) / (pol_length as f32);
let phasor = Complex32::cexp(angle);
output[i + s] = a + phasor * b;
output[i + s + mid_point] = a - phasor * b;
self.output[i + s] = a + phasor * b;
self.output[i + s + mid_point] = a - phasor * b;
}
}
}
}
fn get_output(&self) -> &[Complex32] {
&self.output
}
}
// Utilities