gds-viz: frequency response and time-domain response plots

[rororowyourboat]

Summary

Add frequency.py and response.py to gds_viz/ providing Bode plots, Nyquist plots, Nichols charts, root locus diagrams, and step/impulse response plots. Follows the same lazy-import pattern as phase.py. Behind a new [control] optional extra (matplotlib + numpy only, no gds-continuous needed).

Parent issue: #198 (classical control theory stack)
Depends on: #201 (numerical linear analysis — data source), #202 (response metrics — annotation data)

Motivation

Frequency response and step response plots are the visual language of control engineering. The numerical data will be computed by gds-analysis/linear.py and gds-analysis/response.py — this package provides the matplotlib rendering. Follows the existing phase.py pattern: accept data, plot it, return a figure.

Proposed API

New file: gds_viz/frequency.py

All functions accept plain list[float] / numpy arrays. No imports from gds-domains or gds-analysis.

import matplotlib.pyplot as plt
from matplotlib.figure import Figure

def bode_plot(
    omega: list[float],
    mag_db: list[float],
    phase_deg: list[float],
    *,
    title: str = "Bode Plot",
    gain_margin: tuple[float, float] | None = None,   # (gm_db, gm_freq)
    phase_margin: tuple[float, float] | None = None,   # (pm_deg, pm_freq)
    ax: tuple | None = None,
    figsize: tuple[float, float] = (10, 8),
) -> Figure:
    """Bode magnitude and phase plot.
    
    Two vertically stacked subplots:
    - Top: magnitude (dB) vs. frequency (rad/s, log scale)
    - Bottom: phase (degrees) vs. frequency (rad/s, log scale)
    
    If gain_margin/phase_margin are provided, annotates them
    with vertical lines and text labels.
    """

def nyquist_plot(
    real: list[float],
    imag: list[float],
    *,
    title: str = "Nyquist Plot",
    unit_circle: bool = True,
    critical_point: bool = True,  # mark (-1, 0)
    ax: plt.Axes | None = None,
    figsize: tuple[float, float] = (8, 8),
) -> Figure:
    """Nyquist diagram (polar frequency response).
    
    Plots H(jω) in the complex plane.
    Optional unit circle and critical point (-1+0j) for stability assessment.
    """

def nichols_plot(
    phase_deg: list[float],
    mag_db: list[float],
    *,
    title: str = "Nichols Chart",
    m_circles: bool = True,       # closed-loop magnitude contours
    n_circles: bool = True,       # closed-loop phase contours
    ax: plt.Axes | None = None,
    figsize: tuple[float, float] = (10, 8),
) -> Figure:
    """Nichols chart: open-loop phase (x) vs. open-loop gain (y).
    
    M-circles: contours of constant closed-loop magnitude
    N-circles: contours of constant closed-loop phase
    Both computed analytically from T = L/(1+L).
    """

def root_locus_plot(
    num: list[float],
    den: list[float],
    *,
    gains: list[float] | None = None,
    title: str = "Root Locus",
    mark_poles: bool = True,
    mark_zeros: bool = True,
    ax: plt.Axes | None = None,
    figsize: tuple[float, float] = (8, 8),
) -> Figure:
    """Root locus diagram.
    
    Shows pole migration as gain K varies from 0 to max.
    Open-loop poles (×) and zeros (○) are marked.
    
    If gains is None, auto-selects 500 gains from 0 to a heuristic max.
    Poles at each gain computed via numpy.roots(den + K*num).
    """

New file: gds_viz/response.py

def step_response_plot(
    times: list[float],
    values: list[float],
    *,
    setpoint: float | None = None,
    metrics: "StepMetrics | None" = None,    # from gds_analysis.response
    title: str = "Step Response",
    ax: plt.Axes | None = None,
    figsize: tuple[float, float] = (10, 6),
) -> Figure:
    """Step response plot with optional metric annotations.
    
    If metrics is provided, annotates:
    - Rise time (shaded region from 10% to 90%)
    - Settling time (vertical dashed line)
    - Overshoot (horizontal dashed line at peak)
    - Steady-state error (gap annotation)
    - ±2% settling band (shaded horizontal band)
    """

