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.)
This notebook requires my jlab toolbox to be installed. You will also need to have set up Jupyter Lab to work with Matlab, following these instructions.
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.
set(groot,'defaultfigurepaperposition',[1 1 7 5]) %set default figure size
t=[1:1001]';t=t-mean(t);
g=zeros(size(t));
g(abs(t)<10)=1/19;%this is a 19-point boxcar filter
subplot(2,1,1),plot(t,g),xlabel('Time')
subplot(2,1,2),plot(abs(fft(g))),xlabel('Frequency'),axis tight
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.
x=randn(1001,1);
fx=vfilt(x,ones(19,1)./19,'mirror');
subplot(2,1,1),plot(t,[x fx]),xlabel('Time'),ylim([-3 3])
subplot(2,1,2),plot(abs([fft(x) fft(fx) fft(x).*fft(g)])),
xlabel('Frequency'),axis tight,ylim([0 80])
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 yellow 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.
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.
load m1244
use m1244
dt=num(2)-num(1); %Sample rate in days
cv=cv(:,4);cv=cv-mean(cv);
[psio,lambda]=sleptap(length(cv),16);
[f,spp,snn,spn]=mspec(dt,cv,psio,'cyclic');
subplot(1,2,1),plot(f,snn),flipx
xlabel('Negative Frequency (cycles/day)')
subplot(1,2,2),plot(f,spp)
xlabel('Positive Frequency (cycles/day)')
for i=1:2
subplot(1,2,i)
xlog,ylog,axis tight,ax=axis;
vlines(2,'k:'),vlines(max(f)/8,'C')
end
set(gcf,'paperposition',[0 0 11 5 ])
packfig(1,2,'columns') %pack the column subplots together
SLEPTAP calculating tapers of length 512.
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.
In this next plot, we will zoom on on the region to the high-frequency side of the vertical gray lines in the previous plot.
omega=2*pi*1200/length(cv);
cvo=frac(1,sqrt(2))*conj(wavetrans(conj(cv),{3,100,omega}));
cvr=cv-cvo;%Residual minus filtered version
[f,sppo,snno]=mspec(dt,cvo,psio,'cyclic');%Spectrum of filtered version
[f,sppr,snnr]=mspec(dt,cvr,psio,'cyclic');%Spectrum of residual
%The next line returns the filter used by wavetrans
[psi,psif]=morsewave(length(cv),3,100,omega,'bandpass');
subplot(1,2,1),plot(f,[snn snnr snno])
axis(ax),xlog,ylog,xlim([max(f)/8 max(f)]),flipx
hold on,plot(f,50*squared(psif(1:length(f))),'k') %get positive side only
xlabel('Negative Frequency (cycles/day)'),vlines(2,'k:')
subplot(1,2,2),plot(f,[spp sppr sppo])
axis(ax),xlog,ylog,xlim([max(f)/8 max(f)])
xlabel('Positive Frequency (cycles/day)'),vlines(2,'k:')
set(gcf,'paperposition',[0 0 11 5 ])
packfig(1,2,'columns') %pack the column subplots together
Blue is the original time series, yellow is the bandpassed time series, orange is the residual (original-bandpassed), and black is the wavelet (with a rescaled y-axis for presentational purposes). Note that the original and bandpassed time series are virtually identical outside of the wavelet footprint.
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.
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 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.
subplot(3,1,1),uvplot(num-datenum(1996,12,31),cv),axis tight
subplot(3,1,2),uvplot(num-datenum(1996,12,31),cvo),axis tight
hold on,plot(num-datenum(1996,12,31),[abs(cvo) -abs(cvo)],'k','linewidth',2)
subplot(3,1,3),uvplot(num-datenum(1996,12,31),cvr),axis tight
for i=1:3,subplot(3,1,i),xlim([0, 60]),end
packfig(3,1,'rows')
xlabel('Day of Year 1997')
Rather than just lowpassing, we have selectively removed only the inertial peak. A very useful facet of this analysis is that we now have the amplitude of our extracted signal, the inertial band in this case, as plotted with the heavy curve in the second panel.
We could also now proceed to look at the variance ellipses, statistics, etc, of the inertial band and residual separately. Combining simple statistics with bandpass filtering in this way can be very illuminating.
Next let's remove the tidal peak instead. This will be done by working with the positive and negative sides separately. However, for reference, code to do the same thing using the real and imaginary parts is also included, If one uses the real and imaginary parts of the data, there is no factor of sqrt(2). But since the wavelet is a one-sided bandpass, It will return something complex-valued if we give it a real-valued input. Consequently we need to take the real part of the both output variables.
cvo1=wavetrans(cv,{3,500,2*pi*1200/length(cv)});
cvo2=conj(wavetrans(conj(cv),{3,500,2*pi*1200/length(cv)}));
cvo=frac(1,sqrt(2))*(cvo1+cvo2);
%%This next bit is the same as the above
%[uplus,vplus]=wavetrans(real(cv),imag(cv),{3,500,2*pi*1200/length(cv)});
%cvo=real(uplus)+1i*real(vplus); %Note must take real part!
cvr=cv-cvo;
[f,sppo,snno]=mspec(dt,cvo,psio,'cyclic');
[f,sppr,snnr]=mspec(dt,cvr,psio,'cyclic');
%figure
subplot(1,2,1),plot(f,[snn snnr snno])
axis(ax),xlog,ylog,xlim([max(f)/8 max(f)]),linestyle U T,flipx
xlabel('Positive Frequency (cycles/day)'),vlines(2,'k:')
subplot(1,2,2),plot(f,[spp sppr sppo])
axis(ax),xlog,ylog,xlim([max(f)/8 max(f)]),linestyle U T
xlabel('Positive Frequency (cycles/day)'),vlines(2,'k:')
set(gcf,'paperposition',[0 0 11 5 ])
packfig(1,2,'columns')
We see that the tidal peak has indeed been removed from both sides.
Now that we have a general understanding of the wavelet transform as implementing a one-sided bandpass, we will turn to the wavelet transform itself.
Wavelet analysis is an informative visualization tool as well as the jumping-off point for a number of rigorous feature identification methods. Let's get a feeling for it by returning to the example of the deepest depth at the m1244 mooring. We'll take a look at the wavelet transforms of the alongstream and cross-stream velocity components we worked with in the earlier labs.
use m1244
dt=num(2)-num(1); %Sample rate in days
phi=angle(mean(cv(:,4)));
cv=rot(-phi)*cv(:,4);
gamma=3;beta=2;
fs=morsespace(gamma,beta,length(cv));
[wu,wv]=wavetrans(real(cv),imag(cv),{gamma,beta,fs});
h=wavespecplot(yearfrac(num),vfilt(cv,24),dt*2*pi./fs,wu,wv);
axes(h(1)),ylim([-15 50]),caxis([0 15])
legend('Alongstream Velocity','Cross-Stream Velocity')
for i=2:3
axes(h(i)),ylabel('Period (days)')
hlines(1./(tidefreq('m2')*24/2/pi),'E')
hlines(1./(corfreq(lat)*24/2/pi),'1.5k:')
%width of edge effect region in days
L=dt*2*sqrt(2)*sqrt(gamma*beta)./fs;
plot(yearfrac([num(1)+L/2 num(end)-L/2]),dt*2*pi./fs,'color','w');
end
axes(h(1)),title('Alongstream and Cross-Stream Wavelet Transforms')
axes(h(2)),text(1997,70,'Alongstream Transform','color','k')
axes(h(3)),text(1997,70,'Cross-Stream Transform','color','k')
set(gcf,'paperposition',[0 0 11 8 ])