.. 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_02_wiener_khinchin.py"
.. LINE NUMBERS ARE GIVEN BELOW.

.. only:: html

    .. note::
        :class: sphx-glr-download-link-note

        :ref:`Go to the end <sphx_glr_download_api_gallery_signal_processing_spectral_plot_02_wiener_khinchin.py>`
        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_02_wiener_khinchin.py:


The Wiener-Khinchin theorem (1934)
==================================

Norbert Wiener (1930) and Aleksandr Khinchin (1934) showed that the power
spectrum of a stationary random process is the Fourier transform of its
autocorrelation function. A random signal has no Fourier transform of its
own, but its correlation does, and that gives its spectrum a precise
meaning. For a finite record the discrete version is exact: the DFT of the
biased sample autocorrelation equals the periodogram
(:func:`~engineeringkit.signal_processing.correlation`,
:func:`~engineeringkit.signal_processing.power_spectral_density`).

.. GENERATED FROM PYTHON SOURCE LINES 16-21

.. code-block:: Python

    import matplotlib.pyplot as plt
    import numpy as np

    from engineeringkit.signal_processing import correlation, dft, filter_signal, power_spectral_density








.. GENERATED FROM PYTHON SOURCE LINES 22-28

An exact identity for a finite record
-------------------------------------

``correlation(x, x)`` returns lags :math:`-(N-1)` to :math:`N-1`. Rotating
lag zero to the front and taking a DFT of length :math:`2N-1` gives
:math:`|X|^2`, the squared DFT of ``x`` zero-padded to the same length.

.. GENERATED FROM PYTHON SOURCE LINES 28-38

.. code-block:: Python


    rng = np.random.default_rng(1934)
    x = rng.normal(size=256)
    n = x.size
    r = correlation(x, x) / n  # biased autocorrelation estimate
    from_correlation = dft(np.roll(r, -(n - 1))).real
    from_signal = abs(dft(np.concatenate([x, np.zeros(n - 1)]))) ** 2 / n
    print(f"largest difference: {np.max(abs(from_correlation - from_signal)):.2e}")
    assert np.allclose(from_correlation, from_signal)





.. rst-class:: sphx-glr-script-out

 .. code-block:: none

    largest difference: 3.11e-15




.. GENERATED FROM PYTHON SOURCE LINES 39-47

A coloured process: theory against estimate
-------------------------------------------

White noise through :math:`y_n = a y_{n-1} + e_n` (an AR(1) process) has
autocorrelation :math:`r_k = \sigma^2 a^{|k|}/(1-a^2)`, whose transform is
the spectrum :math:`S(f) = \sigma^2/|1 - a e^{-j2\pi f/f_s}|^2` per unit
of :math:`f_s`. The sample autocorrelation of a long record follows
:math:`a^{|k|}`, and its transform follows :math:`S(f)`.

.. GENERATED FROM PYTHON SOURCE LINES 47-76

.. code-block:: Python


    a, fs = 0.8, 1.0
    y = filter_signal([1.0], [1.0, -a], rng.normal(size=2**16))
    lags = np.arange(-30, 31)
    r_y = correlation(y, y)[y.size - 1 + lags] / y.size
    r_theory = a ** abs(lags) / (1 - a**2)
    assert np.allclose(r_y, r_theory, atol=0.1)

    f, psd = power_spectral_density(y, fs)
    theory = 2 / (fs * abs(1 - a * np.exp(-2j * np.pi * f / fs)) ** 2)  # one-sided
    smoothed = np.convolve(psd, np.ones(1024) / 1024, "same")
    band = (f > 0.05) & (f < 0.45)
    ratio = smoothed[band] / theory[band]
    print(f"smoothed periodogram / theory over the band: {ratio.min():.3f} to {ratio.max():.3f}")
    assert np.all(abs(ratio - 1) < 0.15)

    fig, axes = plt.subplots(1, 2, figsize=(10, 3.8))
    axes[0].stem(lags, r_y, basefmt=" ", label="sample")
    axes[0].plot(lags, r_theory, "r", label=r"$a^{|k|}/(1-a^2)$")
    axes[0].set(xlabel="lag $k$", ylabel="autocorrelation", title="Correlation")
    axes[0].legend()
    axes[1].semilogy(f, psd, lw=0.3, color="0.7", label="periodogram")
    axes[1].semilogy(f, smoothed, label="smoothed")
    axes[1].semilogy(f, theory, "r--", label="transform of $r_k$")
    axes[1].set(xlabel="frequency [cycles/sample]", ylabel="PSD", title="Spectrum")
    axes[1].legend()
    fig.tight_layout()

    plt.show()



.. image-sg:: /api/gallery/signal_processing/spectral/images/sphx_glr_plot_02_wiener_khinchin_001.png
   :alt: Correlation, Spectrum
   :srcset: /api/gallery/signal_processing/spectral/images/sphx_glr_plot_02_wiener_khinchin_001.png
   :class: sphx-glr-single-img


.. rst-class:: sphx-glr-script-out

 .. code-block:: none

    smoothed periodogram / theory over the band: 0.922 to 1.096





.. rst-class:: sphx-glr-timing

   **Total running time of the script:** (0 minutes 0.118 seconds)


.. _sphx_glr_download_api_gallery_signal_processing_spectral_plot_02_wiener_khinchin.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_02_wiener_khinchin.ipynb
        :alt: Launch JupyterLite
        :width: 150 px

    .. container:: sphx-glr-download sphx-glr-download-jupyter

      :download:`Download Jupyter notebook: plot_02_wiener_khinchin.ipynb <plot_02_wiener_khinchin.ipynb>`

    .. container:: sphx-glr-download sphx-glr-download-python

      :download:`Download Python source code: plot_02_wiener_khinchin.py <plot_02_wiener_khinchin.py>`

    .. container:: sphx-glr-download sphx-glr-download-zip

      :download:`Download zipped: plot_02_wiener_khinchin.zip <plot_02_wiener_khinchin.zip>`


.. only:: html

 .. rst-class:: sphx-glr-signature

    `Gallery generated by Sphinx-Gallery <https://sphinx-gallery.github.io>`_