def impulse_response_plot(
    times: list[float],
    values: list[float],
    *,
    title: str = "Impulse Response",
    ax: plt.Axes | None = None,
    figsize: tuple[float, float] = (10, 6),
) -> Figure:
    """Impulse response plot."""

def compare_responses(
    responses: list[tuple[list[float], list[float], str]],
    *,
    title: str = "Response Comparison",
    figsize: tuple[float, float] = (10, 6),
) -> Figure:
    """Overlay multiple step/impulse responses for comparison.
    
    Each tuple is (times, values, label).
    Useful for comparing before/after controller tuning,
    or different operating points in a gain schedule.
    """

Export from __init__.py

Use the same lazy __getattr__ pattern as phase_portrait:

# In gds_viz/__init__.py
_CONTROL_EXPORTS = {
    "bode_plot", "nyquist_plot", "nichols_plot", "root_locus_plot",
    "step_response_plot", "impulse_response_plot", "compare_responses",
}

def __getattr__(name):
    if name in _CONTROL_EXPORTS:
        if name in {"bode_plot", "nyquist_plot", "nichols_plot", "root_locus_plot"}:
            from gds_viz.frequency import ...
        else:
            from gds_viz.response import ...
    ...

Implementation Notes

M-circles and N-circles (Nichols chart)

M-circles are contours of |T(jω)| = M where T = L/(1+L). In Nichols coordinates (open-loop phase θ, open-loop gain G dB):

M-circle center: (phase, gain) = (-180°, 20*log10(M²/(M²-1)))
M-circle radius: 20*log10(M/(M²-1))

Standard M-values: 0.25, 0.5, 1, 2, 3, 6 dB. Drawn as closed contours on the Nichols grid.

Root locus computation

For each gain K in gains:

closed_loop_den = [d + K*n for d, n in zip_longest(den, num, fillvalue=0)]
poles = numpy.roots(closed_loop_den)

Track pole trajectories by sorting poles at each K to minimize distance from previous K's poles (Hungarian algorithm or greedy nearest-neighbor).

Style

Follow existing phase.py conventions:

• Default matplotlib style (no seaborn or custom stylesheets)
• Return Figure object (caller decides plt.show() vs. fig.savefig())
• Accept optional ax parameter for subplot embedding
• Label axes with units (rad/s, dB, degrees)

Dependency

New [control] extra in pyproject.toml:

[project.optional-dependencies]
phase = ["matplotlib>=3.8", "numpy>=1.26", "gds-continuous>=0.1.0"]
control = ["matplotlib>=3.8", "numpy>=1.26"]
plots = ["matplotlib>=3.8", "numpy>=1.26", "gds-continuous>=0.1.0"]  # union

[control] does NOT require gds-continuous — these plots take precomputed data. [plots] is a convenience union of both.

Key Files

• packages/gds-viz/gds_viz/phase.py — reference implementation for matplotlib usage, lazy import, style
• packages/gds-viz/gds_viz/__init__.py — lazy __getattr__ pattern to replicate
• packages/gds-viz/pyproject.toml — add [control] and [plots] extras

Testing

• test_frequency.py:
  • bode_plot() with synthetic data → verify figure has 2 axes, correct labels, semilog x-axis
  • nyquist_plot() with unit circle data → verify critical point marker exists
  • nichols_plot() with simple data → verify M-circle contours rendered
  • root_locus_plot() with 1/(s²+s+1) → verify 2 pole trajectories, correct open-loop pole markers
• test_response_viz.py:
  • step_response_plot() with metrics annotation → verify dashed lines for settling time, rise time band
  • compare_responses() with 3 responses → verify 3 line artists in axes

Note: plot tests validate figure structure (axes count, artist count, labels), not pixel-level rendering.

Concepts Addressed (MATLAB Tech Talks)

• Video 3: Step response visualization with performance metric annotations
• Video 4: Bode plot, Nyquist plot, Nichols chart — the three frequency response representations
• Video 7: Loop shaping visualization via Nichols chart
• Video 10: Notch filter visualization via Bode plot (filter response is just a TF → Bode)
• Video 13: Root locus for non-minimum phase analysis