.. 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_09_bogert_healy_tukey_cepstrum.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_09_bogert_healy_tukey_cepstrum.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_09_bogert_healy_tukey_cepstrum.py:


Bogert, Healy and Tukey's cepstrum (1963)
=========================================

Seismologists wanted to know how deep an earthquake was, from the delay
between a wave and its echo off the surface. An echo multiplies the
spectrum by :math:`1 + \alpha e^{-j\omega d}`: a ripple in the log
spectrum whose "frequency" is the delay. Bruce Bogert, Michael Healy and
John Tukey (1963) took the spectrum of the log spectrum and found the echo
as a peak. They reversed letters to name the new domain: *cepstrum*,
*quefrency*, *liftering*. Logarithms turn the convolution of a source and
a filter into a sum, which is how pitch is separated from the vocal tract
in speech (:func:`~engineeringkit.signal_processing.real_cepstrum`).

.. GENERATED FROM PYTHON SOURCE LINES 17-23

.. code-block:: Python

    import matplotlib.pyplot as plt
    import numpy as np
    from scipy import signal

    from engineeringkit.signal_processing import real_cepstrum








.. GENERATED FROM PYTHON SOURCE LINES 24-29

An echo becomes a peak
----------------------

A coloured pulse and its echo, delayed 37 ms with half the amplitude,
at 1 kHz.

.. GENERATED FROM PYTHON SOURCE LINES 29-49

.. code-block:: Python


    fs, delay, alpha = 1000.0, 37, 0.5
    rng = np.random.default_rng(1963)
    pulse = signal.lfilter([1.0], [1.0, -1.6, 0.8], rng.normal(size=60)) * np.hanning(60)
    x = np.zeros(1024)
    x[100:160] = pulse
    x[100 + delay : 160 + delay] += alpha * pulse
    c = real_cepstrum(x)
    quefrency = np.arange(c.size) / fs
    search = slice(10, 512)
    found = (np.argmax(c[search]) + 10) / fs
    print(f"cepstral peak at {found * 1e3:.0f} ms, height {c[delay]:.3f} (alpha/2 = {alpha / 2})")
    fig, axes = plt.subplots(1, 3, figsize=(13, 3.6))
    spectrum = np.abs(np.fft.rfft(x))
    axes[0].plot(np.fft.rfftfreq(x.size, 1 / fs), 20 * np.log10(spectrum))
    axes[0].set(xlabel="frequency [Hz]", ylabel="log magnitude [dB]", title="Rippled log spectrum")
    axes[1].plot(quefrency[1:200] * 1e3, c[1:200])
    axes[1].axvline(delay, color="C3", ls="--", lw=0.8)
    axes[1].set(xlabel="quefrency [ms]", title="Real cepstrum")




.. image-sg:: /api/gallery/signal_processing/spectral/images/sphx_glr_plot_09_bogert_healy_tukey_cepstrum_001.png
   :alt: Rippled log spectrum, Real cepstrum
   :srcset: /api/gallery/signal_processing/spectral/images/sphx_glr_plot_09_bogert_healy_tukey_cepstrum_001.png
   :class: sphx-glr-single-img


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

 .. code-block:: none

    cepstral peak at 37 ms, height 0.242 (alpha/2 = 0.25)

    [Text(0.5, 10.433333333333328, 'quefrency [ms]'), Text(0.5, 1.0, 'Real cepstrum')]



.. GENERATED FROM PYTHON SOURCE LINES 50-56

Pitch from a pulse train
------------------------

A crude voiced sound: pulses every 8 ms through a resonance. The cepstrum
separates the slowly varying resonance (low quefrency) from the 8 ms
pitch period (a peak at 8 ms).

.. GENERATED FROM PYTHON SOURCE LINES 56-71

.. code-block:: Python


    pitch = 8
    train = np.zeros(1024)
    train[::pitch] = 1.0
    voiced = signal.lfilter([1.0], [1.0, -1.3, 0.9], train) + 1e-3 * rng.normal(size=train.size)
    cv = real_cepstrum(voiced)
    pitch_found = np.argmax(cv[4:100]) + 4
    axes[2].plot(quefrency[1:100] * 1e3, cv[1:100])
    axes[2].set(xlabel="quefrency [ms]", title=f"Voiced sound: pitch {pitch_found} ms")
    fig.tight_layout()
    assert np.isclose(found * fs, delay)
    assert np.isclose(c[delay], alpha / 2, atol=0.05)
    assert pitch_found == pitch

    plt.show()








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

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


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

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

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

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

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

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

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


.. only:: html

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

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