Frequency-Domain Analysis & Stability Margins

freq_domain

Frequency-domain analysis and stability margin calculations for LTI systems.

BodeData dataclass

Frequency response data for Bode diagram analysis.

Parameters:

Name	Type	Description	Default
w	NDArray[float64]	Frequencies in radians per second (rad/s).	required
mag	NDArray[float64]	Linear magnitude |G(jw)|.	required
phase	NDArray[float64]	Unwrapped phase in degrees.	required
mag_db	NDArray[float64]	Magnitude in decibels: 20 * log10(|G(jw)|).	required
phase_rad	NDArray[float64]	Unwrapped phase in radians.	required
response	NDArray[complex128]	Complex frequency response values G(jw).	required

Source code in src/ctrlpy/freq_domain.py

@dataclass
class BodeData:
    """Frequency response data for Bode diagram analysis.

    Parameters
    ----------
    w : NDArray[np.float64]
        Frequencies in radians per second (rad/s).
    mag : NDArray[np.float64]
        Linear magnitude |G(jw)|.
    phase : NDArray[np.float64]
        Unwrapped phase in degrees.
    mag_db : NDArray[np.float64]
        Magnitude in decibels: 20 * log10(|G(jw)|).
    phase_rad : NDArray[np.float64]
        Unwrapped phase in radians.
    response : NDArray[np.complex128]
        Complex frequency response values G(jw).
    """

    w: NDArray[np.float64]
    mag: NDArray[np.float64]
    phase: NDArray[np.float64]
    mag_db: NDArray[np.float64]
    phase_rad: NDArray[np.float64]
    response: NDArray[np.complex128]

    def __iter__(self) -> Iterator[NDArray[np.float64]]:
        """Enable unpacking: w, mag, phase = bode_data(sys).

        Yields
        ------
        NDArray[np.float64]
            Frequency array `w`, linear magnitude array `mag`, and phase in degrees `phase`.
        """
        yield self.w
        yield self.mag
        yield self.phase

    def _repr_latex_(self) -> str:
        """Return a LaTeX representation for Jupyter environments."""
        return (
            rf"$$\text{{BodeData: }} {len(self.w)}\text{{ frequency points from }} "
            rf"\omega = {self.w[0]:.3g}\text{{ to }} {self.w[-1]:.3g}\text{{ rad/s}}$$"
        )

    def _repr_markdown_(self) -> str:
        """Return a Markdown representation for Jupyter environments."""
        return self._repr_latex_()

__iter__

__iter__()

Enable unpacking: w, mag, phase = bode_data(sys).

Yields:

Type	Description
NDArray[float64]	Frequency array w, linear magnitude array mag, and phase in degrees phase.

Source code in src/ctrlpy/freq_domain.py

def __iter__(self) -> Iterator[NDArray[np.float64]]:
    """Enable unpacking: w, mag, phase = bode_data(sys).

    Yields
    ------
    NDArray[np.float64]
        Frequency array `w`, linear magnitude array `mag`, and phase in degrees `phase`.
    """
    yield self.w
    yield self.mag
    yield self.phase

NyquistData dataclass

Frequency response data for Nyquist diagram analysis.

Parameters:

Name	Type	Description	Default
w	NDArray[float64]	Frequencies in radians per second (rad/s).	required
response	NDArray[complex128]	Complex frequency response values G(jw).	required
real	NDArray[float64]	Real part Re(G(jw)).	required
imag	NDArray[float64]	Imaginary part Im(G(jw)).	required
arc_s	NDArray[complex128] | None	Complex frequencies s along the indented arc around s=0.	None
arc_response	NDArray[complex128] | None	Complex frequency response G(s) along the indented arc around s=0.	None

Source code in src/ctrlpy/freq_domain.py

