"""
analysis.py
-----------
Pure numerical analysis functions for extracting resonance frequencies
and linewidths from S-parameter networks.
None of these functions build networks or do optimization — they only
inspect existing rf.Network objects. This makes them easy to test
independently and reuse across models.
"""
from __future__ import annotations
import logging
import numpy as np
import skrf as rf
from scipy.signal import savgol_filter, find_peaks
logger = logging.getLogger(__name__)
def _fit_circle_kasa(z: np.ndarray) -> tuple[complex, float]:
"""
Algebraic Kasa circle fit to complex points z = x + i y.
Returns
-------
center : complex
Fitted circle center.
radius : float
Fitted circle radius.
"""
x = np.real(z)
y = np.imag(z)
A = np.column_stack([x, y, np.ones_like(x)])
b = -(x**2 + y**2)
coeffs, *_ = np.linalg.lstsq(A, b, rcond=None)
a, b_, c = coeffs
xc = -a / 2.0
yc = -b_ / 2.0
r = np.sqrt(max(xc**2 + yc**2 - c, 0.0))
return xc + 1j * yc, r
[docs]
def circle_fit_f0_kappa(
ntwk: rf.Network,
m: int = 0,
n: int = 0,
smooth_window: int | None = 51,
smooth_polyorder: int = 3,
) -> tuple[float, float]:
"""
Extract resonance frequency (f0) and linewidth (kappa) from a circle fit
to S[m,n] in the complex plane (reflection S11/S22).
Returns
-------
f0 : float
Resonance frequency in Hz.
kappa : float
Linewidth (FWHM of the angular velocity peak) in Hz.
"""
f = ntwk.frequency.f
s = ntwk.s[:, m, n]
if len(f) < 5:
raise ValueError("Need ≥5 frequency points for circle-fit.")
center, radius = _fit_circle_kasa(s)
if not np.isfinite(radius) or radius <= 0:
raise ValueError("Circle fit failed: non-physical radius.")
theta = np.unwrap(np.angle(s - center))
sw = smooth_window
if sw is not None:
if sw >= len(theta):
sw = len(theta) - 1 if len(theta) % 2 == 0 else len(theta)
if sw < 5:
sw = None
elif sw % 2 == 0:
sw += 1
if sw is not None:
theta = savgol_filter(theta, window_length=sw, polyorder=smooth_polyorder)
dtheta_df = np.abs(np.gradient(theta, f))
idx_peak = int(np.argmax(dtheta_df))
f0 = float(f[idx_peak])
half_max = dtheta_df[idx_peak] / 2.0
above = dtheta_df > half_max
crossings = np.where(np.diff(above.astype(int)))[0]
left_cross = crossings[crossings < idx_peak]
right_cross = crossings[crossings > idx_peak]
if len(left_cross) == 0 or len(right_cross) == 0:
if idx_peak > 0 and idx_peak < len(f) - 1:
df = f[1] - f[0]
kappa = 2.0 * df
else:
kappa = float("nan")
else:
iL = left_cross[-1]
iR = right_cross[0]
def _interp(i: int) -> float:
f0_, f1_ = f[i], f[i + 1]
y0_, y1_ = dtheta_df[i], dtheta_df[i + 1]
if y1_ == y0_:
return 0.5 * (f0_ + f1_)
return f0_ + (half_max - y0_) * (f1_ - f0_) / (y1_ - y0_)
kappa = float(_interp(iR) - _interp(iL))
return f0, kappa
[docs]
def circle_fit_modes(
ntwk: rf.Network,
m: int = 0,
n: int = 0,
n_modes: int = 3,
prominence: float = None,
min_spacing_hz: float = None,
smooth_window: int | None = 51,
smooth_polyorder: int = 3,
) -> tuple[np.ndarray, np.ndarray]:
"""
Find multiple resonance modes in S[m,n] using the reflection circle fit.
Returns
-------
f0s : np.ndarray
Mode frequencies in Hz, length n_modes (NaN-padded).
kappas : np.ndarray
Mode linewidths in Hz, length n_modes (NaN-padded).
"""
f = ntwk.frequency.f
s = ntwk.s[:, m, n]
df = f[1] - f[0]
if min_spacing_hz is None:
min_spacing_hz = (f[-1] - f[0]) / (10 * n_modes)
min_distance_pts = max(1, int(min_spacing_hz / df))
center, _ = _fit_circle_kasa(s)
theta = np.unwrap(np.angle(s - center))
sw = smooth_window
if sw is not None:
if sw >= len(theta):
sw = len(theta) - 1 if len(theta) % 2 == 0 else len(theta)
if sw < 5:
sw = None
elif sw % 2 == 0:
sw += 1
if sw is not None:
theta = savgol_filter(theta, window_length=sw, polyorder=smooth_polyorder)
trace = np.abs(np.gradient(theta, f))
if prominence is None:
prominence = 0.0005 * (trace.max() - trace.min())
peak_indices, _ = find_peaks(
trace,
prominence=prominence,
distance=min_distance_pts,
)
if len(peak_indices) > n_modes:
heights = trace[peak_indices]
top = np.argsort(heights)[::-1][:n_modes]
peak_indices = np.sort(peak_indices[top])
f0s = np.full(n_modes, np.nan)
kappas = np.full(n_modes, np.nan)
for k, pidx in enumerate(peak_indices):
if k >= n_modes:
break
half_w_pts = max(20, int(0.5 * min_distance_pts))
i_lo = max(0, pidx - 5 * half_w_pts)
i_hi = min(len(f) - 1, pidx + 5 * half_w_pts)
local_freq = rf.Frequency.from_f(f[i_lo:i_hi + 1], unit="Hz")
local_s = s[i_lo:i_hi + 1]
if len(local_s) < 5:
f0s[k] = float(f[pidx])
kappas[k] = float("nan")
continue
local_ntwk = rf.Network(
frequency=local_freq,
s=local_s[:, np.newaxis, np.newaxis],
)
try:
f0_loc, kappa_loc = circle_fit_f0_kappa(
local_ntwk,
m=0,
n=0,
smooth_window=smooth_window,
smooth_polyorder=smooth_polyorder,
)
f0s[k] = f0_loc
kappas[k] = kappa_loc
except Exception:
f0s[k] = float(f[pidx])
kappas[k] = float("nan")
return f0s, kappas
[docs]
def stitch_shifted_freqs(
*arrays: np.ndarray,
dedup_tol_ghz: float = 0.05,
return_sources: bool = False,
):
"""Merge S11/S22 mode frequency lists, deduplicating nearby modes."""
tagged = []
for src_idx, arr in enumerate(arrays):
for freq in arr[~np.isnan(arr)]:
tagged.append((freq, src_idx))
if not tagged:
empty = np.array([])
return (empty, np.array([], dtype=object)) if return_sources else empty
tagged = sorted(tagged, key=lambda x: x[0])
merged_freqs = []
merged_sources = []
cluster_freqs = [tagged[0][0]]
cluster_srcs = [tagged[0][1]]
for freq, src in tagged[1:]:
if freq - cluster_freqs[0] < dedup_tol_ghz:
cluster_freqs.append(freq)
cluster_srcs.append(src)
else:
merged_freqs.append(float(np.mean(cluster_freqs)))
merged_sources.append(max(set(cluster_srcs), key=cluster_srcs.count))
cluster_freqs = [freq]
cluster_srcs = [src]
merged_freqs.append(float(np.mean(cluster_freqs)))
merged_sources.append(max(set(cluster_srcs), key=cluster_srcs.count))
freqs = np.array(merged_freqs)
sources = np.array(["s11" if s == 0 else "s22" for s in merged_sources], dtype=object)
return (freqs, sources) if return_sources else freqs
[docs]
def resonance_circle_fit(
ntwk: rf.Network,
m: int = 0,
n: int = 0,
smooth_window: int | None = 51,
smooth_polyorder: int = 3,
) -> float:
"""
Estimate resonance frequency from a circle fit to S[m,n] in the complex plane.
Strategy
--------
1. Fit a circle to the complex trace S[m,n](f).
2. Compute the angle of each point about the fitted center.
3. Take the frequency where |d(theta)/df| is largest.
This is usually a good reflection-based estimator for S11/S22.
Parameters
----------
ntwk : rf.Network
m, n : int
S-parameter indices.
smooth_window : int or None
Optional Savitzky-Golay smoothing window for theta before differentiation.
Must be odd if provided.
smooth_polyorder : int
Polynomial order for Savitzky-Golay smoothing.
Returns
-------
float
Resonance frequency in Hz.
"""
f0, _ = circle_fit_f0_kappa(
ntwk,
m=m,
n=n,
smooth_window=smooth_window,
smooth_polyorder=smooth_polyorder,
)
return f0
[docs]
def resonance(
ntwk: rf.Network,
m: int = 0,
n: int = 0,
use_max: bool = False,
method: str = "min_re",
) -> float:
"""
Estimate resonance frequency from S[m,n].
Parameters
----------
ntwk : rf.Network
m, n : int
S-parameter indices (0-based).
use_max : bool
Only used for method='min_re'. If True, find frequency of maximum
Re(S) instead of minimum.
method : str
Resonance estimator. Supported:
- 'min_re' : min/max of Re(S)
- 'circle_fit' : circle fit in complex plane
Returns
-------
float
Resonance frequency in Hz.
"""
if method == "circle_fit":
return resonance_circle_fit(ntwk, m=m, n=n)
if method == "min_re":
f = ntwk.frequency.f
re_s = np.real(ntwk.s[:, m, n])
idx = np.argmax(re_s) if use_max else np.argmin(re_s)
return float(f[idx])
raise ValueError(f"Unknown resonance method: {method}")
# def resonance(ntwk: rf.Network, m: int = 0, n: int = 0, use_max: bool = False) -> float:
# """
# Estimate resonance frequency from the real part of S[m,n].
# By default returns the frequency of the minimum of Re(S), which
# corresponds to the resonance dip in S11/S22. Set use_max=True
# for transmission peaks.
# Parameters
# ----------
# ntwk : rf.Network
# m, n : int
# S-parameter indices (0-based).
# use_max : bool
# If True, find frequency of maximum Re(S) instead of minimum.
# Returns
# -------
# float
# Resonance frequency in Hz.
# """
# f = ntwk.frequency.f
# re_s = np.real(ntwk.s[:, m, n])
# idx = np.argmax(re_s) if use_max else np.argmin(re_s)
# return float(f[idx])
[docs]
def resonances(
ntwk: rf.Network,
m: int = 0,
n: int = 0,
n_modes: int = 3,
method: str = "min_re",
prominence: float = None,
min_spacing_hz: float = None,
) -> np.ndarray:
"""
Estimate resonance frequencies for multiple modes from S[m,n].
Parameters
----------
ntwk : rf.Network
m, n : int
S-parameter indices (0-based).
n_modes : int
Number of resonant modes to find.
method : str
Resonance estimator. Supported:
- 'min_re' : finds n_modes minima of Re(S)
- 'circle_fit' : finds n_modes peaks of |dθ/df| from circle fit
prominence : float, optional
Minimum peak prominence passed to scipy.signal.find_peaks.
If None, defaults to 0.1 * (max - min) of the trace being searched.
min_spacing_hz : float, optional
Minimum frequency separation between modes in Hz.
If None, defaults to (f_max - f_min) / (10 * n_modes).
Returns
-------
np.ndarray
Resonance frequencies in Hz, sorted ascending, length n_modes.
If fewer than n_modes peaks are found, the array is padded with NaN.
"""
from scipy.signal import find_peaks
f = ntwk.frequency.f
df = f[1] - f[0] # assumes uniform spacing
if min_spacing_hz is None:
min_spacing_hz = (f[-1] - f[0]) / (10 * n_modes)
min_distance_pts = max(1, int(min_spacing_hz / df))
if method == "min_re":
# Search for minima of Re(S) by inverting the trace
trace = -np.real(ntwk.s[:, m, n])
elif method == "circle_fit":
# Build the angular velocity |dθ/df| trace from the circle fit
s = ntwk.s[:, m, n]
x, y = s.real, s.imag
# Kasa algebraic circle fit
A = np.column_stack([x, y, np.ones(len(x))])
b = -(x**2 + y**2)
coeffs, _, _, _ = np.linalg.lstsq(A, b, rcond=None)
cx, cy = -coeffs[0] / 2.0, -coeffs[1] / 2.0
theta = np.unwrap(np.angle((x - cx) + 1j * (y - cy)))
trace = np.abs(np.gradient(theta, f))
else:
raise ValueError(f"Unknown resonance method: {method!r}")
if prominence is None:
prominence = 0.1 * (trace.max() - trace.min())
peak_indices, _ = find_peaks(
trace,
prominence=prominence,
distance=min_distance_pts,
)
# Sort by peak height descending, take top n_modes, then re-sort by frequency
if len(peak_indices) > n_modes:
heights = trace[peak_indices]
top = np.argsort(heights)[::-1][:n_modes]
peak_indices = np.sort(peak_indices[top])
freqs = f[peak_indices].astype(float)
# Pad with NaN if not enough peaks found
if len(freqs) < n_modes:
freqs = np.concatenate([freqs, np.full(n_modes - len(freqs), np.nan)])
return freqs
[docs]
def resonance_from_s_max(network: rf.Network, m: int = 0, n: int = 0) -> float:
"""
Resonance frequency from the dominant peak in |S[m,n]|.
Finds all peaks in the magnitude and returns the frequency of the
largest one. Useful for transmission parameters (S12/S21) where the
resonance appears as a peak rather than a dip.
Parameters
----------
network : rf.Network
m, n : int
S-parameter indices (0-based).
Returns
-------
float
Dominant resonance frequency in GHz.
Raises
------
ValueError
If no peaks are found in |S[m,n]|.
"""
f = network.frequency.f / 1e9
s_mn = network.s[:, m, n]
peaks, _ = find_peaks(np.abs(s_mn))
if len(peaks) == 0:
raise ValueError("No peaks found in |S[{},{}]|.".format(m, n))
peak_idx = peaks[np.argmax(np.abs(s_mn)[peaks])]
return float(f[peak_idx])
[docs]
def resonances_from_s_max(network: rf.Network, m: int = 0, n: int = 0) -> float:
"""
Resonance frequency from the dominant peak in |S[m,n]|.
Finds all peaks in the magnitude and returns the frequency of the
largest one. Useful for transmission parameters (S12/S21) where the
resonance appears as a peak rather than a dip.
Parameters
----------
network : rf.Network
m, n : int
S-parameter indices (0-based).
Returns
-------
float
Dominant resonance frequency in GHz.
Raises
------
ValueError
If no peaks are found in |S[m,n]|.
"""
f = network.frequency.f / 1e9
s_mn = network.s[:, m, n]
peaks, _ = find_peaks(np.abs(s_mn))
freqs = f[peaks]
logger.debug("Resonance frequencies (GHz): %s", freqs)
if len(peaks) == 0:
raise ValueError("No peaks found in |S[{},{}]|.".format(m, n))
return freqs
[docs]
def resonances_from_s(network: rf.Network, m: int = 0, n: int = 0) -> np.ndarray:
"""
Find all resonance frequencies from peaks in Re(S[m,n]).
Parameters
----------
network : rf.Network
m, n : int
S-parameter indices (0-based).
Returns
-------
np.ndarray
Array of resonance frequencies in GHz.
"""
f = network.frequency.f / 1e9
s_mn = network.s[:, m, n]
peaks, _ = find_peaks(np.real(s_mn))
freqs = f[peaks]
logger.debug("Resonance frequencies (GHz): %s", freqs)
return freqs
[docs]
def fwhm_from_trace_db(
ntwk: rf.Network,
m: int = 0,
n: int = 0,
kind: str = "dip",
smooth: int = None,
) -> float:
"""
Full-Width at Half Maximum (FWHM) from S-parameter magnitude in dB.
Uses crossing-interpolation to find the two frequencies where the
magnitude crosses the half-depth level, giving a more accurate result
than a pure index-based approach.
Parameters
----------
ntwk : rf.Network
m, n : int
S-parameter indices (0-based).
kind : {"dip", "peak"}
"dip" — resonance appears as a downward dip (typical for S11/S22).
"peak" — resonance appears as an upward peak (typical for S12/S21).
smooth : int or None
Optional moving-average window length for noisy traces.
Use an odd integer; None disables smoothing.
Returns
-------
float
FWHM linewidth in Hz.
Raises
------
ValueError
If two crossings at the half-depth level cannot be found.
"""
f = np.asarray(ntwk.frequency.f, dtype=float)
s = ntwk.s[:, m, n]
mag_db = 20 * np.log10(np.maximum(np.abs(s), 1e-300))
if smooth is not None and smooth > 1:
kernel = np.ones(int(smooth)) / int(smooth)
mag_db = np.convolve(mag_db, kernel, mode="same")
if kind == "dip":
i0 = np.argmin(mag_db)
y0 = mag_db[i0]
baseline = np.median(np.r_[
mag_db[:max(3, i0 // 8)],
mag_db[min(len(mag_db) - 3, i0 + len(mag_db) // 8):],
])
level = 0.5 * (baseline + y0)
above = mag_db > level
elif kind == "peak":
i0 = np.argmax(mag_db)
y0 = mag_db[i0]
baseline = np.median(np.r_[
mag_db[:max(3, i0 // 8)],
mag_db[min(len(mag_db) - 3, i0 + len(mag_db) // 8):],
])
level = 0.5 * (baseline + y0)
above = mag_db < level
else:
raise ValueError("kind must be 'dip' or 'peak', got '{}'.".format(kind))
flips = np.where(above[:-1] != above[1:])[0]
if flips.size < 2:
raise ValueError("Could not find two crossings for FWHM.")
left_flips = flips[flips < i0]
right_flips = flips[flips >= i0]
if left_flips.size == 0 or right_flips.size == 0:
raise ValueError("Could not bracket resonance with crossings.")
iL = left_flips[-1]
iR = right_flips[0]
def _interp_crossing(i: int) -> float:
x0, x1 = f[i], f[i + 1]
y0_, y1_ = mag_db[i], mag_db[i + 1]
if y1_ == y0_:
return 0.5 * (x0 + x1)
return x0 + (level - y0_) * (x1 - x0) / (y1_ - y0_)
return float(_interp_crossing(iR) - _interp_crossing(iL))
[docs]
def fwhm_from_res11(ntwk: rf.Network) -> float:
"""
Linewidth from zero crossings of Re(S11).
Re(S11) crosses zero at the two half-power frequencies of the resonance.
Uses linear interpolation between samples for accuracy.
Parameters
----------
ntwk : rf.Network
Returns
-------
float
Linewidth (distance between the two zero crossings) in Hz.
Raises
------
ValueError
If fewer than two zero crossings are found.
"""
f = ntwk.frequency.f
re_s11 = np.real(ntwk.s[:, 0, 0])
sign_changes = np.where(np.diff(np.sign(re_s11)))[0]
if len(sign_changes) < 2:
raise ValueError("Could not determine linewidth: fewer than 2 zero crossings found.")
zeros = []
for idx in sign_changes[:2]:
f1, f2 = f[idx], f[idx + 1]
y1, y2 = re_s11[idx], re_s11[idx + 1]
f_zero = f1 - y1 * (f2 - f1) / (y2 - y1)
zeros.append(f_zero)
return float(zeros[1] - zeros[0])
[docs]
def resonance_from_res11(ntwk: rf.Network) -> float:
"""
Resonance frequency from the minimum of Re(S11).
Thin convenience wrapper around `resonance` for the common single-port
S11 case.
Parameters
----------
ntwk : rf.Network
Returns
-------
float
Resonance frequency in Hz.
"""
return resonance(ntwk, m=0, n=0, use_max=False)
def find_resonant_frequency(network: rf.Network) -> complex:
freqs = network.f
S21 = network.s[:, 1, 0]
mag = np.abs(S21)
peak_idx = np.argmax(mag)
omega_0 = 2 * np.pi * freqs[peak_idx]
half_pow = mag[peak_idx] / np.sqrt(2)
left = np.where(mag[:peak_idx] < half_pow)[0]
right = np.where(mag[peak_idx:] < half_pow)[0]
alpha = (2 * np.pi * freqs[peak_idx + right[0]] - 2 * np.pi * freqs[left[-1]]) / 2
pole = complex(-alpha, omega_0)
f0 = pole.imag / (2 * np.pi)
return f0
[docs]
def find_resonant_frequency(network: rf.Network) -> complex:
freqs = network.f
S21 = network.s[:, 1, 0]
mag = np.abs(S21)
peak_idx = np.argmax(mag)
omega_0 = 2 * np.pi * freqs[peak_idx]
half_pow = mag[peak_idx] / np.sqrt(2)
left = np.where(mag[:peak_idx] < half_pow)[0]
right = np.where(mag[peak_idx:] < half_pow)[0]
alpha = (2 * np.pi * freqs[peak_idx + right[0]] - 2 * np.pi * freqs[left[-1]]) / 2
pole = complex(-alpha, omega_0)
f0 = pole.imag / (2 * np.pi)
return f0
