Remove polynom term without coefficient

8 files changed, 185 insertions, 45 deletions

diff --git a/dpd/main.py b/dpd/main.py
index 560d142..4d2b93d 100755
--- a/dpd/main.py
+++ b/dpd/main.py
@@ -10,10 +10,16 @@
 This engine calculates and updates the parameter of the digital
 predistortion module of ODR-DabMod."""
 
+import datetime
+import os
+
 import logging
+dt = datetime.datetime.now().isoformat()
+logging_path = "/tmp/dpd_{}".format(dt).replace(".","_").replace(":","-")
+os.makedirs(logging_path)
 logging.basicConfig(format='%(asctime)s - %(module)s - %(levelname)s - %(message)s',
                     datefmt='%Y-%m-%d %H:%M:%S',
-                    filename='/tmp/dpd.log',
+                    filename='{}/dpd.log'.format(logging_path),
                     filemode='w',
                     level=logging.DEBUG)
 
@@ -32,7 +38,7 @@ parser.add_argument('--rc-port', default='9400',
 parser.add_argument('--samplerate', default='8192000',
         help='Sample rate',
         required=False)
-parser.add_argument('--coefs', default='dpdpoly.coef',
+parser.add_argument('--coefs', default='poly.coef',
         help='File with DPD coefficients, which will be read by ODR-DabMod',
         required=False)
 parser.add_argument('--samps', default='10240',
@@ -52,8 +58,9 @@ meas = Measure.Measure(samplerate, port, num_req)
 adapt = Adapt.Adapt(port_rc, coef_path)
 coefs = adapt.get_coefs()
 #model = Model.Model(coefs)
-model = Model.Model([0.8, 0, 0, 0, 0])
-adapt.set_txgain(60)
+model = Model.Model([2.2, 0, 0, 0, 0])
+adapt.set_txgain(79)
+adapt.set_rxgain(15+20)
 
 tx_gain   = adapt.get_txgain()
 rx_gain   = adapt.get_rxgain()
@@ -64,7 +71,7 @@ logging.info(
     )
 )
 
-for i in range(10):
+for i in range(500):
     txframe_aligned, tx_ts, rxframe_aligned, rx_ts = meas.get_samples()
     logging.debug("tx_ts {}, rx_ts {}".format(tx_ts, rx_ts))
     coefs = model.get_next_coefs(txframe_aligned, rxframe_aligned)
diff --git a/dpd/poly_am.coef b/dpd/poly_am.coef
index ce11551..ff09e5a 100644
--- a/dpd/poly_am.coef
+++ b/dpd/poly_am.coef
@@ -1,6 +1,6 @@
 5
 2.2
-2.34402815751
--25.0487363081
-93.1884045291
--85.3795993469
+1.16110192544
+-13.7895532562
+55.107515965
+-53.8583673922
diff --git a/dpd/poly_pm.coef b/dpd/poly_pm.coef
index 2a5bf1c..ae91b12 100644
--- a/dpd/poly_pm.coef
+++ b/dpd/poly_pm.coef
@@ -1,6 +1,6 @@
 5
-0.126147521763
--0.507247679586
-0.645526051204
--0.313225683513
-0.0522039980935
+0.141477421551
+-0.654864288614
+1.00568534673
+-0.588530075442
+0.0935391293974
diff --git a/dpd/src/Dab_Util.py b/dpd/src/Dab_Util.py
index 1f88ae4..175b744 100644
--- a/dpd/src/Dab_Util.py
+++ b/dpd/src/Dab_Util.py
@@ -1,15 +1,18 @@
 # -*- coding: utf-8 -*-
 
+import datetime
+import os
+import logging
+logging_path = os.path.dirname(logging.getLoggerClass().root.handlers[0].baseFilename)
+
 import numpy as np
 import scipy
 import matplotlib
 matplotlib.use('agg')
 import matplotlib.pyplot as plt
-import datetime
 import src.subsample_align as sa
 import src.phase_align as pa
 from scipy import signal
-import logging
 
 class Dab_Util:
     """Collection of methods that can be applied to an array
@@ -35,7 +38,7 @@ class Dab_Util:
 
         if logging.getLogger().getEffectiveLevel() == logging.DEBUG:
             dt = datetime.datetime.now().isoformat()
-            corr_path = ("/tmp/" + dt + "_tx_rx_corr.pdf")
+            corr_path = (logging_path + "/" + dt + "_tx_rx_corr.pdf")
             plt.plot(c, label="corr")
             plt.legend()
             plt.savefig(corr_path)
@@ -89,7 +92,7 @@ class Dab_Util:
 
         if logging.getLogger().getEffectiveLevel() == logging.DEBUG:
             dt = datetime.datetime.now().isoformat()
