The post is a brief introduction and overview of wavelets. Wavelets are a powerful tool for time series analysis, de-noising, and data compression, but have recently exploded in fields like hydrology and geophysics. We are often interested in discovering periodic phenomena or underlying frequencies that are characteristic in a signal of interest, such as low-frequency teleconnections. We may also want to better understand very rapidly changing seismic signals during an earthquake or ocean wave dispersion. Wavelets can help us answer many questions in these applications!

### Basic Background

The typical approach to understanding what frequencies are present in a signal involves transforming the time-domain signal into the frequency domain by means of the Fourier Transform. A well-known deficiency of the Fourier Transform is that is provides information about the frequencies that are present, but all information about the time at which these frequencies are present is lost. Thus, signals that are non-stationary, or exhibit changing frequencies over time cannot be analyzed by the Fourier Transform (Deshpande, 2015). One solution that was proposed to address this limitation was the Short-time Fourier Transform (STFT). In a very basic sense, the STFT involves applying a window function to segment the signal in time and then performing a Fourier transform on each segment. Then, the frequencies for each segment can be plotted to better understand the changing spectra (Smith, 2011). A considerable amount of research was spent on determining appropriate windowing functions between the 1940s and 1970s. However, limitations still existed with this approach. The STFT utilizes the same windowing function across the whole time series, which may not be appropriate for all the different frequencies that may be characterize a signal. When a signal has very high frequency content, a narrow window is preferable, but this results in poor frequency resolution. When a signal has lower frequency content, a wider window is preferable, but this results in poor time resolution. This tradeoff is often termed an example of the Heisenberg Uncertainty Principle (Marković et al, 2012).

In order to properly analyze signals that are non-stationary or transient (most signals of interest) and to appropriately address both large and small scale features in the frequency space, we can make use of wavelets! Wavelet transforms are especially useful for analyzing signals which have low frequencies for long durations and short frequencies for short durations. Furthermore, the basis functions are not restricted to only sinusoids, which can make it easier to approximate functions with sharper peaks or corners.

In a continuous wavelet transform, expressed below in Equation (1), the basis function or window is termed the “mother” wavelet, which is designated by Ψ .

This mother wavelet can be scaled and translated with the s and τ parameters, respectively, to produce “daughter” or “child” wavelets, which are simply variations of the mother wavelet. The wavelets are not only translated across the signal but can be dilated or contracted to capture longer or shorter term events.

By definition, a continuous wavelet transform is a convolution of the input signal with the daughter wavelets, which is inherently a measure of similarity. The wavelet is “slid” over the signal of interest and similarities at each time point are measured. Therefore, the wavelet will be correlated with the parts of the series that contain a similar frequency.

The convolution with the selection of wavelets can lead to redundancy of information when the basis functions are not orthogonal. A family of wavelet basis functions have been developed to address this limitation. The simplest orthonormal wavelet is the Haar wavelet, which is a discrete, square shaped wavelet. The Haar wavelet is not continuous or differentiable, however it is particularly useful for approximating the response of systems that may experience a sudden transition. A Morlet wavelet is a complex sinusoid enclosed in a Gaussian envelope and may be more useful in applications such as music, hearing/vision, and for analyzing electrocardiograms.

In a discrete wavelet transform (DWT), the translation and scale parameters, s and τ are discretized in such a way that each successive wavelet is twice the dimension as the one before, to cover all but very low frequencies. This is termed a dyadic filter bank. Each daughter wavelet can therefore be thought of as a bandpass filter that represents a specific frequency of interest or scale.

The correlation across time associated with the signal and each daughter wavelet can be plotted in a scalogram as shown below.

The coefficients of the wavelet indicate the correlation between the wavelet and the signal of interest at a specific time and scale. The amplitude squared of the wavelet coefficient,|W_{i}|^{2} defines the wavelet power and can be used to create a Wavelet Power Spectrum. A larger power corresponds to a higher correlation, therefore, the regions of high power in the spectrum correspond to areas of interest^{5} .

