Diagonal MIMO equals independent SISO

Tutorial goal

Use independent SISO lattice filters and online predictors as sanity checks for diagonal MIMO.

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

Before studying coupled MIMO lattice systems, it helps to inspect the diagonal case. If every matrix Markov parameter or matrix reflection coefficient is diagonal, no output channel depends on any other input channel. The MIMO system is exactly several scalar SISO systems running side by side.

This is the intuition bridge between scalar lattice filters and matrix-valued MIMO filters: scalar reflection coefficients become diagonal entries first, and only later become genuinely coupled matrices.

Key idea and equations

A diagonal MIMO impulse response has Markov matrices

\[M_k = \operatorname{diag}(h^{(1)}_k, \ldots, h^{(p)}_k),\]

and therefore a diagonal transfer matrix

\[H(z)=\operatorname{diag}\bigl(H_1(z),\ldots,H_p(z)\bigr).\]

The convolution separates by channel:

\[y_i[n] = \sum_k h^{(i)}_k x_i[n-k], \qquad y_i \text{ does not depend on } x_j \text{ for } j\ne i.\]

The same diagonal reduction should hold for the online lattice predictor. If the matrix reflection coefficients are diagonal,

\[K_m = \operatorname{diag}(k_{m,1},\ldots,k_{m,p}), \qquad L_m = \operatorname{diag}(\ell_{m,1},\ldots,\ell_{m,p}),\]

then the vector recursion decouples into independent one-channel recursions:

\[\widehat y[n] = \begin{bmatrix} \widehat y_1[n] & \cdots & \widehat y_p[n] \end{bmatrix}^T.\]

This is the algebraic reason the diagonal MIMO result must match running scalar SISO filters or one-channel online lattice predictors independently. Any nonzero off-diagonal Markov block or reflection entry would create true channel coupling.

Causality and data use

The scalar LatticeIIR filters and the online MIMOLatticePredictor comparison are causal. The Markov-convolution part is evaluated on a stored array as a validation check, but each output sample uses only lags x[n-k]. The online predictor part explicitly calls predict() before update(y_n) and therefore uses only previous samples when forming \widehat y[n].

What this example verifies

This is the reduction-to-SISO sanity check. It verifies that when the MIMO Markov matrices or matrix reflections are diagonal, the multichannel object decomposes into independent scalar systems. The expected diagnostic is roundoff-level agreement between the MIMO result and stacked SISO results, including the online predict()/update() predictor path.

How to read the result

The MIMO/SISO differences should be near numerical precision for both the finite impulse-response comparison and the online predict-before-update comparison.

Run command

python examples/mimo_diagonal_equals_independent_siso.py

Run status

Return code: 0

Captured stdout

channels: 5
samples: 700
impulse truncation length: 160
off-diagonal Markov energy: 0.000e+00
diagonal Markov energy: 5.389e+00
max |MIMO output - independent SISO output|: 4.441e-15
RMS error: 7.699e-16
online predictor channels: 3
max |online diagonal MIMO prediction - independent SISO prediction|: 0.000e+00
max |online diagonal MIMO error - independent SISO error|: 0.000e+00

Figures

mimo diagonal error

mimo_diagonal_error.png

mimo diagonal online prediction error

mimo_diagonal_online_prediction_error.png

mimo diagonal outputs

mimo_diagonal_outputs.png

Generated data files

