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

Correlation bfsk fir design

Browse Source
This commit is contained in:
Albin Chaboissier
2025-09-28 20:41:13 +02:00
parent 4a6d79c0b4
commit 1445887f2f
15 changed files with 4032 additions and 137 deletions
Download Patch File Download Diff File Expand all files Collapse all files

52
src/bfsk.rs
View File

@ -3,7 +3,8 @@
use std::f32::consts::PI;
use crate::complex::{Complex, Complex32};
use crate::fft::{self, FFT, FFTDirection, windows};
use crate::fft::{self, FFT, FFTDirection};
use crate::{units, windows};
use crate::map;
use crate::nco::Nco;
@ -54,6 +55,52 @@ impl<'a, T: Iterator<Item = bool>> Iterator for BFSKMod<'a, T> {
}
}
// FSK Demodulator (dumb non coherent + no symbol timing recovery)
// By correlation
pub struct BFSKDem {
samples_per_bit: u32,
deviation: f32,
// State
sample_index: u32,
// Correlators
pos_correlator: Box<[Complex32]>,
neg_correlator: Box<[Complex32]>,
}
impl BFSKDem {
pub fn new(samples_per_bit: u32, deviation: f32) -> Self {
let mut nco_pos = Nco::new(deviation);
let mut nco_neg = Nco::new(-deviation);
BFSKDem {
samples_per_bit,
deviation,
sample_index: 0,
pos_correlator: (0..samples_per_bit).map(|_| {nco_pos.step(); nco_pos.cexp()}).collect(),
neg_correlator: (0..samples_per_bit).map(|_| {nco_neg.step(); nco_neg.cexp()}).collect()
}
}
pub fn demod(&mut self, baseband: &[Complex32]) -> bool {
assert!(baseband.len() >= self.samples_per_bit as usize);
//perform correlation
let positive_energy: Complex32 = self.pos_correlator.iter().zip(baseband.iter())
.enumerate()
.map(|(i, (a, b))| *a * *b * windows::blackmann(i as f32 / self.samples_per_bit as f32)).sum();
let negative_energy: Complex32 = self.neg_correlator.iter().zip(baseband.iter())
.enumerate()
.map(|(i, (a, b))| *a * *b * windows::blackmann(i as f32 / self.samples_per_bit as f32)).sum();
positive_energy.mag().powi(2) < negative_energy.mag().powi(2)
}
}
/*
// FSK Demodulator (dumb non coherent + no symbol timing recovery)
pub struct BFSKDem {
samples_per_bit: u32,
@ -77,7 +124,7 @@ impl BFSKDem {
samples_per_bit,
deviation,
sample_index: 0,
fft: FFT::new(samples_per_bit as usize, windows::bartlett),
fft: FFT::new(samples_per_bit as usize, windows::blackmann),
bin_pos: bin_index as usize,
bin_neg: (samples_per_bit - bin_index - 1) as usize, // -deviation = negative frequency = upper half
}
@ -93,3 +140,4 @@ impl BFSKDem {
positive_energy.mag() > negative_energy.mag()
}
}
*/

8
src/complex.rs
View File

@ -130,11 +130,3 @@ impl Sum for Complex32 {
iter.fold(Complex32::zero(), |acc, el| acc + el)
}
}
impl<T> Deref for Complex<T> {
type Target = Complex<T>;
fn deref(&self) -> &Self::Target {
self
}
}

4
src/fft.rs
View File

@ -7,9 +7,9 @@ pub mod radix2;
use crate::{
complex::Complex32,
fft::{
dft::NaiveDFT, mixed_radix::MixedRadixFFT, rader::RaderFFT, rader2::Rader2FFT,
dft::NaiveDFT, mixed_radix::MixedRadixFFT, rader2::Rader2FFT,
radix2::Radix2FFT,
},
}, windows,
};
pub type FFTWindow = fn(f32) -> f32;

1
src/fft/rader2.rs
View File

