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

I want to die

Browse Source
This commit is contained in:
Albin Chaboissier
2025-09-24 23:10:28 +02:00
parent f62ef05cb8
commit 00b4756138
9 changed files with 4939 additions and 15151 deletions
Download Patch File Download Diff File Expand all files Collapse all files

19800
out.csv
View File

File diff suppressed because it is too large Load Diff

BIN
out.txt
View File

Binary file not shown.

BIN
sine.wav
View File

Binary file not shown.

68
src/bfsk.rs
View File

@ -3,13 +3,13 @@
use std::f32::consts::PI;
use crate::complex::{Complex, Complex32};
use crate::fft::{self, windows};
use crate::fft::{self, windows, FFTDirection, FFT};
use crate::map;
use crate::nco::Nco;
pub struct BFSKMod<'a, T: Iterator<Item = bool>> {
samples_per_bit: u32,
bandwidth: f32,
deviation: f32,
bit_stream: &'a mut T,
// State
@ -21,10 +21,10 @@ impl<'a, T> BFSKMod<'a, T>
where
T: Iterator<Item = bool>,
{
pub fn new(samples_per_bit: u32, bandwidth: f32, bit_stream: &'a mut T) -> Self {
pub fn new(samples_per_bit: u32, deviation: f32, bit_stream: &'a mut T) -> Self {
BFSKMod {
samples_per_bit,
bandwidth,
deviation,
oscillator: Nco::new(0.0),
bit_stream,
sample_index: samples_per_bit,
@ -37,9 +37,9 @@ where
let bit = self.bit_stream.next()?;
let frequency = if bit {
self.bandwidth / 2.0
self.deviation
} else {
-self.bandwidth / 2.0
-self.deviation
};
self.oscillator.set_frequency(frequency);
}
@ -56,62 +56,36 @@ pub struct BFSKDem {
deviation: f32,
// State
sample_index: u32,
fft: FFT,
//fft: Box<dyn DFT>,
bin_pos: usize,
bin_neg: usize,
}
impl BFSKDem {
pub fn new(samples_per_bit: u32, deviation: f32) -> Self {
// Calculate bin locations :
let bin_index = map(deviation, 0., 2. * PI, 0., samples_per_bit as f32).floor() as u32;
println!("bin_index: {bin_index}");
BFSKDem {
samples_per_bit,
deviation,
sample_index: 0,
//fft: fft::create_fft(samples_per_bit as usize, fft::FFTDirection::Forward),
fft: FFT::new(samples_per_bit as usize, windows::rectangular),
bin_pos: bin_index as usize,
bin_neg: (samples_per_bit - bin_index - 1) as usize, // -deviation = negative frequency = upper half
}
}
pub fn demod(&mut self, baseband: &[Complex32]) -> bool {
assert!(baseband.len() >= self.samples_per_bit as usize);
self.fft.execute(baseband);
/*
self.fft
.get_input()
.iter_mut()
.enumerate()
.for_each(|(i, x)| *x = baseband[i]);
self.fft.execute(windows::rectanguar);
*/
let positive_energy = self.fft.get_output()[self.bin_pos];
let negative_energy = self.fft.get_output()[self.bin_neg];
let bin_id = map(
self.deviation,
0.,
PI,
0.,
(self.samples_per_bit / 2) as f32,
)
.floor() as i32;
let bin_width = 5;
/*
let mut positive_energy = 0.0;
for i in (bin_id - bin_width)..(bin_id + bin_width) {
if i >= 0 && i < self.samples_per_bit as i32 {
positive_energy += self.fft.get_output()[i as usize].mag();
}
}
let mut negative_energy = 0.0;
for i in (self.samples_per_bit as i32 - bin_id - bin_width)
..(self.samples_per_bit as i32 - bin_id + bin_width)
{
if i >= 0 && i < self.samples_per_bit as i32 {
negative_energy += self.fft.get_output()[i as usize].mag();
}
}
return positive_energy < negative_energy;
*/
false
positive_energy.mag() < negative_energy.mag()
}
}

12
src/fft.rs
View File

@ -39,23 +39,25 @@ pub trait DFTAlgorithm {
fn get_output(&self) -> &[Complex32];
}
fn create_fft(size: usize, direction: FFTDirection) -> Box<dyn DFTAlgorithm> {
pub fn create_fft(size: usize, direction: FFTDirection) -> Box<dyn DFTAlgorithm> {
if size <= 16 {
//println!("Naive {size}");
println!("Naive {size}");
return Box::new(NaiveDFT::create(size, direction));
}
if size.count_ones() == 1 {
//println!("Radix 2 {size}");
println!("Radix 2 {size}");
return Box::new(Radix2FFT::create(size, direction));
}
if is_prime(size) {
//println!("Prime rader {size}");
println!("Prime rader {size}");
return Box::new(RaderFFT::create(size, direction));
//return Box::new(NaiveDFT::create(size, direction));
}
//println!("Mixed radix {size}");
println!("Mixed radix {size}");
Box::new(MixedRadixFFT::create(size, direction))
//Box::new(NaiveDFT::create(size, direction))
}
pub struct FFT

3
src/fft/mixed_radix.rs
View File

@ -29,6 +29,9 @@ impl DFTAlgorithm for MixedRadixFFT {
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 {
twiddle_factors: compute_twiddle_factors(size, direction),
qfft,

8
src/fft/rader.rs
View File

@ -4,7 +4,7 @@ use std::f32::consts::PI;
use crate::{
complex::Complex32,
fft::{create_fft, is_prime , DFTAlgorithm, FFTDirection},
fft::{create_fft, dft::NaiveDFT, is_prime, DFTAlgorithm, FFTDirection},
};
pub struct RaderFFT {
@ -31,7 +31,8 @@ impl DFTAlgorithm for RaderFFT {
let permutations: Box<[usize]> = (0..(size - 1)).map(|i| exp_mod(g, i + 1, size)).collect();
// Compute fourrier transform of twiddle factors
let mut convolution_fft = create_fft(size - 1, FFTDirection::Forward);
//let mut convolution_fft = create_fft(size - 1, FFTDirection::Forward);
let mut convolution_fft = Box::new(NaiveDFT::create(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>>();
@ -42,7 +43,8 @@ impl DFTAlgorithm for RaderFFT {
permutations,
convolution_operand: convolution_operand.into(),
convolution_ifft: create_fft(size - 1, FFTDirection::Inverse),
//convolution_ifft: create_fft(size - 1, FFTDirection::Inverse),
convolution_ifft: Box::new(NaiveDFT::create(size - 1, FFTDirection::Inverse)),
convolution_fft,
output: vec![Complex32::zero(); size].into(),

141
src/fft/rader2.rs
View File

@ -1,118 +1,91 @@
// Implementation of raders's fft for prime sized ffts
/*
use std::{f32::consts::PI, ops::Deref};
use super::mixed_radix;
use std::f32::consts::PI;
use crate::{
complex::Complex32,
fft::{DFT, FFTDirection, create_fft, dft::NaiveDFT, is_prime, windows},
fft::{create_fft, dft::NaiveDFT, is_prime, DFTAlgorithm, FFTDirection},
};
pub struct Rader2FFT {
input_buffer: Box<[Complex32]>,
output_buffer: Box<[Complex32]>,
pub struct RaderFFT {
permutations: Box<[usize]>,
convolution_operand: Box<[Complex32]>,
convolution_ifft: Box<dyn DFTAlgorithm>,
convolution_fft: Box<dyn DFTAlgorithm>,
output: 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 {
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 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(),
let permutations: Box<[usize]> = (0..(size - 1)).map(|i| exp_mod(g, i + 1, size)).collect();
// Compute fourrier transform of twiddle factors
let mut convolution_fft = create_fft(size - 1, FFTDirection::Forward);
//let mut convolution_fft = Box::new(NaiveDFT::create(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_operand: convolution_operand.into(),
convolution_ifft: create_fft(size - 1, FFTDirection::Inverse),
//convolution_ifft: Box::new(NaiveDFT::create(size - 1, FFTDirection::Inverse)),
convolution_fft,
output: vec![Complex32::zero(); size].into(),
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)
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];
// Using output as staging buffer
self.output[i] = input[k];
}
self.convolution_fft.execute(windows::rectanguar);
self.convolution_fft.execute(&self.output);
// 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];
// Compute convolution by multiplying in freq domain
for i in 0..(self.size - 1) {
// Using output as staging buffer
self.output[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.);
*/
self.convolution_ifft.execute(&self.output);
// Compute ifft to obtain convolution
self.convolution_fft.execute(window);
self.output[0] = input[0];
for i in 0..(self.size - 1) {
self.output_buffer[self.permutation[i]] =
self.convolution_fft.get_output()[i] / self.sub_size as f32;
// Actually compute the output
let k = self.permutations[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];
}
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
&self.output
}
}
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));
@ -123,7 +96,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;
}
@ -158,12 +131,4 @@ pub fn exp_mod(mut n: usize, mut exp: usize, m: usize) -> usize {
r
}
pub fn next_pow2(mut n: usize) -> usize {
let mut pow = 0;
while n > 0 {
n >>= 1;
pow += 1;
}
1 << pow
}
*/

58
src/main.rs
View File

@ -19,7 +19,7 @@ use nco::Nco;
use plotters::prelude::*;
use fft::DFTAlgorithm;
use crate::{bfsk::BFSKDem, fft::{dft::NaiveDFT, mixed_radix::MixedRadixFFT, rader::RaderFFT, radix2::Radix2FFT, windows, FFT}};
use crate::{bfsk::BFSKDem, fft::{create_fft, dft::NaiveDFT, mixed_radix::MixedRadixFFT, rader::RaderFFT, radix2::Radix2FFT, windows, FFTDirection, FFT}};
// Utilities
fn map<T>(input: T, in_min: T, in_max: T, out_min: T, out_max: T) -> T
@ -29,15 +29,46 @@ where
((input - in_min.clone()) / (in_max - in_min)) * (out_max - out_min.clone()) + out_min
}
fn main() {
modulate();
//modulate();
test();
}
fn test()
{
let mut o1 = Nco::new(PI / 2.0);
let mut o2 = Nco::new(PI / 4.0);
let sample_count = 4800;
//let mut fft = FFT::new(sample_count, windows::rectangular);
let mut dft = NaiveDFT::create(sample_count, FFTDirection::Forward);
let mut fft = FFT::new(sample_count, windows::rectangular);
//let mut fft = RaderFFT::create(sample_count, FFTDirection::Forward);
let mut fft_input = vec![Complex32::zero(); sample_count];
for x in fft_input.iter_mut()
{
*x = o1.cexp() + o2.cexp();
o1.step();
o2.step();
}
fft.execute(&fft_input);
dft.execute(&fft_input);
let mut out_file = File::create("out.csv").unwrap();
for (x, y) in fft.get_output().iter().zip(dft.get_output())
{
out_file.write_all(
format!("{},{},\n", x.mag(), y.mag()).as_bytes()
).unwrap();
}
}
fn modulate() {
let sample_rate = 44100;
let mut frequency = 2000.0; //HZ
let mut bandwidth = 500.0; //HZ
let frequency = 2000.0; //HZ
let bandwidth = 1000.0; //HZ
println!("deviation: {}", PI * (bandwidth / sample_rate as f32));
let path = "s.txt";
let file = File::open(path).unwrap();
@ -61,7 +92,7 @@ fn modulate() {
println!("{} samples/bit", sample_rate / baud_rate);
let mut bfsk = BFSKMod::new(
sample_rate / baud_rate,
2. * PI * (bandwidth / sample_rate as f32),
PI * 0.05, //PI * (bandwidth / sample_rate as f32),
&mut bit_stream,
);
@ -90,6 +121,17 @@ fn modulate() {
lo.step();
}
writer.finalize().unwrap();
let mut tfft = FFT::new(44100, windows::rectangular);
tfft.execute(&output_samples);
// Write csv
let mut out_csv = File::create("out.csv").unwrap();
for x in output_samples.iter().take(4400)
{
out_csv.write_all(
format!("{},\n", x.mag()).as_bytes()
).unwrap();
}
let mut of = File::create("out.txt").unwrap();
@ -97,7 +139,7 @@ fn modulate() {
let mut lodem = Nco::new(-2. * PI * (frequency / sample_rate as f32));
let mut demod = BFSKDem::new(
sample_rate / baud_rate,
PI * (bandwidth / sample_rate as f32),
PI * 0.05, //PI * (bandwidth / sample_rate as f32),
);
for chunk in output_samples.chunks((sample_rate / baud_rate) as usize) {
let base_chunk: Vec<Complex32> = chunk
@ -109,7 +151,7 @@ fn modulate() {
.collect();
let bit = demod.demod(base_chunk.as_slice());
bits.push(bit);
println!("{:?}", bit)
//println!("{:?}", bit)
}
for b in bits.chunks(8) {