@dataclass
class NyquistData:
    """Frequency response data for Nyquist diagram analysis.

    Parameters
    ----------
    w : NDArray[np.float64]
        Frequencies in radians per second (rad/s).
    response : NDArray[np.complex128]
        Complex frequency response values G(jw).
    real : NDArray[np.float64]
        Real part Re(G(jw)).
    imag : NDArray[np.float64]
        Imaginary part Im(G(jw)).
    arc_s : NDArray[np.complex128] | None, optional
        Complex frequencies s along the indented arc around s=0.
    arc_response : NDArray[np.complex128] | None, optional
        Complex frequency response G(s) along the indented arc around s=0.
    """

    w: NDArray[np.float64]
    response: NDArray[np.complex128]
    real: NDArray[np.float64]
    imag: NDArray[np.float64]
    arc_s: NDArray[np.complex128] | None = None
    arc_response: NDArray[np.complex128] | None = None

    def __iter__(self) -> Iterator[Any]:
        """Enable unpacking: w, response = nyquist_data(sys).

        Yields
        ------
        NDArray[np.float64] | NDArray[np.complex128]
            Frequency array `w` and complex response array `response`.
        """
        yield self.w
        yield self.response

    def _repr_latex_(self) -> str:
        """Return a LaTeX representation for Jupyter environments."""
        arc_info = (
            rf",\text{{ indented arc: }} {len(self.arc_s)}\text{{ points}}"
            if self.arc_s is not None
            else ""
        )
        return (
            rf"$$\text{{NyquistData: }} {len(self.w)}\text{{ frequency points}}"
            rf"{arc_info}$$"
        )

    def _repr_markdown_(self) -> str:
        """Return a Markdown representation for Jupyter environments."""
        return self._repr_latex_()

__iter__

__iter__()

Enable unpacking: w, response = nyquist_data(sys).

Yields:

Type	Description
NDArray[float64] | NDArray[complex128]	Frequency array w and complex response array response.

Source code in src/ctrlpy/freq_domain.py

def __iter__(self) -> Iterator[Any]:
    """Enable unpacking: w, response = nyquist_data(sys).

    Yields
    ------
    NDArray[np.float64] | NDArray[np.complex128]
        Frequency array `w` and complex response array `response`.
    """
    yield self.w
    yield self.response

RootLocusData dataclass

Root locus pole trajectory data across varying gains.

Parameters:

Name	Type	Description	Default
gains	NDArray[float64]	1D array of non-negative feedback gain values k >= 0.	required
roots	NDArray[complex128]	2D array of closed-loop pole locations of shape (n_gains, n_poles).	required
poles	NDArray[complex128]	Open-loop poles (corresponding to k = 0).	required
zeros	NDArray[complex128]	Open-loop zeros.	required

Source code in src/ctrlpy/freq_domain.py

@dataclass
class RootLocusData:
    """Root locus pole trajectory data across varying gains.

    Parameters
    ----------
    gains : NDArray[np.float64]
        1D array of non-negative feedback gain values k >= 0.
    roots : NDArray[np.complex128]
        2D array of closed-loop pole locations of shape (n_gains, n_poles).
    poles : NDArray[np.complex128]
        Open-loop poles (corresponding to k = 0).
    zeros : NDArray[np.complex128]
        Open-loop zeros.
    """

    gains: NDArray[np.float64]
    roots: NDArray[np.complex128]
    poles: NDArray[np.complex128]
    zeros: NDArray[np.complex128]

    def __iter__(self) -> Iterator[Any]:
        """Enable unpacking: gains, roots = root_locus_data(sys).

        Yields
        ------
        NDArray[np.float64] | NDArray[np.complex128]
            Gains array and roots 2D array.
        """
        yield self.gains
        yield self.roots

    def _repr_latex_(self) -> str:
        """Return a LaTeX representation for Jupyter environments."""
        return (
            rf"$$\text{{RootLocusData: }} {len(self.gains)}\text{{ gains, }} "
            rf"{len(self.poles)}\text{{ open-loop poles, }} "
            rf"{len(self.zeros)}\text{{ open-loop zeros}}$$"
        )

    def _repr_markdown_(self) -> str:
        """Return a Markdown representation for Jupyter environments."""
        return self._repr_latex_()

