charon_vna/charon_vna/vna.py

# %% imports
import copy
from pathlib import Path
from typing import Optional

import adi
import iio
import numpy as np
import skrf as rf
import xarray as xr
from matplotlib import pyplot as plt
from matplotlib.ticker import EngFormatter
from numpy import typing as npt
from scipy import signal
from util import HAM_BANDS, db20, net2s, s2net

dir_ = Path(__file__).parent


# https://wiki.analog.com/resources/tools-software/linux-drivers/iio-transceiver/ad9361-customization


# %% helper functions
def get_config(sdr: adi.ad9361):
    config = dict()
    config["rx_lo"] = sdr.rx_lo
    config["rx_rf_bandwidth"] = sdr.rx_rf_bandwidth
    config["rx_enabled_channels"] = sdr.rx_enabled_channels
    for chan in config["rx_enabled_channels"]:
        config[f"rx_hardwaregain_chan{chan}"] = getattr(sdr, f"rx_hardwaregain_chan{chan}")
        config[f"gain_control_mode_chan{chan}"] = getattr(sdr, f"gain_control_mode_chan{chan}")

    config["tx_lo"] = sdr.tx_lo
    config["tx_rf_bandwidth"] = sdr.tx_rf_bandwidth
    config["tx_cyclic_buffer"] = sdr.tx_cyclic_buffer
    config["tx_enabled_channels"] = sdr.tx_enabled_channels
    for chan in config["tx_enabled_channels"]:
        config[f"tx_hardwaregain_chan{chan}"] = getattr(sdr, f"tx_hardwaregain_chan{chan}")

    config["filter"] = sdr.filter
    config["sample_rate"] = sdr.sample_rate
    config["loopback"] = sdr.loopback

    return config


def generate_tone(f: float, N: int = 1024, fs: Optional[float] = None):
    if fs is None:
        fs = sdr.sample_rate
    fs = int(fs)
    fc = int(f / (fs / N)) * (fs / N)
    ts = 1 / float(fs)
    t = np.arange(0, N * ts, ts)
    i = np.cos(2 * np.pi * t * fc) * 2**14
    q = np.sin(2 * np.pi * t * fc) * 2**14
    iq = i + 1j * q

    return iq


# %% connection
sdr = adi.ad9361(uri="ip:192.168.3.1")

# %% verify device configuration
mode_2r2t = bool(sdr._get_iio_debug_attr("adi,2rx-2tx-mode-enable"))
if not mode_2r2t:
    raise ValueError("'adi,2rx-2tx-mode-enable' is not set in pluto. See README.md for instructions for changing this")
# TODO: it might be possible to change this on the fly. I think we'll actually just fail in __init__ for sdr

# %% switch control outputs
# NOTE: this doesn't appear to work
sdr._set_iio_debug_attr_str("adi,gpo-manual-mode-enable", "1")
sdr._get_iio_debug_attr_str("adi,gpo-manual-mode-enable-mask")
# but direct register access does
# https://ez.analog.com/linux-software-drivers/f/q-a/120853/control-fmcomms3-s-gpo-with-python
ctx = iio.Context("ip:192.168.3.1")
ctrl = ctx.find_device("ad9361-phy")
# https://www.analog.com/media/cn/technical-documentation/user-guides/ad9364_register_map_reference_manual_ug-672.pdf
ctrl.reg_write(0x26, 0x90)  # bit 7: AuxDAC Manual, bit 4: GPO Manual
ctrl.reg_write(0x27, 0x10)  # bits 7-4: GPO3-0


# %% initialization
sdr.rx_lo = int(2.0e9)
sdr.tx_lo = int(2.0e9)
sdr.sample_rate = 30e6
sdr.rx_rf_bandwidth = int(4e6)
sdr.tx_rf_bandwidth = int(4e6)
sdr.rx_destroy_buffer()
sdr.tx_destroy_buffer()
sdr.rx_enabled_channels = [0, 1]
sdr.tx_enabled_channels = [0]
sdr.loopback = 0
sdr.gain_control_mode_chan0 = "manual"
sdr.gain_control_mode_chan1 = "manual"
sdr.rx_hardwaregain_chan0 = 40
sdr.rx_hardwaregain_chan1 = 40
sdr.tx_hardwaregain_chan0 = -10

config = get_config(sdr)
config

# %%
sdr.tx_destroy_buffer()  # must destroy buffer before changing cyclicity
sdr.tx_cyclic_buffer = True
sdr.tx(generate_tone(1e6))

# %%
sdr.rx_destroy_buffer()
data = sdr.rx()

# %% Plot in time
fig, axs = plt.subplots(2, 1, sharex=True, tight_layout=True)
axs[0].plot(np.real(data).T)
axs[1].plot(np.imag(data).T)
axs[0].set_ylabel("Real")
axs[1].set_ylabel("Imag")
axs[-1].set_xlabel("Time")
fig.show()

# %% Plot in frequency
f, Pxx_den = signal.periodogram(data, sdr.sample_rate, axis=-1, return_onesided=False)
plt.figure()
for cc, chan in enumerate(sdr.rx_enabled_channels):
    plt.semilogy(f, Pxx_den[cc], label=f"Channel {chan}")
plt.legend()
plt.ylim([1e-7, 1e2])
plt.xlabel("frequency [Hz]")
plt.ylabel("PSD [V**2/Hz]")
plt.grid(True)
plt.show()


