From c2aa039fe0d0987888568fb03559d3c9e355aade Mon Sep 17 00:00:00 2001 From: Daniel Stephan Date: Fri, 8 Mar 2024 19:51:15 +0100 Subject: [PATCH] Additional documentation for cwt (issue #676) (#678) Closes gh-676 --- .../cwt_wavelet_frequency_bandwidth_demo.py | 60 +++++++++++++ doc/source/pyplots/plot_cwt_scaleogram.py | 62 ++++++++++++++ doc/source/pyplots/plot_wavelets.py | 28 +++++++ doc/source/ref/cwt.rst | 84 ++++++++++++++++++- doc/source/ref/wavelets.rst | 2 + 5 files changed, 233 insertions(+), 3 deletions(-) create mode 100644 doc/source/pyplots/cwt_wavelet_frequency_bandwidth_demo.py create mode 100644 doc/source/pyplots/plot_cwt_scaleogram.py create mode 100644 doc/source/pyplots/plot_wavelets.py diff --git a/doc/source/pyplots/cwt_wavelet_frequency_bandwidth_demo.py b/doc/source/pyplots/cwt_wavelet_frequency_bandwidth_demo.py new file mode 100644 index 000000000..3d44042d9 --- /dev/null +++ b/doc/source/pyplots/cwt_wavelet_frequency_bandwidth_demo.py @@ -0,0 +1,60 @@ +import matplotlib.pyplot as plt +import numpy as np + +import pywt + +# plot complex morlet wavelets with different center frequencies and bandwidths +wavelets = [f"cmor{x:.1f}-{y:.1f}" for x in [0.5, 1.5, 2.5] for y in [0.5, 1.0, 1.5]] +fig, axs = plt.subplots(3, 3, figsize=(10, 10), sharex=True, sharey=True) +for ax, wavelet in zip(axs.flatten(), wavelets): + [psi, x] = pywt.ContinuousWavelet(wavelet).wavefun(10) + ax.plot(x, np.real(psi), label="real") + ax.plot(x, np.imag(psi), label="imag") + ax.set_title(wavelet) + ax.set_xlim([-5, 5]) + ax.set_ylim([-0.8, 1]) +ax.legend() +plt.suptitle("Complex Morlet Wavelets with different center frequencies and bandwidths") +plt.show() + + +def gaussian(x, x0, sigma): + return np.exp(-np.power((x - x0) / sigma, 2.0) / 2.0) + + +def make_chirp(t, t0, a): + frequency = (a * (t + t0)) ** 2 + chirp = np.sin(2 * np.pi * frequency * t) + return chirp, frequency + + +def plot_wavelet(time, data, wavelet, title, ax): + widths = np.geomspace(1, 1024, num=75) + cwtmatr, freqs = pywt.cwt( + data, widths, wavelet, sampling_period=np.diff(time).mean() + ) + cwtmatr = np.abs(cwtmatr[:-1, :-1]) + pcm = ax.pcolormesh(time, freqs, cwtmatr) + ax.set_yscale("log") + ax.set_xlabel("Time (s)") + ax.set_ylabel("Frequency (Hz)") + ax.set_title(title) + plt.colorbar(pcm, ax=ax) + return ax + + +# generate signal +time = np.linspace(0, 1, 1000) +chirp1, frequency1 = make_chirp(time, 0.2, 9) +chirp2, frequency2 = make_chirp(time, 0.1, 5) +chirp = chirp1 + 0.6 * chirp2 +chirp *= gaussian(time, 0.5, 0.2) + +# perform CWT with different wavelets on same signal and plot results +wavelets = [f"cmor{x:.1f}-{y:.1f}" for x in [0.5, 1.5, 2.5] for y in [0.5, 1.0, 1.5]] +fig, axs = plt.subplots(3, 3, figsize=(10, 10), sharex=True) +for ax, wavelet in zip(axs.flatten(), wavelets): + plot_wavelet(time, chirp, wavelet, wavelet, ax) +plt.tight_layout(rect=[0, 0.03, 1, 0.95]) +plt.suptitle("Scaleograms of the same signal with different wavelets") +plt.show() diff --git a/doc/source/pyplots/plot_cwt_scaleogram.py b/doc/source/pyplots/plot_cwt_scaleogram.py new file mode 100644 index 000000000..01c925597 --- /dev/null +++ b/doc/source/pyplots/plot_cwt_scaleogram.py @@ -0,0 +1,62 @@ +import matplotlib.pyplot as plt +import numpy as np + +import pywt + + +def gaussian(x, x0, sigma): + return np.exp(-np.power((x - x0) / sigma, 2.0) / 2.0) + + +def make_chirp(t, t0, a): + frequency = (a * (t + t0)) ** 2 + chirp = np.sin(2 * np.pi * frequency * t) + return chirp, frequency + + +# generate signal +time = np.linspace(0, 1, 2000) +chirp1, frequency1 = make_chirp(time, 0.2, 9) +chirp2, frequency2 = make_chirp(time, 0.1, 5) +chirp = chirp1 + 0.6 * chirp2 +chirp *= gaussian(time, 0.5, 0.2) + +# plot signal +fig, axs = plt.subplots(2, 1, sharex=True) +axs[0].plot(time, chirp) +axs[1].plot(time, frequency1) +axs[1].plot(time, frequency2) +axs[1].set_yscale("log") +axs[1].set_xlabel("Time (s)") +axs[0].set_ylabel("Signal") +axs[1].set_ylabel("True frequency (Hz)") +plt.suptitle("Input signal") +plt.show() + +# perform CWT +wavelet = "cmor1.5-1.0" +# logarithmic scale for scales, as suggested by Torrence & Compo: +widths = np.geomspace(1, 1024, num=100) +sampling_period = np.diff(time).mean() +cwtmatr, freqs = pywt.cwt(chirp, widths, wavelet, sampling_period=sampling_period) +# absolute take absolute value of complex result +cwtmatr = np.abs(cwtmatr[:-1, :-1]) + +# plot result using matplotlib's pcolormesh (image with annoted axes) +fig, axs = plt.subplots(2, 1) +pcm = axs[0].pcolormesh(time, freqs, cwtmatr) +axs[0].set_yscale("log") +axs[0].set_xlabel("Time (s)") +axs[0].set_ylabel("Frequency (Hz)") +axs[0].set_title("Continuous Wavelet Transform (Scaleogram)") +fig.colorbar(pcm, ax=axs[0]) + +# plot fourier transform for comparison +from numpy.fft import rfft, rfftfreq + +yf = rfft(chirp) +xf = rfftfreq(len(chirp), sampling_period) +plt.semilogx(xf, np.abs(yf)) +axs[1].set_xlabel("Frequency (Hz)") +axs[1].set_title("Fourier Transform") +plt.tight_layout() diff --git a/doc/source/pyplots/plot_wavelets.py b/doc/source/pyplots/plot_wavelets.py new file mode 100644 index 000000000..5d057e6df --- /dev/null +++ b/doc/source/pyplots/plot_wavelets.py @@ -0,0 +1,28 @@ +import matplotlib.pyplot as plt +import numpy as np + +import pywt + +wavlist = pywt.wavelist(kind="continuous") +cols = 3 +rows = (len(wavlist) + cols - 1) // cols +fig, axs = plt.subplots(rows, cols, figsize=(10, 10), + sharex=True, sharey=True) +for ax, wavelet in zip(axs.flatten(), wavlist): + # A few wavelet families require parameters in the string name + if wavelet in ['cmor', 'shan']: + wavelet += '1-1' + elif wavelet == 'fbsp': + wavelet += '1-1.5-1.0' + + [psi, x] = pywt.ContinuousWavelet(wavelet).wavefun(10) + ax.plot(x, np.real(psi), label="real") + ax.plot(x, np.imag(psi), label="imag") + ax.set_title(wavelet) + ax.set_xlim([-5, 5]) + ax.set_ylim([-0.8, 1]) + +ax.legend(loc="upper right") +plt.suptitle("Available wavelets for CWT") +plt.tight_layout() +plt.show() diff --git a/doc/source/ref/cwt.rst b/doc/source/ref/cwt.rst index e48e3ae99..731d6705e 100644 --- a/doc/source/ref/cwt.rst +++ b/doc/source/ref/cwt.rst @@ -6,14 +6,82 @@ Continuous Wavelet Transform (CWT) ================================== -This section describes functions used to perform single continuous wavelet -transforms. +This section focuses on the one-dimensional Continuous Wavelet Transform. It +introduces the main function ``cwt`` alongside several helper function, and +also gives an overview over the available wavelets for this transfom. -Single level - ``cwt`` + +Introduction +------------ + +In simple terms, the Continuous Wavelet Transform is an analysis tool similar +to the Fourier Transform, in that it takes a time-domain signal and returns +the signal's components in the frequency domain. However, in contrast to the +Fourier Transform, the Continuous Wavelet Transform returns a two-dimensional +result, providing information in the frequency- as well as in time-domain. +Therefore, it is useful for periodic signals which change over time, such as +audio, seismic signals and many others (see below for examples). + +For more background and an in-depth guide to the application of the Continuous +Wavelet Transform, including topics such as statistical significance, the +following well-known article is highly recommended: + +`C. Torrence and G. Compo: "A Practical Guide to Wavelet Analysis", Bulletin of the American Meteorological Society, vol. 79, no. 1, pp. 61-78, January 1998 `_ + + +The ``cwt`` Function ---------------------- +This is the main function, which calculates the Continuous Wavelet Transform +of a one-dimensional signal. + .. autofunction:: cwt +A comprehensive example of the CWT +---------------------------------- + +Here is a simple end-to-end example of how to calculate the CWT of a simple +signal, and how to plot it using ``matplotlib``. + +First, we generate an artificial signal to be analyzed. We are +using the sum of two sine functions with increasing frequency, known as "chirp". +For reference, we also generate a plot of the signal and the two time-dependent +frequency components it contains. + +We then apply the Continuous Wavelet Transform +using a complex Morlet wavlet with a given center frequency and bandwidth +(namely ``cmor1.5-1.0``). We then plot the so-called "scaleogram", which is the +2D plot of the signal strength vs. time and frequency. + +.. plot:: pyplots/plot_cwt_scaleogram.py + +The Continuous Wavelet Transform can resolve the two frequency components clearly, +which is an obvious advantage over the Fourier Transform in this case. The scales +(widths) are given on a logarithmic scale in the example. The scales determine the +frequency resolution of the scaleogram. However, it is not straightforward to +convert them to frequencies, but luckily, ``cwt`` calculates the correct frequencies +for us. There are also helper functions, that perform this conversion in both ways. +For more information, see :ref:`Choosing scales` and :ref:`Converting frequency`. + +Also note, that the raw output of ``cwt`` is complex if a complex wavelet is used. +For visualization, it is therefore necessary to use the absolute value. + + +Wavelet bandwidth and center frequencies +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +This example shows how the Complex Morlet Wavelet can be configured for optimum +results using the ``center_frequency`` and ``bandwidth_frequency`` parameters, +which can simply be appended to the wavelet's string identifier ``cmor`` for +convenience. It also demonstrates the importance of choosing appropriate values +for the wavelet's center frequency and bandwidth. The right values will depend +on the signal being analyzed. As shown below, bad values may lead to poor +resolution or artifacts. + +.. plot:: pyplots/cwt_wavelet_frequency_bandwidth_demo.py +.. Sphinx seems to take a long time to generate this plot, even though the +.. corresponding script is relatively fast when run on its own. + Continuous Wavelet Families --------------------------- @@ -25,6 +93,12 @@ wavelet names compatible with ``cwt`` can be obtained by: wavlist = pywt.wavelist(kind='continuous') +Here is an overview of all available wavelets for ``cwt``. Note, that they can be +customized by passing parameters such as ``center_frequency`` and ``bandwidth_frequency`` +(see :ref:`ContinuousWavelet` for details). + +.. plot:: pyplots/plot_wavelets.py + Mexican Hat Wavelet ^^^^^^^^^^^^^^^^^^^ @@ -104,6 +178,8 @@ where :math:`M` is the spline order, :math:`B` is the bandwidth and :math:`C` is the center frequency. +.. _Choosing scales: + Choosing the scales for ``cwt`` ------------------------------- @@ -147,6 +223,8 @@ scales. The right column are the corresponding Fourier power spectra of each filter.. For scales 1 and 2 it can be seen that aliasing due to violation of the Nyquist limit occurs. +.. _Converting frequency: + Converting frequency to scale for ``cwt`` ----------------------------------------- diff --git a/doc/source/ref/wavelets.rst b/doc/source/ref/wavelets.rst index 7ff04a0ce..8fb692ae1 100644 --- a/doc/source/ref/wavelets.rst +++ b/doc/source/ref/wavelets.rst @@ -257,6 +257,8 @@ from plain Python lists of filter coefficients and a *filter bank-like* object. >>> myOtherWavelet = pywt.Wavelet(name="myHaarWavelet", filter_bank=filter_bank) +.. _ContinuousWavelet: + ``ContinuousWavelet`` object ----------------------------