Remove polynom term without coefficient

4 files changed, 162 insertions, 28 deletions

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)