Fix lag calculator and add test

4 files changed, 71 insertions, 461 deletions

diff --git a/analyze_dab_nonlinearity.ipynb b/analyze_dab_nonlinearity.ipynb
deleted file mode 100644
index 0bd874a..0000000
--- a/analyze_dab_nonlinearity.ipynb
+++ /dev/null
@@ -1,451 +0,0 @@
-{
- "cells": [
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {
-    "deletable": true,
-    "editable": true
-   },
-   "outputs": [],
-   "source": [
-    "%matplotlib inline\n",
-    "import matplotlib.pyplot as plt\n",
-    "import numpy as np\n",
-    "import src.dab_util as du"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {
-    "deletable": true,
-    "editable": true
-   },
-   "source": [
-    "# Align Original / Recored Files\n",
-    "\n",
-    "Using Txgain 81"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {
-    "deletable": true,
-    "editable": true
-   },
-   "outputs": [],
-   "source": [
-    "#path = \"../test_dat/out.iq\"\n",
-    "path = \"./recored.dat\"\n",
-    "frame = 96*8192\n",
-    "offset = 5 * frame + 8192\n",
-    "length  = 2 * frame"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {
-    "deletable": true,
-    "editable": true
-   },
-   "outputs": [],
-   "source": [
-    "sig_rec = np.fromfile(path, dtype=np.complex64, count=offset + length)[offset:]"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {
-    "deletable": true,
-    "editable": true
-   },
-   "outputs": [],
-   "source": [
-    "#sig_rec = sig_rec[30370+0*frame:30370+1*frame]"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {
-    "deletable": true,
-    "editable": true
-   },
-   "outputs": [],
-   "source": [
-    "i = 0\n",
-    "l = length\n",
-    "plt.plot(np.abs(sig_rec[i:i + l]))"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {
-    "deletable": true,
-    "editable": true
-   },
-   "outputs": [],
-   "source": [
-    "n = sig_rec.shape[0]\n",
-    "max_mean = np.mean(np.sort(np.abs(sig_rec))[int(n*0.90):])"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {
-    "deletable": true,
-    "editable": true
-   },
-   "outputs": [],
-   "source": []
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {
-    "deletable": true,
-    "editable": true
-   },
-   "outputs": [],
-   "source": [
-    "sym = int(frame / 76.0)\n",
-    "mask = np.ones((int(4*sym)))\n",
-    "mask[int(2*sym) : int(3*sym)] = -1"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {
-    "deletable": true,
-    "editable": true
-   },
-   "outputs": [],
-   "source": [
-    "import scipy.signal as ss"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {
-    "deletable": true,
-    "editable": true
-   },
-   "outputs": [],
-   "source": [
-    "c = ss.correlate(np.abs(sig_rec), mask, mode='same')"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {
-    "deletable": true,
-    "editable": true
-   },
-   "outputs": [],
-   "source": [
-    "max1 = np.argmax(c[:frame])\n",
-    "max2 = frame + np.argmax(c[frame:])\n",
-    "max2 - max1 - frame < sym / 100"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {
-    "deletable": true,
-    "editable": true
-   },
-   "outputs": [],
-   "source": [
-    "plt.plot(c)\n",
-    "plt.plot((max1, max1), (0, 1200))\n",
-    "plt.plot((max2, max2), (0, 1200))"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {
-    "deletable": true,
-    "editable": true
-   },
-   "outputs": [],
-   "source": [
-    "i = max2\n",
-    "l = 2*sym\n",
-    "plt.plot(np.abs(sig_rec[i:i + l]))"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {
-    "deletable": true,
-    "editable": true
-   },
-   "outputs": [],
-   "source": [
-    "sig_rec.tofile(\"./orig_frame_8.iq\")"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {
-    "deletable": true,
-    "editable": true
-   },
-   "outputs": [],
-   "source": []
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {
-    "deletable": true,
-    "editable": true
-   },
-   "outputs": [],
-   "source": [
-    "path = \"../test_dat/out.iq\"\n",
-    "frame = 96*8192\n",
-    "offset = 8 * frame + 8192 + 10500 + 72\n",
-    "length  = 1 * frame"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {
-    "deletable": true,
-    "editable": true
-   },
-   "outputs": [],
-   "source": [
-    "sig_orig = np.fromfile(path, dtype=np.complex64, count=offset + length)[offset:]"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {
-    "deletable": true,
-    "editable": true
-   },
-   "outputs": [],
-   "source": [
-    "plt.plot(sig_orig)"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {
-    "deletable": true,
-    "editable": true
-   },
-   "outputs": [],
-   "source": [
-    "sig_orig.tofile(\"./orig_frame_8.iq\")"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {
-    "deletable": true,
-    "editable": true
-   },
-   "outputs": [],
-   "source": []
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {
-    "deletable": true,
-    "editable": true
-   },
-   "outputs": [],
-   "source": [
-    "sig_rec = sig_rec / sig_rec.std()\n",
-    "sig_orig = sig_orig / sig_orig.std()"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {
-    "deletable": true,
-    "editable": true
-   },
-   "outputs": [],
-   "source": [
-    "i = 10000\n",
-    "l = 400\n",
-    "plt.plot(sig_orig[i:i+l], label=\"in\")\n",
-    "plt.plot(sig_rec[i:i+l], label=\"out\")\n",
-    "plt.legend()"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {
-    "deletable": true,
-    "editable": true
-   },
-   "outputs": [],
-   "source": [
-    "from sklearn.linear_model import Ridge\n",
-    "from sklearn.preprocessing import PolynomialFeatures\n",
-    "from sklearn.pipeline import make_pipeline"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {
-    "deletable": true,
-    "editable": true
-   },
-   "source": [
-    "# Approximate Amplifier Function"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {
-    "deletable": true,
-    "editable": true
-   },
-   "outputs": [],
-   "source": [
-    "plt.plot(np.abs(sig_orig), np.abs(sig_rec), 'o')\n",
-    "plt.xlabel(\"Amplitude In\")\n",
-    "plt.ylabel(\"Amplitude Out\")\n",
-    "\n",
-    "deg = 8\n",
-    "t = np.linspace(0,5)\n",
-    "model = make_pipeline(PolynomialFeatures(deg), Ridge())\n",
-    "model.fit(np.abs(sig_orig.reshape(-1,1)), np.abs(sig_rec.reshape(-1,1)))\n",
-    "y = model.predict(t.reshape(-1,1))\n",
-    "plt.plot(t, np.abs(y))"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {
-    "deletable": true,
-    "editable": true
-   },
-   "source": [
-    "# Approximate Inverse Function"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {
-    "deletable": true,
-    "editable": true
-   },
-   "outputs": [],
-   "source": [
-    "plt.plot(np.abs(sig_orig), np.abs(sig_rec), 'o')\n",
-    "plt.xlabel(\"Amplitude In\")\n",
-    "plt.ylabel(\"Amplitude Out\")\n",
-    "\n",
-    "deg = 8\n",
-    "t = np.linspace(0,3)\n",
-    "model = make_pipeline(PolynomialFeatures(deg), Ridge(alpha=100))\n",
-    "model.fit(np.abs(sig_rec.reshape(-1,1)),np.abs(sig_orig.reshape(-1,1)) )\n",
-    "y = model.predict(t.reshape(-1,1))\n",
-    "plt.plot(t, np.abs(y))"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {
-    "deletable": true,
-    "editable": true
-   },
-   "outputs": [],
-   "source": []
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {
-    "deletable": true,
-    "editable": true
-   },
-   "source": [
-    "# Phase Error"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {
-    "deletable": true,
-    "editable": true
-   },
-   "outputs": [],
-   "source": [
-    "l = 100000\n",
-    "plt.plot(np.abs(sig_orig[i:i+l]), (np.angle(sig_rec[i:i+l], deg = 1) - np.angle(sig_orig[i:i+l], deg = 1)) % 180, 'o')\n",
-    "plt.xlabel(\"Amplitude In\")\n",
-    "plt.ylabel(\"Phase difference\")"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {
-    "deletable": true,
-    "editable": true
-   },
-   "outputs": [],
-   "source": []
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {
-    "deletable": true,
-    "editable": true
-   },
-   "outputs": [],
-   "source": []
-  }
- ],
- "metadata": {
-  "kernelspec": {
-   "display_name": "Python 2",
-   "language": "python",
-   "name": "python2"
-  },
-  "language_info": {
-   "codemirror_mode": {
-    "name": "ipython",
-    "version": 2
-   },
-   "file_extension": ".py",
-   "mimetype": "text/x-python",
-   "name": "python",
-   "nbconvert_exporter": "python",
-   "pygments_lexer": "ipython2"
-  }
- },
- "nbformat": 4,
- "nbformat_minor": 2
-}
diff --git a/src/dab_util.py b/src/dab_util.py
index 2b23812..843f8a5 100644
--- a/src/dab_util.py
+++ b/src/dab_util.py
@@ -3,6 +3,7 @@ import scipy
 import matplotlib.pyplot as plt
 import fftconvolve
 import src.dabconst as dabconst
+from scipy import signal
 
 c = {}
 c["bw"]=1536000
@@ -47,27 +48,45 @@ def crop_signal(signal, n_window = 1000, n_zeros = 120000, debug = False):
     signal = signal[max(0,idx_start - n_zeros): min(idx_end + n_zeros, signal.shape[0] -1)]
     return signal
 
-#def fftlag(signal_original, signal_rec):
+#def fftlag(sig_orig, sig_rec):
 #    """
 #    Efficient way to find lag between two signals
 #    Args:
-#        signal_original: The signal that has been sent
-#        signal_rec: The signal that has been recored
+#        sig_orig: The signal that has been sent
+#        sig_rec: The signal that has been recored
 #    """
-#    c = np.flipud(scipy.signal.fftconvolve(signal_original,np.flipud(signal_rec)))
+#    c = np.flipud(scipy.signal.fftconvolve(sig_orig,np.flipud(sig_rec)))
 #    #plt.plot(c)
-#    return np.argmax(c) - signal_original.shape[0] + 1
+#    return np.argmax(c) - sig_orig.shape[0] + 1
 
-def fftlag(signal_original, signal_rec, n_upsampling = 1):
+def lag(sig_orig, sig_rec):
+    """
+    Find lag between two signals
+    Args:
+        sig_orig: The signal that has been sent
+        sig_rec: The signal that has been recored
+    """
+    off = sig_rec.shape[0]
+    c = signal.correlate(sig_orig, sig_rec)
+    return np.argmax(c) - off + 1
+
+def lag_upsampling(sig_orig, sig_rec, n_up):
+    sig_orig_up = signal.resample(sig_orig, sig_orig.shape[0] * n_up)
+    sig_rec_up  = signal.resample(sig_rec, sig_rec.shape[0] * n_up)
+    l = lag(sig_orig_up, sig_rec_up)
+    l_orig = float(l) / n_up
+    return l_orig
+
+def fftlag(sig_orig, sig_rec, n_upsampling = 1):
     """
     Efficient way to find lag between two signals
     Args:
-        signal_original: The signal that has been sent
-        signal_rec: The signal that has been recored
+        sig_orig: The signal that has been sent
+        sig_rec: The signal that has been recored
     """
-    c = np.flipud(fftconvolve.fftconvolve(signal_original,np.flipud(signal_rec), n_upsampling))
+    c = np.flipud(fftconvolve.fftconvolve(sig_orig,np.flipud(sig_rec), n_upsampling))
     #plt.plot(c)
-    return (np.argmax(c) - signal_original.shape[0] + 1)
+    return (np.argmax(c) - sig_orig.shape[0] + 1)
 
 def get_amp_ratio(ampl_1, ampl_2, a_out_abs, a_in_abs):
     idxs = (a_in_abs > ampl_1) & (a_in_abs < ampl_2)
diff --git a/src/dab_util_test.py b/src/dab_util_test.py
new file mode 100644
index 0000000..83813a9
--- /dev/null
+++ b/src/dab_util_test.py
@@ -0,0 +1,36 @@
+from scipy import signal
+import numpy as np
+import src.gen_source as gs
+reload(gs)
+import src.dab_util as du
+reload(du)
+
+def gen_test_signals(oversampling=4, sample_offset_float=0):
+    off = int(sample_offset_float)
+    phi_samples = sample_offset_float - off
+    phi = phi_samples*360/oversampling
+
+    s1 = np.zeros((1024))
+    s1[256:768] = gs.gen_sin(512, oversampling, 0)
+    s2 = np.zeros((1024))
+    s2[256+off:768+off] = gs.gen_sin(512, oversampling, phi)
+
+    return s1, s2
+
+def test_phase_offset(lag_function, tol):
+    def r():
+        return np.random.rand(1)*100-50
+    res = []
+    for i in range(100):
+        off = r()
+        s1, s2 = gen_test_signals(
+            oversampling=4, sample_offset_float=off)
+
+        off_meas = lag_function(s2, s1)
+        res.append(np.abs(off-off_meas)<tol)
+    return np.mean(res)
+
+for n_up in [1, 2, 3, 4, 7, 8, 16]:
+   correct_ratio = test_phase_offset(lambda x,y: du.lag_upsampling(x,y,n_up), tol=1./n_up)
+   print("%.1f%% of the tested offsets were measured within tolerance %.4f for n_up = %d" % (correct_ratio * 100, 1./n_up, n_up))
+
diff --git a/src/gen_source.py b/src/gen_source.py
index 7620bc3..c8509ed 100644
--- a/src/gen_source.py
+++ b/src/gen_source.py
@@ -97,3 +97,9 @@ def gen_file_i(frequency_0, frequency_1, x1 = 0, x2 = 0, x3 = 0, x4 = 0, samp_ra
     assert(np.isclose(a_load, two_tone).all()), "Inconsistent stored file"
 
     return path
+
+def gen_sin(samples, oversampling, phi):
+    t = np.arange(samples, dtype=np.float)
+    sig = np.sin(((2*np.pi)/oversampling) * t - np.pi*phi/180.)
+    return sig
+