path: root/src/main.rs

/*
 * Copyright (C) 2022 by Matthias P. Braendli <matthias.braendli@mpb.li>
 *
 * Based on previous work by
 * Copyright (C) 2022 by Felix Erckenbrecht <eligs@eligs.de>
 * Copyright (C) 2016-2018 by Steve Markgraf <steve@steve-m.de>
 * Copyright (C) 2009 by Bartek Kania <mbk@gnarf.org>
 *
 * 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 <http://www.gnu.org/licenses/>.
 */

use std::{env, thread};
use std::io::{prelude::*, BufReader, BufWriter};
use std::fs::File;
use std::sync::mpsc;
use getopts::Options;

//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<i16>,
    trig_table_inphase : Vec<i16>,

    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<String> = 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::<Vec<f32>>(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::<i16>()
                )
            };

            source_file.read_exact(&mut buf_u8).expect("Read from source file");

            let buf : Vec<f32> = 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");
}