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

Raders algorithm

Browse Source
This commit is contained in:
Albin Chaboissier
2025-09-21 19:42:51 +02:00
parent 79f03a071a
commit 01e3657b55
5 changed files with 4975 additions and 4968 deletions
Download Patch File Download Diff File Expand all files Collapse all files

9598
out.csv
View File

File diff suppressed because it is too large Load Diff

58
src/fft.rs
View File

@ -1,16 +1,66 @@
pub mod mixed_radix;
pub mod radix2;
pub mod dft;
pub mod mixed_radix;
pub mod rader;
pub mod radix2;
use crate::complex::Complex32;
use crate::{
complex::Complex32,
fft::{dft::NaiveDFT, mixed_radix::MixedRadixFFT, rader::RaderFFT, radix2::Radix2FFT},
};
pub trait DFT {
fn create(size: usize) -> Self
where Self: Sized;
where
Self: Sized;
fn get_input(&mut self) -> &mut [Complex32];
fn get_output(&self) -> &[Complex32];
fn execute(&mut self);
}
pub fn create_fft(size: usize) -> Box<dyn DFT> {
if size == 1 {
return Box::new(NaiveDFT::create(size));
}
if size.count_ones() == 1 {
// TODO: Return hardcoded fft for small sized
return Box::new(Radix2FFT::create(size));
}
// Get factors
if is_prime(size) {
return Box::new(RaderFFT::create(size));
}
return Box::new(MixedRadixFFT::create(size));
}
// Utilities
fn prime_factors(n: usize) -> Vec<usize> {
let mut factors = vec![];
let mut num = n;
// Divide num successively
while num != 1 {
// Try divisors from 2 up to n (included)
for i in 2..n + 1 {
// if i divides num, it is a prime factor (if it wasn't, then i would have prime
// factors that would divide into num before i)
if num % i == 0 {
factors.push(i);
num /= i;
}
}
}
// If n = 1 then it does not have any prime factors
// The prime factor decomposition theorem states that any integer
// greater than TWO has a unique decomposition
factors
}
pub fn is_prime(n: usize) -> bool {
prime_factors(n).len() == 1
}

78
src/fft/mixed_radix.rs
View File