-            fig_path = "/tmp/" + dt + "_sync_raw.pdf"
+            fig_path = logging_path + "/" + dt + "_sync_raw.pdf"
 
             fig, axs = plt.subplots(2)
             axs[0].plot(np.abs(sig1[:128]), label="TX Frame")
@@ -127,7 +130,7 @@ class Dab_Util:
 
         if logging.getLogger().getEffectiveLevel() == logging.DEBUG:
             dt = datetime.datetime.now().isoformat()
-            fig_path = "/tmp/" + dt + "_sync_sample_aligned.pdf"
+            fig_path = logging_path + "/" + dt + "_sync_sample_aligned.pdf"
 
             fig, axs = plt.subplots(2)
             axs[0].plot(np.abs(sig1[:128]), label="TX Frame")
@@ -152,7 +155,7 @@ class Dab_Util:
 
         if logging.getLogger().getEffectiveLevel() == logging.DEBUG:
             dt = datetime.datetime.now().isoformat()
-            fig_path = "/tmp/" + dt + "_sync_subsample_aligned.pdf"
+            fig_path = logging_path + "/" + dt + "_sync_subsample_aligned.pdf"
 
             fig, axs = plt.subplots(2)
             axs[0].plot(np.abs(sig1[:128]), label="TX Frame")
@@ -176,7 +179,7 @@ class Dab_Util:
 
         if logging.getLogger().getEffectiveLevel() == logging.DEBUG:
             dt = datetime.datetime.now().isoformat()
-            fig_path = "/tmp/" + dt + "_sync_phase_aligned.pdf"
+            fig_path = logging_path + "/" + dt + "_sync_phase_aligned.pdf"
 
             fig, axs = plt.subplots(2)
             axs[0].plot(np.abs(sig1[:128]), label="TX Frame")
diff --git a/dpd/src/Model.py b/dpd/src/Model.py
index f66ba8f..8df3925 100644
--- a/dpd/src/Model.py
+++ b/dpd/src/Model.py
@@ -1,11 +1,16 @@
 # -*- coding: utf-8 -*-
 
-import numpy as np
 import datetime
+import os
 import logging
+logging_path = os.path.dirname(logging.getLoggerClass().root.handlers[0].baseFilename)
+
+from pynverse import inversefunc
+import numpy as np
 import matplotlib
 matplotlib.use('agg')
 import matplotlib.pyplot as plt
