Simplify DDS, add plots

5 files changed, 178 insertions, 65 deletions

diff --git a/.gitignore b/.gitignore
index ea8c4bf..cb8a96f 100644
--- a/.gitignore
+++ b/.gitignore
@@ -1 +1,6 @@
 /target
+debug-in.f32
+debug-out.i8
+debug-pd.csv
+sine.dat
+debug-dds.csv
diff --git a/gensine.py b/gensine.py
index 872337e..f8b5161 100755
--- a/gensine.py
+++ b/gensine.py
@@ -1,13 +1,17 @@
+#!/usr/bin/env python
 import matplotlib.pyplot as plt
 import numpy as np
 
+amplitude = 0.8
+
+dac_rate = 10000 #Hz
 sample_rate = 1000 #Hz
 freq = 5 #Hz
-duration = 2# seconds
+duration = 2 #seconds
 
 num_samples = sample_rate * duration
 
-sine_table = (32767 * np.sin(2.0 * np.pi * freq * np.arange(num_samples) / sample_rate)).astype("i2")
+sine_table = (amplitude * 32767 * np.sin(2.0 * np.pi * freq * np.arange(num_samples) / sample_rate)).astype("i2")
 
 sine_table.tofile("sine.dat")
 
@@ -15,4 +19,6 @@ plt.plot(sine_table)
 plt.show()
 
 print("Now run:")
