r"""
Qfactor (:mod:`skrf.qfactor`)
========================================
Module for estimating the Quality (Q) factor(s) from S-parameters.
This class implements methods for determining *loaded* and *unloaded* Q-factor
from frequency-domain S-parameters, that can be applied to measurements
of transmission or reflection.
Documentation and implementation are adapted from [MAT58]_
Q-factor
--------
.. autosummary::
:toctree: generated/
:nosignatures:
:recursive:
Qfactor
.. currentmodule:: skrf.qfactor
Loaded and Unloaded Q-factor
----------------------------
The Quality factor (Q-factor) of a resonator is defined by [Pozar]_:
.. math::
Q = \frac{2 \pi U}{\Delta U}
where :math:`U` is the average energy stored by the resonator and
:math:`\Delta U` is the decrease in the average stored energy per wave cycle
at the resonant frequency, that is, the average power loss.
The *loaded* Q-factor, :math:`Q_L`, describes energy dissipation within the
entire resonant system comprising of the resonator itself and the instrument
used for observing resonances. The term loading refers to the effect that the
external circuit has on measured quantities.
The external circuit consists of the measuring instrument and uncalibrated lines,
but not the couplings of microwave resonators. Loading by an instrument that
has 50 Ohm impedance, such as a VNA, causes :math:`Q_L` to be reduced substantially
if strong coupling is used.
For many applications the quantity that is desired is the *unloaded* Q-factor :math:`Q_0`,
which is determined by energy dissipation associated with the resonator only
and therefore gives the best description of the resonant mode(s).
In other words, :math:`Q_0` is the Q-factor of the uncoupled resonator. The value of
:math:`Q_0` can be estimated from measurements of :math:`Q_L`, but cannot be measured directly.
:math:`Q_0` is largely governed by ohmic loss arising from surface currents
in the metal conductors (walls and loop couplings), and from dielectric loss
in any insulating materials that may be present.
Relationships between Loaded and Unloaded Q-factors
---------------------------------------------------
Energy dissipation in the external circuit is characterised by the *external* Q-factor,
:math:`Q_e`. For both series and parallel equivalent circuits, the three
Q-factors are related by [Pozar]_:
.. math::
\frac{1}{Q_L} = \frac{1}{Q_0} + \frac{1}{Q_e}
The coupling factor :math:`\beta` is defined for each port as:
.. math::
\beta = \frac{Q_0}{Q_e}
Finding the unloaded Q from measured S-parameters consists in first finding
the coupling factor, then measure :math:`Q_L` from the 3 dB bandwidth
and using the relationships above.
Fortunately, scikit-rf implements methods for determining loaded and
unloaded Q-factors from frequency-domain S-parameters. The implemented methods
are described in detail in [MAT58]_, and can be applied to measurements of
transmission or reflection.
Q-factor determination through equivalent-circuit models
--------------------------------------------------------
Characterisation of resonances from measurements in the frequency-domain
can be achieved through equivalent-circuit models [MAT58]_. Resonators can be
modelled as an ideal RLC resonator connected to an external circuit,
incorporating elements to account for a lossy coupling and coupling reactances.
For high Q-factor resonators (in practice, :math:`Q_L` > 100), the S-parameter
response of a resonator measured in a calibrated system with reference planes
at the resonator couplings can be expressed like [MAT58]_, [Galwas]_ :
.. math::
S = S_D + d \frac{e^{−2j\delta}}{1 + j Q_L t}
where :math:`S_D` is the detuned S-parameter measured at frequencies far above or below
resonance, :math:`d` is The diameter of the Q-circle, :math:`\delta` is a
real-valued constant that defines the orientation of the Q-circle, and :math:`t`
is the fractional offset frequency given by:
.. math::
t = \frac{f}{f_L} - \frac{f_L}{f} \approx 2 \frac{f − f_L}{f_L}
where :math:`f_L` is the loaded resonant frequency and :math:`f` the frequency
at which S is measured. This equation can be applied to measurements
by transmission (S21 or S12) or reflection (S11 or S22).
The S-parameters are fitted against a modified expression of the above equations
to deduce the resonant frequency, loaded and unloaded Q-factors and other properties.
References
----------
.. [MAT58] "Q-factor Measurement by using a Vector Network Analyser",
A. P. Gregory, National Physical Laboratory Report MAT 58 (2021)
https://eprintspublications.npl.co.uk/9304/
.. [Pozar] D. M. Pozar, Microwave engineering, 4th ed. J. Wiley, 2012.
.. [Galwas] B. A. Galwas, ‘Scattering Matrix Description of Microwave
Resonators’, IEEE Trans. Microwave Theory Techn., vol. 31, no. 8,
pp. 669–671, Aug. 1983, doi: 10.1109/TMTT.1983.1131566.
"""
from __future__ import annotations
from warnings import warn
import numpy as np
from .constants import NumberLike
from .frequency import Frequency
from .network import Network
# Available resonance types
RESONANCE_TYPES = ['reflection', 'reflection_method2',
'transmission', 'absorption']
class OptimizedResult(dict):
"""Represent Q-factor optimisation result.
Attributes
----------
Q_L : float
Loaded Quality factor
f_L : float
Resonance frequency [Hz]
success: bool
Is the fit method has been successfully performed
method: str
Fitting method used.
m1, m2, m3, ...: float
Coefficients described in [MAT58]_
number_iterations: int
Number of iterations performed.
Notes
-----
`OptimizedResult` may have additional attributes not listed here depending
on the specific fitting method used. Since this class is essentially a
subclass of dict with attribute accessors, one can see which
attributes are available using the `OptimizedResult.keys` method.
References
----------
.. [MAT58] "Q-factor Measurement by using a Vector Network Analyser",
A. P. Gregory, National Physical Laboratory Report MAT 58 (2021)
https://eprintspublications.npl.co.uk/9304/
"""
def __getattr__(self, name):
try:
return self[name]
except KeyError as e:
raise AttributeError(name) from e
__setattr__ = dict.__setitem__
__delattr__ = dict.__delitem__
def __repr__(self):
if self.keys():
m = max(map(len, list(self.keys()))) + 1
return '\n'.join([k.rjust(m) + ': ' + repr(v)
for k, v in sorted(self.items())])
else:
return self.__class__.__name__ + "()"
def __dir__(self):
return list(self.keys())
[docs]
class Qfactor:
"""
Q-factor calculation class.
This class implements methods for determining *loaded* and *unloaded* Q-factor
from frequency-domain S-parameters, that can be applied to measurements
of transmission or reflection.
Parameters
----------
ntwk : :class:`~skrf.network.Network` object
A 1-port scikit-rf Network.
If your device is a N-port, pass the desired sub-S-parameters to fit
the data from, like `ntwk.s21`.
res_type : str
Specifies the resonance type: 'reflection', 'transmission',
'reflection_method2' or 'absorption':
'reflection' is generally suited for undercoupled resonators,
while 'reflection_method2' is favoured for coupling with large loop.
Q_L0 : float, optional. Default is None.
Estimated loaded Q-factor, used to improve fitting.
f_L0 : float, optional. Default is None.
Estimated loaded resonant frequency, used to improve fitting [Hz]
If None, automatically search for the min or max,
depending on the resonance type defined by `res_type`.
verbose : bool, optional. Default is False.
Boolean flag controlling output of information to the console.
Raises
------
ValueError
If the passed Network is not 1-port or if the resonance type is unknown.
Notes
-----
Uncalibrated line should be de-embedded (if it has significant
length) from the S-parameter data before calling the functions
in this module to get best results.
"""
[docs]
def __init__(self,
ntwk: Network,
res_type: str,
Q_L0: None | float = None,
f_L0: None | float = None,
verbose: bool = False):
"""Q-factor initializer."""
# check ntwk is a 1-port
if ntwk.nports != 1:
raise ValueError('The Network is not a 1-port Network.')
if res_type not in RESONANCE_TYPES:
raise ValueError(f'res_type must be in: {RESONANCE_TYPES}.')
self.s = ntwk.s
self.f = ntwk.f
self.f_scaled = ntwk.frequency.f_scaled
self.f_multiplier = ntwk.frequency.multiplier
self.f_unit = ntwk.frequency.unit
self._ntwk = ntwk
self.res_type = res_type
self.tol = 1.0e-5
self.verbose = verbose
self.fitted = False
self.opt_res = None
self.N = len(self.f)
# step 1: initial_fit. Deduce premilinary values for Q_L and f_L.
self._initial_fit(self.N, Q_L0, f_L0)
def __str__(self) -> str:
if self.fitted:
status = f"fitted: f_L={float(self.f_L/self.f_multiplier):.3f}{self.f_unit}, Q_L={float(self.Q_L):.3f}"
else:
status = 'not fitted'
_str = f"Q-factor of Network {self._ntwk.name}. ({status})"
return _str
def __repr__(self) -> str:
return self.__str__()
[docs]
def fit(self,
method: str = "NLQFIT6",
loop_plan: str = 'fwfwc'
) -> OptimizedResult:
"""Fit Q-factor from S-parameter data.
Fitting overwrites the parameters `Q_L` and `f_L`.
Parameters
----------
method : str, optional
Fitting method : 'NLQFIT6' (default), 'NLQFIT7', 'NLQFIT8':
- 'NLQFIT6': Least Square Minimum of Eq.21 [MAT58]_ with 6 coeffcients.
- 'NLQFIT7': Least Square Minimum of Eq.26 [MAT58]_ with 7 coeffcients,
including one that characterize the trans. line length.
- 'NLQFIT8': Least Square Minimum of Eq.43 [MAT58]_ with 8 coeffcients,
A model for frequency-dependent leakage.
loop_plan : str, optional
Defines order of steps used by the fitting process.
The convergence algorithm uses a number of steps set by loop_plan,
a string of characters as follows:
- 'f': fit once without testing for convergence
- 'c': repeated fit, iterating until convergence is obtained
- 'w': re-calculate weighting factors on basis of previous fit
- Initially the weighting factors are all unity.
- The first character in `loop_plan` must not be 'w'.
e.g.: 'fwfwc' (default).
Returns
-------
result : :class:`~skrf.qfactor.OptimizedResult`
References
----------
.. [MAT58] "Q-factor Measurement by using a Vector Network Analyser",
A. P. Gregory, National Physical Laboratory Report MAT 58 (2021)
https://eprintspublications.npl.co.uk/9304/
"""
for op in loop_plan:
if op not in ['f', 'c', 'w']:
raise ValueError("Unexpected character in loop_plan")
if loop_plan[-1] == "w":
raise ValueError("Last item in loop_plan must not be w (weight calculation)")
if loop_plan[0] == "w":
raise ValueError("First item in loop_plan must not be w (weight calculation)")
if loop_plan[-1] != "c":
warn("Last item in loop_plan is not c so convergence not tested!", stacklevel=2)
self.method = method
self.loop_plan = loop_plan
# step 2: least square fitting
if method == 'NLQFIT6':
result = self._optimise_fit6(self.N)
elif method == 'NLQFIT7':
result = self._optimise_fit7(self.N)
elif method == 'NLQFIT8':
result = self._optimise_fit8(self.N)
# overwrite results in self
self.Q_L = result.Q_L
self.f_L = result.f_L
self.fitted = True
self.opt_res = result
if result.Q_L < 0:
warn('Negative Q_L, fitting may be inaccurate.', stacklevel=2)
return result
[docs]
@staticmethod
def angular_weights(f: NumberLike,
f_L: NumberLike,
Q_L: NumberLike
) -> NumberLike:
r"""Diagonal elements W_i of weights matrix.
.. math::
W_i = \frac{1}{\left[ \frac{2 Q_L (f_i - f_L)}{f_L} \right]^2 + 1}
The weights are needed when frequencies are equally-spaced
(rather than points equally spaced around the Q-circle), and help
reducing systematic error [MAT58]_.
Parameters
----------
f : np.ndarray
Frequency values array.
f_L : float
Loaded resonant frequency.
Q_L : float
Loaded Q-factor.
Returns
-------
W_i : np.ndarray
Weighting factors in proportion to the rate of change of angle
with frequency relative to the centre of the Q-circle.
References
----------
.. [MAT58] "Q-factor Measurement by using a Vector Network Analyser",
A. P. Gregory, National Physical Laboratory Report MAT 58 (2021)
https://eprintspublications.npl.co.uk/9304/
section 2.4, eqn. (28).
"""
ptmp = 2 * Q_L * (f - f_L) / f_L
W_i = 1 / (ptmp ** 2 + 1)
return W_i
def _initial_fit(self,
N: int,
Q_L0: float | None = None,
f_L0: None | float = None
):
"""Initial Linear least squares Q-factor fit.
As this is not optimised in this function (use `fit`), the solution
will only be approximate. This method is called during the
initialization of the `Qfactor` class. Note that a reasonable estimate
for the resonant frequency should be supplied is multiple resonances
are present.
Also calculate the internal parameters a, b and QL from ([MAT58]_, eqn. 17)
Parameters
----------
N : int
Number of points.
Q_L0 : float, optional. Default is None.
Estimated loaded Q-factor (will be improved by fitting).
f_L0 : float, optional. Default is None.
Estimated loaded resonant frequency, used to improve fitting [Hz]
Returns
-------
None
References
----------
.. [MAT58] "Q-factor Measurement by using a Vector Network Analyser",
A. P. Gregory, National Physical Laboratory Report MAT 58 (2021)
https://eprintspublications.npl.co.uk/9304/
section 2.1, eqn. (17).
"""
if f_L0 is None:
# search for the initial value of the resonance frequency
if self.res_type in ['reflection', 'reflection_method2', 'absorption']:
# Find minimum in |S11|
index_min = np.argmin(np.abs(self.s))
f_L0 = self.f[index_min]
else:
# Find peak in |S21|
index_max = np.argmax(np.abs(self.s))
f_L0 = self.f[index_max]
# Q_L0 : An order-of-magnitude estimate for Q_L-factor
if Q_L0 is None:
# The value 5.0 should work well
# for initial and optimised fits (Section 2.6).
mult = 5.0
Q_L0 = mult * f_L0/(self.f[-1] - self.f[0])
if self.verbose:
print(f'Initial estimation: Q_L0={Q_L0}, f_L0={f_L0}')
N2 = 2 * N
M = np.zeros([N2, 5])
G = np.zeros(N2)[:, np.newaxis]
for i in range(N):
i2 = i + N
t = 2.0 * (self.f[i] / f_L0 - 1.0)
y = 1.0 / complex(1.0, Q_L0 * t)
v = t * y
v1 = y * self.s[i, 0, 0]
G[i] = v1.real
G[i2] = v1.imag
v2 = v1 * t
M[i, :] = np.array([v.real, -v.imag, y.real, -y.imag, v2.imag])
M[i2, :] = np.array([v.imag, v.real, y.imag, y.real, -v2.real])
T = M.transpose() # unweighted
C = T @ M
q = T @ G
sv = np.linalg.solve(C, q)
a_re, a_im, b_re, b_im, Q_L = sv
self._a = (a_re + 1j*a_im)[0]
self._b = (b_re + 1j*b_im)[0]
self.Q_L = Q_L[0]
self.f_L = f_L0
if self.verbose:
print(f'Preliminary estimation: Q_L={self.Q_L}, f_L={self.f_L}')
def _optimise_fit6(self, N: int):
"""Iterative non-linear fit, NLQFIT6 Step (2).
Optimised fit of Q-factor (Q_L) and resonant frequency (f_L)
by the gradient-descent method [MAT58]_.
Uses the results of the initial fit as the starting
values for the iteration.
Parameters
----------
N : int
Number of points.
Returns
-------
res : OptimizedResult
Fitted values.
References
----------
.. [MAT58] MAT 58, section 2.2, eqn. (22).
"""
N2 = N * 2
iterations = 0
PV = np.ones(N) # default weights vector
PV2 = np.ones(N2)
m1 = self._a.imag / self.Q_L # a''/QL
m2 = -self._a.real / self.Q_L
m3 = self._b.real - m1
m4 = self._b.imag - m2
m5 = self.Q_L
Flwst = self.f[0] # lowest freq. is a convenient normalisation factor.
m6 = Flwst * m5 / self.f_L
last_op = "n"
weighting_ratio = None
number_iterations = 0
## Loop through all of the operations specified in loop_plan
for op in self.loop_plan:
if op == "w":
PV = self.angular_weights(self.f, Flwst * float(m5) / float(m6), float(m5))
# PV = self.angular_weights(m5)
weighting_ratio = max(PV) / min(PV)
PV2 = np.concatenate((PV, PV))
if self.verbose:
print("Op w, Calculate weights")
last_op = "n"
elif op == "c":
seek_convergence = True
elif op == "f":
seek_convergence = False
TerminationConditionMet = False
RMS_Error = None
while not (TerminationConditionMet):
number_iterations += 1
M = np.zeros([N2, 6]) # X is the transpose of M
G = np.zeros(N2)[:, np.newaxis]
c1 = complex(-m4, m3)
c2 = complex(m1, m2)
c3 = complex(m3, m4)
for i in range(N):
i2 = i + N
y = 1.0 / complex(1.0, 2 * (m6 * self.f[i] / Flwst - m5))
u = c1 * y * y * 2
u2 = -u * self.f[i] / Flwst
M[i, :] = 1.0, 0.0, y.real, -y.imag, u.real, u2.real
M[i2, :] = 0.0, 1.0, y.imag, y.real, u.imag, u2.imag
v = c2 + c3 * y
r = self.s[i] - v # residual
G[i] = r.real
G[i2] = r.imag
X = M.transpose()
T = np.multiply(X, PV2)
C = np.dot(T, M)
q = np.dot(T, G)
dm = np.linalg.solve(C, q)[:,0]
m1 += dm[0]
m2 += dm[1]
m3 += dm[2]
m4 += dm[3]
m5 += dm[4]
m6 += dm[5]
del G, X, T, C, dm
iterations = iterations + 1
if RMS_Error is not None:
Last_RMS_Error = RMS_Error
else:
Last_RMS_Error = None
SumNum = 0.0
SumDen = 0.0
for i in range(N):
den = complex(1.0, 2 * (m6 * self.f[i] / Flwst - m5))
ip = PV[i]
E = self.s[i] - complex(m1, m2) - complex(m3, m4) / den
SumNum = SumNum + ip * (E.real * E.real + E.imag * E.imag)
SumDen = SumDen + ip
RMS_Error = np.sqrt(SumNum / SumDen)
if self.verbose:
if last_op == "c":
print(f"Iteration {iterations}, RMS Error: {RMS_Error}")
else:
print(f"op {op}, Iteration {iterations}, RMS Error: {RMS_Error}")
last_op = op
if seek_convergence:
if Last_RMS_Error is not None:
delta_S = abs(RMS_Error - Last_RMS_Error)
TerminationConditionMet = delta_S < self.tol
else:
TerminationConditionMet = True
# After last operation, we end up here ...
if self.verbose:
print("Optimization done.")
return OptimizedResult({
'success': TerminationConditionMet,
'm1': m1, 'm2': m2, 'm3': m3, 'm4': m4,
'Q_L': m5,
'f_L': m5 * Flwst / m6,
'weighting_ratio': weighting_ratio,
'number_iterations': number_iterations,
'RMS_Error': RMS_Error,
'method': self.method,
})
def _optimise_fit7(self, N):
"""Iterative non-linear fit, NLQFIT7 Step (2).
Optimised fit of Q-factor (QL) and resonant frequency (FL)
by the gradient-descent method.
Uses the results of the initial fit (sv) as the starting
values for the iteration.
Parameters
----------
F : np.ndarray
Frequency points
S : np.ndarray
Complex data S-parameter to be fitted.
N : int
Number of points.
Fseed : float
Estimated resonant frequency.
sv : list
Initial solution (numpy vector or a list) found with initial_fit.
loop_plan : str
Characters which defines order of steps used by the fitting process
e.g. 'fwfwc':
- 'f': fit once without testing for convergence.
- 'c': repeated fit, iterating until convergence is obtained.
- 'w': re-calculate weighting factors on basis of previous fit.
Tol : float
Criterion for the convergence test.
Recommend using 1.0E-5 for reflection or max(abs(S))*1.0E-5
for transmission.
quiet : bool
Boolean flag controlling output of information to the console.
Returns
-------
list
list of fitted parameters: [m1, m2, m3, m4, m5, m5 * Flwst / m6, m7 / Flwst]
weighting_ratio : float
number_iterations : int
RMS_Error : float
References
----------
.. [MAT58] MAT 58, section 2.3, eqn. (26).
"""
N2 = N * 2
iterations = 0
PV = np.ones(N) # default weights vector
PV2 = np.ones(N2)
m1 = self._a.imag / self.Q_L # a''/QL
m2 = -self._a.real / self.Q_L
m3 = self._b.real - m1
m4 = self._b.imag - m2
m5 = self.Q_L
Flwst = self.f[0] # lowest freq. is a convenient normalisation factor.
m6 = Flwst * self.Q_L / self.f_L
m7 = 0.0
last_op = "n"
weighting_ratio = None
number_iterations = 0
## Loop through all of the operations specified in loop_plan
for op in self.loop_plan:
if op == "w":
PV = self.angular_weights(self.f, Flwst * m5 / m6, m5)
weighting_ratio = max(PV) / min(PV)
PV2 = np.concatenate((PV, PV))
if self.verbose:
print("Op w, Calculate weights")
last_op = "n"
continue
if op == "c":
seek_convergence = True
elif op == "f":
seek_convergence = False
else:
assert 0, "Unexpected character in loop_plan"
TerminationConditionMet = False
RMS_Error = None
while not (TerminationConditionMet):
number_iterations += 1
M = np.zeros([N2, 7])
G = np.zeros(N2)[:, np.newaxis]
c1 = complex(-m4, m3)
c2 = complex(m1, m2)
c3 = complex(m3, m4)
for i in range(N):
i2 = i + N
y = 1.0 / complex(1.0, 2 * (m6 * self.f[i] / Flwst - m5))
fdn = self.f[i] / Flwst - m5 / m6
pj = complex(0.0, m7 * fdn)
expm7 = np.exp(pj)
ym = y * expm7
u = c1 * y * ym * 2
u2 = -u * self.f[i] / Flwst
v = (c2 + y * c3) * expm7
u3 = v * fdn
M[i, :] = np.array(
[expm7.real,
-expm7.imag,
ym.real,
-ym.imag,
u.real,
u2.real,
-u3.imag]
)
M[i2, :] = np.array(
[expm7.imag,
expm7.real,
ym.imag,
ym.real,
u.imag,
u2.imag,
u3.real]
)
r = self.s[i] - v # residual
G[i] = r.real
G[i2] = r.imag
X = M.transpose()
T = np.multiply(X, PV2)
C = np.dot(T, M)
q = np.dot(T, G)
dm = np.linalg.solve(C, q)[:,0]
m1 += dm[0]
m2 += dm[1]
m3 += dm[2]
m4 += dm[3]
m5 += dm[4]
m6 += dm[5]
m7 += dm[6]
del G, X, T, C, dm
iterations = iterations + 1
if RMS_Error is not None:
Last_RMS_Error = RMS_Error
else:
Last_RMS_Error = None
SumNum = 0.0
SumDen = 0.0
for i in range(N):
fdn = self.f[i] / Flwst - m5 / m6
den = complex(1.0, 2 * (m6 * self.f[i] / Flwst - m5))
pj = complex(0.0, m7 * fdn)
E = self.s[i] - (c2 + c3 / den) * np.exp(pj)
ip = PV[i]
SumNum = SumNum + ip * (E.real * E.real + E.imag * E.imag)
SumDen = SumDen + ip
RMS_Error = np.sqrt(SumNum / SumDen)
if self.verbose:
if last_op == "c":
print(f"Iteration {iterations}, RMS Error: {RMS_Error}")
else:
print(f"op {op}, Iteration {iterations}, RMS Error: {RMS_Error}")
last_op = op
if seek_convergence:
if Last_RMS_Error is not None:
delta_S = abs(RMS_Error - Last_RMS_Error)
TerminationConditionMet = delta_S < self.tol
else:
TerminationConditionMet = True
# After last operation, we end up here ...
if self.verbose:
print("Optimization done.")
return OptimizedResult({
'success': TerminationConditionMet,
'm1': m1, 'm2': m2, 'm3': m3, 'm4': m4,
'Q_L': m5,
'f_L': m5 * Flwst / m6,
'm7a' : m7 / Flwst,
'weighting_ratio': weighting_ratio,
'number_iterations': number_iterations,
'RMS_Error': RMS_Error,
'method': self.method,
})
def _optimise_fit8(self, N):
"""Iterative non-linear fit, NLQFIT8 Step (2).
Optimised fit of Q-factor (QL) and resonant frequency (FL)
by the gradient-descent method.
Uses the results of the initial fit (sv) as the starting
values for the iteration.
Parameters
----------
F : np.ndarray
Frequency points
S : np.ndarray
Complex data S-parameter to be fitted.
N : int
Number of points.
Fseed : float
Estimated resonant frequency.
sv : list
Initial solution (numpy vector or a list) found with initial_fit.
loop_plan : str
Characters which defines order of steps used by the fitting process
e.g. 'fwfwc':
'f' - fit once without testing for convergence.
'c' - repeated fit, iterating until convergence is obtained.
'w' - re-calculate weighting factors on basis of previous fit.
Tol : float
Criterion for the convergence test.
Recommend using 1.0E-5 for reflection or max(abs(S))*1.0E-5
for transmission.
quiet : bool
Boolean flag controlling output of information to the console.
Returns
-------
list
list of fitted parameters: [m1, m2, m3, m4, m5, m5 * Flwst / m6, m7 / Flwst]
weighting_ratio : float
number_iterations : int
RMS_Error : float
References
----------
.. [MAT58] MAT 58, sec 4.5, eqn. (43).
"""
N2 = N * 2
iterations = 0
PV = np.ones(N) # default weights vector
PV2 = np.ones(N2)
m1 = self._a.imag / self.Q_L # a''/QL
m2 = -self._a.real / self.Q_L
m3 = self._b.real - m1
m4 = self._b.imag - m2
m5 = self.Q_L
Flwst = self.f[0] # lowest freq. is a convenient normalisation factor.
m6 = Flwst * self.Q_L / self.f_L
m8 = 0.0
m9 = 0.0
last_op = "n"
weighting_ratio = None
number_iterations = 0
## Loop through all of the operations specified in loop_plan
for op in self.loop_plan:
if op == "w": # Fr QL
PV = self.angular_weights(self.f, Flwst * float(m5) / float(m6), float(m5))
weighting_ratio = max(PV) / min(PV)
PV2 = np.concatenate((PV, PV))
if self.verbose:
print("Op w, Calculate weights")
last_op = "n"
continue
if op == "c":
seek_convergence = True
elif op == "f":
seek_convergence = False
else:
assert 0, "Unexpected character in loop_plan"
TerminationConditionMet = False
RMS_Error = None
while not (TerminationConditionMet):
number_iterations += 1
M = np.zeros([N2, 8])
G = np.zeros(N2)[:, np.newaxis]
c1 = complex(-m4, m3)
c2 = complex(m1, m2)
c3 = complex(m3, m4)
for i in range(N):
i2 = i + N
y = 1.0 / complex(1.0, 2 * (m6 * self.f[i] / Flwst - m5))
u = c1 * y * y * 2
u2 = -u * self.f[i] / Flwst
FL = Flwst * m5 / m6
t = 2 * (self.f[i] - FL) / FL
M[i, :] = np.array([1.0, 0.0, y.real, -y.imag, u.real, u2.real, t, 0.0])
M[i2, :] = np.array([0.0, 1.0, y.imag, y.real, u.imag, u2.imag, 0.0, t])
v = c2 + c3 * y + (m8 + 1j * m9) * t
r = self.s[i] - v # residual
G[i] = r.real
G[i2] = r.imag
X = M.transpose()
T = np.multiply(X, PV2)
C = np.dot(T, M)
q = np.dot(T, G)
dm = np.linalg.solve(C, q)[:,0]
m1 += dm[0]
m2 += dm[1]
m3 += dm[2]
m4 += dm[3]
m5 += dm[4]
m6 += dm[5]
m8 += dm[6]
m9 += dm[7]
del G, X, T, C, dm
iterations = iterations + 1
if RMS_Error is not None:
Last_RMS_Error = RMS_Error
else:
Last_RMS_Error = None
SumNum = 0.0
SumDen = 0.0
for i in range(N):
den = complex(1.0, 2 * (m6 * self.f[i] / Flwst - m5))
FL = Flwst * m5 / m6
t = 2 * (self.f[i] - FL) / FL
ip = PV[i]
E = (
self.s[i]
- complex(m1, m2)
- complex(m8, m9) * t
- complex(m3, m4) / den
)
SumNum = SumNum + ip * (E.real * E.real + E.imag * E.imag)
SumDen = SumDen + ip
RMS_Error = np.sqrt(SumNum / SumDen)
if self.verbose:
if last_op == "c":
print(f"Iteration {iterations}, RMS Error: {RMS_Error}")
else:
print(f"{op}, Iteration {iterations}, RMS Error: {RMS_Error}")
last_op = op
if seek_convergence:
if Last_RMS_Error is not None:
delta_S = abs(RMS_Error - Last_RMS_Error)
TerminationConditionMet = delta_S < self.tol
else:
TerminationConditionMet = True
# After last operation, we end up here ...
if self.verbose:
print("Optimization done.")
return OptimizedResult({
'success': TerminationConditionMet,
'm1': m1, 'm2': m2, 'm3': m3, 'm4': m4,
'Q_L': m5,
'f_L': m5 * Flwst / m6,
'weighting_ratio': weighting_ratio,
'number_iterations': number_iterations,
'RMS_Error': RMS_Error,
'method': self.method,
})
[docs]
def Q_circle(self,
opt_res: None | OptimizedResult = None,
A: None | float = None
) -> list:
r"""Q-circle diameter.
The diameter of the Q-circle (as displayed in a VNA) provides a visual
indication of whether the coupling is strong or weak [MAT58]_.
Parameters
----------
opt_res : None or :class:`~skrf.qfactor.OptimizedResult`. Default is None.
Solution produced by the :meth:`~skrf.qfactor.Qfactor.fit` method.
If None, uses the solution previously calculated, if performed.
A : None of float. Default is None.
Scaling factor as defined in MAT 58 [MAT58]_.
For `reflection` resonance type, can be set as None
to use the magnitude of the fitted detuned reflection coefficient S_V
Returns
-------
diam : float
Q-circle diameter d.
S_V : complex
Off-resonance Reflection Coefficient.
S_T : complex
Tuned Reflection Coefficient.
References
----------
.. [MAT58] "Q-factor Measurement by using a Vector Network Analyser",
A. P. Gregory, National Physical Laboratory Report MAT 58 (2021)
https://eprintspublications.npl.co.uk/9304/,
section 2.5, eqn. (31).
"""
# if no solution passed, use internal solution if exist
if opt_res is None:
if self.opt_res is None:
raise ValueError('No solution found or passed.')
else:
opt_res = self.opt_res
# m1 : real part of cal_gamma_V
# m2 : imag part of cal_gamma_V
# m3 : real part of b + j a/Q_L
# m4 : imag part of b + j a/Q_L
m1, m2, m3, m4 = (opt_res[key] for key in ['m1', 'm2', 'm3', 'm4'])
if A is None:
A = 1.0 / abs(complex(m1, m2)) # scale to S_V
elif not isinstance(A, (int, float)):
raise ValueError("A should be a float or None")
aqratio = complex(m1, m2)
b = complex(m1 + m3, m2 + m4)
diam = abs(b - aqratio) * A
S_V = complex(m1, m2) * A
S_T = b * A
return diam, S_V, S_T
[docs]
def Q_unloaded(self,
opt_res: None| OptimizedResult = None,
A: None | float = None
) -> float:
"""Unloaded Q-factor Q0.
The value of the unloaded Q-factor Q0 cannot be measured directly but
can be estimated from the measurement of the loaded Q-factor Q_L [MAT58]_ .
Parameters
----------
opt_res : None or :class:`~skrf.qfactor.OptimizedResult`. Default is None.
Solution produced by the :meth:`~skrf.qfactor.Qfactor.fit` method.
If None, uses the solution previously calculated, if performed.
A : float or None. Default is None.
Scaling factor as defined in MAT 58 [MAT58]_.
For `reflection` resonance type, can be set as None
to use the magnitude of the fitted detuned reflection coefficient S_V
Returns
-------
Q0 : float
Unloaded Q-factor.
References
----------
.. [MAT58] "Q-factor Measurement by using a Vector Network Analyser",
A. P. Gregory, National Physical Laboratory Report MAT 58 (2021)
https://eprintspublications.npl.co.uk/9304/
"""
# if no solution passed, use internal solution if exist
if opt_res is None:
if self.opt_res is None:
raise ValueError('No solution found or passed.')
else:
opt_res = self.opt_res
if A is None:
auto_flag = True
elif isinstance(A, (int, float)):
auto_flag = False
else:
raise ValueError("Illegal Scaling factor; should be a float or None")
m1, m2, m3, m4, m5 = (opt_res[key] for key in ['m1', 'm2', 'm3', 'm4', 'Q_L'])
if self.res_type == "transmission":
if auto_flag:
raise ValueError('Scaling factor must be defined for transmission case')
cal_diam, cal_gamma_V, cal_gamma_T = self.Q_circle(opt_res, A)
if cal_diam == 1.0:
raise ZeroDivisionError("Divide by zero forestalled in calculation of Q0")
Q0 = m5 / (1.0 - cal_diam)
elif self.res_type == "reflection":
if auto_flag:
if self.verbose:
print('A is undefined: using fitted data to estimate it')
A = 1.0 / abs(complex(m1, m2)) # scale to S_V if A not defined
cal_diam, S_V, S_T = self.Q_circle(opt_res, A)
cal_touching_circle_diam = 2.0
if self.verbose:
print(f"Q-circle diam = {cal_diam}, touching_circle_diam = {cal_touching_circle_diam}")
den = cal_touching_circle_diam / cal_diam - 1.0
Q0 = m5 * (1.0 + 1.0 / den)
elif self.res_type == "reflection_method2":
if auto_flag:
raise ValueError('Scaling factor must be defined for Method 2')
cal_diam, S_V, S_T = self.Q_circle(opt_res, A)
gv = abs(S_V)
gv2 = gv * gv
mb = abs(S_T)
cosphi = (gv2 + cal_diam * cal_diam - mb * mb) / (
2.0 * gv * cal_diam
) # Cosine rule
cal_touching_circle_diam = (1.0 - gv2) / (1.0 - gv * cosphi)
if self.verbose:
print(f"Q-circle diam = {cal_diam}, touching_circle_diam = {cal_touching_circle_diam}")
den = cal_touching_circle_diam / cal_diam - 1.0
Q0 = m5 * (1.0 + 1.0 / den)
elif self.res_type == "notch" or self.res_type == "absorption": # By transmission
if auto_flag:
if self.verbose:
print(
'Notch/absorption Qo calculation: using fitted data to estimate scaling factor'
)
# scale to S_V if A is undefined
A = 1.0 / abs(complex(m1, m2))
cal_diam, S_V, S_T = self.Q_circle(opt_res, A)
if self.verbose:
print(f"Q-circle diam = {cal_diam}")
if cal_diam == 1.0:
raise ZeroDivisionError("Divide by zero forestalled in calculation of Qo")
den = 1.0 / cal_diam - 1.0 # Gao thesis (2008) 4.35 and 4.40
Q0 = m5 * (
1.0 + 1.0 / den
) # https://resolver.caltech.edu/CaltechETD:etd-06092008-235549
# For this type of resonator, critical coupling occurs for cal_diam = 0.5.
else:
raise ValueError("Unknown resonance type {self.res_type}")
return Q0
[docs]
def fitted_s(self,
opt_res: None | OptimizedResult = None,
f: None | np.ndarray = None
) -> np.ndarray:
r"""S-parameter response of an equivalent circuit model resonator.
The approximate solution and estimated is obtained from [MAT58]_:
.. math::
S = m_1 + j m_2 + \frac{m_3 + j m_4}{1 + j Q_L t}
where the m coefficients come from the fitted solution and
.. math::
t = \frac{f}{f_L} - \frac{f_L}{f} \approx 2 \frac{f − f_L}{f_L}
Parameters
----------
opt_res : None or :class:`~skrf.qfactor.OptimizedResult`. Default is None.
Solution produced by the :meth:`~skrf.qfactor.Qfactor.fit` method.
If None, uses the solution previously calculated, if performed.
f : None or np.ndarray. Default is None.
frequency array [Hz]. If None, use the self frequencies.
Returns
-------
s : np.ndarray
S-parameter response.
References
----------
.. [MAT58] "Q-factor Measurement by using a Vector Network Analyser",
A. P. Gregory, National Physical Laboratory Report MAT 58 (2021)
https://eprintspublications.npl.co.uk/9304/,
section 2.2, eqn. (21).
"""
# if no solution passed, use internal solution if exist
if opt_res is None:
if self.opt_res is None:
raise ValueError('No solution found or passed.')
else:
opt_res = self.opt_res
if f is None:
f = self.f
# fractional offset frequency
t = f/opt_res.f_L - opt_res.f_L/f
y = 1/(1 + 1j*opt_res.Q_L*t)
s = opt_res.m1 +1j*opt_res.m2 + (opt_res.m3 + 1j*opt_res.m4) * y
return s
[docs]
def fitted_network(self,
opt_res: None | OptimizedResult = None,
frequency: None | Frequency = None,
) -> Network:
"""Fitted Network.
Return the Network corresponding to the fitted response.
Parameters
----------
opt_res : None or :class:`~skrf.qfactor.OptimizedResult`. Default is None.
Solution produced by the :meth:`~skrf.qfactor.Qfactor.fit` method.
If None, uses the solution previously calculated, if performed.
frequency : None or :class:`~skrf.frequency.Frequency`. Default is None.
Frequency for the fitted Network. If None, use the same
Frequency than the one used to create the QFactor.
Returns
-------
ntwk : :class:`~skrf.network.Network`
Fitted Network for the passed Frequency.
"""
# if no solution passed, use internal solution if exist
if opt_res is None:
if self.opt_res is None:
raise ValueError('No solution found or passed.')
else:
opt_res = self.opt_res
if frequency is None:
frequency = self._ntwk.frequency
s = self.fitted_s(opt_res, f=frequency.f)
ntwk = Network(s=s, frequency=frequency)
return ntwk
@property
def f_L_scaled(self) -> float:
"""
Resonant Frequency in the frequency unit..
Returns
-------
float
Resonant frequency in the frequency unit.
See Also
--------
f_L : Resonant Frequency in Hz.
"""
if not self.fitted:
warn('Q-factor not fitted, result may be inaccurate. Use the .fit() method before.', stacklevel=2)
return self.f_L/self.f_multiplier
@property
def BW(self) -> float:
r"""3-dB Bandwidth.
Return the half-power fractional bandwidth (aka 3-dB bandwidth)
defined as:
.. math::
BW = \frac{f_L}{Q_L}
Returns
-------
float
3-dB Bandwidth in Hz.
See Also
--------
BW_scaled : 3-dB Bandwidth scaled to the frequency unit.
"""
if not self.fitted:
warn('Q-factor not fitted, result may be inaccurate. Use the .fit() method before.', stacklevel=2)
return self.f_L/self.Q_L
@property
def BW_scaled(self) -> float:
r"""3-dB Bandwidth scaled to the frequency unit.
Returns
-------
float
3-dB Bandwidth in the frequency unit.
See Also
--------
BW : 3-dB Bandwidth in Hz.
"""
if not self.fitted:
warn('Q-factor not fitted, result may be inaccurate. Use the .fit() method before.', stacklevel=2)
return self.BW/self.f_multiplier