# %% TX helper functions
def set_output_power(power: float):
    if power == -5:
        # FIXME: this is a hack because I don't want to go through re-calibration
        tx_gain = -8
    else:
        raise NotImplementedError()
        # # TODO: correct over frequency
        # tx_gain_idx = np.abs(pout.sel(tx_channel=0) - power).argmin(dim="tx_gain")
        # tx_gain = pout.coords["tx_gain"][tx_gain_idx]
    sdr.tx_hardwaregain_chan0 = float(tx_gain)


def set_output(frequency: float, power: float, offset_frequency: float = 1e6):
    sdr.tx_destroy_buffer()
    set_output_power(power)
    sdr.tx_lo = int(frequency - offset_frequency)
    offset_frequency = frequency - sdr.tx_lo
    sdr.tx_cyclic_buffer = True
    sdr.tx(generate_tone(offset_frequency))


# %%
def vna_capture(frequency: npt.ArrayLike):
    s = xr.DataArray(
        np.empty(len(frequency), dtype=np.complex128),
        dims=["frequency"],
        coords=dict(
            frequency=frequency,
        ),
    )
    for freq in s.frequency.data:
        set_output(frequency=freq, power=-5)
        sdr.rx_destroy_buffer()
        sdr.rx_lo = int(freq)
        sdr.rx_enabled_channels = [0, 1]
        sdr.gain_control_mode_chan0 = "manual"
        sdr.gain_control_mode_chan1 = "manual"
        sdr.rx_hardwaregain_chan0 = 40
        sdr.rx_hardwaregain_chan1 = 40
        rx = sdr.rx()
        s.loc[dict(frequency=freq)] = np.mean(rx[1] / rx[0])

    return s


# %%
s = vna_capture(frequency=np.linspace(70e6, 200e6, 101))

# %% Plot Logmag
fig, axs = plt.subplots(2, 1, sharex=True, tight_layout=True)

axs[0].plot(s.frequency, db20(s), label="Measured")
axs[1].plot(s.frequency, np.rad2deg(np.angle((s))), label="Measured")

axs[0].grid(True)
axs[1].grid(True)

axs[0].set_ylim(-80, 0)
axs[1].set_ylim(-200, 200)
axs[1].set_xlim(np.min(s.frequency), np.max(s.frequency))
axs[1].xaxis.set_major_formatter(EngFormatter(places=1))
axs[1].set_xlabel("Frequency")

axs[0].set_ylabel("|S11| [dB]")
axs[1].set_ylabel("∠S11 [deg]")

reference_sparams = None
reference_sparams = dir_ / "RBP-135+_Plus25degC.s2p"
if reference_sparams is not None:
    ref = rf.Network(reference_sparams)
    rbp135 = net2s(ref)

    axs[0].plot(rbp135.frequency, db20(rbp135.sel(m=1, n=1)), label="Datasheet")
    axs[1].plot(rbp135.frequency, np.rad2deg(np.angle(rbp135.sel(m=2, n=1))), label="Datasheet")
    axs[0].legend()
    axs[1].legend()

plt.show()


# %% SOL calibration
cal_frequency = np.linspace(70e6, 600e6, 101)
ideal_cal_frequency = rf.Frequency(np.min(cal_frequency), np.max(cal_frequency), len(cal_frequency))
input("Connect SHORT and press ENTER...")
short = vna_capture(frequency=cal_frequency)
input("Connect OPEN and press ENTER...")
open = vna_capture(frequency=cal_frequency)
input("Connect LOAD and press ENTER...")
load = vna_capture(frequency=cal_frequency)

short_net = s2net(short)
open_net = s2net(open)
load_net = s2net(load)

cal_ideal = rf.media.DefinedGammaZ0(frequency=ideal_cal_frequency)
calibration = rf.calibration.OnePort(
    [short_net, open_net, load_net],
    [cal_ideal.short(), cal_ideal.open(), cal_ideal.load(0)],
)


# %%
s = vna_capture(frequency=cal_frequency)

# %%
s_calibrated = calibration.apply_cal(s2net(s))

plt.figure()
s_calibrated.plot_s_smith()
# ref.plot_s_smith(m=1, n=1)
plt.show()

plt.figure()
for start, stop in HAM_BANDS:
    plt.axvspan(start, stop, alpha=0.1, color="k")
s_calibrated.plot_s_db()
# ref.plot_s_db(m=1, n=1)
plt.gca().xaxis.set_major_formatter(EngFormatter())
plt.grid(True)
plt.xlim(s_calibrated.f[0], s_calibrated.f[-1])
plt.show()

plt.figure()
for start, stop in HAM_BANDS:
    plt.axvspan(start, stop, alpha=0.1, color="k")
# s_calibrated.plot_s_vswr()
# drop invalid points
vswr = copy.deepcopy(s_calibrated.s_vswr[:, 0, 0])
vswr[vswr < 1] = np.nan
plt.plot(s_calibrated.f, vswr)
plt.axhline(1, color="k", linestyle="--")
plt.ylabel("VSWR")
plt.xlabel("Frequency [Hz]")
# ref.plot_s_vswr(m=1, n=1)
plt.gca().xaxis.set_major_formatter(EngFormatter())
plt.grid(True)
plt.ylim(0, 10)
plt.xlim(s_calibrated.f[0], s_calibrated.f[-1])
plt.show()

# %%