Add iterative dpd

2 files changed, 126 insertions, 38 deletions

diff --git a/dpd/main.py b/dpd/main.py
index 72026dc..b0b7a87 100755
--- a/dpd/main.py
+++ b/dpd/main.py
@@ -59,7 +59,7 @@ adapt = Adapt.Adapt(port_rc, coef_path)
 coefs = adapt.get_coefs()
 #model = Model.Model(coefs)
 model = Model.Model([2.2, 0, 0, 0, 0])
-adapt.set_txgain(82)
+adapt.set_txgain(84)
 
 tx_gain   = adapt.get_txgain()
 rx_gain   = adapt.get_rxgain()
@@ -70,7 +70,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/src/Model.py b/dpd/src/Model.py
index 7664dd3..aef112c 100644
--- a/dpd/src/Model.py
+++ b/dpd/src/Model.py
@@ -5,6 +5,7 @@ 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')
@@ -17,32 +18,71 @@ 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):
         rx_abs = np.abs(rxframe_aligned)
-        A = np.vstack([rx_abs,
+        rx_A = np.vstack([rx_abs,
                        rx_abs**3,
                        rx_abs**5,
                        rx_abs**7,
                        rx_abs**9,
                        ]).T
-        y = np.abs(txframe_aligned)
-
-        clf = Ridge(alpha=10)
-        clf.fit(A, y)
-        sol = clf.coef_
-
-        rx_range = np.linspace(0,1,50)
-        A_range = np.vstack([
-                       rx_range,
-                       rx_range**3,
-                       rx_range**5,
-                       rx_range**7,
-                       rx_range**9,
-                       ]).T
-        y_est = np.sum(A_range * sol, axis=1)
-
-        logging.debug("New coefficents {}".format(sol))
+        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
+        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" %
@@ -60,19 +100,33 @@ class Model:
             dt = datetime.datetime.now().isoformat()
             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")
@@ -82,13 +136,7 @@ class Model:
             ax.scatter(
                 np.abs(txframe_aligned[:1024]),
                 np.abs(rxframe_aligned[:1024]),
-                s = 0.1
-            )
-            ax.plot(
-                y_est,
-                rx_range,
-                linewidth=0.25
-            )
+                s = 0.1)
             ax.set_title("Amplifier Characteristic")
             ax.set_xlabel("TX Amplitude")
             ax.set_ylabel("RX Amplitude")
@@ -109,15 +157,55 @@ 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))
-
-        sol = sol * 1.7/sol[0]
-        return sol
+        self.coefs = new_coefs
+        self.coefs_history.append(self.coefs)
+        return self.coefs
 
 # The MIT License (MIT)
 #