.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "auto_examples/grp2_PyScriptExamples/plot_e00_SOL.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note :ref:`Go to the end ` to download the full example code. .. rst-class:: sphx-glr-example-title .. _sphx_glr_auto_examples_grp2_PyScriptExamples_plot_e00_SOL.py: SOL Calibration in Python ========================= This example a performs a short, open, and load (SOL) calibration with a linear sensitivity analysis. In this scenario, we have some raw S-parameter measurements stored in tab separated text files of our calibration standards and of our DUT - which is a load. We also have some data-defined models of our calibration standards stored in the HDF5 format. These models contain some linear uncertainty mechanisms which carry some category information describing their origin. This tutorial will walk through defining file reader functions, analysis functions, and performing the analysis with and without uncertainties. The data for this example is stored in th 'sol_demo_data' folder of the source code repository hosted on github. The data set is intended to act as an example use case of the software library only. .. GENERATED FROM PYTHON SOURCE LINES 20-28 .. code-block:: Python from pathlib import Path import h5py import xarray as xr import numpy as np import matplotlib.pyplot as plt import rmellipse as rme .. GENERATED FROM PYTHON SOURCE LINES 29-34 Demo Data --------- The sample data are store inside the `rmellipse.sol_demo` submodule, and the paths to the files can be imported as well. .. GENERATED FROM PYTHON SOURCE LINES 34-47 .. code-block:: Python # text files paths as Path objects # replace these with the correct data paths # after downloading the sample data short_raw_path = sample_data_dir / 'Short.s1p' open_raw_path = sample_data_dir / 'Open.s1p' load_raw_path = sample_data_dir / 'Load.s1p' dut_raw_path = sample_data_dir / 'Dut.s1p' # HDF5 file containing the standard models defs_path = sample_data_dir / 'definitions.hdf5' .. GENERATED FROM PYTHON SOURCE LINES 58-64 Writing Functions ----------------- The first thing we need to do is write some functions. These include the text file reader and analysis functions to calculate our error box and correct our raw data. .. GENERATED FROM PYTHON SOURCE LINES 67-73 Text File Reader ^^^^^^^^^^^^^^^^ First we define a function that can read our text files into a DataArray. We assign a dimension called Frequency (GHz), where we store the frequencies corresponding to each S-parameter. We also store our data as a complex-valued array. .. GENERATED FROM PYTHON SOURCE LINES 73-84 .. code-block:: Python def from_s1p(path) -> xr.DataArray: arr = np.loadtxt(path, float, delimiter='\t') # define a coordinate set coords = {'Frequency (GHz)': arr[:, 0]} # convert to 1 port complex data values = arr[:, 1] + 1.0j * arr[:, 2] # create an xarray data set out = xr.DataArray(values, coords=coords, dims=('Frequency (GHz)')) return out .. GENERATED FROM PYTHON SOURCE LINES 85-86 Lets confirm this works by reading in our raw measurements .. GENERATED FROM PYTHON SOURCE LINES 86-92 .. code-block:: Python raw_dut = from_s1p(dut_raw_path) raw_short = from_s1p(short_raw_path) raw_open = from_s1p(open_raw_path) raw_load = from_s1p(load_raw_path) .. GENERATED FROM PYTHON SOURCE LINES 93-104 Calculating the Error Box ^^^^^^^^^^^^^^^^^^^^^^^^^ We want two analysis functions. The first takes in our definitions and raw measurements, and calculates the error-box terms of the 1-port calibration. Importantly, all the import objects that might be ``RMEMeas`` objects are exposed as input arguments. The ``**`` packs any keyword arguments ino a dictionary. This function matches any inputs by matching device_names, and assuming keywords are formatted by 'def_' and 'meas_' only. Also, note that the function signature says that the function is expecting xr.DataArray objects. This is what will be passed in when the function is wrapped in a propagator. .. GENERATED FROM PYTHON SOURCE LINES 104-167 .. code-block:: Python def SOL_cal(**stds: xr.DataArray) -> xr.DataArray: # get list of standards and definitions, ordered to correspond defs = [] ms = [] for k, v in stds.items(): if 'def' in k: def_key = k m_key = def_key.replace('def_', 'raw_') defs.append(stds[def_key]) ms.append(stds[m_key]) N = len(defs) frequencies = ms[-1]['Frequency (GHz)'] # output has the same shape as the inputs except an # additional dimensions to hold the error terms output_shape = list(ms[0].shape) + [3] output_dims = list(ms[0].dims) + ['errterm'] output_coords = dict(ms[0].coords) output_coords.update({'errterm': ['e00', 'e11', 'delta']}) # pre allocate output xarray # 3 output has 3 result = np.zeros(output_shape, complex) result = xr.DataArray(result, dims=output_dims, coords=output_coords) # pre allocated temporary arrays # that will be used to solve the set of equations mshape = list(result.shape)[:-1] + [N, 3] M = np.zeros(mshape, complex) yshape = list(result.shape)[:-1] + [N, 1] y = np.zeros(yshape, complex) # I am going to work with the underlying numpy arrays # here because it is convenient for linear algebra # for each device add a row to the regressor matrix for i in range(N): S11_meas = ms[i].values S11_def = defs[i].values M[..., i, 0] = 1 M[..., i, 1] = S11_meas * S11_def M[..., i, 2] = -S11_def y[..., i, 0] = S11_meas # this transposes M along last 2 dimensions n_dims = len(M.shape) transpose_dims = np.arange(n_dims) transpose_dims[[n_dims - 1, n_dims - 2]] = transpose_dims[[n_dims - 2, n_dims - 1]] Mt = np.transpose(M, transpose_dims) # do least squares coeff = np.linalg.inv(Mt @ M) @ Mt @ y # reassign values to output result.loc[..., 'e00'] = coeff[..., 0, 0] result.loc[..., 'e11'] = coeff[..., 1, 0] result.loc[..., 'delta'] = coeff[..., 2, 0] return result .. GENERATED FROM PYTHON SOURCE LINES 168-171 Lets use the function without any uncertainty propagation to verify it works. I will do this by reading in the definitions and grabbing just a view into the underlying nominal DataArray. .. GENERATED FROM PYTHON SOURCE LINES 171-189 .. code-block:: Python with h5py.File(defs_path, 'r') as f: def_short = rme.load_object(f['Short']).nom def_open = rme.load_object(f['Open']).nom def_load = rme.load_object(f['Load']).nom flist = raw_short['Frequency (GHz)'] errbox = SOL_cal( def_short=def_short, def_open=def_open, def_load=def_load, raw_short=raw_short, raw_open=raw_open, raw_load=raw_load, ) print(errbox) .. GENERATED FROM PYTHON SOURCE LINES 190-195 Correcting Raw Data ^^^^^^^^^^^^^^^^^^^ Finally, lets write a function that takes in our error box and corrects a raw measurement. .. GENERATED FROM PYTHON SOURCE LINES 195-211 .. code-block:: Python def SOL_correct(errorbox: xr.DataArray, device: xr.DataArray) -> xr.DataArray: S11 = device e00 = errorbox.sel(errterm='e00') e11 = errorbox.sel(errterm='e11') delta = errorbox.sel(errterm='delta') corrected = (-e00 + S11) / (-delta + e11 * S11) # return the corrected result r = xr.zeros_like(device) r.loc[...] = corrected return r .. GENERATED FROM PYTHON SOURCE LINES 212-213 And lets use it to correct our raw device. .. GENERATED FROM PYTHON SOURCE LINES 213-216 .. code-block:: Python dut = SOL_correct(errbox, raw_dut) .. GENERATED FROM PYTHON SOURCE LINES 217-219 We can plot this to check that it matches our expectations. This DUT is a load so we expect a magnitude close to 0. .. GENERATED FROM PYTHON SOURCE LINES 219-231 .. code-block:: Python plt.close('all') fig, ax = plt.subplots(2, 1) fig.suptitle('DUT Corrected Data') ax[0].plot(dut['Frequency (GHz)'], np.abs(dut)) ax[0].set_ylabel('Linear Magnitude') ax[1].plot(dut['Frequency (GHz)'], np.angle(dut, deg=True)) ax[1].set_xlabel('Frequency (GHz)') ax[1].set_ylabel('Phase (deg)') .. GENERATED FROM PYTHON SOURCE LINES 232-236 Linear Uncertainty Analysis --------------------------- In this step, we will add the ability to propagate uncertainties associated with the calibration standard definitions. .. GENERATED FROM PYTHON SOURCE LINES 239-246 Import and Wrap Functions ^^^^^^^^^^^^^^^^^^^^^^^^^ The first step is to define a propagator and wrap our analysis functions in it. We are going to turn on the sensitivity analysis only. When we wrap our functions in the propagator, any RMEMeas objects passed through the function will be run through the linear sensitivity analysis. .. GENERATED FROM PYTHON SOURCE LINES 246-259 .. code-block:: Python myprop = rme.RMEProp(sensitivity=True) SOL_cal = myprop.propagate(SOL_cal) SOL_correct = myprop.propagate(SOL_correct) # turn xarray's into uncertainty objects with RMEMeas raw_dut = rme.RMEMeas.from_nom('DUT', raw_dut) raw_short = rme.RMEMeas.from_nom('DUT', raw_short) raw_open = rme.RMEMeas.from_nom('DUT', raw_open) raw_load = rme.RMEMeas.from_nom('DUT', raw_load) .. GENERATED FROM PYTHON SOURCE LINES 260-261 Lets also grab our definitions with the full uncertainty information .. GENERATED FROM PYTHON SOURCE LINES 261-266 .. code-block:: Python with h5py.File(defs_path, 'r') as f: def_short = rme.load_object(f['Short']) def_open = rme.load_object(f['Open']) def_load = rme.load_object(f['Load']) .. GENERATED FROM PYTHON SOURCE LINES 267-272 Propagate Functions ^^^^^^^^^^^^^^^^^^^ Now we can use our analysis functions that were wrapped in the propagator to do our analysis. .. GENERATED FROM PYTHON SOURCE LINES 272-285 .. code-block:: Python errbox = SOL_cal( def_short=def_short, def_open=def_open, def_load=def_load, raw_short=raw_short, raw_open=raw_open, raw_load=raw_load, ) dut = SOL_correct(errbox, raw_dut) .. GENERATED FROM PYTHON SOURCE LINES 286-293 Plot Results ------------ Lets plot our corrected device like last time, and lets instead inspect the uncertainty measurements. We would like to look at the magnitude and phase, so lets define some functions to propagate those conversions. .. GENERATED FROM PYTHON SOURCE LINES 295-297 Nominal Values ^^^^^^^^^^^^^^ .. GENERATED FROM PYTHON SOURCE LINES 297-353 .. code-block:: Python @myprop.propagate def calc_mag(arr): out = xr.zeros_like(arr, dtype=float) out.values = np.abs(arr) return out @myprop.propagate def calc_phase(arr): out = xr.zeros_like(arr, dtype=float) out.values = np.angle(arr, deg=True) return out mag = calc_mag(dut) phase = calc_phase(dut) k = 2 fig, ax = plt.subplots(2, 1) mag_lower = mag.uncbounds(k=k).cov mag_upper = mag.uncbounds(k=-k).cov phs_lower = phase.uncbounds(k=k, deg=True).cov phs_upper = phase.uncbounds(k=-k, deg=True).cov ax[0].fill_between( dut.nom['Frequency (GHz)'], y1=mag_lower, y2=mag_upper, color='k', alpha=0.5, label=f'k = {k} Uncertainty', ) ax[1].fill_between( dut.nom['Frequency (GHz)'], y1=phs_lower, y2=phs_upper, alpha=0.5, color='k', label=f'k = {k} Uncertainty', ) ax[0].plot(mag.nom['Frequency (GHz)'], mag.nom, label='nominal') ax[0].set_ylabel('Linear Magnitude') ax[1].plot(phase.nom['Frequency (GHz)'], phase.nom, label='nominal') ax[1].set_xlabel('Frequency (GHz)') ax[1].set_ylabel('Phase (deg)') handles, labels = ax[0].get_legend_handles_labels() fig.legend(handles, labels, loc='upper center', ncols=3) .. GENERATED FROM PYTHON SOURCE LINES 354-361 Uncertainties Only ^^^^^^^^^^^^^^^^^^ Lets look at the k=2 uncertainties only. Notice how the phase uncertainty increases as the magnitude approaches zero. This is natural, as phase is not well defined if your magnitude is on the origin. .. GENERATED FROM PYTHON SOURCE LINES 361-371 .. code-block:: Python fig, ax = plt.subplots(2, 1) ax[0].plot( mag.nom['Frequency (GHz)'], mag.stdunc(k=k).cov, ) ax[0].set_ylabel(f'Lin Magn k={k} Uncertainty') ax[1].plot(phase.nom['Frequency (GHz)'], phase.stdunc(k=k).cov, label='nominal') ax[1].set_xlabel('Frequency (GHz)') ax[1].set_ylabel(f'Phase k={k} Uncertainty (deg)') .. _sphx_glr_download_auto_examples_grp2_PyScriptExamples_plot_e00_SOL.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: plot_e00_SOL.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: plot_e00_SOL.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: plot_e00_SOL.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_