@ -2,7 +2,10 @@
use std::f32::consts::PI;
use crate::{complex::Complex32, fft::{dft::NaiveDFT, DFT}};
use crate::{
complex::Complex32,
fft::{DFT, create_fft, dft::NaiveDFT, prime_factors},
};
pub struct MixedRadixFFT {
input_buffer: Box<[Complex32]>,
@ -16,23 +19,23 @@ pub struct MixedRadixFFT {
qfft: Box<dyn DFT>,
pfft: Box<dyn DFT>,
staging_buffer: Box<[Complex32]>
staging_buffer: Box<[Complex32]>,
}
impl DFT for MixedRadixFFT {
fn create(size: usize) -> Self {
let q = decide_radix_factor(size);
let p = size / q;
let qfft = NaiveDFT::create(q);
let pfft = NaiveDFT::create(p);
let qfft = create_fft(q);
let pfft = create_fft(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),
qfft: Box::new(qfft),
pfft: Box::new(pfft),
qfft,
pfft,
staging_buffer: vec![Complex32::zero(); size].into_boxed_slice(),
p,
@ -40,43 +43,36 @@ impl DFT for MixedRadixFFT {
}
}
fn execute(&mut self)
{
fn execute(&mut self) {
// Perform p ffts of size q
for k0 in 0..self.p
{
for k0 in 0..self.p {
// Copy samples into input buffer
for k1 in 0..self.q
{
for k1 in 0..self.q {
self.qfft.get_input()[k1] = self.input_buffer[k1 * self.p + k0];
}
self.qfft.execute();
for j0 in 0..self.q
{
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[k0 * self.q + j0] =
self.qfft.get_output()[j0] * self.twiddle_factors[k0 * self.q + j0];
}
}
// Perform q ffts of size p
for j0 in 0..self.q
{
for j0 in 0..self.q {
// Copy input
for k0 in 0..self.p
{
for k0 in 0..self.p {
self.pfft.get_input()[k0] = self.staging_buffer[k0 * self.q + j0];
}
self.pfft.execute();
for j1 in 0..self.p
{
for j1 in 0..self.p {
self.output_buffer[j1 * self.q + j0] = self.pfft.get_output()[j1];
}
}
}
fn get_input(&mut self) -> &mut [Complex32] {
@ -88,15 +84,12 @@ impl DFT for MixedRadixFFT {
}
}
fn compute_twiddle_factors(q: usize, p: usize) -> Box<[Complex32]>
{
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;
for i in 0..q
{
for j in 0..p
{
for i in 0..q {
for j in 0..p {
factors[i * p + j] = Complex32::cexp(2. * PI / (n as f32));
}
}
@ -116,32 +109,3 @@ fn decide_radix_factor(n: usize) -> usize {
// Otherwise take next big prime
return *factors.iter().skip(two_count).next().unwrap();
}
// Utilities
fn prime_factors(n: usize) -> Vec<usize> {
let mut factors = vec![];
let mut num = n;
// Divide num successively
while num != 1 {
// Try divisors from 2 up to n (included)
for i in 2..n + 1 {
// if i divides num, it is a prime factor (if it wasn't, then i would have prime
// factors that would divide into num before i)
if num % i == 0 {
factors.push(i);
num /= i;
}
}
}
// If n = 1 then it does not have any prime factors
// The prime factor decomposition theorem states that any integer
// greater than TWO has a unique decomposition
factors
}
pub fn is_prime(n: usize) -> bool {
prime_factors(n).len() == 1
}

173
src/fft/rader.rs
View File

@ -1,78 +1,74 @@
// Implementation of raders's fft for prime sized ffts
use std::f32::consts::PI;
use std::{f32::consts::PI, ops::Deref};
use crate::{complex::Complex32, fft::{dft::NaiveDFT, mixed_radix::is_prime, DFT}};
use super::mixed_radix;
use crate::{
complex::Complex32,
fft::{DFT, create_fft, dft::NaiveDFT, is_prime},
};
pub struct RaderFFT {
input_buffer: Box<[Complex32]>,
output_buffer: Box<[Complex32]>,
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
convolution_fft: Box<dyn DFT>, // TODO: Use fft
permutation: Box<[usize]>,
}
impl DFT for RaderFFT
{
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(),
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();
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),
}
size,
g,
convolution_operand: compute_convolution_operand(size, &permutation),
convolution_fft: Box::new(NaiveDFT::create(size - 1)),
permutation,
}
}
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)];
for i in 0..(self.size - 1) {
self.convolution_fft.get_input()[i] =
self.input_buffer[self.permutation[self.size - 1 - i - 1]]
}
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];
// Use output buffer as staging buffer
for i in 0..(self.size - 1) {
self.output_buffer[i] =
self.convolution_fft.get_output()[i] * self.convolution_operand[i];
}
for i in 0..(self.size - 1) {
self.convolution_fft.get_input()[i] = self.output_buffer[self.size - 1 - i - 1];
}
self.convolution_fft.get_input()[0] =
self.convolution_fft.get_input()[0] + self.input_buffer[0];
// Compute ifft to obtain convolution
self.convolution_fft.execute();
for i in 0..(self.size - 1) {
self.output_buffer[self.permutation[i]] =
self.convolution_fft.get_output()[i] / (self.size - 1) as f32;
}
self.output_buffer[0] = self.input_buffer.iter().copied().sum();
}
@ -83,46 +79,39 @@ impl DFT for RaderFFT
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
pub fn compute_convolution_operand(n: usize, permutation: &[usize]) -> Box<[Complex32]> {
//let mut fft = create_fft(n - 1);
let mut fft = NaiveDFT::create(n - 1);
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.get_input()
.iter_mut()
.enumerate()
.for_each(|(i, x)| *x = Complex32::cexp(-2. * PI * (permutation[i] as f32) / (n as f32)));
fft.execute();
fft.get_output().iter().map(|x| *x).collect::<Vec<_>>().into_boxed_slice()
fft.get_output().iter().copied().collect()
}
pub fn compute_prime_primitive_root(n: usize) -> usize
{
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
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
{
for j in 0..n {
if val == 1 {
break;
}
val = (val * i) % n;
order += 1;
}
if order == phi
{
if order == phi {
return i;
}
}
@ -130,31 +119,21 @@ pub fn compute_prime_primitive_root(n: usize) -> usize
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;
pub fn exp_mod(mut n: usize, mut exp: usize, m: usize) -> usize {
if m == 1 {
return 0;
}
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;
n %= m;
let mut r = 1;
while exp > 0 {
if exp % 2 == 1 {
r = (r * n) % m;
}
n = (n * n) % m;
exp >>= 1;
}
(num * acc) % m
r
}

36
src/main.rs
View File

@ -12,11 +12,17 @@ mod nco;
use bfsk::BFSKMod;
use complex::Complex;
use fft::rader;
use complex::Complex32;
use fft::rader;
use nco::Nco;
use crate::fft::{dft::NaiveDFT, mixed_radix::MixedRadixFFT, rader::{compute_prime_primitive_root, RaderFFT}, radix2::Radix2FFT, DFT};
use crate::fft::{
DFT, create_fft,
dft::NaiveDFT,
mixed_radix::MixedRadixFFT,
rader::{RaderFFT, compute_prime_primitive_root, exp_mod},
radix2::Radix2FFT,
};
// Utilities
fn map<T>(input: T, in_min: T, in_max: T, out_min: T, out_max: T) -> T
@ -38,16 +44,18 @@ fn test() {
let freq1 = 2. * PI / 4.0;
let freq2 = 2. * PI / 8.0;
let sample_count = 4799;
let mut o1 = Nco::new(freq1);
let mut o2 = Nco::new(freq2);
let mut fft = RaderFFT::create(4799);
let mut dft = NaiveDFT::create(4799);
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, y) in vals.iter_mut().zip(vals_dft.iter_mut()) {
//let vals_dft = dft.get_input();
for x in vals.iter_mut() {
*x = o1.cexp() + o2.cexp();
*y = *x;
//*y = *x;
//*x = o2.cexp(); //+ o2.cexp();
//*x = *x * (1. / x.mag());
o1.step();
@ -55,14 +63,20 @@ fn test() {
}
fft.execute();
//dft.execute();
let output = fft.get_output();
let mut f = File::create("out.csv").unwrap();
for (i, (v, v2)) in output.iter().zip(dft.get_output()).enumerate() {
for (i, v) in output.iter().enumerate() {
f.write_all(
format!("{},{},{},\n", i as f32 / 8192., v.mag(), v2.mag())
.to_string()
.as_bytes(),
format!(
"{},{},\n",
i as f32 / sample_count as f32,
v.mag(),
//v2.mag()
)
.to_string()
.as_bytes(),
)
.unwrap();
}