__iter__

__iter__()

Enable unpacking: gains, roots = root_locus_data(sys).

Yields:

Type	Description
NDArray[float64] | NDArray[complex128]	Gains array and roots 2D array.

Source code in src/ctrlpy/freq_domain.py

def __iter__(self) -> Iterator[Any]:
    """Enable unpacking: gains, roots = root_locus_data(sys).

    Yields
    ------
    NDArray[np.float64] | NDArray[np.complex128]
        Gains array and roots 2D array.
    """
    yield self.gains
    yield self.roots

StabilityMargins dataclass

Container for classical stability margins and crossover frequencies.

Parameters:

Name	Type	Description	Default
gm_db	float	Gain margin in decibels (dB). Inf if no phase crossover occurs.	required
pm_deg	float	Phase margin in degrees. Inf if no gain crossover occurs.	required
wcg	float	Gain crossover frequency in rad/s (where |G(jw)| = 1 / 0 dB). NaN if no crossover.	required
wcp	float	Phase crossover frequency in rad/s (where phase crosses -180 deg). NaN if no crossover.	required

Source code in src/ctrlpy/freq_domain.py

@dataclass
class StabilityMargins:
    """Container for classical stability margins and crossover frequencies.

    Parameters
    ----------
    gm_db : float
        Gain margin in decibels (dB). Inf if no phase crossover occurs.
    pm_deg : float
        Phase margin in degrees. Inf if no gain crossover occurs.
    wcg : float
        Gain crossover frequency in rad/s (where |G(jw)| = 1 / 0 dB). NaN if no crossover.
    wcp : float
        Phase crossover frequency in rad/s (where phase crosses -180 deg). NaN if no crossover.
    """

    gm_db: float
    pm_deg: float
    wcg: float
    wcp: float

    @property
    def gm(self) -> float:
        """Gain margin in dB."""
        return self.gm_db

    @property
    def pm(self) -> float:
        """Phase margin in degrees."""
        return self.pm_deg

    @property
    def Wcg(self) -> float:
        """Gain crossover frequency in rad/s."""
        return self.wcg

    @property
    def Wcp(self) -> float:
        """Phase crossover frequency in rad/s."""
        return self.wcp

    def __iter__(self) -> Iterator[float]:
        """Enable unpacking: gm_db, pm_deg, wcg, wcp = margin(sys).

        Yields
        ------
        float
            gm_db, pm_deg, wcg, wcp.
        """
        yield self.gm_db
        yield self.pm_deg
        yield self.wcg
        yield self.wcp

    def _repr_latex_(self) -> str:
        """Return a LaTeX representation for Jupyter environments."""
        gm_str = (
            f"{self.gm_db:.2f}\\text{{ dB}}" if not np.isinf(self.gm_db) else r"\infty\text{ dB}"
        )
        pm_str = (
            f"{self.pm_deg:.2f}^\\circ"
            if (not np.isinf(self.pm_deg) and not np.isnan(self.pm_deg))
            else (r"\infty" if np.isinf(self.pm_deg) else r"\text{N/A}")
        )
        wcg_str = f"{self.wcg:.4g}\\text{{ rad/s}}" if not np.isnan(self.wcg) else r"\text{N/A}"
        wcp_str = f"{self.wcp:.4g}\\text{{ rad/s}}" if not np.isnan(self.wcp) else r"\text{N/A}"

        lines = [
            r"$$\begin{array}{|l|c|}",
            r"\hline",
            r"\textbf{Stability Margin / Crossover Frequency} & \textbf{Value} \\",
            r"\hline",
            rf"\text{{Gain Margin }} (\mathrm{{GM}}) & {gm_str} \\",
            rf"\text{{Phase Margin }} (\mathrm{{PM}}) & {pm_str} \\",
            rf"\text{{Gain Crossover Frequency }} (\omega_{{cg}}) & {wcg_str} \\",
            rf"\text{{Phase Crossover Frequency }} (\omega_{{cp}}) & {wcp_str} \\",
            r"\hline",
            r"\end{array}$$",
        ]
        return "\n".join(lines)

    def _repr_markdown_(self) -> str:
        """Return a Markdown representation for Jupyter environments."""
        return self._repr_latex_()