### WaveletComp

We will use the library WaveletComp in R to demonstrate how to find the Wavelet Power Spectrum. Many wavelet libraries exist in Python and MATLAB work equivalently and just as simply, but I like WaveletComp due the huge supplement in the package repository that contains many examples on how to use the functions. For this post, I took quarterly El Niño 3 Region (NINO3) SST surface anomalies from NOAA recorded over 1871-1996 and applied a single function from the package to create the spectrum.

library(waveletComp) my.w = analyze.wavelet(mydata, "Index", loess.span = 0, dt = 0.25, dj = 1/250, lowerPeriod = 0.25, upperPeriod = 32, make.pval = TRUE, n.sim = 30) wt.image(my.w, n.levels = 250, legend.params = list(lab = "wavelet power levels") )

Here we specify the name of the dataframe that contains the data (mydata), the column that contains the index (“Index”), the number of observations per unit time (dt=0.25 due to the quarterly resolution) and the upper and lower period that bounds the search area. The range of periods is generally represented as a series of octaves of 2^{j }and each octave is divided into 250 sub-octaves based on the *dj* term. The *“make.pvalue=TRUE” *argument draws a white line around the areas that are deemed significant. Finally, we plot the wavelet objecive using *wt.image*.

As seen in the power spectrum, the highest powers across the time series are recorded within the 2-7 year frequency bands, which matches with the period that El Niño events tend to occur. Interestingly, the El Niño signal seems strongest at the earliest and later parts of the decade and is notably less prominent between the years of 1920-1960, which has been observed in other work (Torrence & Webster, 1997). The wavelet spectrum allows us to see how the periodicity of the El Niño signal has changed over time, with longer periods being observed around 1920 and shorter periods closer to the beginning and end of the century.

Because wavelets are shifted across each time-point in the convolution operation, the coefficients at both edges (half of the wavelet duration at each frequency) are not accurate. The smaller frequencies utilize smaller wavelet durations, creating a “cone of influence” that is shown by the white shading on the edges of the plot. The cone of influence designates the areas of the plot that should be disregarded.

Coherence is a measure of the common oscillations that two signals share. Generally, coherence or cross-correlation is used to assess similarity in the time or frequency domain. However, if the two signals being compared are non-stationary, the correlation can change over time. Therefore, the coherence must be represented in a way to show changes across frequency and time. We can create cross-correlation plots in WaveletComp as well. I have extracted some data supplied by MathWorks which contains a monthly Niño3 index along with an average All-India Rainfall Index^{6}. In order to assess the how these two time series are linked, we use the analyze.coherency function in WaveletComp.

my.wc = analyze.coherency(data, my.pair = c("Nino_Index", "Rainfall"), loess.span = 0, dt = 1/12, dj = 1/50, lowerPeriod = 0.25, upperPeriod = 32, make.pval = TRUE, n.sim = 10) wc.image(my.wc, n.levels = 250, legend.params = list(lab = "cross-wavelet power levels"), color.key = "interval",show.date=TRUE)

The maximum correlation aligns well within the 2-7 year bands that were observed in the above plot. The orientation of the arrows in the figure signify the delay between the two signals at that time period, The vertical to horizontal arrows denote a ¼ – ½ cycle delay within the significant areas or ½-3.5 years of delay between the El Niño SST’s observed off the coast of South America to influence rainfall in the Indian subcontinent. Pretty cool!

There is quite a bit of solid math and proofs that one should go through to truly understand the theory behind wavelets, but hopefully this post serves as a good introduction to show how wavelets can be useful for your respective analysis. When I first learned about wavelets and tried out WaveletComp, I immediately began conducting a wavelet analysis on every time series that I had on hand to look for hidden frequencies! I hope that this tutorial motivates you to explore wavelets for analyzing your non-stationary signals.

### Sources:

Baker, J. W. “Quantitative Classification of Near-Fault Ground Motions Using Wavelet Analysis.” *Bulletin of the Seismological Society of America*, vol. 97, no. 5, 2007, pp. 1486–1501., doi:10.1785/0120060255.

