Wavelets#
Overview#
Computes the wavelet of time-series data
Background
Wavelet Inputs
Mother Wavelets
Choosing Wavelet Scales
Wavelets Example in Python using PyWavelets
Power Spectrum
Phase Spectrum
Background#
Wavelets are a powerful tool for analyzing time-series data. Time-series data is data that depends on time as well as frequency. For example, modeling changes in temperature across atmospheric data over time or analyzing a short audio sample where a limited number of notes can be repeated in any order.
Fourier Transform is the most common form of signal processing. Fourier Transforms can return precise information about the frequencies present in a signal (known as the frequency domain). However, due to the Uncertainty Principle, it is impossible to know the exact frequency (frequency domain) and exact time (time domain) when a frequency occurs.
Wavelet transforms provide a potential solution and can be used to return an analysis with both time and frequency at a cost of reduced precision in frequency. Wavelets work by squishing and stretching a known wave (known as a Mother Wavelet) to match a range of possible frequencies over the short length of a signal. The modified squished/stretched wavelet (known as a Daughter Wavelet) will be shifted along the entire signal to compare the known wavelet frequency with the signal. Because the comparison occurs over a known time range with a known wavelet frequency, wavelet analysis can return data both about the frequency of the signal and the time in which it occurs.
Wavelet Inputs#
x: Input time-series data, for example: time and temperature in nino3_sst
wavelet: mother wavelet name
dt: sampling period (time between each y-value)
s0: smallest scale
dj: spacing between each discrete scales
jtot: largest scale
Mother Wavelets#
There are many different kinds of wavelets to choose from. The three most common (and used by NCL) are Morlet, Paul, and DOG. During wavelet analysis, a Mother wavelet is squished/stretched into many different Daughter wavelets based on the scales to fit different shaped waves. A Mother wavelet is chosen based on the characteristics of the data and the features that are being investigated
Some wavelets include: Mexican Hat, Morlet, Complex Morlet, Gaussian Derivative, Complex Gaussian Derivative, Shannon, etc..
# Available Mother Wavelets within the PyWavelets package
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()
---------------------------------------------------------------------------
ModuleNotFoundError Traceback (most recent call last)
Cell In[1], line 6
3 import matplotlib.pyplot as plt
4 import numpy as np
----> 6 import pywt
8 wavlist = pywt.wavelist(kind="continuous")
9 cols = 3
ModuleNotFoundError: No module named 'pywt'
Choosing Wavelet Scales#
Scales determine the range that the wavelet analysis is senstitive to where:
Large scales correspond to stretched wavelets that are sensitive to low frequencies
Small scales correspond to squished wavelets that are sensitive to high frequencies
The smallest scale is most commonly equal to:
s0 = 2 * sampling period
s0 = 2 * dt
Note: NCL sets smallest scale (s0) to sampling period when using Morlet wavelet (s0=dt)
Note: NCL sets smallest scale (s0) to 1/4th of the sampling period when using Paul wavelet (s0=dt/4)
The largest scale is most commonly equal to:
jtot = 1 + floattointeger(((log10(N*dt/s0))/dj)/log10(2.))
jtot = 1 + int(((log10(length of input data * sampling period/smallest scale))/spacing between each discrete scales)/log10(2))
The spacing between scales is commonly equal to 0.25 (dj=0.25) where a smallest spacing will produce a greater resolution of data but will run slower
Wavelets Example in Python using PyWavelets#
PyWavelets: pywt.cwt#
Wavelet analysis is accomplished in Python via external package. PyWavelets is an open source Python package for applying wavelet analysis.
coeffs, freqs = pywt.cwt(data, scales, wavelet, sampling_period)
Input Values#
data: input data as a array_like
scales: array_like collection of the scales to use (np.arange(s0, jtot, dj))
wavelet: name of Mother wavelet
sampling_period: optional sampling period for frequencies output
Return Value#
coefs: array_like collection of complex number output for wavelet coefficients
freqs: array_like collection of frequencies
# Download nino3 data
import geocat.datafiles as gcd
nino3_data = gcd.get('ascii_files/sst_nino3.dat')
sst_data = np.loadtxt(nino3_data)
PyWavelets: Morlet Wavelet#
The Morlet wavelet is a complex wavelet, with both a real and imaginary component. In PyWavelets this is known as a Complex Morlet ("cmor"). PyWavelets functions are derived from “Computational Signal Processing with Wavelets” (Teolis) to be compatible with how Matlab is defined.
import numpy as np # access to complex numbers (real and imaginary)
import matplotlib.pyplot as plt # plot data
import pywt # PyWavelets
Define the Complex Morlet#
PyWavelets defines a Complex Morlet Mother wavelet as cmorB-C where B is the bandwidth and C is the center frequency. For example:
wavelet_mother = "cmor1.5-1"
# PyWavelet Input Values
wavelet_mother = "cmor1.5-1"
dt = 0.25 # sampling period (time between each y-value)
s0 = 0.25 # smallest scale
dj = 0.25 # spacing between each discrete scales
jtot = 44 # largest scale
scales = np.arange(1, jtot + 1, dj)
# PyWavelets
wavelet_coeffs, freqs = pywt.cwt(
data=sst_data, scales=scales, wavelet=wavelet_mother, sampling_period=dt
)
Power Spectrum#
wavelet_coeffs is a complex number with a real and an imaginary number (1 + 2i). The power spectrum plots the real component of the complex number
real_component = abs(wavelet_coeffs)
The real component represents the magntiude of the wavelet coefficient displayed as the absolute value of the coefficients squared
# compare the power spectrum (absolute value squared)
power = np.power((abs(wavelet_coeffs)), 2)
fig, ax = plt.subplots(figsize=(10, 10))
# Convert Y-Axis from default to symmetrical log (symlog) with labels
ax.set_yscale("symlog")
ax.invert_yaxis()
ax.set_yticks([10, 20, 30, 40, 50])
ax.set_yticklabels([10, 20, 30, 40, 50])
# Plot Scalogram
plt.imshow(power, vmax=(power).max(), vmin=(power).min(), aspect="auto")
# Convert default X-axis from time steps of 0-504 (0-len(sst_data)) to Years
start_year = 1871
end_year = 1871 + (len(sst_data) * dt)
x_tickrange = np.arange(start_year, end_year, dt)
start = int(9 / dt) # 36, starts the x-axis label at 1880 (9 years after start of data)
display_nth = int(20 / dt) # 80, display x-axis label every 20 years
plt.xticks(range(len(x_tickrange))[start::display_nth], x_tickrange[start::display_nth])
plt.title("Wavelet Power Spectrum")
plt.xlabel("Year")
plt.ylabel("Scale")
plt.colorbar()
plt.show()
Phase Spectrum#
wavelet_coeffs is a complex number with a real and an imaginary number (1 + 2i). While less commonly used in frequency and signal analysis, the phase spectrum plots the imaginary component of the complex number
import numpy as np
# note: np.angle() returns the angle of a complex number
imaginary_component = np.angle(wavelet_coeffs)
The imaginary component represents the direction of the wavelet coefficient
# compare the phase spectrum
phase = np.angle(wavelet_coeffs)
fig, ax = plt.subplots(figsize=(10, 10))
# Convert Y-Axis from default to symmetrical log (symlog) with labels
ax.set_yscale("symlog")
ax.invert_yaxis()
ax.set_yticks([10, 20, 30, 40, 50])
ax.set_yticklabels([10, 20, 30, 40, 50])
# Plot scalogram
plt.imshow(phase, vmax=(phase).max(), vmin=(phase).min(), aspect="auto")
# Convert default X-axis from time steps of 0-504 (0-len(sst_data)) to Years
start_year = 1871
end_year = 1871 + (len(sst_data) * dt)
x_tickrange = np.arange(start_year, end_year, dt)
start = int(9 / dt) # 36, starts the x-axis label at 1880 (9 years after start of data)
display_nth = int(20 / dt) # 80, display x-axis label every 20 years
plt.xticks(range(len(x_tickrange))[start::display_nth], x_tickrange[start::display_nth])
plt.title("Wavelet Phase Spectrum")
plt.xlabel("Year")
plt.ylabel("Scale")
plt.colorbar()
plt.show()
