From ba0aeee005d5fdaaab59cd7c099a237f51eddc86 Mon Sep 17 00:00:00 2001 From: "Matthias P. Braendli" Date: Sat, 31 Dec 2022 14:36:18 +0100 Subject: Create project --- src/main.rs | 346 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 346 insertions(+) create mode 100644 src/main.rs (limited to 'src/main.rs') diff --git a/src/main.rs b/src/main.rs new file mode 100644 index 0000000..4bdb879 --- /dev/null +++ b/src/main.rs @@ -0,0 +1,346 @@ +/* + * Copyright (C) 2022 by Matthias P. Braendli + * + * Based on previous work by + * Copyright (C) 2022 by Felix Erckenbrecht + * Copyright (C) 2016-2018 by Steve Markgraf + * Copyright (C) 2009 by Bartek Kania + * + * SPDX-License-Identifier: GPL-2.0+ + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 2 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see . + */ + +use std::{env, thread}; +use std::io::{prelude::*, BufReader, BufWriter}; +use std::fs::File; +use std::sync::mpsc; +use getopts::Options; +use num_complex::Complex; + +//mod fl2k; + +const TRIG_TABLE_ORDER : usize = 8; +const TRIG_TABLE_SHIFT : usize = 32 - TRIG_TABLE_ORDER; +const TRIG_TABLE_LEN : usize = 1 << TRIG_TABLE_ORDER; + +const PI : f32 = std::f32::consts::PI; +const INT32_MAX_AS_FLOAT : f32 = 0x1_0000_0000u64 as f32; +const ANG_INCR : f32 = INT32_MAX_AS_FLOAT / (2.0 * PI); + +enum Waveform { Sine, Rect } + +enum Output { Debug } + +struct DDS { + trig_table_quadrature : Vec, + trig_table_inphase : Vec, + + samp_rate : f32, + freq : f32, + + /* instantaneous phase */ + phase : u32, + /* phase increment */ + phase_step : u32, + + /* for phase modulation */ + phase_delta : i32, + phase_slope : i32, + + amplitude : f32, +} + +impl DDS { + fn init(samp_rate : f32, freq : f32, phase : f32, amp : f32, waveform : Waveform) -> Self { + let mut trig_table_inphase = Vec::with_capacity(TRIG_TABLE_LEN); + let mut trig_table_quadrature = Vec::with_capacity(TRIG_TABLE_LEN); + + let incr = 1.0f32 / TRIG_TABLE_LEN as f32; + for i in 0..TRIG_TABLE_LEN { + let i = f32::cos(incr * i as f32 * 2.0 * PI) * 32767.0; + let q = f32::sin(incr * i as f32 * 2.0 * PI) * 32767.8; + + match waveform { + Waveform::Sine => { + trig_table_inphase.push(f32::round(i) as i16); + trig_table_quadrature.push(f32::round(q) as i16); + } + Waveform::Rect => { + trig_table_inphase.push(if i >= 0.0 { 32767 } else { -32767 }); + trig_table_quadrature.push(if q >= 0.0 { 32767 } else { -32767 }); + } + } + } + + let phase_step = (freq / samp_rate) * 2.0 * PI * ANG_INCR; + + DDS { + trig_table_quadrature, + trig_table_inphase, + samp_rate, + freq, + phase: f32::round(phase * ANG_INCR) as u32, + phase_step: f32::round(phase_step) as u32, + phase_delta: 0, + phase_slope: 0, + amplitude: amp, + } + } + + fn set_phase(&mut self, phase_delta : i32, phase_slope : i32) { + self.phase_delta = phase_delta; + self.phase_slope = phase_slope; + } + + fn generate_iq(&mut self, len : usize) -> Vec<(i8, i8)> { + let mut out = Vec::with_capacity(len); + for _ in 0..len { + let phase = self.phase as i32; + // get current carrier phase, add phase mod, calculate table index + let phase_idx_i = (phase - self.phase_delta) >> TRIG_TABLE_SHIFT; + let phase_idx_q = (phase + self.phase_delta) >> TRIG_TABLE_SHIFT; + + if phase_idx_q > 255 || phase_idx_i > 255 { + panic!("Phase IDX out of bounds"); + } + + self.phase = phase as u32 + self.phase_step; + + let amp = (self.amplitude * 32767.0) as i32; // 0..15 + let amp_i = amp * self.trig_table_inphase[phase_idx_i as usize] as i32; // 0..31 + let amp_q = amp * self.trig_table_quadrature[phase_idx_q as usize] as i32; // 0..31 + // + let i = (amp_i >> 24) as i8; // 0..31 >> 24 => 0..8 + let q = (amp_q >> 24) as i8; // 0..31 >> 24 => 0..8 + out.push((i, q)); + + self.phase_delta += self.phase_slope; + } + out + } +} + +fn print_usage(program: &str, opts: Options) { + let brief = format!("Usage: {} FILE [options]", program); + eprint!("{}", opts.usage(&brief)); +} + +fn main() { + let args: Vec = env::args().collect(); + let program = args[0].clone(); + + let mut opts = Options::new(); + opts.optopt("f", "file", "Input file, containing signed 16-bit samples.", "FILE"); + opts.optopt("d", "device-index", "Select device index", "DEVINDEX"); + opts.optopt("c", "center-freq", "Center frequency in Hz (default: 1440 kHz)", "FREQ"); + opts.optopt("s", "samplerate", "Samplerate in Hz (default: 96 MS/s)", "RATE"); + opts.optopt("m", "mod-index", "Modulation index (default: 0.25)]", "FACTOR"); + opts.optopt("i", "input-rate", "Input baseband sample rate (default: 48000 Hz)", "RATE"); + opts.optopt("w", "waveform", "(sine|rect) default: rect", "WAVEFORM"); + opts.optflag("D", "debug", "Write to debug files instead of FL2K"); + opts.optflag("h", "help", "print this help menu"); + let cli_args = match opts.parse(&args[1..]) { + Ok(m) => { m } + Err(f) => { panic!("{}", f.to_string()) } + }; + if cli_args.opt_present("h") { + print_usage(&program, opts); + std::process::exit(1); + } + + let output = if cli_args.opt_str("D").is_none() { + eprintln!("Only debug supported for now"); + std::process::exit(1); + } + else { + Output::Debug + }; + + let samp_rate : u32 = match cli_args.opt_str("s") { + Some(s) => s.parse().expect("floating point value"), + None => 96_000_000, + }; + + let base_freq = match cli_args.opt_str("c") { + Some(s) => s.parse().expect("floating point value"), + None => 1_440_000.0, + }; + + let input_freq : u32 = match cli_args.opt_str("i") { + Some(s) => s.parse().expect("floating point value"), + None => 48_000, + }; + + let modulation_index = match cli_args.opt_str("i") { + Some(s) => s.parse().expect("floating point value"), + None => 0.25, + }; + + let waveform = match cli_args.opt_str("w") { + None => Waveform::Rect, + Some(w) if w == "sine" => Waveform::Sine, + Some(w) if w == "rect" => Waveform::Rect, + _ => { + eprintln!("Waveform must be 'sine' or 'rect'"); + print_usage(&program, opts); + std::process::exit(1); + } + }; + + //eprintln!("Device count {}", fl2k::get_device_count()); + + let source_file_name = match cli_args.opt_str("f") { + Some(f) => f, + None => { + eprintln!("Specify input file!"); + print_usage(&program, opts); + std::process::exit(1); + } + }; + + if samp_rate % input_freq != 0 { + eprintln!("WARNING: input freq does not divide sample rate."); + } + let rf_to_baseband_sample_ratio = samp_rate / input_freq; + + eprintln!("Samplerate: {} MHz", (samp_rate as f32)/1000000.0); + eprintln!("Center frequency: {} kHz", base_freq/1000.0); + + + let (input_samples_tx, input_samples_rx) = mpsc::sync_channel::>(2); + let (pd_tx, pd_rx) = mpsc::sync_channel(2); + let (iq_tx, iq_rx) = mpsc::sync_channel(2); + + let source_file = File::open(source_file_name).expect("open file"); + let mut source_file = BufReader::new(source_file); + + const BASEBAND_BUF_SAMPS : usize = 1024; + + // Read file and convert samples + thread::spawn(move || { + loop { + let mut buf = Vec::with_capacity(BASEBAND_BUF_SAMPS); + buf.resize(BASEBAND_BUF_SAMPS, 0i16); + + let mut buf_u8: &mut [u8] = unsafe { + std::slice::from_raw_parts_mut( + buf.as_mut_ptr() as *mut u8, + buf.len() * std::mem::size_of::() + ) + }; + + source_file.read_exact(&mut buf_u8).expect("Read from source file"); + + let buf : Vec = buf + .iter() + .map(|v| (v/2 + i16::MAX) as f32 / 32768.0) + .collect(); + + if let Err(_) = input_samples_tx.send(buf) { + eprintln!("Quit read thread"); + break; + } + } + }); + + // Read samples, calculate PD and PDSLOPE + thread::spawn(move || { + let mut lastamp = 0f32; + loop { + let Ok(buf) = input_samples_rx.recv() else { break }; + + let mut pd_buf = Vec::with_capacity(buf.len()); + + /* What we do here is calculate a linear "slope" from + the previous sample to this one. This is then used by + the modulator to gently increase/decrease the phase + with each sample without the need to recalculate + the dds parameters. In fact this gives us a very + efficient and pretty good interpolation filter. */ + + for sample in buf { + let slope = sample - lastamp; + let slope = slope * 1.0 / rf_to_baseband_sample_ratio as f32; + + let pd = lastamp * modulation_index * INT32_MAX_AS_FLOAT; + + const MIN_VAL : f32 = std::i32::MIN as f32; + const MAX_VAL : f32 = std::i32::MAX as f32; + + if pd < MIN_VAL || pd > MAX_VAL { + panic!("pd out of bounds {}", pd); + } + + let pdslope = slope * modulation_index * INT32_MAX_AS_FLOAT; + if pdslope < MIN_VAL || pdslope > MAX_VAL { + panic!("pdslope out of bounds {}", pdslope); + } + + pd_buf.push((pd as i32, pdslope as i32)); + + lastamp = sample; + } + + if let Err(_) = pd_tx.send(pd_buf) { + eprintln!("Quit pd thread"); + break; + } + } + }); + + // Read PD and PDSLOPE, interpolate to higher rate + thread::spawn(move || { + let mut dds = DDS::init(samp_rate as f32, base_freq, 0.0, 1.0, waveform); + loop { + let Ok(buf) = pd_rx.recv() else { break }; + + for (pd, pdslope) in buf { + dds.set_phase(pd, pdslope); + + let iq_buf = dds.generate_iq(rf_to_baseband_sample_ratio as usize); + + if let Err(_) = iq_tx.send(iq_buf) { + eprintln!("Quit dds thread"); + break; + } + } + } + }); + + // Main thread, output to file/device + match output { + Output::Debug => { + let out_file = File::create("debug-out.i8").expect("create file"); + let mut out_file = BufWriter::new(out_file); + loop { + let Ok(buf) = iq_rx.recv() else { break }; + + let buf_u8: &[u8] = unsafe { + std::slice::from_raw_parts( + buf.as_ptr() as *const u8, + buf.len() + ) + }; + + if let Err(e) = out_file.write_all(buf_u8) { + eprintln!("Write output error: {}", e); + break; + } + } + } + } + + eprintln!("Leaving main thread"); +} -- cgit v1.2.3