gm property

gm

Gain margin in dB.

pm property

pm

Phase margin in degrees.

Wcg property

Wcg

Gain crossover frequency in rad/s.

Wcp property

Wcp

Phase crossover frequency in rad/s.

__iter__

__iter__()

Enable unpacking: gm_db, pm_deg, wcg, wcp = margin(sys).

Yields:

Type	Description
float	gm_db, pm_deg, wcg, wcp.

Source code in src/ctrlpy/freq_domain.py

def __iter__(self) -> Iterator[float]:
    """Enable unpacking: gm_db, pm_deg, wcg, wcp = margin(sys).

    Yields
    ------
    float
        gm_db, pm_deg, wcg, wcp.
    """
    yield self.gm_db
    yield self.pm_deg
    yield self.wcg
    yield self.wcp

bode_data

bode_data(sys, omega=None)

Calculate frequency response data for Bode diagram analysis.

Evaluates the frequency response:

.. math::

G(j\omega) = |G(j\omega)| e^{j \angle G(j\omega)}

where magnitude in decibels is:

.. math::

|G(j\omega)|_{\text{dB}} = 20 \log_{10} |G(j\omega)|

Parameters:

Name	Type	Description	Default
sys	LinearTimeInvariant	Linear Time-Invariant system (TransferFunction or StateSpace).	required
omega	Sequence[float] | NDArray[floating] | None	Frequencies \(\omega\) in rad/s. If None, auto-generated based on poles and zeros.	None

Returns:

Type	Description
BodeData	Dataclass containing frequencies \(\omega\), linear magnitude \(|G(j\omega)|\), unwrapped phase in degrees, magnitude in dB, phase in radians, and complex response \(G(j\omega)\).

Source code in src/ctrlpy/freq_domain.py

def bode_data(
    sys: LinearTimeInvariant,
    omega: Sequence[float] | NDArray[np.floating] | None = None,
) -> BodeData:
    r"""Calculate frequency response data for Bode diagram analysis.

    Evaluates the frequency response:

    .. math::

        G(j\omega) = |G(j\omega)| e^{j \angle G(j\omega)}

    where magnitude in decibels is:

    .. math::

        |G(j\omega)|_{\text{dB}} = 20 \log_{10} |G(j\omega)|

    Parameters
    ----------
    sys : LinearTimeInvariant
        Linear Time-Invariant system (TransferFunction or StateSpace).
    omega : Sequence[float] | NDArray[np.floating] | None, optional
        Frequencies $\omega$ in rad/s. If None, auto-generated based on poles and zeros.

    Returns
    -------
    BodeData
        Dataclass containing frequencies $\omega$, linear magnitude $|G(j\omega)|$,
        unwrapped phase in degrees, magnitude in dB, phase in radians,
        and complex response $G(j\omega)$.
    """
    w = _generate_omega(sys, omega)
    _, resp = _evaluate_siso_freqresp(sys, w)

    mag = np.abs(resp).astype(np.float64)
    with np.errstate(divide="ignore"):
        mag_db = np.where(mag > 0.0, 20.0 * np.log10(mag), -np.inf)

    phase_rad = np.unwrap(np.angle(resp)).astype(np.float64)
    phase_deg = np.rad2deg(phase_rad)

    return BodeData(
        w=w,
        mag=mag,
        phase=phase_deg,
        mag_db=mag_db,
        phase_rad=phase_rad,
        response=resp,
    )