+from sklearn.linear_model import Ridge
 
 class Model:
     """Calculates new coefficients using the measurement and the old
@@ -13,8 +18,76 @@ class Model:
 
     def __init__(self, coefs):
         self.coefs = coefs
+        self.coefs_history = [coefs,]
+        self.mses = [0,]
+        self.errs = [0,]
 
     def get_next_coefs(self, txframe_aligned, rxframe_aligned):
+        dt = datetime.datetime.now().isoformat()
+        txframe_aligned.tofile(logging_path + "/txframe_" + dt + ".iq")
+        rxframe_aligned.tofile(logging_path + "/rxframe_" + dt + ".iq")
+        rx_abs = np.abs(rxframe_aligned)
+        rx_A = np.vstack([rx_abs,
+                       rx_abs**3,
+                       rx_abs**5,
+                       rx_abs**7,
+                       rx_abs**9,
+                       ]).T
+        rx_dpd = np.sum(rx_A * self.coefs, axis=1)
+        rx_dpd = rx_dpd * (
+            np.median(np.abs(txframe_aligned)) / np.median(np.abs(rx_dpd)))
+
+        err = rx_dpd - np.abs(txframe_aligned)
+        self.errs.append(np.mean(np.abs(err**2)))
+
+        a_delta = np.linalg.lstsq(rx_A, err)[0]
+        new_coefs = self.coefs - 0.1 * a_delta
+        new_coefs = new_coefs * (self.coefs[0] / new_coefs[0])
+        logging.debug("a_delta {}".format(a_delta))
+        logging.debug("new coefs {}".format(new_coefs))
+
+        tx_abs = np.abs(rxframe_aligned)
+        tx_A = np.vstack([tx_abs,
+                          tx_abs**3,
+                          tx_abs**5,
+                          tx_abs**7,
+                          tx_abs**9,
+                          ]).T
+        tx_dpd = np.sum(tx_A * new_coefs, axis=1)
+
+        tx_dpd_norm = tx_dpd * (
+            np.median(np.abs(txframe_aligned)) / np.median(np.abs(tx_dpd)))
+
+        rx_A_complex = np.vstack([rxframe_aligned,
+                          rxframe_aligned * rx_abs**2,
+                          rxframe_aligned * rx_abs**4,
+                          rxframe_aligned * rx_abs**6,
+                          rxframe_aligned * rx_abs**8,
+                          ]).T
+        rx_post_distored = np.sum(rx_A_complex * self.coefs, axis=1)
+        rx_post_distored = rx_post_distored * (
+            np.median(np.abs(txframe_aligned)) /
+            np.median(np.abs(rx_post_distored)))
+        mse = np.mean(np.abs((txframe_aligned - rx_post_distored)**2))
+        logging.debug("MSE: {}".format(mse))
+        self.mses.append(mse)
+
+        def dpd(tx):
+            tx_abs = np.abs(tx)
+            tx_A_complex = np.vstack([tx,
+                                  tx * tx_abs**2,
+                                  tx * tx_abs**4,
+                                  tx * tx_abs**6,
+                                  tx * tx_abs**8,
+                                  ]).T
+            tx_dpd = np.sum(tx_A_complex * self.coefs, axis=1)
+            return tx_dpd
+        tx_inverse_dpd = inversefunc(dpd, y_values=txframe_aligned[:128])
+        tx_inverse_dpd = tx_inverse_dpd * (
+            np.median(np.abs(txframe_aligned)) /
+            np.median(np.abs(tx_inverse_dpd))
+        )
+
         if logging.getLogger().getEffectiveLevel() == logging.DEBUG:
             logging.debug("txframe: min %f, max %f, median %f" %
                           (np.min(np.abs(txframe_aligned)),
@@ -29,21 +102,35 @@ class Model:
                            ))
 
             dt = datetime.datetime.now().isoformat()
-            fig_path = "/tmp/" + dt + "_Model.pdf"
+            fig_path = logging_path + "/" + dt + "_Model.pdf"
 
-            fig, axs = plt.subplots(4, figsize=(6,2*6))
+            fig, axs = plt.subplots(7, figsize=(6,3*6))
 
             ax = axs[0]
-            ax.plot(np.abs(txframe_aligned[:128]), label="TX Frame")
-            ax.plot(np.abs(rxframe_aligned[:128]), label="RX Frame")
-            ax.set_title("Synchronized Signals")
+            ax.plot(np.abs(txframe_aligned[:128]),
+                    label="TX sent",
+                    linestyle=":")
+            ax.plot(np.abs(tx_inverse_dpd[:128]),
+                    label="TX inverse dpd",
+                    color="green")
+            ax.plot(np.abs(rxframe_aligned[:128]),
+                    label="RX received",
+                    color="red")
+            ax.set_title("Synchronized Signals of Iteration {}".format(len(self.coefs_history)))
             ax.set_xlabel("Samples")
             ax.set_ylabel("Amplitude")
             ax.legend(loc=4)
 
             ax = axs[1]
-            ax.plot(np.real(txframe_aligned[:128]), label="TX Frame")
-            ax.plot(np.real(rxframe_aligned[:128]), label="RX Frame")
+            ax.plot(np.real(txframe_aligned[:128]),
+                    label="TX sent",
+                    linestyle=":")
+            ax.plot(np.real(tx_inverse_dpd[:128]),
+                    label="TX inverse dpd",
+                    color="green")
+            ax.plot(np.real(rxframe_aligned[:128]),
+                    label="RX received",
+                    color="red")
             ax.set_title("Synchronized Signals")
             ax.set_xlabel("Samples")
             ax.set_ylabel("Real Part")
@@ -53,8 +140,7 @@ class Model:
             ax.scatter(
                 np.abs(txframe_aligned[:1024]),
                 np.abs(rxframe_aligned[:1024]),
-                s = 0.1
-            )
+                s = 0.1)
             ax.set_title("Amplifier Characteristic")
             ax.set_xlabel("TX Amplitude")
             ax.set_ylabel("RX Amplitude")
@@ -75,13 +161,54 @@ class Model:
             ax.set_xlabel("TX Amplitude")
             ax.set_ylabel("Phase Difference [deg]")
 
+            ax = axs[4]
+            ax.plot(np.abs(txframe_aligned[:128]),
+                     label="TX Frame",
+                     linestyle=":",
+                     linewidth=0.5)
+            ax.plot(np.abs(rxframe_aligned[:128]),
+                     label="RX Frame",
+                     linestyle="--",
+                     linewidth=0.5)
+            ax.plot(np.abs(rx_dpd[:128]),
+                     label="RX DPD Frame",
+                     linestyle="-.",
+                     linewidth=0.5)
+            ax.plot(np.abs(tx_dpd_norm[:128]),
+                     label="TX DPD Frame Norm",
+                     linestyle="-.",
+                     linewidth=0.5)
+            ax.legend(loc=4)
+            ax.set_title("RX DPD")
+            ax.set_xlabel("Samples")
+            ax.set_ylabel("Amplitude")
+
+            ax = axs[5]
+            coefs_history = np.array(self.coefs_history)
+            for idx, coef_hist in enumerate(coefs_history.T):
+                ax.plot(coef_hist,
+                     label="Coef {}".format(idx),
+                     linewidth=0.5)
+            ax.legend(loc=4)
+            ax.set_title("Coefficient History")
+            ax.set_xlabel("Iterations")
+            ax.set_ylabel("Coefficient Value")
+
+            ax = axs[6]
+            coefs_history = np.array(self.coefs_history)
+            ax.plot(self.mses, label="MSE")
+            ax.plot(self.errs, label="ERR")
+            ax.legend(loc=4)
+            ax.set_title("MSE History")
+            ax.set_xlabel("Iterations")
+            ax.set_ylabel("MSE")
+
             fig.tight_layout()
             fig.savefig(fig_path)
             fig.clf()
 
-        mse = np.mean(np.abs(np.square(txframe_aligned[:1024] - rxframe_aligned[:1024])))
-        logging.debug("MSE: {}".format(mse))
-
+        self.coefs = new_coefs
+        self.coefs_history.append(self.coefs)
         return self.coefs
 
 # The MIT License (MIT)
diff --git a/dpd/src/phase_align.py b/dpd/src/phase_align.py
index bea0b82..f03184b 100644
--- a/dpd/src/phase_align.py
+++ b/dpd/src/phase_align.py
@@ -1,9 +1,12 @@
+import datetime
+import os
+import logging
+logging_path = os.path.dirname(logging.getLoggerClass().root.handlers[0].baseFilename)
+
 import numpy as np
 from scipy import signal, optimize
 import sys
 import matplotlib.pyplot as plt
-import datetime
-import logging
 
 
 def phase_align(sig, ref_sig):
@@ -19,7 +22,7 @@ def phase_align(sig, ref_sig):
 
     if logging.getLogger().getEffectiveLevel() == logging.DEBUG:
         dt = datetime.datetime.now().isoformat()
-        fig_path = "/tmp/" + dt + "_phase_align.pdf"
+        fig_path = logging_path + "/" + dt + "_phase_align.pdf"
 
         plt.subplot(511)
         plt.hist(angle_diff, bins=60, label="Angle Diff")
diff --git a/dpd/src/subsample_align.py b/dpd/src/subsample_align.py
index 0dc78c1..0a51593 100755
--- a/dpd/src/subsample_align.py
+++ b/dpd/src/subsample_align.py
@@ -1,8 +1,11 @@
+import datetime
+import os
+import logging
+logging_path = os.path.dirname(logging.getLoggerClass().root.handlers[0].baseFilename)
+
 import numpy as np
 from scipy import signal, optimize
-import sys
 import matplotlib.pyplot as plt
-import datetime
 
 def gen_omega(length):
     if (length % 2) == 1:
@@ -59,15 +62,13 @@ def subsample_align(sig, ref_sig):
     if optim_result.success:
         best_tau = optim_result.x
 
-        #print("Found subsample delay = {}".format(best_tau))
-
         if 1:
             corr = np.vectorize(correlate_for_delay)
             ixs = np.linspace(-1, 1, 100)
             taus = corr(ixs)
 
             dt = datetime.datetime.now().isoformat()
-            tau_path = ("/tmp/" + dt + "_tau.pdf")
+            tau_path = (logging_path + "/" + dt + "_tau.pdf")
             plt.plot(ixs, taus)
             plt.title("Subsample correlation, minimum is best: {}".format(best_tau))
             plt.savefig(tau_path)
diff --git a/src/MemlessPoly.cpp b/src/MemlessPoly.cpp
index d665839..e103b73 100644
--- a/src/MemlessPoly.cpp
+++ b/src/MemlessPoly.cpp
@@ -175,7 +175,6 @@ void MemlessPoly::load_coefficients_pm(const std::string &coefFile_pm)
 /* The restrict keyword is C99, g++ and clang++ however support __restrict
  * instead, and this allows the compiler to auto-vectorize the loop.
  */
-
 static void apply_coeff(
         const vector<float> &coefs_am,
         const vector<float> &coefs_pm,
@@ -191,14 +190,14 @@ static void apply_coeff(
               ( coefs_am[1] + in_mag_sq *
                 ( coefs_am[2] + in_mag_sq *
                   ( coefs_am[3] + in_mag_sq *
-                    ( coefs_am[4] + in_mag_sq )))));
+                    coefs_am[4]))));
 
         float phase_correction = -1 *
             ( coefs_pm[0] + in_mag_sq *
               ( coefs_pm[1] + in_mag_sq *
                 ( coefs_pm[2] + in_mag_sq *
                   ( coefs_pm[3] + in_mag_sq *
-                    ( coefs_pm[4] + in_mag_sq )))));
+                    coefs_pm[4]))));
 
         float phase_correction_sq = phase_correction * phase_correction;