add some analysis of generated signal

README.md

@@ -1,9 +1,9 @@
 # Charon VNA
+
 <!-- ![PyPi Downloads](https://img.shields.io/pypi/dm/charon-vna) -->
 <!-- ![Last Commit](https://img.shields.io/gitea/last-commit/brendanhaines/charon-vna?gitea_url=https%3A%2F%2Fgit.brendanhaines.com) -->
 <!-- ![Workflow Status](https://git.brendanhaines.com/brendanhaines/charon-vna/actions/workflows/python_publish.yml/badge.svg) -->
 
-
 Named after [Pluto's moon](https://en.wikipedia.org/wiki/Charon_(moon)), Charon uses the [ADI Pluto SDR](https://www.analog.com/en/resources/evaluation-hardware-and-software/evaluation-boards-kits/adalm-pluto.html) as a vector network analyzer. The basic usage is as a 1 port VNA but this can be extended to arbitrarily many ports with the addition of a couple RF switches.
 
 ## Installation
@@ -22,13 +22,14 @@ You need a few things:
   - Note that you _must_ have two receive ports which means revision C or later of the basic Pluto
 - Directional couplers (1 per port up to 4 ports)
   - I have been using [AAMCS-UDC-0.5G-18G-10dB-Sf](http://www.aa-mcs.com/wp-content/uploads/documents/AAMCS-UDC-0.5G-18G-10dB-Sf.pdf)
-- Charon switch board  - coming soon.
+- [Cerberus RF switch](https://git.brendanhaines.com/brendanhaines/cerberus_sp4t) + [Pluto IO Shield](https://git.brendanhaines.com/brendanhaines/pluto_io_shield)
   - Optional. Without this you'll be limited to S11 and uncalibrated S21 measurements (with required re-cabling)
   - There's nothing special about this particular board, if you want more than 4 ports you can make your own pretty easily. You just need 3 SPxT switches. Note that these switches will see tons of cycles so avoid mechanical switches
 - SMA cables
 - Calibration standard
   - Ideally something with s-parameters measured on a better VNA
   - I have used a basic SMA load and two modified SMA jacks with decent results
+    ![calibration standard](img/calibration_standard.jpg)
 
 ### Pluto Configuration
 
@@ -53,6 +54,7 @@ It will also be accessible over a socket to enable test automation with external
 TBD
 
 ### Power Calibration
+
 I include a default output power lookup table. This is derived from two TX channels of two Pluto SDRs and does not include any of the loss of a coupler or Charon switch board.
 
 Absolute output power is generally not well calibrated for VNAs anyway and has negligible impact on most measurements so this is probably sufficient for most users. If you're trying to run a power sweep this may be insufficient.
@@ -64,11 +66,13 @@ Note that unlike the main calibration, power calibration frequencies do not need
 ## References
 
 #### Pluto Default Connection Settings
-user: `root`
-password: `analog`
-ip: `192.168.2.1`
+
+- user: `root`
+- password: `analog`
+- ip: `192.168.2.1`
 
 ## Alternatives
+
 - [NanoVNA](https://nanovna.com/). 2-ports. 50 kHz - 2.7 GHz. Degraded performance above 1.5 GHz. S11 and S21 only.
 - [pluto-network-analyzer](https://github.com/fromconcepttocircuit/pluto-network-analyzer). 2-ports. 100 MHz - 3 GHz. S11 and S21 only. Uses a [wideband RF bridge](https://www.60dbm.com/product/rf-bridge-1-3000-mhz/) instead of a coupler
 - [LibreVNA](https://github.com/jankae/LibreVNA). 2-ports. 100 KHz - 6 GHz. I've never used this but it is almost certainly faster than Charon. Not sure how the performance compares. $700 on [AliExpress](https://www.aliexpress.us/item/3256802242049773.html?gatewayAdapt=glo2usa4itemAdapt)

pyproject.toml

@@ -9,7 +9,7 @@ description = "RF Network Analyzer based on the Pluto SDR"
 readme = "README.md"
 requires-python = ">=3"
 # keywords = ["one", "two"]
-license = { text = "MIT License" }
+license = "MIT"
 classifiers = ["Programming Language :: Python :: 3"]
 dependencies = [
   "numpy",
@@ -34,7 +34,7 @@ charon-cli = "charon_vna.cli:main"
 charon-gui = "charon_vna.gui:main"
 
 [tool.setuptools_scm]
-version_file = "charon_vna/_version.py"
+version_file = "src/charon_vna/_version.py"
 
 [tool.black]
 line-length = 120

src/charon_vna/gui.py

@@ -56,6 +56,9 @@ class MainWindow(QMainWindow):
             vna_kwargs["ip"] = ip
         self.vna = Charon(**vna_kwargs)
 
+        self.active_port = 0
+        self.vna.set_switches(a=self.active_port, b=self.active_port)
+
         mpl.use("QtAgg")
 
         self.setWindowTitle("Charon VNA")
@@ -90,6 +93,9 @@ class MainWindow(QMainWindow):
         action_trigger.triggered.connect(self.capture)
         action_trigger.setShortcut("Ctrl+T")
         menu_stimulus.addAction(action_trigger)
+        action_p0 = QAction("Switch &Port", self)
+        action_p0.triggered.connect(self.toggle_port)
+        menu_stimulus.addAction(action_p0)
 
         menu_calibration = QMenu("&Calibration")
         menubar.addMenu(menu_calibration)
@@ -139,6 +145,11 @@ class MainWindow(QMainWindow):
         widget.setLayout(window_layout)
         self.setCentralWidget(widget)
 
+    def toggle_port(self):
+        self.active_port = int(not self.active_port)
+        print(f"Activating port {self.active_port}")
+        self.vna.set_switches(a=self.active_port, b=self.active_port)
+
     def saveas_config(self) -> None:
         print("Prompting for save path...")
         dialog = QFileDialog(self)
@@ -263,7 +274,7 @@ def main() -> None:
                 "Pluto IP Address",
                 "Enter Pluto IP Address",
                 QLineEdit.Normal,
-                "192.168.2.1",
+                Charon.DEFAULT_IP,
             )
             match = re.match(r"(\d{1,3}\.){3}\d{1,3}", text)
             if not match:

src/charon_vna/vna.py

@@ -1,8 +1,9 @@
 # %% imports
 import copy
+import pickle
 from enum import IntEnum, unique
 from pathlib import Path
-from typing import Any, Callable, Dict, Literal, Tuple
+from typing import Any, Callable, Dict, List, Literal, Tuple
 
 import adi
 import iio
@@ -82,7 +83,7 @@ class AD9361DacStepFactor(IntEnum):
 
 class Charon:
     FREQUENCY_OFFSET = 1e6
-    DEFAULT_IP = "192.168.2.1"
+    DEFAULT_IP = "192.168.3.1"
 
     calibration: rf.calibration.Calibration | None = None
 
@@ -90,8 +91,11 @@ class Charon:
         self,
         ip: str = DEFAULT_IP,
         frequency: npt.ArrayLike = np.linspace(1e9, 2e9, 3),
-        ports: Tuple[int] = (1,),
+        ports: Tuple[int] | int = 1,
     ):
+        if isinstance(ports, int):
+            ports = (np.arange(ports) + 1).tolist()
+        ports = tuple(ports)
         self.ports = ports
         self.frequency = frequency
 
@@ -128,18 +132,23 @@ class Charon:
         self.sdr.rx_hardwaregain_chan1 = 10
         self.sdr.tx_hardwaregain_chan0 = -10
 
-        # # switch control
+        # switch control
         ctx = iio.Context(uri)
         self.ctrl = ctx.find_device("ad9361-phy")
         # raw ad9361 register accesss:
         # https://ez.analog.com/linux-software-drivers/f/q-a/120853/control-fmcomms3-s-gpo-with-python
-        # https://www.analog.com/media/cn/technical-documentation/user-guides/ad9364_register_map_reference_manual_ug-672.pdf  # noqa: E501
+        # https://www.analog.com/media/cn/technical-documentation/user-guides/ad9364_register_map_reference_manual_ug-672.pdf
         self.ctrl.reg_write(AD9361Register.EXTERNAL_LNA_CONTROL, 0x90)  # bit 7: AuxDAC Manual, bit 4: GPO Manual
         self.ctrl.reg_write(AD9361Register.AUXDAC_ENABLE_CONTROL, 0x3F)
+
+        # initialize switch control outputs
         self._set_gpo(0b0000)
         self._set_dac_code(value=0, channel=1)
         self._set_dac_code(value=0, channel=2)
 
+        # set default switch state
+        self.set_switches(a=self.ports[0] - 1, b=self.ports[0] - 1)
+
     def get_config(self) -> Dict[str, Any]:
         config = dict()
         config["rx_lo"] = self.sdr.rx_lo
@@ -168,8 +177,15 @@ class Charon:
     def _set_gpo(self, value: int) -> None:
         self.ctrl.reg_write(AD9361Register.GPO_FORCE_AND_INIT, (value & 0x0F) << 4)  # bits 7-4: GPO3-0
 
-    def _set_dac_voltage(self, voltage: float, channel: Literal[1, 2]):
-        raise NotImplementedError()
+    def _get_dac_code(self, channel: Literal[1, 2]) -> Tuple[float, AD9361DacVref, AD9361DacStepFactor]:
+        word = self.ctrl.reg_read(AD9361Register.__getitem__(f"AUXDAC{channel}_WORD"))
+        config = self.ctrl.reg_read(AD9361Register.__getitem__(f"AUXDAC{channel}_CONFIG"))
+
+        value = (word << 2) + (config & 0x3)
+        vref = AD9361DacVref((config >> 2) & 0x3)
+        step_factor = AD9361DacStepFactor((config >> 4) & 0x1)
+
+        return (value, vref, step_factor)
 
     def _set_dac_code(
         self,
@@ -189,7 +205,12 @@ class Charon:
 
         # https://www.analog.com/media/cn/technical-documentation/user-guides/ad9364_register_map_reference_manual_ug-672.pdf
         # page 13
-        # vout = 0.97 * vref + (0.000738 + 9e-6 * (vref * 1.6 - 2)) * auxdac_word[9:0] * step_factor - 0.3572 * step_factor + 0.05
+        # vout = (
+        #     0.97 * vref
+        #     + (0.000738 + 9e-6 * (vref * 1.6 - 2)) * auxdac_word[9:0] * step_factor
+        #     - 0.3572 * step_factor
+        #     + 0.05
+        # )
         # vout ~= (vref - 0.3572 * step_factor) + 0.000738 * auxdac_word[9:0] * step_factor
         # which gives a 1.5V swing with step_factor == 2 and 0.75V swing with step_factor == 1
         # vref basically just changes the minimum voltage with negligible impact on output scaling
@@ -203,18 +224,17 @@ class Charon:
             (value & 0x3) | (vref.value << 2) | (step_factor << 4),
         )
 
-    def _get_dac_code(self, channel: Literal[1, 2]) -> Tuple[float, AD9361DacVref, AD9361DacStepFactor]:
-        word = self.ctrl.reg_read(AD9361Register.__getitem__(f"AUXDAC{channel}_WORD"))
-        config = self.ctrl.reg_read(AD9361Register.__getitem__(f"AUXDAC{channel}_CONFIG"))
+    def set_switches(self, b: int, a: int, excitation: int | None = None):
+        if excitation is None:
+            excitation = a
 
-        value = (word << 2) + (config & 0x3)
-        vref = AD9361DacVref((config >> 2) & 0x3)
-        step_factor = AD9361DacStepFactor((config >> 4) & 0x1)
+        val = 0
 
-        return (value, vref, step_factor)
+        val |= int(bool(excitation)) << 0  # exc = GPO0
+        val |= int(bool(a)) << 2  # a = GPO2
+        val |= int(bool(b)) << 1  # b = GPO1
 
-    def _get_dac_voltage(self) -> float:
-        raise NotImplementedError()
+        self._set_gpo(val)
 
     def set_output_power(self, power: float):
         pout = xr.DataArray(
@@ -291,50 +311,147 @@ class Charon:
 
         return np.mean(data[1] / data[0])
 
-    def sweep_b_over_a(self):
+    def capture(
+        self,
+        callback: Callable[[int, int], None] | None = None,
+        *,
+        measurements: List[Tuple[int, int]] = None,
+    ):
+        if measurements is None:
+            measurements = [(m, n) for n in self.ports for m in self.ports]
+        measurements = list(measurements)
+
         s = xr.DataArray(
             np.zeros(
-                len(self.frequency),
+                [len(self.frequency), len(self.ports), len(self.ports)],
                 dtype=np.complex128,
             ),
-            dims=["frequency"],
+            dims=["frequency", "m", "n"],
             coords=dict(
                 frequency=self.frequency,
+                m=list(self.ports),
+                n=list(self.ports),
             ),
         )
-        for frequency in self.frequency:
-            s.loc[dict(frequency=frequency)] = self.get_b_over_a(frequency=frequency)
-        return s
 
-    def vna_capture(self, frequency: npt.ArrayLike, callback: Callable[int, int] | None):
-        s = xr.DataArray(
-            np.empty(len(frequency), dtype=np.complex128),
-            dims=["frequency"],
-            coords=dict(
-                frequency=frequency,
-            ),
-        )
-        for ff, freq in enumerate(s.frequency.data):
-            if callback is not None:
-                # report progress during sweep
-                callback(ff, len(s.frequency))
-            self.set_output(frequency=freq, power=-5)
-            self.sdr.rx_destroy_buffer()
-            self.sdr.rx_lo = int(freq)
-            self.sdr.rx_enabled_channels = [0, 1]
-            self.sdr.gain_control_mode_chan0 = "manual"
-            self.sdr.gain_control_mode_chan1 = "manual"
-            self.sdr.rx_hardwaregain_chan0 = 40
-            self.sdr.rx_hardwaregain_chan1 = 40
-            rx = self.sdr.rx()
-            s.loc[dict(frequency=freq)] = np.mean(rx[1] / rx[0])
+        total_count = len(measurements) * len(s.frequency)
+        count = 0
+
+        for m in s.m.data:
+            for n in s.n.data:
+                if (m, n) in measurements:
+                    self.set_switches(b=m - 1, a=n - 1)
+
+                    for ff, freq in enumerate(s.frequency.data):
+                        if callback is not None:
+                            # report progress during sweep
+                            callback(count, total_count)
+
+                        self.set_output(frequency=freq, power=-5)
+                        self.sdr.rx_destroy_buffer()
+                        self.sdr.rx_lo = int(freq)
+                        self.sdr.rx_enabled_channels = [0, 1]
+                        self.sdr.gain_control_mode_chan0 = "manual"
+                        self.sdr.gain_control_mode_chan1 = "manual"
+                        self.sdr.rx_hardwaregain_chan0 = 40
+                        self.sdr.rx_hardwaregain_chan1 = 40
+                        rx = self.sdr.rx()
+                        s.loc[dict(frequency=freq, m=m, n=n)] = np.mean(rx[1] / rx[0])
+
+                        count += 1
 
         if callback is not None:
             # mark capture as complete
-            callback(len(s.frequency), len(s.frequency))
+            callback(total_count, total_count)
 
         return s
 
+    def calibrate_sol(self, prompt: Callable[[str], None] | None = None, **kwargs) -> None:
+        if len(self.ports) != 1:
+            raise ValueError(
+                f"SOL calibration needs only one port but {len(self.ports)} ports are enabled. "
+                "Did you mean to use SOLT?"
+            )
+
+        if prompt is None:
+            prompt = lambda s: input(f"{s}\nENTER to continue...")
+
+        ideal = rf.media.DefinedGammaZ0(frequency=rf.media.Frequency.from_f(self.frequency, unit="Hz"))
+        ideals = [ideal.short(), ideal.open(), ideal.load(0)]
+
+        names = ["short", "open", "load"]
+
+        measured = list()
+        for name in names:
+            prompt(f"Connect standard {name} to port {self.ports[0]}")
+            measured.append(self.capture(**kwargs))
+
+        cal = rf.OnePort(measured=[s2net(m) for m in measured], ideals=ideals)
+
+        self.calibration = cal
+
+    def calibrate_solt(self, prompt: Callable[[str], None] | None = None, **kwargs) -> None:
+        if len(self.ports) < 2:
+            raise ValueError(
+                f"SOLT calibration needs at least two ports but {len(self.ports)} ports are enabled. "
+                "Did you mean to use SOL?"
+            )
+
+        if len(self.ports) > 2:
+            raise NotImplementedError("SOLT calibration with more than two ports not yet supported")
+
+        if prompt is None:
+            prompt = lambda s: input(f"{s}\nENTER to continue...")
+
+        ideal = rf.media.DefinedGammaZ0(frequency=rf.media.Frequency.from_f(self.frequency, unit="Hz"))
+        ideals = [ideal.short(), ideal.open(), ideal.load(0)]
+        ideals = [rf.two_port_reflect(id, id) for id in ideals]
+
+        thru = np.zeros((len(self.frequency), 2, 2), dtype=np.complex128)
+        thru[:, 0, 1] = 1
+        thru[:, 1, 0] = 1
+        thru = rf.Network(frequency=self.frequency, f_unit="Hz", s=thru)
+
+        ideals.append(thru)
+
+        names_1p = ["short", "open", "load"]
+        names_2p = ["thru"]
+
+        measured = list()
+        for name in names_1p:
+            measured_param = list()
+            for port in self.ports:
+                prompt(f"Connect standard {name} to port {port}")
+                measured_param.append(self.capture(measurements=[(port, port)], **kwargs).sel(m=port, n=port))
+            measured.append(rf.two_port_reflect(*[s2net(m) for m in measured_param]))
+
+        for name in names_2p:
+            prompt(f"Connect standard {name} between ports {self.ports[0]} and {self.ports[1]}")
+            measured.append(s2net(self.capture(**kwargs)))
+
+        cal = rf.SOLT(measured=measured, ideals=ideals)
+
+        self.calibration = cal
+
+    def save_calibration(self, path: Path | str):
+        path = Path(path)
+        if path.suffix.lower() == ".pkl":
+            with open(str(path), "wb") as f:
+                pickle.dump(self.calibration, f)
+        else:
+            raise NotImplementedError(f"Unknown calibration file extension: {path.suffix}")
+
+    def load_calibration(self, path: Path | str):
+        path = Path(path)
+        if path.suffix.lower() == ".pkl":
+            with open(str(path), "rb") as f:
+                cal = pickle.load(f)
+                if not isinstance(cal, rf.calibration.Calibration):
+                    raise ValueError(f"Expected {rf.calibration.Calibration}, got {type(cal)}")
+                self.calibration = cal
+        else:
+            raise NotImplementedError(f"Unknown calibration file extension: {path.suffix}")
+
 
 # %%
 if __name__ == "__main__":

src/charon_vna/vna_dev.py

@@ -0,0 +1,81 @@
+# %% imports
+import numpy as np
+from matplotlib import pyplot as plt
+from matplotlib.ticker import EngFormatter
+
+from charon_vna.util import db20, net2s, s2net
+from charon_vna.vna import Charon
+
+# %%
+frequency = np.linspace(80e6, 280e6, 301)
+
+# %%
+vna = Charon(frequency=frequency, ports=2)
+
+# %%
+s = vna.capture()
+
+# %%
+for m in s.m.data:
+    for n in s.n.data:
+        plt.plot(s.frequency, db20(s.sel(m=m, n=n)), label="$S_{" + str(m) + str(n) + "}$")
+plt.grid(True)
+plt.legend()
+plt.show()
+
+
+# %%
+vna.calibrate_sol()
+
+# %%
+vna.calibrate_solt()
+
+# %%
+vna.save_calibration("./calibration.pkl")
+
+# %%
+vna.load_calibration("./calibration.pkl")
+
+# %%
+s2 = net2s(vna.calibration.apply_cal(s2net(s)))
+# s2.coords["m"] = s.m
+# s2.coords["n"] = s.n
+
+for m in s.m.data:
+    for n in s.n.data:
+        plt.plot(s.frequency, db20(s.sel(m=m, n=n)), label="$S_{" + str(m) + str(n) + "}$ (uncalibrated)")
+        plt.plot(s2.frequency, db20(s2.sel(m=m, n=n)), label="$S_{" + str(m) + str(n) + "}$ (calibrated)")
+plt.grid(True)
+plt.legend()
+plt.xlabel("Frequency [Hz]")
+plt.ylabel("Magnitude [dB]")
+# plt.ylim(-30, 5)
+plt.ylim(-25, 5)
+plt.xlim(s.frequency[0], s.frequency[-1])
+plt.gca().xaxis.set_major_formatter(EngFormatter())
+plt.tight_layout()
+plt.show()
+
+for m in s.m.data:
+    for n in s.n.data:
+        if m != n:
+            plt.plot(
+                s.frequency,
+                np.angle(s.sel(m=m, n=n), deg=True),
+                label="$S_{" + str(m) + str(n) + "}$ (uncalibrated)",
+            )
+            plt.plot(
+                s2.frequency,
+                np.angle(s2.sel(m=m, n=n), deg=True),
+                label="$S_{" + str(m) + str(n) + "}$ (calibrated)",
+            )
+plt.grid(True)
+plt.legend()
+plt.ylabel("Phase [deg]")
+plt.xlabel("Frequency [Hz]")
+plt.xlim(s.frequency[0], s.frequency[-1])
+plt.gca().xaxis.set_major_formatter(EngFormatter())
+plt.tight_layout()
+plt.show()
+
+# %%