nyquist_data

nyquist_data(sys, omega=None, n_arc_points=100)

Calculate complex frequency response data for Nyquist diagram analysis.

Evaluates \(G(j\omega)\) along the imaginary axis and optionally around indented arcs \(s = \epsilon e^{j\theta}\) (\(\theta \in [-\pi/2, \pi/2]\)) for poles at \(s=0\).

Parameters:

Name	Type	Description	Default
sys	LinearTimeInvariant	Linear Time-Invariant system (TransferFunction or StateSpace).	required
omega	Sequence[float] | NDArray[floating] | None	Frequencies \(\omega\) in rad/s. If None, auto-generated based on poles and zeros.	None
n_arc_points	int	Number of points to evaluate along the indented D-contour around \(s=0\) when an origin pole is present, by default 100.	100

Returns:

Type	Description
NyquistData	Dataclass containing frequencies \(\omega\), complex frequency response \(G(j\omega)\), real part \(\mathrm{Re}(G(j\omega))\), imaginary part \(\mathrm{Im}(G(j\omega))\), and indented arc data if applicable.

Source code in src/ctrlpy/freq_domain.py

def nyquist_data(
    sys: LinearTimeInvariant,
    omega: Sequence[float] | NDArray[np.floating] | None = None,
    n_arc_points: int = 100,
) -> NyquistData:
    r"""Calculate complex frequency response data for Nyquist diagram analysis.

    Evaluates $G(j\omega)$ along the imaginary axis and optionally around
    indented arcs $s = \epsilon e^{j\theta}$ ($\theta \in [-\pi/2, \pi/2]$) for poles at $s=0$.

    Parameters
    ----------
    sys : LinearTimeInvariant
        Linear Time-Invariant system (TransferFunction or StateSpace).
    omega : Sequence[float] | NDArray[np.floating] | None, optional
        Frequencies $\omega$ in rad/s. If None, auto-generated based on poles and zeros.
    n_arc_points : int, optional
        Number of points to evaluate along the indented D-contour around $s=0$
        when an origin pole is present, by default 100.

    Returns
    -------
    NyquistData
        Dataclass containing frequencies $\omega$, complex frequency response $G(j\omega)$,
        real part $\mathrm{Re}(G(j\omega))$, imaginary part $\mathrm{Im}(G(j\omega))$, and indented arc data if applicable.
    """
    w = _generate_omega(sys, omega)
    tf_sys, resp = _evaluate_siso_freqresp(sys, w)

    poles = tf_sys.poles()
    has_origin_pole = np.any(np.isclose(poles, 0.0, atol=1e-5))

    arc_s: NDArray[np.complex128] | None = None
    arc_resp: NDArray[np.complex128] | None = None

    if has_origin_pole:
        eps = float(w[0])
        # Indented Nyquist D-contour around s=0: evaluate s = eps * exp(j * theta)
        # for theta in [-pi/2, +pi/2] connecting omega = 0^- (s = -j*eps) to omega = 0^+ (s = +j*eps).
        theta = np.linspace(-np.pi / 2.0, np.pi / 2.0, n_arc_points, dtype=np.float64)
        s_arc: NDArray[np.complex128] = (eps * np.exp(1j * theta)).astype(np.complex128)
        num_arc = np.polyval(tf_sys.num, s_arc)
        den_arc = np.polyval(tf_sys.den, s_arc)
        arc_s = s_arc
        arc_resp = np.asarray(num_arc / den_arc, dtype=np.complex128)

    return NyquistData(
        w=w,
        response=resp,
        real=np.real(resp).astype(np.float64),
        imag=np.imag(resp).astype(np.float64),
        arc_s=arc_s,
        arc_response=arc_resp,
    )

root_locus_data

root_locus_data(sys, gains=None)

Calculate closed-loop root locus trajectories for varying gains \(k \ge 0\).

Computes the roots of the characteristic equation:

.. math::