Source code

  1"""Tutorial: diagonal MIMO equals independent SISO filters.
  2
  3A useful sanity check for MIMO DSP is the diagonal case.  If every matrix
  4Markov parameter or reflection matrix is diagonal, the MIMO system does not mix
  5channels; it is just several scalar SISO systems running side by side.
  6
  7This tutorial performs two checks:
  8
  9* five stable SISO lattice IIR filters are placed on the diagonal of a MIMO
 10  Markov sequence and compared with independent SISO filtering; and
 11* a three-channel online MIMO lattice predictor with diagonal reflection
 12  matrices is compared sample-by-sample with three independent one-channel
 13  online lattice predictors.
 14"""
 15
 16from __future__ import annotations
 17
 18import csv
 19import os
 20from pathlib import Path
 21
 22import numpy as np
 23
 24import lattice_dsp as ld
 25
 26
 27def artifact_dir() -> Path:
 28    path = Path(os.environ.get("LATTICE_DSP_ARTIFACT_DIR", "reports/example-artifacts"))
 29    path.mkdir(parents=True, exist_ok=True)
 30    return path
 31
 32
 33def diagonal_markov_from_siso(
 34    filters: list[tuple[np.ndarray, np.ndarray]], n_impulse: int
 35) -> np.ndarray:
 36    channels = len(filters)
 37    markov = np.zeros((n_impulse, channels, channels), dtype=float)
 38    for ch, (reflection, numerator) in enumerate(filters):
 39        denominator = ld.reflection_to_denominator(reflection.tolist())
 40        h = ld.iir_impulse_response(denominator, numerator.tolist(), n_impulse)
 41        markov[:, ch, ch] = np.asarray(h, dtype=float)
 42    return markov
 43
 44
 45def mimo_convolve(markov: np.ndarray, x: np.ndarray) -> np.ndarray:
 46    samples, inputs = x.shape
 47    _, outputs, markov_inputs = markov.shape
 48    if inputs != markov_inputs:
 49        raise ValueError("input dimension does not match Markov parameters")
 50    y = np.zeros((samples, outputs), dtype=float)
 51    horizon = markov.shape[0]
 52    for n in range(samples):
 53        max_lag = min(horizon, n + 1)
 54        for lag in range(max_lag):
 55            y[n] += markov[lag] @ x[n - lag]
 56    return y
 57
 58
 59def online_mimo_vs_independent_siso(
 60    x: np.ndarray,
 61    forward_scalars: np.ndarray,
 62    backward_scalars: np.ndarray,
 63) -> tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray]:
 64    """Compare a diagonal online MIMO lattice with independent one-channel predictors."""
 65
 66    kf = np.asarray([np.diag(row) for row in forward_scalars], dtype=float)
 67    kb = np.asarray([np.diag(row) for row in backward_scalars], dtype=float)
 68    mimo = ld.MIMOLatticePredictor(kf, kb)
 69    siso = [
 70        ld.MIMOLatticePredictor(
 71            forward_scalars[:, ch].reshape(-1, 1, 1),
 72            backward_scalars[:, ch].reshape(-1, 1, 1),
 73        )
 74        for ch in range(x.shape[1])
 75    ]
 76
 77    prediction_mimo = np.empty_like(x, dtype=float)
 78    prediction_siso = np.empty_like(x, dtype=float)
 79    error_mimo = np.empty_like(x, dtype=float)
 80    error_siso = np.empty_like(x, dtype=float)
 81
 82    for n, sample in enumerate(x):
 83        # The prediction is intentionally requested before the current sample is
 84        # passed into update().  This is the causal one-step prediction contract.
 85        prediction_mimo[n] = mimo.predict().real
 86        error_mimo[n] = mimo.update(sample).real
 87        for ch, predictor in enumerate(siso):
 88            prediction_siso[n, ch] = predictor.predict()[0].real
 89            error_siso[n, ch] = predictor.update(np.array([sample[ch]]))[0].real
 90
 91    return prediction_mimo, prediction_siso, error_mimo, error_siso
 92
 93
 94def main() -> None:
 95    out_dir = artifact_dir()
 96    rng = np.random.default_rng(2026)
 97
 98    channels = 5
 99    samples = 700
