Wavelet analysis. Part 2

Introduction

This publication deals with wavelet - time series analysis. The basic idea of ​​the wavelet transform meets the specifics of many time series, demonstrating the evolution in time of its main characteristics - the average value, dispersion, periods, amplitudes and phases of harmonic components. The overwhelming majority of the processes studied in various fields of knowledge have the above listed features.

The purpose of this publication is to describe the continuous wavelet transform time series technique using the PyWavelets library.

A bit of history
')
Geophysicist D. Morlet in the late 70s of XX century. faced the problem of analyzing signals from seismic sensors, which contained a high-frequency component (seismic activity) for a short period of time and low-frequency components (a quiet state of the earth's crust) - for a long period. The window Fourier transform allows you to analyze either the high-frequency component or the low-frequency component, but not both components at once.

Therefore, an analysis method was proposed in which the width of the window function increased for low frequencies and decreased for high frequencies. The new window transformation was obtained as a result of stretching (compressing) and shifting in time one generating (the so-called scaling function - scaling function, scalet) function. This generating function was called D. Morlet wavelet.

Continuous Wavelet Transform (CWT)

Single-level wavelet transform:

coefs, frequencies = pywt.cwt (data, scales, wavelet)

Where:

data: (array_like) -Input signal;

scales (array_like) : - The scale of the wavelet to use. You can use the ratio f = scale2frequency (scale, wavelet) / sampling_period and determine which physical frequency f. Here, f is in hertz when sampling_period in seconds;

wavelet (Wavelet object or name) : - Wavelet to use;

sampling_period (float): - The sampling period for output frequencies (optional). The values ​​calculated for the coefs do not depend on the choice of sampling_period (i.e. scales are not scaled by the period sample.);

coefs (array_like) : - Continuous wavelet - conversion of the input signal for a given scale and wavelet;

frequencies (array_like) : - If the unit of the sampling period is a second and is specified, then the frequencies are output in hertz. Otherwise, a sample period of 1 is assumed.

Note: The size of the coefficient arrays depends on the length of the input array and the lengths of the specified scales.

Get the list of available wavelet names compatible with cwt:

 >>> import pywt >>> wavlist = pywt.wavelist(kind='continuous') >>> wavlist 

['cgau1', 'cgau2', 'cgau3', 'cgau4', 'cgau5', 'cgau6', 'cgau7', 'cgau8', 'cmor', 'fbsp', 'gaus1', 'gaus2', ' gaus3 ',' gaus4 ',' gaus5 ',' gaus6 ',' gaus7 ',' gaus8 ',' mexh ',' morl ',' shan ']

For wavelet analysis, the choice of maternal wavelet is one of the key tasks. Therefore, before each choice of the wavelet it is simply necessary to use the following program, which allows to determine the dynamic properties of the wavelet:


Further we will consider the wavelet functions and their properties using the above program:

Mexican hat wavelet "mexh":

$$$\ varphi \ left (t \ right) = \ frac {2} {\ sqrt {3} \ sqrt [4] {\ pi}} \ cdot exp ^ {- \ frac {t ^ {2}} {2} } \ cdot \ left (1-t ^ {2} \ right)$



Morle Morlet Wavelet:

$$$\ varphi \ left (t \ right) = exp ^ {- \ frac {t ^ {2}} {2}} \ cdot cos (5t)$

Complex Morlet wavelet cmorB-C

$$$\ varphi \ left (t \ right) = \ frac {1} {\ sqrt {\ pi} B} exp ^ {- \ frac {t ^ {2}} {B}} \ cdot exp ^ {j2 \ pi C t}$

where B is the throughput; C is the center frequency.

Hereinafter (without special instructions) the values ​​of B, C are set with a floating point.



Gaussian wavelets "gausP"

$$$\ varphi \ left (t \ right) = C \ cdot exp ^ {- t ^ {2}}$

where: P is an integer from 1 to 8 is inserted into the name of the wavelet, for example "gaus8"; C is a built-in normalization constant.



Complex Gaussian wavelets "cgauP"

$$$\ varphi \ left (t \ right) = C \ cdot exp ^ {- t ^ {2}} \ cdot exp ^ {- jt}$

where: P - an integer from 1 to 8 is inserted into the name of the wavelet, for example, “gaus8” and correspond to the order of the derivative of the wavelet function; C is a built-in normalization constant.



Shannon wavelets "shanB-C"

$$$\ varphi \ left (t \ right) = \ sqrt {B} \ cdot \ frac {sin (\ pi B t)} {\ pi B T} \ cdot exp ^ {j 2 pi C t}$

where: B is the width of the strip; C is the center frequency.



CWT in PyWavelets is applied to discrete data - convolutions with samples of the wavelet integral. If scale is too small, then we get a discrete filter with inadequate sampling, which leads to smoothing, as shown in the graph for the 'cmor1.5-1.0' wavelet.

The left column on the graphs shows the discrete filters used in convolution with different scales. The right column is the corresponding Fourier power spectra of each filter. With scales 1 and 2 for the 'cmor1.5-1.0' graph, it can be seen that anti-aliasing occurs due to violation of the Nyquist constraint. This phenomenon may also occur in other wavelets when choosing C and B, therefore, when working with wavelets, you must use the program listing 1.

To associate a given scale with a specific signal frequency, the signal sampling period must be known. Using the pywt.scale2frequency () function, you can convert the list of scales to the corresponding frequencies. The correct choice of scales depends on the selected wavelet, so pywt.scale2frequency () should be used to get an idea of ​​the corresponding frequency range of the signal.

 >>> import pywt >>> dt = 0.01 # 100 Hz sampling >>> frequencies = pywt.scale2frequency('cmor1.5-1.0', [1, 2, 3, 4]) / dt >>> frequencies 

array ([100., 50., 33.33333333, 25.])

Wavelet - time series analysis using CWT PyWavelets

The el-Nino dataset is a time series dataset used to track El Nino and contains quarterly measurements of sea surface temperature from 1871 to 1997. To understand the power of a scalegram, let's visualize it for the el-Nino dataset along with the original data of the time series and its Fourier transform.

For the wavelet analysis of time series, you must perform the following steps for the items:

1. Choose maternal wavelet: Select the complex Morlet wavelet “cmorB-C”:

$$$\ varphi \ left (t \ right) = \ frac {1} {\ sqrt {\ pi} B} exp ^ {- \ frac {t ^ {2}} {B}} \ cdot exp ^ {j2 \ pi C t}$

Bandwidth - B and center frequency - From which we will determine in the next stage.

2. Determine the center frequency, taking dt = 0.25 for our time series:

 import pywt dt = 0.25 # 4 Hz sampling scale=range(1,4) frequencies = pywt.scale2frequency('cmor1.0-0.5', scale) / dt print(frequencies) 

[2. 1. 0.66666667]

3. Find the Fourier transform of the mother wavelet cmor1.0-0.5 (we use Listing 1):
Continuous wavelet will be evaluated in the whole range [-8.0, 8.0]



4. Select the data for the time series:

 pd.read_csv("http://paos.colorado.edu/research/wavelets/wave_idl/sst_nino3.dat", sep = "|") 

The data is a quarterly measurement of the temperature of the sea surface from 1871 to 1997. For this data: t0 = 1871 dt = 0.25

5. Build a time series with a signal and its moving average on one chart:




6. We carry out the Fourier transform and the modality of the spectrum from the time series:




7. Determine the scale: scales = arange (1, 128); levels = [2 ** - 4, 2 ** - 3, 2 ** - 2, 2 ** - 1, 2 ** 0, 2 ** 1, 2 ** 2, 2 ** 3].

8. Build a scalegram:




The scale diagram shows that most of the spectrum power is concentrated over a period of 2-8 years, this corresponds to a frequency of 0.125 - 0.5 Hz. In the Fourier spectrum, fashion also concentrates around these frequency values. However, the wavelet transform also gives us temporal information, but the Fourier transform does not.

For example, in the scale diagram we see that before 1920 there were many fluctuations, while between 1960 and 1990 there were not so many. We can also see that over time there is a shift from shorter periods to longer ones.

Using the scaleogram module

Scaleogram is a handy tool for analyzing 1D data with a continuous wavelet transform built on the basis of the PyWavelets library. The module is designed to provide a reliable tool for quick data analysis.

Scaleogram has the following features:

• simple call sign for beginners
• readable axes and matplotlib net integration
• Many options for zooming in, filtering the spectrum, embedding colorbar, etc ...
• support of frequency and frequency units in accordance with the marking
• speed, uses N * log (N) algorithm for conversions
• portability: tested with python2.7 and python3.7
• detailed error messages and documentation with examples

However, in the usage examples, the data access is not correctly recorded:

 nino3 = "http://paos.colorado.edu/research/wavelets/wave_idl/sst_nino3.dat" 

The following warning appears:

nino3 = pd.read_table (nino3)
which leads to a warning:
Warning (from warnings module):
File "C: /Users/User/Desktop/wavelet/pr.py", line 7
nino3 = pd.read_table (nino3)
FutureWarning: read_table is deprecated, use read_csv instead, passing sep = '\ t'

Access to data should be written like this:

 url="http://paos.colorado.edu/research/wavelets/wave_idl/sst_nino3.dat" nino3 = pd.read_csv(url, sep = "|") 

To show the use of the scalar and compare the results with the results of the previous example, let us use the same data - for the quarterly measurements of the sea surface temperature from 1871 to 1997. For this data: t0 = 1871 dt = 0.25



If we compare the scalegram with the received scangram, then, with the exception of the color palette, they are identical, and therefore show the same dynamics of the time series.
Listing 4 is simpler than Listing 3 and is faster.

When the frequency spectrum and the Fourier power spectrum do not provide additional information on the dynamics of the time series, then the use of a scalar is preferable.

findings

An example of a wavelet analysis of a time series using the PyWavelets library is given.