1 + k \cdot G(s) = 0 \iff D(s) + k \cdot N(s) = 0

Parameters:

Name	Type	Description	Default
sys	LinearTimeInvariant	Open-loop SISO LTI system \(G(s) = N(s)/D(s)\).	required
gains	Sequence[float] | NDArray[floating] | None	1D array of non-negative gain values \(k \ge 0\). If None, automatically generated.	None

Returns:

Type	Description
RootLocusData	Dataclass containing gains, 2D roots array of shape (n_gains, n_poles), open-loop poles, and open-loop zeros.

Raises:

Type	Description
ValueError	If gains array is empty or contains negative values.

Source code in src/ctrlpy/freq_domain.py

def root_locus_data(
    sys: LinearTimeInvariant,
    gains: Sequence[float] | NDArray[np.floating] | None = None,
) -> RootLocusData:
    r"""Calculate closed-loop root locus trajectories for varying gains $k \ge 0$.

    Computes the roots of the characteristic equation:

    .. math::

        1 + k \cdot G(s) = 0 \iff D(s) + k \cdot N(s) = 0

    Parameters
    ----------
    sys : LinearTimeInvariant
        Open-loop SISO LTI system $G(s) = N(s)/D(s)$.
    gains : Sequence[float] | NDArray[np.floating] | None, optional
        1D array of non-negative gain values $k \ge 0$. If None, automatically generated.

    Returns
    -------
    RootLocusData
        Dataclass containing gains, 2D roots array of shape `(n_gains, n_poles)`,
        open-loop poles, and open-loop zeros.

    Raises
    ------
    ValueError
        If gains array is empty or contains negative values.
    """
    if not isinstance(sys, LinearTimeInvariant):
        raise TypeError(f"Expected LinearTimeInvariant instance, got {type(sys).__name__}.")
    if not sys.is_siso:
        raise ValueError(
            f"Root locus requires a SISO system, got ({sys.inputs} inputs, {sys.outputs} outputs)."
        )

    if isinstance(sys, TransferFunction):
        tf_sys = sys
    elif isinstance(sys, StateSpace):
        tf_sys = sys.to_tf()
    else:
        raise TypeError(f"Unsupported system type: {type(sys).__name__}")

    ol_poles = tf_sys.poles()
    ol_zeros = tf_sys.zeros()

    if gains is not None:
        gains_arr = np.asarray(gains, dtype=np.float64).ravel()
        if gains_arr.size == 0:
            raise ValueError("gains array cannot be empty.")
        if np.any(gains_arr < 0.0):
            raise ValueError("All gains in root locus must be non-negative (k >= 0).")
    else:
        gains_arr = np.concatenate(([0.0], np.geomspace(1e-4, 1e4, num=500, dtype=np.float64)))

    num = tf_sys.num
    den = tf_sys.den

    n_poles = len(den) - 1
    if n_poles <= 0:
        roots_arr = np.empty((len(gains_arr), 0), dtype=np.complex128)
        return RootLocusData(
            gains=gains_arr,
            roots=roots_arr,
            poles=ol_poles,
            zeros=ol_zeros,
        )

    deg_d = len(den)
    deg_n = len(num)
    max_deg = max(deg_d, deg_n)

    num_pad = np.pad(num, (max_deg - deg_n, 0), mode="constant")
    den_pad = np.pad(den, (max_deg - deg_d, 0), mode="constant")

    n_gains = len(gains_arr)
    roots_matrix = np.zeros((n_gains, n_poles), dtype=np.complex128)

    p0 = np.roots(den)
    idx0 = np.lexsort((p0.real, p0.imag))
    prev_roots = p0[idx0].astype(np.complex128)
    roots_matrix[0] = prev_roots

    for i in range(1, n_gains):
        k = gains_arr[i]
        poly_k = den_pad + k * num_pad
        curr_roots = np.roots(poly_k).astype(np.complex128)

        if curr_roots.size == n_poles:
            dist_matrix = np.abs(curr_roots[:, np.newaxis] - prev_roots[np.newaxis, :])
            row_ind, col_ind = linear_sum_assignment(dist_matrix)
            matched_roots = np.empty(n_poles, dtype=np.complex128)
            matched_roots[col_ind] = curr_roots[row_ind]
            roots_matrix[i] = matched_roots
            prev_roots = matched_roots
        else:
            padded = np.full(n_poles, np.nan + 1j * np.nan, dtype=np.complex128)
            padded[: curr_roots.size] = curr_roots
            roots_matrix[i] = padded

    return RootLocusData(
        gains=gains_arr,
        roots=roots_matrix,
        poles=ol_poles,
        zeros=ol_zeros,
    )

