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

Still not working

Browse Source
This commit is contained in:
Albin Chaboissier
2025-09-21 22:43:24 +02:00
parent 01e3657b55
commit a300386f7f
8 changed files with 5695 additions and 4848 deletions
Download Patch File Download Diff File Expand all files Collapse all files

10
src/fft.rs
View File

@ -3,6 +3,8 @@ pub mod mixed_radix;
pub mod rader;
pub mod radix2;
use std::iter::Map;
use crate::{
complex::Complex32,
fft::{dft::NaiveDFT, mixed_radix::MixedRadixFFT, rader::RaderFFT, radix2::Radix2FFT},
@ -28,16 +30,15 @@ pub fn create_fft(size: usize) -> Box<dyn DFT> {
return Box::new(Radix2FFT::create(size));
}
// Get factors
if is_prime(size) {
return Box::new(RaderFFT::create(size));
return Box::new(NaiveDFT::create(size));
}
return Box::new(MixedRadixFFT::create(size));
Box::new(MixedRadixFFT::create(size))
}
// Utilities
fn prime_factors(n: usize) -> Vec<usize> {
pub fn prime_factors(n: usize) -> Vec<usize> {
let mut factors = vec![];
let mut num = n;
@ -50,6 +51,7 @@ fn prime_factors(n: usize) -> Vec<usize> {
if num % i == 0 {
factors.push(i);
num /= i;
break;
}
}
}

27
src/fft/mixed_radix.rs
View File

@ -26,14 +26,19 @@ impl DFT for MixedRadixFFT {
fn create(size: usize) -> Self {
let q = decide_radix_factor(size);
let p = size / q;
// TODO: Figure out why it does not work in the other direction ...
let (p, q) = (q, p);
let qfft = create_fft(q);
let pfft = create_fft(p);
//let qfft = Box::new(NaiveDFT::create(q));
//let pfft = Box::new(NaiveDFT::create(p));
MixedRadixFFT {
input_buffer: vec![Complex32::zero(); size].into_boxed_slice(),
output_buffer: vec![Complex32::zero(); size].into_boxed_slice(),
size,
twiddle_factors: compute_twiddle_factors(q, p),
twiddle_factors: compute_twiddle_factors(size),
qfft,
pfft,
@ -55,8 +60,8 @@ impl DFT for MixedRadixFFT {
for j0 in 0..self.q {
// "Store j0'th of k0'th fft into staging buffer"
self.staging_buffer[k0 * self.q + j0] =
self.qfft.get_output()[j0] * self.twiddle_factors[k0 * self.q + j0];
self.staging_buffer[j0 * self.p + k0] =
self.qfft.get_output()[j0] * self.twiddle_factors[j0 * k0];
}
}
@ -64,7 +69,7 @@ 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[k0 * self.q + j0];
self.pfft.get_input()[k0] = self.staging_buffer[j0 * self.p + k0];
}
self.pfft.execute();
@ -84,14 +89,11 @@ impl DFT for MixedRadixFFT {
}
}
fn compute_twiddle_factors(q: usize, p: usize) -> Box<[Complex32]> {
let mut factors = vec![Complex32::zero(); q * p].into_boxed_slice();
let n = p * q;
fn compute_twiddle_factors(size: usize) -> Box<[Complex32]> {
let mut factors = vec![Complex32::zero(); size].into_boxed_slice();
for i in 0..q {
for j in 0..p {
factors[i * p + j] = Complex32::cexp(2. * PI / (n as f32));
}
for i in 0..size {
factors[i] = Complex32::cexp(2. * PI * i as f32 / (size as f32));
}
factors
}
@ -99,6 +101,7 @@ fn compute_twiddle_factors(q: usize, p: usize) -> Box<[Complex32]> {
// This will decide on a good factor to use for the mixed radix fft
fn decide_radix_factor(n: usize) -> usize {
let factors = prime_factors(n);
//factors.iter().for_each(|a| print!(" {a} "));
let two_count = factors.iter().take_while(|i| **i == 2).count();
// If there is a lot of two, we can use them as q factor to be able to use radix2 later on
@ -107,5 +110,5 @@ fn decide_radix_factor(n: usize) -> usize {
}
// Otherwise take next big prime
return *factors.iter().skip(two_count).next().unwrap();
*factors.get(two_count).unwrap()
}

4
src/fft/rader.rs
View File

@ -13,7 +13,6 @@ pub struct RaderFFT {
output_buffer: Box<[Complex32]>,
size: usize,
g: usize,
// Fourrier transform of the exponential terms
convolution_operand: Box<[Complex32]>,
@ -33,10 +32,9 @@ impl DFT for RaderFFT {
output_buffer: vec![Complex32::zero(); size].into_boxed_slice(),
size,
g,
convolution_operand: compute_convolution_operand(size, &permutation),
convolution_fft: Box::new(NaiveDFT::create(size - 1)),
convolution_fft: create_fft(size - 1),
permutation,
}
}

66
src/main.rs
View File

@ -15,13 +15,11 @@ use complex::Complex;
use complex::Complex32;
use fft::rader;
use nco::Nco;
use plotters::prelude::*;
use crate::fft::{
DFT, create_fft,
dft::NaiveDFT,
mixed_radix::MixedRadixFFT,
rader::{RaderFFT, compute_prime_primitive_root, exp_mod},
radix2::Radix2FFT,
create_fft, dft::NaiveDFT, mixed_radix::MixedRadixFFT, prime_factors, rader::{compute_prime_primitive_root, exp_mod, RaderFFT}, radix2::Radix2FFT, DFT
};
// Utilities
@ -44,42 +42,50 @@ fn test() {
let freq1 = 2. * PI / 4.0;
let freq2 = 2. * PI / 8.0;
let sample_count = 4799;
//let sample_count = 71*71;
//let sample_count = 71*71;
let sample_count = 4800;
let mut o1 = Nco::new(freq1);
let mut o2 = Nco::new(freq2);
let mut fft = RaderFFT::create(sample_count);
//let mut dft = NaiveDFT::create(sample_count);
let vals = fft.get_input();
//let vals_dft = dft.get_input();
for x in vals.iter_mut() {
let mut fft = create_fft(sample_count);
//let mut fft = create_fft(sample_count);
for x in fft.get_input().iter_mut() {
*x = o1.cexp() + o2.cexp();
//*y = *x;
//*x = o2.cexp(); //+ o2.cexp();
//*x = *x * (1. / x.mag());
o1.step();
o2.step();
}
fft.execute();
//dft.execute();
let output = fft.get_output();
let mut f = File::create("out.csv").unwrap();
for (i, v) in output.iter().enumerate() {
f.write_all(
format!(
"{},{},\n",
i as f32 / sample_count as f32,
v.mag(),
//v2.mag()
)
.to_string()
.as_bytes(),
)
.unwrap();
}
let root = BitMapBackend::new("out.png", (640, 480)).into_drawing_area();
root.fill(&WHITE).unwrap();
let mut chart = ChartBuilder::on(&root)
.caption("fft", ("sans-serif", 50).into_font())
.margin(5)
.x_label_area_size(30)
.y_label_area_size(30)
.build_cartesian_2d(0f32..(sample_count as f32), 0f32..(sample_count as f32)).unwrap();
//chart.configure_mesh().draw()?;
chart
.draw_series(LineSeries::new(
(0..sample_count).zip(fft.get_output().iter()).map(|(x, y)| (x as f32, y.mag())),
&RED,
)).unwrap()
.legend(|(x, y)| PathElement::new(vec![(x, y), (x + 20, y)], &RED));
chart
.configure_series_labels()
.background_style(&WHITE.mix(0.8))
.border_style(&BLACK)
.draw().unwrap();
root.present().unwrap();
}
fn modulate() {