Apply asin to input samplesHEADmain

2 files changed, 44 insertions, 8 deletions

diff --git a/plot.py b/plot.py
index 21219aa..9ff3b2a 100755
--- a/plot.py
+++ b/plot.py
@@ -1,6 +1,26 @@
 #!/usr/bin/env python
+import sys
 import matplotlib.pyplot as plt
 import numpy as np
+from scipy.signal import butter, lfilter, freqz
+import matplotlib.pyplot as plt
+
+samp_rate = 10000
+input_rate = 1000
+assert(samp_rate % input_rate == 0)
+rf_to_baseband_sample_ratio = samp_rate // input_rate;
+
+L = 200
+
+def butter_lowpass(cutoff, fs, order=5):
+    nyq = 0.5 * fs
+    normal_cutoff = cutoff / nyq
+    b, a = butter(order, normal_cutoff, btype='low', analog=False)
+    return b, a
+
+def butter_lowpass_filter(data, cutoff, fs, order=6):
+    b, a = butter_lowpass(cutoff, fs, order=order)
+    return lfilter(b, a, data)
 
 debug_pd = np.loadtxt("debug-pd.csv", delimiter=",")
 sample, slope, pd, pdslope = np.split(debug_pd, 4, 1)
@@ -11,12 +31,6 @@ dds_ix, dds_phase, dds_phase_delta, dds_phase_idx_i, dds_phase_idx_q = np.split(
 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
 
@@ -52,4 +66,16 @@ plt.title("output")
 plt.plot(out_r[0:L_out])
 plt.plot(out_g[0:L_out])
 
+plt.figure()
+plt.subplot(2, 1, 1)
+plt.title("out sum")
+out_sum = out_r[0:L_out].astype(float) + out_g[0:L_out].astype(float)
+
+plt.plot(out_sum)
+
+plt.subplot(2, 1, 2)
+plt.title("out sum filt")
+
+plt.plot(butter_lowpass_filter(out_sum, cutoff=10, fs=2000))
+
 plt.show()
diff --git a/src/main.rs b/src/main.rs
index e35062e..8519158 100644
--- a/src/main.rs
+++ b/src/main.rs
@@ -270,10 +270,20 @@ fn main() {
                 }
             }
 
-            // Downmix stereo to mono, apply (v/2 + i16::MAX/2), and normalise to [-1.0; 1.0]
+            // Downmix stereo to mono, apply (v/2 + i16::MAX/2),
+            // normalise to [0.0; 1.0],
+            // take the arcsine to compensate for the phasing,
+            // and divide by PI/2 to normalise to [0.0; 1.0]
+            use std::f32::consts::FRAC_2_PI;
             let buf: Vec<f32> = buf
                 .chunks_exact(2)
-                .map(|v| ((v[0]/2 + v[1]/2)/2 + i16::MAX/2) as f32 / 32768.0)
+                .map(|v| {
+                    let samp = ((v[0]/2 + v[1]/2)/2 + i16::MAX/2) as f32 / 32768.0;
+                    if samp < -1.0 || samp > 1.0 {
+                        panic!("sample normalistion broken")
+                    }
+                    f32::asin(samp) * FRAC_2_PI
+                })
                 .collect();
 
             if let Some(w) = &mut in_debug {