Add control for dpd phase coefficient

3 files changed, 158 insertions, 96 deletions

diff --git a/dpd/main.py b/dpd/main.py
index 4d2b93d..7e1041b 100755
--- a/dpd/main.py
+++ b/dpd/main.py
@@ -38,9 +38,12 @@ parser.add_argument('--rc-port', default='9400',
 parser.add_argument('--samplerate', default='8192000',
         help='Sample rate',
         required=False)
-parser.add_argument('--coefs', default='poly.coef',
-        help='File with DPD coefficients, which will be read by ODR-DabMod',
+parser.add_argument('--coefs_am', default='poly_am.coef',
+        help='File with DPD Amplitude coefficients, which will be read by ODR-DabMod',
         required=False)
+parser.add_argument('--coefs_pm', default='poly_am.coef',
+                    help='File with DPD Phase coefficients, which will be read by ODR-DabMod',
+                    required=False)
 parser.add_argument('--samps', default='10240',
         help='Number of samples to request from ODR-DabMod',
         required=False)
@@ -49,33 +52,37 @@ cli_args = parser.parse_args()
 
 port = int(cli_args.port)
 port_rc = int(cli_args.rc_port)
-coef_path = cli_args.coefs
+coef_am_path = cli_args.coefs_am
+coef_pm_path = cli_args.coefs_pm
 num_req = int(cli_args.samps)
 samplerate = int(cli_args.samplerate)
 
 meas = Measure.Measure(samplerate, port, num_req)
 
-adapt = Adapt.Adapt(port_rc, coef_path)
-coefs = adapt.get_coefs()
+adapt = Adapt.Adapt(port_rc, coef_am_path, coef_pm_path)
+coefs_am = adapt.get_coefs_am()
+coefs_pm = adapt.get_coefs_pm()
 #model = Model.Model(coefs)
-model = Model.Model([2.2, 0, 0, 0, 0])
-adapt.set_txgain(79)
-adapt.set_rxgain(15+20)
+model = Model.Model([2.2, 0, 0, 0, 0], [0, 0, 0, 0, 0])
+adapt.set_txgain(70)
+adapt.set_rxgain(30)
 
 tx_gain   = adapt.get_txgain()
 rx_gain   = adapt.get_rxgain()
-dpd_coefs = adapt.get_coefs()
+dpd_coefs_am = adapt.get_coefs_am()
+dpd_coefs_pm = adapt.get_coefs_pm()
 logging.info(
-    "TX gain {}, RX gain {}, dpd_coefs {}".format(
-        tx_gain, rx_gain, dpd_coefs
+    "TX gain {}, RX gain {}, dpd_coefs_am {}, dpd_coefs_pm {}".format(
+        tx_gain, rx_gain, dpd_coefs_am, dpd_coefs_pm
     )
 )
 
 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)
-    adapt.set_coefs(coefs)
+    coefs_am, coefs_pm = model.get_next_coefs(txframe_aligned, rxframe_aligned)
+    adapt.set_coefs_am(coefs_am)
+    adapt.set_coefs_pm(coefs_pm)
 
 # The MIT License (MIT)
 #
diff --git a/dpd/src/Adapt.py b/dpd/src/Adapt.py
index 4117856..cc4c8c7 100644
--- a/dpd/src/Adapt.py
+++ b/dpd/src/Adapt.py
@@ -23,10 +23,11 @@ class Adapt:
         ZMQ remote control.
     """
 
-    def __init__(self, port, coef_path):
+    def __init__(self, port, coef_am_path, coef_pm_path):
         logging.info("Instantiate Adapt object")
         self.port = port
-        self.coef_path = coef_path
+        self.coef_am_path = coef_am_path
+        self.coef_pm_path = coef_pm_path
         self.host = "localhost"
         self._context = zmq.Context()
 
@@ -110,10 +111,10 @@ class Adapt:
         # TODO handle failure
         return self.send_receive("get uhd rxgain")
 
-    def _read_coef_file(self):
+    def _read_coef_file(self, path):
         """Load the coefficients from the file in the format given in the README"""
         coefs_out = []
-        f = open(self.coef_path, 'r')
+        f = open(path, 'r')
         lines = f.readlines()
         n_coefs = int(lines[0])
         coefs = [float(l) for l in lines[1:]]
@@ -124,24 +125,31 @@ class Adapt:
             else:
                 raise ValueError(
                     "Incorrect coef file format: too many coefficients in {}, should be {}, coefs are {}"
-                        .format(self.coef_path, n_coefs, coefs))
+                        .format(path, n_coefs, coefs))
             i += 1
         f.close()
         return coefs_out
 
-    def get_coefs(self):
-        return self._read_coef_file()
+    def get_coefs_am(self):
+        return self._read_coef_file(self.coef_am_path)
 
-    def _write_coef_file(self, coefs):
-        f = open(self.coef_path, 'w')
+    def get_coefs_pm(self):
+        return self._read_coef_file(self.coef_pm_path)
+
+    def _write_coef_file(self, coefs, path):
+        f = open(path, 'w')
         f.write("{}\n".format(len(coefs)))
         for coef in coefs:
             f.write("{}\n".format(coef))
         f.close()
 
-    def set_coefs(self, coefs):
-        self._write_coef_file(coefs)
-        self.send_receive("set memlesspoly coeffile {}".format(self.coef_path))
+    def set_coefs_am(self, coefs):
+        self._write_coef_file(coefs, self.coef_am_path)
+        self.send_receive("set memlesspoly coeffile_am {}".format(self.coef_am_path))
+
+    def set_coefs_pm(self, coefs):
+        self._write_coef_file(coefs, self.coef_pm_path)
+        self.send_receive("set memlesspoly coeffile_pm {}".format(self.coef_pm_path))
 
 # The MIT License (MIT)
 #
diff --git a/dpd/src/Model.py b/dpd/src/Model.py
index 8df3925..b962140 100644
--- a/dpd/src/Model.py
+++ b/dpd/src/Model.py
@@ -16,42 +16,86 @@ class Model:
     """Calculates new coefficients using the measurement and the old
     coefficients"""
 
-    def __init__(self, coefs):
-        self.coefs = coefs
-        self.coefs_history = [coefs,]
-        self.mses = [0,]
-        self.errs = [0,]
+    def __init__(self, coefs_am, coefs_pm):
+        self.coefs_am = coefs_am
+        self.coefs_history = [coefs_am, ]
+        self.mses = [0, ]
+        self.errs = [0, ]
+
+        self.coefs_pm = coefs_pm
+        self.coefs_pm_history = [coefs_pm, ]
+        self.errs_phase = [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")
+
+        # Calculate new coefficients for AM/AM correction
         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_abs ** 3,
+                          rx_abs ** 5,
+                          rx_abs ** 7,
+                          rx_abs ** 9,
+                          ]).T
+        rx_dpd = np.sum(rx_A * self.coefs_am, 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)))
+        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])
+        new_coefs = self.coefs_am - 0.1 * a_delta
+        new_coefs = new_coefs * (self.coefs_am[0] / new_coefs[0])
+        logging.debug("a_delta {}".format(a_delta))
+        logging.debug("new coefs_am {}".format(new_coefs))
+
+        # Calculate new coefficients for AM/PM correction
+        phase_diff_rad = ((
+                              (np.angle(txframe_aligned) -
+                               np.angle(rxframe_aligned) +
+                               np.pi) % (2 * np.pi)) -
+                          np.pi
+                          )
+
+        tx_abs = np.abs(txframe_aligned)
+        tx_abs_A = np.vstack([tx_abs,
+                             tx_abs ** 2,
+                             tx_abs ** 3,
+                             tx_abs ** 4,
+                             tx_abs ** 5,
+                             ]).T
+        phase_dpd = np.sum(tx_abs_A * self.coefs_pm, axis=1)
+
+        err_phase = phase_dpd - phase_diff_rad
+        self.errs_phase.append(np.mean(np.abs(err_phase ** 2)))
+        a_delta = np.linalg.lstsq(tx_abs_A, err_phase)[0]
+        new_coefs_pm = self.coefs_pm - 0.1 * a_delta
         logging.debug("a_delta {}".format(a_delta))
-        logging.debug("new coefs {}".format(new_coefs))
+        logging.debug("new new_coefs_pm {}".format(new_coefs_pm))
+
+        def dpd_phase(tx):
+            tx_abs = np.abs(tx)
+            tx_A_complex = np.vstack([tx,
+                                      tx * tx_abs ** 1,
+                                      tx * tx_abs ** 2,
+                                      tx * tx_abs ** 3,
+                                      tx * tx_abs ** 4,
+                                      ]).T
+            tx_dpd = np.sum(tx_A_complex * self.coefs_pm, axis=1)
+            return tx_dpd
+
+        tx_range = np.linspace(0, 2)
+        phase_range_dpd = dpd_phase(tx_range)
 
         tx_abs = np.abs(rxframe_aligned)
         tx_A = np.vstack([tx_abs,
-                          tx_abs**3,
-                          tx_abs**5,
-                          tx_abs**7,
-                          tx_abs**9,
+                          tx_abs ** 3,
+                          tx_abs ** 5,
+                          tx_abs ** 7,
+                          tx_abs ** 9,
                           ]).T
         tx_dpd = np.sum(tx_A * new_coefs, axis=1)
 
@@ -59,34 +103,34 @@ class Model:
             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)
+                                  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_am, 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))
+        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)
+                                      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_am, 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))
-        )
+
+        rx_range = np.linspace(0, 1, num=100)
+        rx_range_dpd = dpd(rx_range)
+        rx_range = rx_range[(rx_range_dpd > 0) & (rx_range_dpd < 2)]
+        rx_range_dpd = rx_range_dpd[(rx_range_dpd > 0) & (rx_range_dpd < 2)]
 
         if logging.getLogger().getEffectiveLevel() == logging.DEBUG:
             logging.debug("txframe: min %f, max %f, median %f" %
@@ -104,15 +148,12 @@ class Model:
             dt = datetime.datetime.now().isoformat()
             fig_path = logging_path + "/" + dt + "_Model.pdf"
 
-            fig, axs = plt.subplots(7, figsize=(6,3*6))
+            fig, axs = plt.subplots(8, figsize=(6, 4 * 6))
 
             ax = axs[0]
             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")
@@ -125,9 +166,6 @@ class Model:
             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")
@@ -140,44 +178,40 @@ class Model:
             ax.scatter(
                 np.abs(txframe_aligned[:1024]),
                 np.abs(rxframe_aligned[:1024]),
-                s = 0.1)
+                s=0.1)
+            ax.plot(rx_range_dpd / self.coefs_am[0], rx_range, linewidth=0.25)
             ax.set_title("Amplifier Characteristic")
             ax.set_xlabel("TX Amplitude")
             ax.set_ylabel("RX Amplitude")
 
             ax = axs[3]
-            angle_diff_rad = ((
-                             (np.angle(txframe_aligned[:1024]) -
-                              np.angle(rxframe_aligned[:1024]) +
-                              np.pi) % (2 * np.pi)) -
-                          np.pi
-                          )
             ax.scatter(
                 np.abs(txframe_aligned[:1024]),
-                angle_diff_rad * 180 / np.pi,
-                s = 0.1
+                phase_diff_rad[:1024] * 180 / np.pi,
+                s=0.1
             )
+            ax.plot(tx_range, phase_range_dpd * 180 / np.pi, linewidth=0.25)
             ax.set_title("Amplifier Characteristic")
             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)
+                    label="TX Frame",
+                    linestyle=":",
+                    linewidth=0.5)
             ax.plot(np.abs(rxframe_aligned[:128]),
-                     label="RX Frame",
-                     linestyle="--",
-                     linewidth=0.5)
+                    label="RX Frame",
+                    linestyle="--",
+                    linewidth=0.5)
             ax.plot(np.abs(rx_dpd[:128]),
-                     label="RX DPD Frame",
-                     linestyle="-.",
-                     linewidth=0.5)
+                    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)
+                    label="TX DPD Frame Norm",
+                    linestyle="-.",
+                    linewidth=0.5)
             ax.legend(loc=4)
             ax.set_title("RX DPD")
             ax.set_xlabel("Samples")
@@ -187,14 +221,25 @@ class Model:
             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)
+                        label="Coef {}".format(idx),
+                        linewidth=0.5)
             ax.legend(loc=4)
-            ax.set_title("Coefficient History")
+            ax.set_title("AM/AM Coefficient History")
             ax.set_xlabel("Iterations")
             ax.set_ylabel("Coefficient Value")
 
             ax = axs[6]
+            coefs_history = np.array(self.coefs_pm_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("AM/PM Coefficient History")
+            ax.set_xlabel("Iterations")
+            ax.set_ylabel("Coefficient Value")
+
+            ax = axs[7]
             coefs_history = np.array(self.coefs_history)
             ax.plot(self.mses, label="MSE")
             ax.plot(self.errs, label="ERR")
@@ -207,9 +252,11 @@ class Model:
             fig.savefig(fig_path)
             fig.clf()
 
-        self.coefs = new_coefs
-        self.coefs_history.append(self.coefs)
-        return self.coefs
+        self.coefs_am = new_coefs
+        self.coefs_history.append(self.coefs_am)
+        self.coefs_pm = new_coefs_pm
+        self.coefs_pm_history.append(self.coefs_pm)
+        return self.coefs_am, self.coefs_pm
 
 # The MIT License (MIT)
 #