Wavelet Analysis

In this demo we will learn to apply wavelet analysis. This is a great way to visualize the variability in your time series.

We'll be working with data from the m1244 mooring in the Labrador Sea, which you can get in NetCDF form here. (This is included in the full distribution of the course materials.)

Many thanks to Jan-Adrian Kallmyr for helping with translating the Matlab tutorial into Python.

As usual, we have to import a few packages.

We'll need to install analytic_wavelet, a reimplementation of my wavelet code into Python by Dan Schwartz.

https://github.com/jonathanlilly/analytic_wavelet

First, if you're using an environment to work in these notebooks, activate that environemnt.

Next you'll need to install pip and git, if they are not installed already, with

conda install -c conda-forge git pip

Then the analytic_wavelet package is installed direclty from GitHub with

pip install git+git://github.com/danrsc/analytic_wavelet.git

Simple Smoothing Revisited

First we will review simple smoothing of our data.

Generally, the goal of smoothing is to remove high frequencies. This is why we often use the term lowpass, because we want to design a filter through which will let the low frequencies pass through, but which will reject the high ones.

Let's see how the boxcar or rectangle function performs as a lowpass filter. Whereas, as we have seen, the Fourier transform of a Gaussian is also a Gaussian, the Fourier transform of a boxcar is a very bouncy function.

Notice that rather than decaying to zero, the magnitude of the Fourier transform hits zero and bounces back repeatedly. This has undesireable consequences, as we will see.

Now we will use that boxcar to smooth a dataset. For simplicity, we will use a dataset consisting of white noise. We will look at what happens in time and also in frequency.

In the time domain, we see that the orange curve is a smoothed version of the blue curve, as expected. However, in the frequency domain, we see something quite odd.

The very lowest frequencies do indeed pass through, as desired. Yet one also sees bumps of energy extending to high frequencies that also passed through. Evidently the boxcar is not a very good lowpass filter! If we are just interested in quickly smoothing the data, it works fine, but if we are interested in looking at the frequency, we see it leaves a strange imprint.

Our time spent studying Fourier theory has an explanation for what is happening:

"Convolution in the time domain is multiplication in the frequency domain!"

When we smooth our time series by the boxcar function, this is accomplished through the convolution operation. Convolving two time seres, as we have stressed repeatedly, is equivalent to multiplying their Fourier transforms.

In fact we have plotted three curves in the lower panel of this plot. The orange curve is the Fourier transform of our smoothed noise time series. The green curve is the Fourier transform of the original noise time series multiplied by the Fourier transform of the boxcar. As you can see, they are virtually identical, as expected from the convolution theorem.

The bumps of the Fourier transform of the boxcar are called “sidelobes” or a “ringing” effect. Such ringing is due to the discontinuous nature of the boxcar. If we want to get rid of this ringing, we need a function that is more smooth in the time domain. This why we use the will typically use the Hanning window for simple smoothing, as it has a nicer Fourier transform.

Bandpass Filtering

Let's turn our attention now to a new application, the construction of a bandpass rather than lowpass filter. For this we will again use the deepest velocity at the m1244 mooring. First recall the rotary spectra of this time series.

We will construct a narrowband filter to remove the whole inertial peak (which includes a contribution from the semidiurnal tide at this latitude). Let's just run the code and then understand it afterward.

Blue is the original time series, orange is the residual (original-bandpassed), green is the bandpassed signal, and black is the wavelet (with a rescaled y-axis for presentational purposes).

As you can see, the negatively-rotating inertial peak has been removed, but the positive spectrum has not been touched at all. This means that the filtered version cvo has effectively isolated the inertial signal.

The wavelet transform essentially convolves a filter, the wavelet, with the original signal. This means we multiply the wavelet’s Fourier transform by the Fourier transform of the data. As the wavelet is a narrowband signal, this essentially shuts off all frequencies except those inside the band. Subtracting gives us the original signal with that band removed.

Filter Parameter Choices

Let's investigate further the choices of parameters.

The morse wavelets, constructed by the lines

morse = analytic_wavelet.GeneralizedMorseWavelet(gamma, beta)

psi, psif = morse.make_wavelet(len(cv_centered), omega)

are controlled by two parameters, gamma and beta, together with the filter frequency omega.

Gamma is a parameter controlling the filter shape, beta is a parameter controlling the filter width in frequency, and omega is the filter center frequency in radians per sample interval. We will understand all of these more later.

Note, a very important feature of the wavelet filter is that it is one-sided. That is, if we tell it to implement a bandpass filter centred around frequency 10, we mean frequency +10 and not frequency -10. Only positively phasors will be affected and not negatively rotating phasors.

To cause the wavelet transform to affect negative frequencies, as we wish to do here for the inertial filtering, we conjugate the input signal---which we saw will conjugate its Fourier transform but also flip positive and negative frequencies. Then to put the frequencies back, we also conjugate the output.

Notice also the factor of 1/sqrt(2) in the line

cvo = 1/np.sqrt(2)*np.conj(analytic_wavelet.analytic_wavelet_transform...)

This factor arises because there is a factor of sqrt(2) ambiguity in defining the amplitude of a real-valued versus a complex-valued sinusoid. Essentially we have defined a real-valued sinusoid to have amplitude one, with a consequence being that when we use wavetrans to form a bandpassed version of a complex valued signal, we need to divide by sqrt(2).

I chose the gamma and beta parameters as follows. Gamma is generally set to 3 for reasons we will see later. I found the frequency by noticing that the tidal peak was centered on about the 1200th Fourier frequency after calling plot(spp). The I picked beta by increasing it until it subjectively looked right.

Here is what the original, bandpassed, and residual time series look like, zoomed in to a 60 day period. Note the different y-limits for the middle subplot.