@ -34,7 +34,6 @@ impl DFTAlgorithm for Rader2FFT {
let g = compute_prime_primitive_root(size);
let permutations: Box<[usize]> = (0..(size - 1)).map(|i| exp_mod(g, i + 1, size)).collect();
let sub_size = next_pow2(2 * size - 3);
println!("{}", sub_size);
// Compute fourrier transform of twiddle factors
let twiddle_factors = (0..sub_size)

4
src/filtering.rs
View File

@ -1,2 +1,2 @@
mod fir;
mod impulse_response;
pub mod fir;
pub mod impulse_response;

46
src/filtering/fir.rs
View File

@ -4,52 +4,36 @@ use std::collections::VecDeque;
use crate::complex::Complex32;
pub struct FIRFilter<T>
where
T: Iterator<Item = Complex32>,
pub struct FIRFilter
{
size: usize,
normalization: f32,
impulse_response: Box<[Complex32]>,
taps: VecDeque<Complex32>,
input: T,
}
impl<T> FIRFilter<T>
where
T: Iterator<Item = Complex32>,
impl FIRFilter
{
pub fn new(impulse_response: &[Complex32], input: T) -> Self {
pub fn new(impulse_response: &[Complex32]) -> Self {
FIRFilter {
size: impulse_response.len(),
impulse_response: impulse_response.iter().copied().collect(),
//normalization: impulse_response.iter().map(|x| x.mag()).reduce(|acc, e| acc.max(e)).unwrap(),
normalization: impulse_response.iter().copied().sum::<Complex32>().mag(), // DC normalization
// TODO: Maybe we'd want other types of normalization (per frequency)
taps: VecDeque::from(vec![Complex32::zero(); impulse_response.len()]),
input,
}
}
}
impl<T> Iterator for FIRFilter<T>
where
T: Iterator<Item = Complex32>,
{
type Item = Complex32;
fn next(&mut self) -> Option<Self::Item> {
pub fn next(&mut self, next: Complex32) -> Complex32 {
let _ = self.taps.pop_front();
let sample = self.input.next();
if let None = sample {
self.taps.push_back(Complex32::zero());
None
} else {
self.taps.push_back(sample.unwrap());
Some(
self.taps
.iter()
.zip(self.impulse_response.iter())
.map(|(tap, coef)| *tap * *coef)
.sum(),
)
}
self.taps.push_back(next);
self.taps
.iter()
.zip(self.impulse_response.iter())
.map(|(tap, coef)| *tap * *coef)
.sum::<Complex32>() / self.normalization
}
}

51
src/filtering/impulse_response.rs
View File

@ -1,8 +1,51 @@
// Utilities for impulse response design
pub struct WindowIRDesigner {
window: fn(f32) -> f32,
use crate::complex::Complex32;
// Ideal requested transfer function
transfer_function: Box<[Complex32]>,
pub mod design
{
use crate::fft::FFT;
use crate::windows::{self, Window};
use crate::complex::Complex32;
///Designs a impulse response from a desired transfer function using windowing technique
pub fn ir_from_transfer_function(transfer_function: &[Complex32], ir_length: usize, window: Window) -> Vec<Complex32>
{
let tf_len = transfer_function.len();
let mut ifft = FFT::new_inv(tf_len);
// Compute ideal convolution kernel/impulse response
ifft.execute(transfer_function);
// Shorten and window
let mut ir = vec![];
for i in 0..ir_length
{
// Get value within ifft result (centering/trimming)
let k = (tf_len - (ir_length / 2) + i) % tf_len;
// Windowing
ir.push(
ifft.get_output()[k] * window(i as f32 / ir_length as f32) / tf_len as f32
);
}
ir
}
pub fn frequency_response(impulse_response: &[Complex32]) -> Vec<Complex32>
{
let len = impulse_response.len();
let mut fft = FFT::new(len, windows::rectangular);
// Recenter impulse response
let mut centered_ir = vec![];
for i in 0..len
{
let k = (len / 2) + i;
centered_ir.push(impulse_response[k % len]);
}
fft.execute(&centered_ir);
Vec::from(fft.get_output())
}
}

61
src/iq.rs
View File

@ -1,37 +1,50 @@
use crate::{complex::Complex32, nco::Nco};
use std::f32::consts::PI;
pub struct IQSampler<T>
where
T: Iterator<Item = f32>,
use crate::{complex::Complex32, filtering::fir::FIRFilter, map, nco::Nco, windows};
pub struct IQSampler
{
passband: T,
local_oscillator: Nco,
low_pass_i: FIRFilter,
low_pass_q: FIRFilter,
}
impl<T> IQSampler<T>
where
T: Iterator<Item = f32>,
impl IQSampler
{
fn new(passband: T, center_freq: f32) -> Self {
pub fn new(center_freq: f32) -> Self {
// Design a lowpass filter that cuts off at the center freq
// Estimate FIR length :
let fir_length = 50;
// Ideal transfer function :
let mut transfer_function = vec![Complex32::zero(); fir_length];
let bin_id = map(center_freq, 0., PI, 0., transfer_function.len() as f32 / 2.).floor() as usize;
for i in 0..bin_id
{
transfer_function[i] = Complex32::new(1., 0.);
transfer_function[fir_length - 1 - i] = Complex32::new(1., 0.);
}
let ir = crate::filtering::impulse_response::design::ir_from_transfer_function(&transfer_function, fir_length, windows::blackmann);
IQSampler {
passband,
local_oscillator: Nco::new(center_freq),
low_pass_i: FIRFilter::new(&ir),
low_pass_q: FIRFilter::new(&ir)
}
}
fn new_from_lo(passband: T, local_oscillator: Nco) -> Self {
IQSampler {
passband,
local_oscillator,
}
pub fn sample(&mut self, input_sample: f32) -> Complex32
{
let i_mixed = self.local_oscillator.cexp().re * input_sample;
let q_mixed = self.local_oscillator.cexp().im * input_sample;
self.local_oscillator.step();
// TODO: Could use one filter for both I and Q
Complex32::new(
self.low_pass_i.next(Complex32::new(i_mixed, 0.)).re,
self.low_pass_q.next(Complex32::new(q_mixed, 0.)).re
) * 2.
}
}
impl<T> Iterator for IQSampler<T>
where
T: Iterator<Item = f32>,
{
type Item = Complex32;
fn next(&mut self) -> Option<Self::Item> {}
}

236
src/main.rs
View File

@ -3,6 +3,7 @@
use std::{
f32::consts::PI,
fs::File,
i16,
io::{Read, Write},
ops::{Add, Div, Mul, Sub},
};
@ -21,33 +22,22 @@ use complex::Complex;
use complex::Complex32;
use fft::DFTAlgorithm;
use nco::Nco;
use plotters::prelude::*;
use eframe::egui::{self, Color32, debug_text::print};
use egui_plot::{self, Bar, BarChart, Legend, Line, LineStyle, Plot, PlotPoints};
use plotters::style::Color;
use crate::{
bfsk::BFSKDem,
fft::{
FFT, FFTDirection, create_fft,
dft::NaiveDFT,
mixed_radix::MixedRadixFFT,
prime_factors,
rader::{RaderFFT, compute_prime_primitive_root, exp_mod},
rader2::{Rader2FFT, next_pow2},
radix2::Radix2FFT,
windows,
fft::{FFT, dft::NaiveDFT},
filtering::{
fir::FIRFilter,
impulse_response::design::{self, frequency_response, ir_from_transfer_function},
},
iq::IQSampler,
units::frequency::{self, hz_to_rad_per_sample},
};
struct QuickLCG(i32);
impl QuickLCG {
pub fn seed(val: i32) -> QuickLCG {
QuickLCG(val % 10)
}
pub fn next(&mut self) -> i32 {
self.0 = self.0.overflowing_mul(9321).0.overflowing_add(5672).0 % 10;
self.0
}
}
// Utilities
fn map<T>(input: T, in_min: T, in_max: T, out_min: T, out_max: T) -> T
where
@ -58,17 +48,185 @@ where
fn main() {
modulate();
//œtest();
println!("Demodulating");
//demodulate();
let native_options = eframe::NativeOptions::default();
let _ = eframe::run_native("Egui", native_options, Box::new(|cc| Ok(Box::new(EguiApp::new(cc)))));
}
#[derive(Default)]
struct EguiApp {
samples: Vec<f32>,
iq: Vec<Complex32>,
dem: Vec<f32>,
}
const BAUD_RATE: u32 = 1000;
impl EguiApp {
fn new(cc: &eframe::CreationContext<'_>) -> Self {
let mut reader = hound::WavReader::open("audio/modulated.wav").unwrap();
let samples = reader.samples::<i16>();
let input_test = samples
.take(10_000)
.map(|x| x.unwrap() as f32 / i16::MAX as f32)
.collect::<Vec<_>>();
// Modulation parameters
let frequency = 1700.;
let deviation = 500.;
// Data parameters
let sample_rate = 48000;
let baud_rate = BAUD_RATE;
let mut iq_sampler = IQSampler::new(hz_to_rad_per_sample(frequency, sample_rate as f32));
let iq = input_test
.iter()
.map(|x| iq_sampler.sample(*x))
.collect::<Vec<_>>();
let mut dem = BFSKDem::new(
sample_rate / baud_rate,
hz_to_rad_per_sample(deviation, sample_rate as f32),
);
let mut demodulated = vec![];
for x in iq.chunks((sample_rate / baud_rate) as usize) {
let samples = x
.iter()
.copied()
.chain(std::iter::repeat(Complex32::zero()))
.take((sample_rate / baud_rate) as usize)
.collect::<Vec<_>>();
demodulated.push(dem.demod(&samples));
}
Self {
samples: input_test.clone(),
iq,
dem: demodulated
.iter()
.map(|b| if *b { 1. } else { 0. })
.collect(),
}
}
}
impl eframe::App for EguiApp {
fn update(&mut self, ctx: &egui::Context, frame: &mut eframe::Frame) {
egui::CentralPanel::default().show(ctx, |ui| {
ui.heading("Hello World!");
Plot::new("Fourrier transform")
.legend(Legend::default())
.show(ui, |plot_ui| {
plot_ui.line(
Line::new(
"Transfer function",
self.samples
.iter()
.enumerate()
.map(|(i, x)| [i as f64 / self.samples.len() as f64, *x as f64])
.collect::<Vec<_>>(),
)
.width(2.),
);
plot_ui.line(
Line::new(
"RE",
self.iq
.iter()
.enumerate()
.map(|(i, x)| [i as f64 / self.iq.len() as f64, x.re as f64])
.collect::<Vec<_>>(),
)
.width(2.),
);
plot_ui.line(
Line::new(
"IM",
self.iq
.iter()
.enumerate()
.map(|(i, x)| [i as f64 / self.iq.len() as f64, x.im as f64])
.collect::<Vec<_>>(),
)
.width(2.),
);
plot_ui.line(
Line::new(
"Bits",
self.iq
.iter()
.enumerate()
.map(|(i, x)| {
[
i as f64 / self.iq.len() as f64,
self.dem[(self.dem.len() as f32 * i as f32
/ self.iq.len() as f32)
.floor()
as usize]
as f64,
]
})
.collect::<Vec<_>>(),
)
.width(2.),
);
})
});
}
}
fn demodulate() {
let mut reader = hound::WavReader::open("audio/modulated.wav").unwrap();
let samples = reader.samples::<i16>();
// Modulation parameters
let frequency = 1700.;
let deviation = 500.;
// Data parameter
let sample_rate = 48000;
let baud_rate = BAUD_RATE;
let mut iq_sampler = IQSampler::new(hz_to_rad_per_sample(frequency, sample_rate as f32));
let iq_samples = samples
.map(|x| iq_sampler.sample(x.unwrap() as f32 / i16::MAX as f32))
.collect::<Vec<_>>();
let mut dem = BFSKDem::new(
sample_rate / baud_rate,
hz_to_rad_per_sample(deviation, sample_rate as f32),
);
let mut bits = vec![];
for x in iq_samples.chunks((sample_rate / baud_rate) as usize) {
// Zero pad
//let zero_padded = x.iter().copied().chain(std::iter::repeat(Complex32::zero())).take(sample_rate as usize / baud_rate as usize).collect::<Vec<_>>();
//bits.push(dem.demod(&zero_padded));
bits.push(dem.demod(x));
}
assert!(bits.len() % 8 == 0);
let mut out_file = File::create("out.txt").unwrap();
for x in bits.chunks(8) {
out_file.write_all(&[!bits_to_byte(x)]).unwrap();
}
}
fn modulate() {
// Modulation parameters
let frequency = 2000.;
let frequency = 1700.;
let deviation = 500.;
// Data parameters
let sample_rate = 44100;
let baud_rate = 400;
// Data parameter
let sample_rate = 48000;
let baud_rate = BAUD_RATE;
// File to modulate
let f = File::open("s.txt").unwrap();
@ -92,7 +250,7 @@ fn modulate() {
sample_format: hound::SampleFormat::Int,
};
let mut writer = hound::WavWriter::create("sine.wav", spec).unwrap();
let mut writer = hound::WavWriter::create("audio/modulated.wav", spec).unwrap();
for (s, up) in modulator.zip(lo) {
let sample = (s * up).re; // Project to I coords
let amplitude = i16::MAX as f32;
@ -113,3 +271,27 @@ fn byte_to_bits(byte: u8) -> Vec<bool> {
(byte >> 7) & 1 == 1,
]
}
/*
fn bits_to_byte(bits: &[bool]) -> u8 {
bits[7] as u8 |
bits[6] as u8 >> 1 |
bits[5] as u8 >> 2 |
bits[4] as u8 >> 3 |
bits[3] as u8 >> 4 |
bits[2] as u8 >> 5 |
bits[1] as u8 >> 6 |
bits[0] as u8 >> 7
}
*/
fn bits_to_byte(bits: &[bool]) -> u8 {
bits[0] as u8
| (bits[1] as u8) << 1
| (bits[2] as u8) << 2
| (bits[3] as u8) << 3
| (bits[4] as u8) << 4
| (bits[5] as u8) << 5
| (bits[6] as u8) << 6
| (bits[7] as u8) << 7
}

20
src/windows.rs
View File

@ -1,3 +1,7 @@
use std::f32::consts::PI;
pub type Window = fn(f32) -> f32;
pub fn rectangular(_: f32) -> f32 {
1.
}
@ -5,3 +9,19 @@ pub fn rectangular(_: f32) -> f32 {
pub fn bartlett(t: f32) -> f32 {
if t < 0.5 { 2. * t } else { 2. - 2. * t }
}
pub fn hann(t: f32) -> f32
{
0.5 - 0.5 * (2. * PI * t).cos()
}
pub fn hamming(t: f32) -> f32
{
0.54 - 0.46 * (2. * PI * t).cos()
}
pub fn blackmann(t: f32) -> f32
{
let x = 2. * PI * t;
0.45 - 0.5 * x.cos() + 0.08 * (2. * x).cos()
}