.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "gallery/frequency_change_detection.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 .. rst-class:: sphx-glr-example-title .. _sphx_glr_gallery_frequency_change_detection.py: Frequency Change Detection =============================== This example considers a signal which is a sinusoid whose frequency changes abruptly at a specific point in time. Looking at the waveform of the signal, it is very difficult to identify the time where the frequency change occurs. However, when we perform the wavelet decomposition of the signal, we can clearly see the location at which the discontinuity in frequency occurs. We will do the following: - Synthetically construct a sinusoid with two different frequencies at different times. - Examine its Fourier spectrum which cannot reveal this discontinuity. - Perform multilevel wavelet decomposition. - Visually inspect the detail coefficients at multiple levels to locate the discontinuity - Algorithmically locate the discontinuity by scanning the detail coefficients. .. GENERATED FROM PYTHON SOURCE LINES 24-28 .. code-block:: default # Configure JAX to work with 64-bit floating point precision. from jax.config import config config.update("jax_enable_x64", True) .. GENERATED FROM PYTHON SOURCE LINES 29-30 Let's import necessary libraries .. GENERATED FROM PYTHON SOURCE LINES 30-39 .. code-block:: default import jax.numpy as jnp # CR-Suite libraries from cr.nimble.dsp import hard_threshold_sorted import cr.wavelets as wt # Utility functions to construct sinusoids import cr.nimble.dsp.signals as signals # Plotting import matplotlib.pyplot as plt .. GENERATED FROM PYTHON SOURCE LINES 40-42 Test signal generation ------------------------------ .. GENERATED FROM PYTHON SOURCE LINES 42-72 .. code-block:: default # Sampling frequency in Hz fs = 100 # Signal duration in seconds T = 20 # Frequency of first part of signal in Hz f = 1/8 # Start and end times of first part of signal start_time = 0 end_time = 12 # Construct the first part of signal t, x1 = signals.transient_sine_wave(fs, T, f, start_time, end_time) # Frequency of second part of signal in Hz f = 1/6 # Start and end times of first part of signal start_time = end_time end_time = T # Adjust the initial phase of second part of signal to be in continuity with the first part initial_phase=jnp.pi # Construct the second part of signal t, x2 = signals.transient_sine_wave(fs, T, f, start_time, end_time, initial_phase=initial_phase) # Combine the first and second parts of signal x = x1 + x2 # Overall signal length n = len(x) # Plot the parts and combined signal fig, axs = plt.subplots(3, figsize=(12,12)) axs[0].plot(t,x1) axs[1].plot(t,x2) axs[2].plot(t,x) .. image-sg:: /gallery/images/sphx_glr_frequency_change_detection_001.png :alt: frequency change detection :srcset: /gallery/images/sphx_glr_frequency_change_detection_001.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none [] .. GENERATED FROM PYTHON SOURCE LINES 73-76 The last plot shows the combined signal with discontinuity at t=12 sec when the frequency changes from 1/8 Hz to 1/6 Hz. The change is very difficult to notice as well as locate in the time domain plot of the signal. .. GENERATED FROM PYTHON SOURCE LINES 78-80 Frequency spectrum ------------------------------ .. GENERATED FROM PYTHON SOURCE LINES 80-86 .. code-block:: default # Compute the FFT f = jnp.fft.fftshift(jnp.fft.fft(x)) # Plot the magnitude fig, axs = plt.subplots(1, figsize=(12,4)) plt.plot(jnp.abs(f[750:1250])) .. image-sg:: /gallery/images/sphx_glr_frequency_change_detection_002.png :alt: frequency change detection :srcset: /gallery/images/sphx_glr_frequency_change_detection_002.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none [] .. GENERATED FROM PYTHON SOURCE LINES 87-90 The frequencies 1/6 and 1/8 Hz are so close that it is difficult to distinguish them in the frequency spectrum. Of course, location of the discontinuity cannot be identified in the frequency spectrum. .. GENERATED FROM PYTHON SOURCE LINES 93-95 Multilevel wavelet decomposition ---------------------------------- .. GENERATED FROM PYTHON SOURCE LINES 95-102 .. code-block:: default # Compute the multilevel wavelet decomposition coeffs = wt.wavedec(x, 'db4') # Total number of decomposition levels levels = len(coeffs) -1 print(levels) .. rst-class:: sphx-glr-script-out .. code-block:: none 8 .. GENERATED FROM PYTHON SOURCE LINES 103-104 First level detail coefficients .. GENERATED FROM PYTHON SOURCE LINES 104-107 .. code-block:: default cd1 = coeffs[-1] # Dyadic upsample them to align with the time values cd1 = wt.dyadup_in(cd1) .. GENERATED FROM PYTHON SOURCE LINES 108-109 Second level detail coefficients .. GENERATED FROM PYTHON SOURCE LINES 109-116 .. code-block:: default cd2 = coeffs[-2] # Dyadic upsample them to align with the time values cd2 = wt.dyadup_in(wt.dyadup_in(cd2)) # Plot the first and second level detail coefficients fig, axs = plt.subplots(2, figsize=(12,12)) axs[0].plot(t, cd1[:n]) axs[1].plot(t, cd2[:n]) .. image-sg:: /gallery/images/sphx_glr_frequency_change_detection_003.png :alt: frequency change detection :srcset: /gallery/images/sphx_glr_frequency_change_detection_003.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none [] .. GENERATED FROM PYTHON SOURCE LINES 117-122 The discontinuities are clearly visible in the plots of detail coefficients both at first level and second level. The detail coefficients have been aligned with the time values. Both plots should large coefficients around t=12 sec. Large coefficients at the boundary are due to boundary effects in DWT and can be safely ignored. .. GENERATED FROM PYTHON SOURCE LINES 124-126 Locate the indices of largest entries (by magnitude) in the detail coefficients at first level .. GENERATED FROM PYTHON SOURCE LINES 126-128 .. code-block:: default idx, values = hard_threshold_sorted(cd1, 8) print(idx) .. rst-class:: sphx-glr-script-out .. code-block:: none [ 0 2 4 1200 1202 1204 2002 2004] .. GENERATED FROM PYTHON SOURCE LINES 129-132 After ignoring the first couple of entries for the boundary effects at the beginning of the data, we can see that the discontinuity has been correctly identified at sample 1200 which happens to correspond to t=12 sec. .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 2.207 seconds) .. _sphx_glr_download_gallery_frequency_change_detection.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: frequency_change_detection.py ` .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: frequency_change_detection.ipynb ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_