-print("target/debug/fl2k_ampliphase -D -c 2000 -s 10000 -i 1000 -f sine.dat")
+print(f"target/debug/fl2k_ampliphase -D -c 2000 -s {dac_rate} -i {sample_rate} -f sine.dat -w sine")
+print("And then")
+print("./plot.py")
diff --git a/plot.py b/plot.py
new file mode 100755
index 0000000..21219aa
--- /dev/null
+++ b/plot.py
@@ -0,0 +1,55 @@
+#!/usr/bin/env python
+import matplotlib.pyplot as plt
+import numpy as np
+
+debug_pd = np.loadtxt("debug-pd.csv", delimiter=",")
+sample, slope, pd, pdslope = np.split(debug_pd, 4, 1)
+
+debug_dds = np.loadtxt("debug-dds.csv", delimiter=",")
+dds_ix, dds_phase, dds_phase_delta, dds_phase_idx_i, dds_phase_idx_q = np.split(debug_dds, 5, 1)
+
+out = np.fromfile("debug-out.i8", dtype="i1")
+out_r, out_g = np.split(np.reshape(out, newshape=(out.shape[0]//2, 2)), 2, 1)
+
+samp_rate = 10000
+input_rate = 1000
+assert(samp_rate % input_rate == 0)
+rf_to_baseband_sample_ratio = samp_rate // input_rate;
+
+L = 200
+
+L_out = L * rf_to_baseband_sample_ratio
+
+plt.figure()
+plt.subplot(3, 1, 1)
+plt.title("sample")
+plt.plot(sample[0:L])
+
+plt.subplot(3, 1, 2)
+plt.title("pd")
+plt.plot(pd[0:L])
+
+plt.subplot(3, 1, 3)
+plt.title("pdslope")
+plt.plot(pdslope[0:L])
+
+
+plt.figure()
+plt.subplot(4, 1, 1)
+plt.title("dds ix")
+plt.plot(dds_ix[0:L_out])
+
+plt.subplot(4, 1, 2)
+plt.title("dds phase")
+plt.plot(dds_phase[0:L_out])
+
+plt.subplot(4, 1, 3)
+plt.title("dds phase_delta")
+plt.plot(dds_phase_delta[0:L_out])
+
+plt.subplot(4, 1, 4)
+plt.title("output")
+plt.plot(out_r[0:L_out])
+plt.plot(out_g[0:L_out])
+
+plt.show()
diff --git a/src/fl2k.rs b/src/fl2k.rs
index f2e197b..74357b3 100644
--- a/src/fl2k.rs
+++ b/src/fl2k.rs
@@ -69,6 +69,7 @@ extern "C" fn tx_callback(data_info: *mut fl2k_ampliphase::fl2k_data_info_t) {
                 },
                 Err(_) => {
                     (*ctx).buffer_underflows += 1;
+                    eprintln!("FL2K buffers not ready on callback");
                 },
             }
         }
diff --git a/src/main.rs b/src/main.rs
index 7d052bc..1c6f520 100644
--- a/src/main.rs
+++ b/src/main.rs
@@ -41,11 +41,35 @@ const ANG_INCR: f32 = INT32_MAX_AS_FLOAT / (2.0 * PI);
 
 enum Waveform { Sine, Rect }
 
+#[derive(Clone, Copy)]
 enum Output { Debug, FL2K }
 
+struct BufferDumper<T> {
+    phantom : std::marker::PhantomData<T>,
+    writer : BufWriter<File>,
+}
+
+impl<T> BufferDumper<T> {
+    fn new(filename: &str) -> Self {
+        let writer = BufWriter::new(File::create(filename).expect("create file"));
+        BufferDumper { phantom: std::marker::PhantomData, writer }
+    }
+
+    fn write_buf(&mut self, buf: &[T]) -> std::io::Result<()> {
+        let buf_u8: &[u8] = unsafe {
+            std::slice::from_raw_parts(
+                buf.as_ptr() as *const u8,
+                buf.len() * std::mem::size_of::<T>()
+            )
+        };
+
+        self.writer.write_all(buf_u8)
+    }
+}
+
 struct DDS {
-    trig_table_quadrature: Vec<i16>,
-    trig_table_inphase: Vec<i16>,
+    trig_table_quadrature: Vec<i8>,
+    trig_table_inphase: Vec<i8>,
 
     /* instantaneous phase */
     phase: u32,
@@ -55,28 +79,26 @@ struct DDS {
     /* 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 {
+    fn init(samp_rate: f32, freq: f32, phase: 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;
+            let i = f32::cos(incr * i as f32 * 2.0 * PI) * 255.0;
+            let q = f32::sin(incr * i as f32 * 2.0 * PI) * 255.0;
 
             match waveform {
                 Waveform::Sine => {
-                    trig_table_inphase.push(f32::round(i) as i16);
-                    trig_table_quadrature.push(f32::round(q) as i16);
+                    trig_table_inphase.push(f32::round(i) as i8);
+                    trig_table_quadrature.push(f32::round(q) as i8);
                 }
                 Waveform::Rect => {
-                    trig_table_inphase.push(if i >= 0.0 { 32767 } else { -32767 });
-                    trig_table_quadrature.push(if q >= 0.0 { 32767 } else { -32767 });
+                    trig_table_inphase.push(if i >= 0.0 { 127 } else { -127 });
+                    trig_table_quadrature.push(if q >= 0.0 { 127 } else { -127 });
                 }
             }
         }
@@ -90,41 +112,7 @@ impl DDS {
             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>, Vec<i8>) {
-        let mut out_i = Vec::with_capacity(len);
-        let mut out_q = Vec::with_capacity(len);
-        for _ in 0..len {
-            // get current carrier phase, add phase mod, calculate table index
-            let phase_idx_i = self.phase.overflowing_sub(self.phase_delta as u32).0 >> TRIG_TABLE_SHIFT;
-            let phase_idx_q = self.phase.overflowing_add(self.phase_delta as u32).0 >> TRIG_TABLE_SHIFT;
-
-            if phase_idx_q > 255 || phase_idx_i > 255 {
-                panic!("Phase IDX out of bounds");
-            }
-
-            self.phase = self.phase.overflowing_add(self.phase_step).0;
-
-            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_i.push(i);
-            out_q.push(q);
-
-            self.phase_delta += self.phase_slope;
         }
-        (out_i, out_q)
     }
 }
 
@@ -184,7 +172,7 @@ fn main() {
         None => 1_440_000.0,
     };
 
-    let input_freq: u32 = match cli_args.opt_str("i") {
+    let input_rate: u32 = match cli_args.opt_str("i") {
         Some(s) => s.parse().expect("integer value"),
         None => 48_000,
     };
@@ -214,13 +202,14 @@ fn main() {
         }
     };
 
-    if samp_rate % input_freq != 0 {
-        eprintln!("WARNING: input freq does not divide sample rate.");
+    if samp_rate % input_rate != 0 {
+        eprintln!("WARNING: input_rate freq does not divide sample rate.");
     }
-    let rf_to_baseband_sample_ratio = samp_rate / input_freq;
+    let rf_to_baseband_sample_ratio = samp_rate / input_rate;
 
-    eprintln!("Samplerate:       {} MHz", (samp_rate as f32)/1000000.0);
-    eprintln!("Center frequency: {} kHz", base_freq/1000.0);
+    eprintln!("Input rate:       {} kHz", (input_rate as f32)/1e3);
+    eprintln!("Samplerate:       {} MHz", (samp_rate as f32)/1e6);
+    eprintln!("Center frequency: {} kHz", base_freq/1e3);
 
     let running = Arc::new(AtomicBool::new(true));
 
@@ -240,6 +229,11 @@ fn main() {
 
     // Read file and convert samples
     thread::spawn(move || {
+        let mut in_debug = match output {
+            Output::Debug => Some(BufferDumper::<f32>::new("debug-in.f32")),
+            Output::FL2K => None,
+        };
+
         while running.load(Ordering::SeqCst) {
             let mut buf = Vec::with_capacity(BASEBAND_BUF_SAMPS);
             buf.resize(BASEBAND_BUF_SAMPS, 0i16);
@@ -251,13 +245,27 @@ fn main() {
                 )
             };
 
-            source_file.read_exact(&mut buf_u8).expect("Read from source file");
+            if let Err(e) = source_file.read_exact(&mut buf_u8) {
+                if e.kind() != std::io::ErrorKind::UnexpectedEof {
+                    eprintln!("Input file read stopped {:?}", e);
+                }
+                else {
+                    eprintln!("Input file EOF reached");
+                }
+                break;
+            }
 
             let buf: Vec<f32> = buf
                 .iter()
                 .map(|v| (v/2 + i16::MAX/2) as f32 / 32768.0)
                 .collect();
 
+            if let Some(w) = &mut in_debug {
+                if let Err(e) = w.write_buf(&buf) {
+                    eprintln!("Error writing debug-in.f32: {}", e);
+                }
+            }
+
             if let Err(_) = input_samples_tx.send(buf) {
                 eprintln!("Quit read thread");
                 break;
@@ -268,6 +276,12 @@ fn main() {
     // Read samples, calculate PD and PDSLOPE
     thread::spawn(move || {
         let mut lastamp = 0f32;
+
+        let mut debug_writer = match output {
+            Output::Debug => Some(BufWriter::new(File::create("debug-pd.csv").expect("create file"))),
+            Output::FL2K => None,
+        };
+
         loop {
             let Ok(buf) = input_samples_rx.recv() else { break };
 
@@ -281,8 +295,7 @@ fn main() {
                efficient and pretty good interpolation filter. */
 
             for sample in buf {
-                let slope = sample - lastamp;
-                let slope = slope / rf_to_baseband_sample_ratio as f32;
+                let slope = (sample - lastamp) / rf_to_baseband_sample_ratio as f32;
 
                 let pd = lastamp * modulation_index * INT32_MAX_AS_FLOAT;
 
@@ -298,6 +311,11 @@ fn main() {
                     panic!("pdslope out of bounds {}", pdslope);
                 }
 
+                if let Some(w) = &mut debug_writer {
+                    writeln!(w, "{},{},{},{}", sample, slope, pd, pdslope)
+                        .expect("write debug-pd.csv");
+                }
+
                 pd_buf.push((pd as i32, pdslope as i32));
 
                 lastamp = sample;
@@ -312,16 +330,45 @@ fn main() {
 
     // 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);
+        let mut dds = DDS::init(samp_rate as f32, base_freq, 0.0, waveform);
+
+        let mut debug_writer = match output {
+            Output::Debug => Some(BufWriter::new(File::create("debug-dds.csv").expect("create file"))),
+            Output::FL2K => None,
+        };
+
         loop {
             let Ok(pd_buf) = pd_rx.recv() else { break };
 
             for (pd, pdslope) in pd_buf {
-                dds.set_phase(pd, pdslope);
-
-                let iq_bufs = dds.generate_iq(rf_to_baseband_sample_ratio as usize);
+                dds.phase_delta = pd;
+                dds.phase_slope = pdslope;
+
+                let len = rf_to_baseband_sample_ratio as usize;
+                let mut out_i = Vec::with_capacity(len);
+                let mut out_q = Vec::with_capacity(len);
+                for ix in 0..len {
+                    // get current carrier phase, add phase mod, calculate table index
+                    let phase_idx_i = dds.phase.overflowing_sub(dds.phase_delta as u32).0 >> TRIG_TABLE_SHIFT;
+                    let phase_idx_q = dds.phase.overflowing_add(dds.phase_delta as u32).0 >> TRIG_TABLE_SHIFT;
+
+                    if phase_idx_q > 255 || phase_idx_i > 255 {
+                        panic!("Phase IDX out of bounds");
+                    }
+
+                    out_i.push(dds.trig_table_inphase[phase_idx_i as usize]);
+                    out_q.push(dds.trig_table_quadrature[phase_idx_q as usize]);
+
+                    if let Some(w) = &mut debug_writer {
+                        writeln!(w, "{},{},{},{},{}", ix, dds.phase, dds.phase_delta, phase_idx_i, phase_idx_q)
+                            .expect("write debug-dds.csv");
+                    }
+
+                    dds.phase = dds.phase.overflowing_add(dds.phase_step).0;
+                    dds.phase_delta += dds.phase_slope;
+                }
 
-                if let Err(_) = iq_tx.send(iq_bufs) {
+                if let Err(_) = iq_tx.send((out_i, out_q)) {
                     eprintln!("Quit dds thread");
                     break;
                 }
@@ -347,8 +394,7 @@ fn main() {
             fl2k.stop_tx().expect("stop tx");
         }
         Output::Debug => {
-            let out_file = File::create("debug-out.i8").expect("create file");
-            let mut out_file = BufWriter::new(out_file);
+            let mut out_file = BufferDumper::new("debug-out.i8");
             loop {
                 let Ok((i_buf, q_buf)) = iq_rx.recv() else { break };
 
@@ -369,7 +415,7 @@ fn main() {
                     )
                 };
 
-                if let Err(e) = out_file.write_all(buf_u8) {
+                if let Err(e) = out_file.write_buf(buf_u8) {
                     eprintln!("Write output error: {}", e);
                     break;
                 }