Deshpande, Jaidev. “Limitations of the Fourier Transform: Need For a Data Driven Approach.” *Limitations of the Fourier Transform: Need For a Data Driven Approach – Pyhht 0.0.1 Documentation*, 2015, pyhht.readthedocs.io/en/latest/tutorials/limitations_fourier.html.

Marković, Dejan, et al. “Time-Frequency Analysis: FFT and Wavelets.” *DSP Architecture Design Essentials*, 2012, pp. 145–170., doi:10.1007/978-1-4419-9660-2_8.

Smith, Julius O. *Spectral Audio Signal Processing*. W3K, 2011.

Torrence, Christopher, and Peter J. Webster. “Interdecadal Changes in the ENSO–Monsoon System.” *Journal of Climate*, vol. 12, no. 8, 1999, pp. 2679–2690., doi:10.1175/1520-0442(1999)0122.0.co;2.

[5] “Wavelet Power Spectrum.” *Wavelet Toolkit Architecture*, Dartmouth, 16 June 2005, northstar-www.dartmouth.edu/doc/idl/html_6.2/Wavelet_Power_Spectrum.html.

[6] “Wcoherence.” *MATLAB & Simulink Example*, The MathWorks, Inc., 2020, http://www.mathworks.com/help/wavelet/examples/compare-time-frequency-content-in-signals-with-wavelet-coherence.html.

Hi Rohini,

Thank you very much for providing this well-explained wavelet introduction post. Pretty cool! Recently I also started to learn wavelet method. I have a question about wavelet coherence. Using your analysis (“Nino_Index” & “Rainfall”) as an example, because of their cross-correlation over time, are we able to extract an estimated rainfall signal series that is contributed by Nino_Index? I’m not expecting a nicely fitted rainfall series but expecting this extracted series can explain 10-30% variance of the observed rainfall series (the explanatory power depends on the input-output data). Not sure if this is doable using this wavelet coherence method.

I have looked through the WaveletComp guided tour document but did not find a solution for this purpose – derive one series based on the other. The wavelet reconstruction function does not work for this purpose.

I would appreciate if you have any suggestions. Thank you!

Gavin

Hi Gavin, I’m so sorry I missed your comment! This is a really great question, but unfortunately, there aren’t too many ways that you can actually make use of the coherence information for reconstruction purposes. Once you find the times and frequencies where you have maximized correlation between the two time series, you can potentially create something like a “mask” where you zero out all other time periods and frequencies that are not meaningful and try to reconstruct a signal for each time series using the wavelets during those times. It would essentially be a very disjoint and not meaningful set of reconstructed signals. I’m curious if you were able to find a different method to do this, instead of using wavelets.

Hi Rohini, Thanks for your comments. The other possible method on top of my mind is something as “lagged regression”, where regression coefficients are optimized globally between 2 time series. The nice feature of wavelet coherence is that coefficients are optimized locally, reflecting the localized time and frequency information. Due to the nonstationarity of time series (e.g., involving relation between 2 time series), it looks to me that regression based on wavelet coherence has great merits compared to the lagged regression. That’s where my question comes from.

Hi Rohini, thanks so much for providing such a clear introductory to wavelet transform! I’m learning wavelets recently and have some questions about wavelet power. I find in Torrence 1998, the wavelet power is presented by the amplitude squared, |wi|^2, however, this will lead to systematically reduction in power of high frequency.(check the frequently asked questions: https://paos.colorado.edu/research/wavelets/faq.html#bias). And later, Liu et al. (2007) has corrected the wavelet power by dividing by its scale. |wi|^2/scale. I’m not sure if I am understanding it correctly, but just leave the comment here in case it is of help to anyone else! Thank you!

Liu, Y., X. San Liang, and R. H. Weisberg. 2007. Rectification of the Bias in the Wavelet Power

Spectrum. Journal of Atmospheric and Oceanic Technology 24:2093-2102.