margin

margin(sys)

Compute gain margin, phase margin, and crossover frequencies for a SISO system.

.. math::

\text{GM} = -20 \log_{10} |G(j\omega_{cp})| \text{ dB}, \quad \angle G(j\omega_{cp}) = -180^\circ

.. math::

\text{PM} = 180^\circ + \angle G(j\omega_{cg}), \quad |G(j\omega_{cg})| = 1 \text{ (0 dB)}

Parameters:

Name	Type	Description	Default
sys	LinearTimeInvariant	Open-loop SISO LTI system \(G(s)\).	required

Returns:

Type

Description

StabilityMargins

Dataclass containing: - gm_db : Gain Margin \(\mathrm{GM}\) in dB (\(\infty\) if phase never crosses \(-180^\circ\)) - pm_deg : Phase Margin \(\mathrm{PM}\) in degrees (\(\infty\) if magnitude never crosses \(0\text{ dB}\)) - wcg : Gain Crossover Frequency \(\omega_{cg}\) in rad/s (where \(|G(j\omega)| = 1\)) - wcp : Phase Crossover Frequency \(\omega_{cp}\) in rad/s (where \(\angle G(j\omega) = -180^\circ\))

Source code in src/ctrlpy/freq_domain.py

def margin(sys: LinearTimeInvariant) -> StabilityMargins:
    r"""Compute gain margin, phase margin, and crossover frequencies for a SISO system.

    .. math::

        \text{GM} = -20 \log_{10} |G(j\omega_{cp})| \text{ dB}, \quad \angle G(j\omega_{cp}) = -180^\circ

    .. math::

        \text{PM} = 180^\circ + \angle G(j\omega_{cg}), \quad |G(j\omega_{cg})| = 1 \text{ (0 dB)}

    Parameters
    ----------
    sys : LinearTimeInvariant
        Open-loop SISO LTI system $G(s)$.

    Returns
    -------
    StabilityMargins
        Dataclass containing:
        - `gm_db` : Gain Margin $\mathrm{GM}$ in dB ($\infty$ if phase never crosses $-180^\circ$)
        - `pm_deg` : Phase Margin $\mathrm{PM}$ in degrees ($\infty$ if magnitude never crosses $0\text{ dB}$)
        - `wcg` : Gain Crossover Frequency $\omega_{cg}$ in rad/s (where $|G(j\omega)| = 1$)
        - `wcp` : Phase Crossover Frequency $\omega_{cp}$ in rad/s (where $\angle G(j\omega) = -180^\circ$)
    """
    if not isinstance(sys, LinearTimeInvariant):
        raise TypeError(f"Expected LinearTimeInvariant instance, got {type(sys).__name__}.")
    if not sys.is_siso:
        raise ValueError(
            f"Stability margin analysis requires a SISO system, got "
            f"({sys.inputs} inputs, {sys.outputs} outputs)."
        )

    if isinstance(sys, TransferFunction):
        tf_sys = sys
    elif isinstance(sys, StateSpace):
        tf_sys = sys.to_tf()
    else:
        raise TypeError(f"Unsupported system type: {type(sys).__name__}")

    # Generate a wide, dense grid for initial bracket detection
    w_grid = _generate_omega(sys, n_points=4000)
    min_log = min(np.log10(w_grid[0]), -4.0)
    max_log = max(np.log10(w_grid[-1]), 5.0)
    w_dense = np.logspace(min_log, max_log, num=5000, dtype=np.float64)

    _, resp_dense = _evaluate_siso_freqresp(tf_sys, w_dense)
    mag_dense = np.abs(resp_dense)
    phase_rad_dense = np.unwrap(np.angle(resp_dense))
    phase_deg_dense = np.rad2deg(phase_rad_dense)

    with np.errstate(divide="ignore"):
        mag_db_dense = np.where(mag_dense > 0.0, 20.0 * np.log10(mag_dense), -np.inf)

    # 1. Gain Crossover Frequency (Wcg) where |G(jw)| = 1 (0 dB)
    wcg_candidates: list[float] = []
    sign_mag = np.sign(mag_db_dense)

    def mag_objective(w: float) -> float:
        _, r = _evaluate_siso_freqresp(tf_sys, np.array([w], dtype=np.float64))
        m = float(np.abs(r[0]))
        return float(20.0 * np.log10(m)) if m > 0 else -1000.0

    for i in range(len(w_dense) - 1):
        if sign_mag[i] * sign_mag[i + 1] <= 0 and sign_mag[i] != sign_mag[i + 1]:
            try:
                sol = root_scalar(
                    mag_objective,
                    bracket=[w_dense[i], w_dense[i + 1]],
                    method="brentq",
                )
                if sol.converged:
                    wcg_candidates.append(float(sol.root))
            except (ValueError, RuntimeError):
                pass

    if wcg_candidates:
        wcg = float(wcg_candidates[0])
        _, r_wcg = _evaluate_siso_freqresp(tf_sys, np.array([wcg], dtype=np.float64))
        raw_phi_deg = float(np.rad2deg(np.angle(r_wcg[0])))
        pm_val = (180.0 + raw_phi_deg + 180.0) % 360.0 - 180.0
        pm_deg = float(pm_val)
    else:
        wcg = float("nan")
        pm_deg = float("inf")

    # 2. Phase Crossover Frequency (Wcp) where phase crosses -180 deg
    wcp_candidates: list[float] = []

    def phase_objective(w: float) -> float:
        _, r = _evaluate_siso_freqresp(tf_sys, np.array([w], dtype=np.float64))
        raw_phi = float(np.rad2deg(np.angle(r[0])))
        err = (raw_phi + 180.0 + 180.0) % 360.0 - 180.0
        return err

    phase_err_dense = (phase_deg_dense + 180.0 + 180.0) % 360.0 - 180.0
    sign_phase = np.sign(phase_err_dense)

    for i in range(len(w_dense) - 1):
        if (
            sign_phase[i] * sign_phase[i + 1] <= 0
            and sign_phase[i] != sign_phase[i + 1]
            and abs(phase_err_dense[i] - phase_err_dense[i + 1]) < 180.0
        ):
            try:
                sol = root_scalar(
                    phase_objective,
                    bracket=[w_dense[i], w_dense[i + 1]],
                    method="brentq",
                )
                if sol.converged:
                    wcp_candidates.append(float(sol.root))
            except (ValueError, RuntimeError):
                pass

    if wcp_candidates:
        wcp = float(wcp_candidates[0])
        _, r_wcp = _evaluate_siso_freqresp(tf_sys, np.array([wcp], dtype=np.float64))
        mag_wcp = float(np.abs(r_wcp[0]))
        if mag_wcp > 0.0:
            gm_db = float(-20.0 * np.log10(mag_wcp))
        else:
            gm_db = float("inf")
    else:
        wcp = float("nan")
        gm_db = float("inf")

    return StabilityMargins(
        gm_db=gm_db,
        pm_deg=pm_deg,
        wcg=wcg,
        wcp=wcp,
    )