.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "api/gallery/signal_processing/spectral/plot_08_levinson_durbin.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note :ref:`Go to the end ` to download the full example code or to run this example in your browser via JupyterLite. .. rst-class:: sphx-glr-example-title .. _sphx_glr_api_gallery_signal_processing_spectral_plot_08_levinson_durbin.py: Durbin's recursion for autoregressive fitting (1960) ==================================================== The Yule-Walker equations have a Toeplitz matrix, and Norman Levinson (1947) had solved such systems, for Wiener's filter, by a recursion on the order. James Durbin (1960) recast it for fitting autoregressions: the order-:math:`m` model is the order-:math:`(m-1)` one plus a correction weighted by a *reflection coefficient* :math:`k_m`. It costs :math:`p^2` operations instead of :math:`p^3`, gives every order on the way, so the prediction error can be watched as the order grows, and :math:`|k_m| < 1` certifies a stable model (:func:`~engineeringkit.signal_processing.levinson_durbin`). Linear predictive coding of speech is built on it. .. GENERATED FROM PYTHON SOURCE LINES 18-24 .. code-block:: Python import matplotlib.pyplot as plt import numpy as np from scipy import linalg, signal from engineeringkit.signal_processing import levinson_durbin .. GENERATED FROM PYTHON SOURCE LINES 25-30 A fourth-order process ---------------------- Two resonances. The biased sample autocorrelation of 20 000 samples feeds the recursion up to order 10. .. GENERATED FROM PYTHON SOURCE LINES 30-39 .. code-block:: Python poles = [0.95 * np.exp(1j * 0.6), 0.85 * np.exp(1j * 1.9)] a_true = np.real(np.poly(poles + [p.conjugate() for p in poles])) rng = np.random.default_rng(1960) x = signal.lfilter([1.0], a_true, rng.normal(size=20000)) n = x.size r = np.array([x[: n - k] @ x[k:] for k in range(11)]) / n model = levinson_durbin(r, 10) .. GENERATED FROM PYTHON SOURCE LINES 40-42 The recursion solves the Toeplitz system ---------------------------------------- .. GENERATED FROM PYTHON SOURCE LINES 42-48 .. code-block:: Python direct = linalg.solve(linalg.toeplitz(r[:10]), -r[1:11]) assert np.allclose(model.coefficients[1:], direct) print("order 4 from the recursion:", np.round(levinson_durbin(r, 4).coefficients, 3)) print("true: ", np.round(a_true, 3)) .. rst-class:: sphx-glr-script-out .. code-block:: none order 4 from the recursion: [ 1. -1.023 0.761 -0.626 0.646] true: [ 1. -1.019 0.763 -0.637 0.652] .. GENERATED FROM PYTHON SOURCE LINES 49-54 Prediction error against order ------------------------------ :math:`E_m = (1 - k_m^2)E_{m-1}` falls until the true order and then stays flat, because the extra reflection coefficients are nearly zero. .. GENERATED FROM PYTHON SOURCE LINES 54-66 .. code-block:: Python k = model.reflection_coefficients E = model.prediction_errors fig, axes = plt.subplots(1, 3, figsize=(13, 3.6)) axes[0].plot(np.arange(11), E / E[0], "o-") axes[0].axvline(4, color="0.6", ls="--") axes[0].set(xlabel="order m", ylabel="$E_m/E_0$", yscale="log", title="Prediction error") axes[1].stem(np.arange(1, 11), k, basefmt=" ") axes[1].axhspan(-1, 1, color="0.95") axes[1].set(xlabel="m", ylabel="$k_m$", ylim=(-1.1, 1.1), title="Reflection coefficients") print("reflection coefficients:", np.round(k, 3)) .. image-sg:: /api/gallery/signal_processing/spectral/images/sphx_glr_plot_08_levinson_durbin_001.png :alt: Prediction error, Reflection coefficients :srcset: /api/gallery/signal_processing/spectral/images/sphx_glr_plot_08_levinson_durbin_001.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none reflection coefficients: [-0.716 0.528 0.06 0.646 0.009 0.006 0.002 -0.01 -0.006 -0.008] .. GENERATED FROM PYTHON SOURCE LINES 67-72 Cost ---- The recursion takes about :math:`p^2` multiplications, Gaussian elimination about :math:`p^3/3`. .. GENERATED FROM PYTHON SOURCE LINES 72-84 .. code-block:: Python orders = np.array([4, 16, 64, 256]) axes[2].loglog(orders, orders**2, "o-", label="Levinson-Durbin, p²") axes[2].loglog(orders, orders**3 / 3, "s-", label="elimination, p³/3") axes[2].set(xlabel="order p", ylabel="multiplications", title="Cost") axes[2].legend(fontsize=8) fig.tight_layout() assert np.all(np.abs(k) < 1) assert np.abs(k[4:]).max() < 0.03 and abs(E[4] / E[-1] - 1) < 0.01 assert np.allclose(levinson_durbin(r, 4).coefficients, a_true, atol=0.02) plt.show() .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 0.296 seconds) .. _sphx_glr_download_api_gallery_signal_processing_spectral_plot_08_levinson_durbin.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: lite-badge .. image:: images/jupyterlite_badge_logo.svg :target: ../../../../lite/lab/index.html?path=api/gallery/signal_processing/spectral/plot_08_levinson_durbin.ipynb :alt: Launch JupyterLite :width: 150 px .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: plot_08_levinson_durbin.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: plot_08_levinson_durbin.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: plot_08_levinson_durbin.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_