100    n_impulse = 160
101
102    # Five different stable SISO lattice IIR filters.  Each channel has its own
103    # reflection coefficients and numerator, but there is no cross-channel mixing.
104    filters: list[tuple[np.ndarray, np.ndarray]] = []
105    for ch in range(channels):
106        scale = 0.58 - 0.04 * ch
107        reflection = scale * np.array([0.55, -0.42, 0.28, -0.16], dtype=float)
108        numerator = np.array([1.0, 0.08 * (ch + 1), -0.04, 0.015, 0.0], dtype=float)
109        filters.append((reflection, numerator))
110
111    x = rng.normal(size=(samples, channels))
112    siso_y = np.zeros_like(x)
113    for ch, (reflection, numerator) in enumerate(filters):
114        filt = ld.LatticeIIR(reflection.tolist(), numerator.tolist())
115        siso_y[:, ch] = filt.process(x[:, ch])
116
117    markov = diagonal_markov_from_siso(filters, n_impulse=n_impulse)
118    mimo_y = mimo_convolve(markov, x)
119
120    # Ignore the very end of the impulse truncation error by using a long enough
121    # impulse response.  With these stable filters, the residual is near roundoff.
122    error = mimo_y - siso_y
123    max_abs_error = float(np.max(np.abs(error)))
124    rms_error = float(np.sqrt(np.mean(error**2)))
125    off_diagonal_energy = float(np.sum(markov[:, ~np.eye(channels, dtype=bool)] ** 2))
126    diagonal_energy = float(np.sum(markov[:, np.eye(channels, dtype=bool)] ** 2))
127
128    rows = []
129    for ch in range(channels):
130        rows.append(
131            {
132                "channel": ch,
133                "max_abs_error": float(np.max(np.abs(error[:, ch]))),
134                "rms_error": float(np.sqrt(np.mean(error[:, ch] ** 2))),
135                "max_abs_reflection": float(np.max(np.abs(filters[ch][0]))),
136            }
137        )
138
139    csv_path = out_dir / "mimo_diagonal_equals_independent_siso_summary.csv"
140    with csv_path.open("w", newline="", encoding="utf-8") as f:
141        writer = csv.DictWriter(f, fieldnames=list(rows[0]))
142        writer.writeheader()
143        writer.writerows(rows)
144
145    online_channels = 3
146    online_x = rng.normal(size=(samples, online_channels))
147    forward_scalars = np.array(
148        [
149            [0.18, -0.12, 0.07],
150            [-0.05, 0.03, -0.02],
151            [0.015, -0.01, 0.012],
152        ],
153        dtype=float,
154    )
155    backward_scalars = np.array(
156        [
157            [0.16, -0.10, 0.06],
158            [-0.04, 0.02, -0.01],
159            [0.012, -0.008, 0.010],
160        ],
161        dtype=float,
162    )
163    pred_mimo, pred_siso, err_mimo, err_siso = online_mimo_vs_independent_siso(
164        online_x,
165        forward_scalars,
166        backward_scalars,
167    )
168    online_prediction_diff = pred_mimo - pred_siso
169    online_error_diff = err_mimo - err_siso
170    max_online_prediction_diff = float(np.max(np.abs(online_prediction_diff)))
171    max_online_error_diff = float(np.max(np.abs(online_error_diff)))
172
173    online_rows = []
174    for ch in range(online_channels):
175        online_rows.append(
176            {
177                "channel": ch,
178                "max_abs_prediction_difference": float(
179                    np.max(np.abs(online_prediction_diff[:, ch]))
180                ),
181                "max_abs_error_difference": float(np.max(np.abs(online_error_diff[:, ch]))),
182                "max_abs_forward_reflection": float(np.max(np.abs(forward_scalars[:, ch]))),
183                "max_abs_backward_reflection": float(np.max(np.abs(backward_scalars[:, ch]))),
184            }
185        )
186
187    online_csv_path = out_dir / "mimo_diagonal_online_predictor_summary.csv"
188    with online_csv_path.open("w", newline="", encoding="utf-8") as f:
189        writer = csv.DictWriter(f, fieldnames=list(online_rows[0]))
190        writer.writeheader()
191        writer.writerows(online_rows)
192
193    print("channels:", channels)
194    print("samples:", samples)
195    print("impulse truncation length:", n_impulse)
196    print("off-diagonal Markov energy:", f"{off_diagonal_energy:.3e}")
197    print("diagonal Markov energy:", f"{diagonal_energy:.3e}")
198    print("max |MIMO output - independent SISO output|:", f"{max_abs_error:.3e}")
199    print("RMS error:", f"{rms_error:.3e}")
200    print("online predictor channels:", online_channels)
201    print(
202        "max |online diagonal MIMO prediction - independent SISO prediction|:",
203        f"{max_online_prediction_diff:.3e}",
204    )
205    print(
206        "max |online diagonal MIMO error - independent SISO error|:", f"{max_online_error_diff:.3e}"
207    )
208    print(f"wrote {csv_path}")
209    print(f"wrote {online_csv_path}")
210
211    try:
212        import matplotlib.pyplot as plt
213    except Exception:
214        print("matplotlib is not installed; skipped figures")
215        return
216
217    t = np.arange(160)
218    fig, axes = plt.subplots(channels, 1, figsize=(9, 8), sharex=True)
219    for ch, ax in enumerate(axes):
220        ax.plot(t, siso_y[: t.size, ch], label="independent SISO", linewidth=1.7)
221        ax.plot(t, mimo_y[: t.size, ch], "--", label="diagonal MIMO", linewidth=1.2)
222        ax.set_ylabel(f"ch {ch}")
223        ax.grid(True, alpha=0.25)
224    axes[0].set_title("Diagonal MIMO reproduces five independent SISO lattice IIR outputs")
225    axes[-1].set_xlabel("sample")
226    axes[0].legend(loc="upper right")
227    fig.tight_layout()
228    fig_path = out_dir / "mimo_diagonal_outputs.png"
229    fig.savefig(fig_path, dpi=160)
230    print(f"wrote {fig_path}")
231
232    fig2, ax2 = plt.subplots(figsize=(8.5, 4.4))
233    for ch in range(channels):
234        ax2.semilogy(np.abs(error[:, ch]) + 1e-18, label=f"channel {ch}")
235    ax2.set_title("Absolute difference between diagonal MIMO and independent SISO outputs")
236    ax2.set_xlabel("sample")
237    ax2.set_ylabel("absolute error")
238    ax2.grid(True, alpha=0.3)
239    ax2.legend(ncol=2)
240    fig2.tight_layout()
241    fig2_path = out_dir / "mimo_diagonal_error.png"
242    fig2.savefig(fig2_path, dpi=160)
243    print(f"wrote {fig2_path}")
244
245    fig3, ax3 = plt.subplots(figsize=(8.5, 4.5))
246    for ch in range(online_channels):
247        ax3.semilogy(np.abs(online_prediction_diff[:, ch]) + 1e-18, label=f"channel {ch}")
248    ax3.set_title("Online diagonal MIMO predictor equals independent one-channel predictors")
249    ax3.set_xlabel("sample")
250    ax3.set_ylabel("absolute prediction difference")
251    ax3.grid(True, alpha=0.3)
252    ax3.legend(ncol=online_channels)
253    fig3.tight_layout()
254    fig3_path = out_dir / "mimo_diagonal_online_prediction_error.png"
255    fig3.savefig(fig3_path, dpi=160)
256    print(f"wrote {fig3_path}")
257
258
259if __name__ == "__main__":
260    main()