Reflection coefficients as a stability coordinate system

Tutorial goal

Show why bounded scalar reflection/PARCOR coefficients are a practical stability coordinate system.

Note

New to the terminology? See the lattice DSP concept map and the causality/data-use guide for how online, offline, block, and MIMO examples should be read.

Context

Scalar IIR stability is hard to read from direct denominator coefficients. In lattice coordinates, the all-pole stability condition is exposed stage by stage. This tutorial compares reflection-parameterized denominators with a few direct denominator examples near the unit-circle boundary.

Key idea and equations

For a scalar all-pole denominator built by the Schur step-up recursion,

\[|k_i| < 1 \quad \text{for every stage}\]

keeps the poles inside the unit disk. Direct denominator coefficients do not provide an equivalent per-coefficient test.

How to read the result

Compare the maximum reflection magnitude with the maximum pole radius. The direct denominator rows show why pole or reflection diagnostics are still needed outside lattice coordinates.

Run command

python examples/reflection_coefficients_stability_demo.py

Run status

Return code: 0

Captured stdout

reflection coefficients and scalar IIR stability
====================================================
mild reflection vector:
  type: reflection
  parameters: [0.2, -0.3, 0.1]
  max |reflection|: 0.3
  max pole radius: 0.71056
  stable: True
moderate reflection vector:
  type: reflection
  parameters: [0.75, -0.45, 0.25]
  max |reflection|: 0.75
  max pole radius: 0.946016
  stable: True
near-boundary reflection vector:
  type: reflection
  parameters: [0.95, -0.8, 0.65]
  max |reflection|: 0.95
  max pole radius: 0.998825
  stable: True
stable direct denominator:
  type: direct denominator
  parameters: [1.0, -1.85, 0.95]
  max |reflection|: 0.95
  max pole radius: 0.974679
  stable: True
nearby unstable direct denominator:
  type: direct denominator
  parameters: [1.0, -2.05, 1.05]
  max |reflection|: nan
  max pole radius: 1.05
  stable: False
another unstable direct denominator:
  type: direct denominator
  parameters: [1.0, 0.0, 1.05]
  max |reflection|: nan
  max pole radius: 1.0247
  stable: False

Captured stderr

Matplotlib is building the font cache; this may take a moment.

Figures

reflection coefficients stability pole radius

reflection_coefficients_stability_pole_radius.png

Generated data files

Source code

  1"""Reflection coefficients as a stability coordinate system.
  2
  3This example complements ``stability_vs_direct_iir.py``.  Instead of an adaptive
  4failure case, it shows the static parameterization: scalar all-pole denominators
  5built from reflection/PARCOR coefficients with ``|k_i| < 1`` have poles inside
  6the unit disk, while direct denominator coefficients do not expose such a simple
  7per-coefficient stability rule.
  8"""
  9
 10from __future__ import annotations
 11
 12import csv
 13import os
 14from pathlib import Path
 15
 16import numpy as np
 17
 18import lattice_dsp as ld
 19
 20
 21def artifact_dir() -> Path:
 22    path = Path(os.environ.get("LATTICE_DSP_ARTIFACT_DIR", "reports/example-artifacts"))
 23    path.mkdir(parents=True, exist_ok=True)
 24    return path
 25
 26
 27def max_pole_radius(denominator: np.ndarray) -> float:
 28    roots = np.roots(np.asarray(denominator, dtype=float))
 29    return float(np.max(np.abs(roots))) if roots.size else 0.0
 30
 31
 32def write_rows(path: Path, rows: list[dict[str, object]]) -> None:
 33    with path.open("w", newline="", encoding="utf-8") as f:
 34        writer = csv.DictWriter(
 35            f,
 36            fieldnames=[
 37                "label",
 38                "parameter_type",
 39                "parameters",
 40                "max_abs_reflection",
 41                "max_pole_radius",
 42                "stable",
 43            ],
 44        )
 45        writer.writeheader()
 46        writer.writerows(rows)
 47
 48
 49def main() -> None:
 50    out_dir = artifact_dir()
 51
 52    reflection_cases = {
 53        "mild reflection vector": np.array([0.20, -0.30, 0.10]),
 54        "moderate reflection vector": np.array([0.75, -0.45, 0.25]),
 55        "near-boundary reflection vector": np.array([0.95, -0.80, 0.65]),
 56    }
 57    direct_denominators = {
 58        "stable direct denominator": np.array([1.0, -1.85, 0.95]),
 59        "nearby unstable direct denominator": np.array([1.0, -2.05, 1.05]),
 60        "another unstable direct denominator": np.array([1.0, 0.0, 1.05]),
 61    }
 62
 63    rows: list[dict[str, object]] = []
 64    for label, k in reflection_cases.items():
 65        den = np.asarray(ld.reflection_to_denominator(k), dtype=float)
 66        radius = max_pole_radius(den)
 67        rows.append(
 68            {
 69                "label": label,
 70                "parameter_type": "reflection",
 71                "parameters": np.round(k, 5).tolist(),
 72                "max_abs_reflection": float(np.max(np.abs(k))),
 73                "max_pole_radius": radius,
 74                "stable": radius < 1.0,
 75            }
 76        )
 77
 78    for label, den in direct_denominators.items():
 79        radius = max_pole_radius(den)
 80        try:
 81            reflection = np.asarray(ld.denominator_to_reflection(den), dtype=float)
 82            max_abs_reflection = float(np.max(np.abs(reflection)))
 83        except Exception:
 84            max_abs_reflection = float("nan")
 85        rows.append(
 86            {
 87                "label": label,
 88                "parameter_type": "direct denominator",
 89                "parameters": np.round(den, 5).tolist(),
 90                "max_abs_reflection": max_abs_reflection,
 91                "max_pole_radius": radius,
 92                "stable": radius < 1.0,
 93            }
 94        )
 95
 96    print("reflection coefficients and scalar IIR stability")
 97    print("=" * 52)
 98    for row in rows:
 99        print(f"{row['label']}:")
100        print(f"  type: {row['parameter_type']}")
101        print(f"  parameters: {row['parameters']}")
102        print(f"  max |reflection|: {row['max_abs_reflection']:.6g}")
103        print(f"  max pole radius: {row['max_pole_radius']:.6g}")
104        print(f"  stable: {row['stable']}")
105
106    csv_path = out_dir / "reflection_coefficients_stability_summary.csv"
107    write_rows(csv_path, rows)
108    print(f"\nwrote {csv_path}")
109
110    try:
111        import matplotlib.pyplot as plt
112    except Exception:
113        print("matplotlib is not installed; skipped figures")
114        return
115
116    labels = [str(row["label"]) for row in rows]
117    pole_radii = [float(row["max_pole_radius"]) for row in rows]
118    x = np.arange(len(labels))
119    fig, ax = plt.subplots(figsize=(9, 4.8))
120    ax.bar(x, pole_radii)
121    ax.axhline(1.0, linestyle="--", linewidth=1.0)
122    ax.set_title("Pole radius from reflection versus direct denominator examples")
123    ax.set_ylabel("max pole radius")
124    ax.set_xticks(x)
125    ax.set_xticklabels(labels, rotation=30, ha="right")
126    fig.tight_layout()
127    fig_path = out_dir / "reflection_coefficients_stability_pole_radius.png"
128    fig.savefig(fig_path, dpi=160)
129    print(f"wrote {fig_path}")
130
131
132if __name__ == "__main__":
133    main()