diff --git a/src/bfsk.rs b/src/bfsk.rs index c288f66..18c11e6 100644 --- a/src/bfsk.rs +++ b/src/bfsk.rs @@ -3,7 +3,7 @@ use std::f32::consts::PI; use crate::complex::{Complex, Complex32}; -use crate::fft::{self, DFT, windows}; +use crate::fft::{self, windows}; use crate::map; use crate::nco::Nco; @@ -57,7 +57,7 @@ pub struct BFSKDem { // State sample_index: u32, - fft: Box, + //fft: Box, } impl BFSKDem { @@ -66,19 +66,21 @@ impl BFSKDem { samples_per_bit, deviation, sample_index: 0, - fft: fft::create_fft(samples_per_bit as usize, fft::FFTDirection::Forward), + //fft: fft::create_fft(samples_per_bit as usize, fft::FFTDirection::Forward), } } pub fn demod(&mut self, baseband: &[Complex32]) -> bool { assert!(baseband.len() >= self.samples_per_bit as usize); + /* self.fft .get_input() .iter_mut() .enumerate() .for_each(|(i, x)| *x = baseband[i]); self.fft.execute(windows::rectanguar); + */ let bin_id = map( self.deviation, @@ -91,6 +93,7 @@ impl BFSKDem { 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 { @@ -108,5 +111,7 @@ impl BFSKDem { } return positive_energy < negative_energy; + */ + false } } diff --git a/src/fft.rs b/src/fft.rs index e808e8a..b9aa676 100644 --- a/src/fft.rs +++ b/src/fft.rs @@ -12,6 +12,8 @@ use crate::{ fft::{dft::NaiveDFT, mixed_radix::MixedRadixFFT, rader::RaderFFT, radix2::Radix2FFT}, }; +pub type FFTWindow = fn(f32) -> f32; + #[derive(Copy, Clone)] pub enum FFTDirection { Forward, @@ -27,19 +29,17 @@ impl FFTDirection { } } -pub trait DFT { +pub trait DFTAlgorithm { fn create(size: usize, direction: FFTDirection) -> Self where Self: Sized; - fn execute(&mut self, input: &[Complex32], output: &mut [Complex32], window: fn(f32) -> f32); + fn execute(&mut self, input: &[Complex32]); + + fn get_output(&self) -> &[Complex32]; } -pub trait DFTWindow { - fn eval(t: f32) -> f32; -} - -pub fn create_fft(size: usize, direction: FFTDirection) -> Box { +fn create_fft(size: usize, direction: FFTDirection) -> Box { if size <= 16 { //println!("Naive {size}"); return Box::new(NaiveDFT::create(size, direction)); @@ -58,6 +58,56 @@ pub fn create_fft(size: usize, direction: FFTDirection) -> Box { Box::new(MixedRadixFFT::create(size, direction)) } +pub struct FFT +{ + fft: Box, + size: usize, + window: FFTWindow, + input_buffer: Box<[Complex32]> +} + +impl FFT +{ + pub fn new(size: usize, window: FFTWindow) -> FFT + { + FFT + { + fft: create_fft(size, FFTDirection::Forward), + window, + size, + input_buffer: vec![Complex32::zero(); size].into(), + } + } + + pub fn new_inv(size: usize) -> FFT + { + FFT + { + fft: create_fft(size, FFTDirection::Inverse), + window: windows::rectangular, + size, + input_buffer: vec![Complex32::zero(); size].into(), + } + } + + pub fn execute(&mut self, input: &[Complex32]) + { + self.input_buffer.iter_mut().zip(input.iter()) + .enumerate() + .for_each(|(i, (x, y))| + { + *x = *y * (self.window)(i as f32 / self.size as f32) + }); + + self.fft.execute(&self.input_buffer); + } + + pub fn get_output(&self) -> &[Complex32] + { + self.fft.get_output() + } +} + // Utilities pub fn prime_factors(n: usize) -> Vec { let mut factors = vec![]; diff --git a/src/fft/dft.rs b/src/fft/dft.rs index 6965caf..5c58cd3 100644 --- a/src/fft/dft.rs +++ b/src/fft/dft.rs @@ -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 + } } diff --git a/src/fft/mixed_radix.rs b/src/fft/mixed_radix.rs index dece321..dc42f13 100644 --- a/src/fft/mixed_radix.rs +++ b/src/fft/mixed_radix.rs @@ -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, - pfft: Box, + qfft: Box, + pfft: Box, 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 } } diff --git a/src/fft/rader.rs b/src/fft/rader.rs index 0ffd6c8..e6bce41 100644 --- a/src/fft/rader.rs +++ b/src/fft/rader.rs @@ -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, - conv_fft: Box, + convolution_operand: Box<[Complex32]>, + + convolution_ifft: Box, + convolution_fft: Box, + + 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 = (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::>(); + 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; } diff --git a/src/fft/radix2.rs b/src/fft/radix2.rs index 6570fca..c900008 100644 --- a/src/fft/radix2.rs +++ b/src/fft/radix2.rs @@ -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 diff --git a/src/main.rs b/src/main.rs index 1ab4f8b..4d8e985 100644 --- a/src/main.rs +++ b/src/main.rs @@ -1,3 +1,5 @@ +#![allow(dead_code)] + use std::{ f32::consts::PI, fs::File, @@ -13,11 +15,11 @@ mod nco; use bfsk::BFSKMod; use complex::Complex; use complex::Complex32; -use fft::rader; use nco::Nco; use plotters::prelude::*; +use fft::DFTAlgorithm; -use crate::bfsk::BFSKDem; +use crate::{bfsk::BFSKDem, fft::{dft::NaiveDFT, mixed_radix::MixedRadixFFT, rader::RaderFFT, radix2::Radix2FFT, windows, FFT}}; // Utilities fn map(input: T, in_min: T, in_max: T, out_min: T, out_max: T) -> T @@ -27,6 +29,7 @@ where ((input - in_min.clone()) / (in_max - in_min)) * (out_max - out_min.clone()) + out_min } + fn main() { modulate(); } diff --git a/src/nco.rs b/src/nco.rs index 5c8856a..7ad39f4 100644 --- a/src/nco.rs +++ b/src/nco.rs @@ -1,6 +1,6 @@ // Numerically controlled oscillator -use crate::complex::Complex; +use crate::complex::{Complex, Complex32}; use std::f32::consts::PI; use std::ops::{Add, Div, Mul, Sub}; @@ -97,3 +97,14 @@ impl Nco { )) } } + +impl Iterator for Nco +{ + type Item = Complex32; + + fn next(&mut self) -> Option { + let val = self.cexp(); + self.step(); + Some(val) + } +}