Additional documentation for cwt (issue #676) (#678)

Closes gh-676

Author: drs251

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()

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()

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()

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 <https://paos.colorado.edu/research/wavelets/bams_79_01_0061.pdf>`_
+
+
+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``
 -----------------------------------------
 

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
 ----------------------------