Island-trapped wave analysis
Here we take a closer look at the current meter record from the Alenuihaha Channel, on the Hawaii side, that we used to introduce EOF analysis. We start with some of the code in that notebook so that we can rotate the currents into approximate along-shore and cross-shore coordinates. We will then look at the spectra of the two components, followed by the rotary spectra. We will isolate the island-trapped wave signal via direct band-pass filtering and via complex demodulation.
%matplotlib inline
import os
from six.moves import urllib
import numpy as np
import matplotlib.pyplot as plt
import matplotlib as mpl
from pycurrents.file import matfile
from pycurrents.num import eof, spectra
from pycurrents.system import Bunch
# Colors from R's ggplot:
colors = ["#E24A33", "#348ABD", "#988ED5", "#777777",
"#FBC15E", "#8EBA42", "#FFB5B8"]
mpl.rcParams['axes.color_cycle'] = colors
mpl.rcParams['figure.figsize'] = (7, 5)
mpl.rcParams['figure.dpi'] = 90
mpl.rcParams['savefig.dpi'] = 90
fname = 'uv_alenuihaha.mat'
url = 'http://currents.soest.hawaii.edu/ocn_data_analysis_files/data/' + fname
if not os.path.exists(fname):
urllib.request.urlretrieve(url, filename=fname)
mcm = matfile.loadmatbunch(fname)
print(matfile.showmatbunch(mcm))
print(mcm.t.min(), mcm.t.max())
print(mcm.u.min(), mcm.u.max())
I don't know much about this data set, but evidently time is in days from some origin and current components are in cm/s. Let's plot the time series to get a feel for the record.
fig, axs = plt.subplots(nrows=2, sharey=True)
ax = axs[0]
ax.plot(mcm.t, mcm.u, 'r', label='u', alpha=0.5)
ax.plot(mcm.t, mcm.v, 'b', label='v', alpha=0.5)
ax.legend(loc='upper right')
ax = axs[1]
ax.plot(mcm.t, mcm.u, 'r', label='u', alpha=0.5)
ax.plot(mcm.t, mcm.v, 'b', label='v', alpha=0.5)
ax.set_xlim(380, 400)
We can see from the raw u, v time series that the two velocity components are correlated, and have roughly similar standard deviations. Another way to look at the correlation is with a scatter plot.
fig, ax = plt.subplots()
ax.plot(mcm.u, mcm.v, '.')
ax.set_xlabel('u (cm/s)')
ax.set_ylabel('v (cm/s)')
ax.set_aspect('equal')
ax.axhline(0, color='0.5')
ax.axvline(0, color='0.5')
Rotate to principle axes: approximately along-shore and cross-shore
The east-north coordinate system is arbitrary. It looks like a more natural coordinate system would be rotated about the mean value, with the first axis aligned with the long axis of the scatter cloud. To look at this more closely, we will glue the u and v components together as columns in an (n,2) matrix, and then use that to calculate the variance-covariance matrix. If we call this $U$, then the variance-covariance matrix will be $U^TU$.
mcm.udm = mcm.u - mcm.u.mean()
mcm.vdm = mcm.v - mcm.v.mean()
uv = np.vstack((mcm.udm, mcm.vdm)).T
print(uv.shape)
nt = len(mcm.u)
The variance-covariance matrix shows strong covariance between u and v.
vcv = np.dot(uv.T, uv) / nt
print(vcv)
vals, vecs = np.linalg.eig(vcv)
print(vals)
print(vecs)
The second eigenvalue is the big one, so that is the one corresponding to the eigenvector that is approximately along-shore. We will reverse the direction so that it is in the first quadrant instead of the third, and calculate the angle by which the velocities must be back-rotated to make the new u component be the along-shore component.
# arguments are the y and x components
# (first index is 1, and then 0)
# of the second column (second index is 1),
# with minus signs to select the first quadrant.
theta = np.arctan2(-vecs[1, 1], -vecs[0, 1])
print(np.rad2deg(theta))
Expressing 2-D vectors as complex numbers and multiplying by a complex exponential provides a simple method of rotating the vector, or equivalently, of transforming to a rotated coordinate system. Here, the semi-major axis is $\theta = 54^\circ$ CCW from the $x$-axis. To put all the velocities in a coordinate system in which this is the $x^\prime$-axis, we rotate them CW by this amount. In the complex plane this is a multiplication by $e^{-i\theta}$.
uvc = uv[:, 0] + 1j * uv[:, 1]
uvc_rot = uvc * np.exp(-1j * theta)
uv_rot = np.vstack((uvc_rot.real, uvc_rot.imag)).T
vcv_rot = np.dot(uv_rot.T, uv_rot) / nt
print(vcv_rot)
u_along = uvc_rot.real
u_cross = uvc_rot.imag
Above, we verified that the variance-covariance matrix of the rotated velocity matrix is diagonal, to within numerical precision. Transformation to the new coordinate system has diagonalized the matrix. The first component of uv_rot is now the along-shore component, and the second component is off-shore.
fig, axs = plt.subplots(nrows=2, sharey=True)
ax = axs[0]
ax.plot(mcm.t, u_along, 'r', label='along', alpha=0.5)
ax.plot(mcm.t, u_cross, 'b', label='offshore', alpha=0.5)
ax.legend(loc='upper right', ncol=2)
ax = axs[1]
ax.plot(mcm.t, u_along, 'r', label='along', alpha=0.5)
ax.plot(mcm.t, u_cross, 'b', label='offshore', alpha=0.5)
ax.set_xlim(380, 400)
Spectra
dt = mcm.t[1] - mcm.t[0]
speckw = dict(dt=dt, nfft=None, window='quadratic', smooth=5)
sr = spectra.spectrum(uvc, **speckw)
s_along = spectra.spectrum(u_along, **speckw)
s_cross = spectra.spectrum(u_cross, **speckw)
print(sr.keys())
print(s_along.keys())
fig, axs = plt.subplots(nrows=2, sharex=True, sharey=True)
ax = axs[1]
ax.loglog(sr.ccwfreqs, sr.psd[sr.ccwtup], 'r', label='CCW')
ax.loglog(sr.cwfreqs, sr.psd[sr.cwtup], 'b', label='CW')
ax = axs[0]
ax.loglog(s_along.freqs, s_along.psd, 'r', label='along')
ax.loglog(s_cross.freqs, s_cross.psd, 'b', label='cross')
ax.set_xlim(1e-2, 20)
ax.set_xlabel('cycles per day')
for ax in axs:
ax.legend(loc='lower left')
The big signals are the island-trapped wave centered near 0.38 CPD and the semi-diurnal tide. A diurnal tide is also visible, as is the first harmonic of the semi-diurnal tide.
The rotary spectra are almost identical; the flow is strongly controlled by the topography, so there is little systematic rotation of the velocity vector. This is confirmed by the spectra in the rotated coordinate system. The island-trapped wave signal is absent from the cross-shore component. There is a little bit of cross-shore velocity in the semidiurnal tide and possibly in its first harmonic, however. It looks like there might be a little energy at the sum of the island-trapped wave and the semi-diurnal tide frequencies--the first little peak to the right of the semi-diurnal. If this is real, it would suggest a nonlinear interaction between these two types of motion.
Complex demodulation
Let's use complex demodulation to look at the two biggest spectral peaks: first the semi-diurnal tide, and then the island-trapped wave. The functions are taken from the complex demodulation notebook, with minor modification.
def bl_filt(y, half_width):
"""
Simple Blackman filter.
"""
nf = half_width * 2 + 1
x = np.linspace(-1, 1, nf, endpoint=True)
x = x[1:-1] # chop off the useless endpoints with zero weight
w = 0.42 + 0.5 * np.cos(x * np.pi) + 0.08 * np.cos(x * 2 * np.pi)
ytop = np.convolve(y, w, mode='same')
ybot = np.convolve(np.ones_like(y), w, mode='same')
return ytop / ybot
def complex_demod(t, x, central_period, hwidth = 2):
"""
Complex demodulation of a real or complex series, *x*
of samples at times *t*, assumed to be uniformly spaced.
*central_period* is the period of the central frequency
for the demodulation. It should be positive for real
signals. For complex signals, a positive value will
return the CCW rotary component, and a negative value
will return the CW component (negative frequency).
Period is in the same time units as are used for *t*.
*hwidth* is the Blackman filter half-width in units of the
*central_period*. For example, the default value of 2
makes the Blackman half-width equal to twice the
central period.
Returns a Bunch; look at the code to see what it contains.
"""
rotary = x.dtype.kind == 'c' # complex input
# Make the complex exponential for demodulation:
c = np.exp(-1j * 2 * np.pi * t / central_period)
product = x * c
# filter half-width number of points
dt = t[1] - t[0]
hwpts = int(round(hwidth * abs(central_period) / dt))
demod = bl_filt(product, hwpts)
if not rotary:
# The factor of 2 below comes from fact that the
# mean value of a squared unit sinusoid is 0.5.
demod *= 2
reconstructed = (demod * np.conj(c))
if not rotary:
reconstructed = reconstructed.real
if np.sign(central_period) < 0:
demod = np.conj(demod)
# This is to make the phase increase in time
# for both positive and negative demod frequency
# when the frequency of the signal exceeds the
# frequency of the demodulation.
return Bunch(t=t,
signal=x,
hwpts=hwpts,
demod=demod,
reconstructed=reconstructed,
resid=(x - reconstructed))
def plot_demod(dm, difference=True):
fig, axs = plt.subplots(3, sharex=True)
if dm.signal.dtype.kind == 'c':
axs[0].plot(dm.t, dm.signal.real, label='signal.real')
axs[0].plot(dm.t, dm.signal.imag, label='signal.imag')
if difference:
axs[0].plot(dm.t, dm.resid.real, label='difference real')
axs[0].plot(dm.t, dm.resid.imag, label='difference imag')
else:
axs[0].plot(dm.t, dm.signal, label='signal')
axs[0].plot(dm.t, dm.reconstructed, label='reconstructed')
if difference:
axs[0].plot(dm.t, dm.resid, label='difference')
axs[0].legend(loc='upper right', fontsize='small', ncol=3)
axs[1].plot(dm.t, np.abs(dm.demod), label='amplitude', color=colors[3])
axs[1].legend(loc='upper right', fontsize='small')
axs[2].plot(dm.t, np.angle(dm.demod, deg=True), '.', label='phase',
color=colors[4])
axs[2].set_ylim(-180, 180)
axs[2].legend(loc='upper right', fontsize='small')
for ax in axs:
ax.locator_params(axis='y', nbins=5)
return fig, axs
Start with the semi-diurnal tide; demodulate at frequency of the strongest constituent, M2.
dm_M2 = complex_demod(mcm.t, u_along, 12.42/24)
plot_demod(dm_M2)
The spring-neap cycle is evident; it is quite regular in the second half of the record, but less so in the first half. In the middle we see the ITW signal in the tidal amplitude.
Merely removing the semidiurnal tide from the record is enough to leave the ITW as the most visually obvious feature. Now, try demodulating at its frequency:
dm_ITW = complex_demod(mcm.t, u_along, 1/0.39)
plot_demod(dm_ITW, difference=False)
A little bit of experimentation showed that a frequency of 0.39 CPD gives a more nearly constant phase during the biggest event than the original estimate of 0.38 CPD. We can also see that during the event the phase first decreases and then increases with time, so the frequency starts a little bit lower than 0.39 CPD and then increases.
We can also experiment with the bandwidth of the demodulation filter. First, double the bandwidth by cutting the filter half-width from 2 periods down to 1:
dm_hw1 = complex_demod(mcm.t, u_along, 1/0.39, hwidth=1)
plot_demod(dm_hw1, difference=False)
It still works--the ITW signal is still shown clearly--but the amplitude and phase is looking a little jittery, and this jitter is probably not useful information here. Now try going the other direction, with a filter half-width of 4 times the period:
dm = complex_demod(mcm.t, u_along, 1/0.39, hwidth=4)
plot_demod(dm, difference=False)
That's also a plausible picture, but comparing it to the case with the half-width of 2 periods, I think the latter is a better representation of the signal. This last plot appears to be spreading the energy out in time more than it should. Let's try one more thing: we will redo the demodulation with a half-width of 2, but it will be applied to the signal with the semidiurnal band removed:
u_detide = u_along - dm_M2.reconstructed
dm = complex_demod(mcm.t, u_detide, 1/0.39, hwidth=2)
fig, axs = plot_demod(dm, difference=False)
axs[0].set_ylim(-200, 200)
As expected, this makes little difference to the result, but it makes it easier to see how the ITW band is being extracted. Notice the lower frequency oscillation at the end of the record. It looks like the period might be 15-20 days. What is it? There is no hint of a corresponding peak in the spectra. Is it nothing more than the manifestation of a broad spectrum at a particular time, or does it correspond to an event with an identifiable dynamic balance?
Direct filtering
Complex demodulation accomplishes bandpass filtering as a side effect; if a bandpass filter is what you need, it is often a good choice. More generally, however, one might need a lowpass filter or a highpass filter instead of a bandpass. In fact, complex demodulation itself is accomplished with a lowpass filter.
There are many kinds of filter, but they all involve tradeoffs: a sharp cutoff comes at the expense of ringing and/or non-localized response to transients. You will often see references in the literature to a Butterworth filter. This is a digital filter with an Infinite Impulse Response (IIR), meaning that the effect of an impulse in the input is not completely localized, but can linger indefinitely. Digital filters like this are run through the data in one direction to yield a particular spectral transfer function, but this distorts the phase. The solution is to run the filter through twice, first going forward and then going back. This doubles the order of the filter, and the corresponding steepness of the response cutoff, and it restores the phase.
I have never used a Butterworth filter myself, but let's give it a try. I will specify a filter of order 2 so that after using filtfilt to run it twice we will have a filter of order 4, which seems to be most commonly used in the oceanographic literature. The response rolloff will therefore be 24 db/octave.
from scipy.signal import butter, filtfilt
# I'm picking cutoff frequencies by eye from the spectrum,
# roughly matching the bounds of the ITW peak.
cutoffs_cpd = np.array([2.5e-1, 6.5e-1])
nyquist = 12 # cycles per day
# Notice how the frequencies are specified: scaled
# by the nyquist.
b, a = butter(2, cutoffs_cpd / nyquist, btype='bandpass')
uf = filtfilt(b, a, u_along)
fig, ax = plt.subplots()
ax.plot(dm.t, uf, label='Butterworth')
ax.plot(dm.t, dm_ITW.reconstructed, label='complex demod')
# Uncomment to compare to complex demod with broader BW.
#ax.plot(dm.t, dm_hw1.reconstructed, label='c.d., hw=1')
ax.legend()
It worked! The result of the Butterworth filter is quite similar to our original complex demodulation, but there are also differences. Over most of the series, the Butterworth result has higher amplitude.
We can also compare the two filter methods on the series with the semidiurnal tide filtered out.
uf = filtfilt(b, a, u_detide)
fig, ax = plt.subplots()
ax.plot(dm.t, uf, label='Butterworth')
ax.plot(dm.t, dm.reconstructed, label='complex demod')
ax.legend()
The only obvious difference is in the last cycle on the right; the complex demodulation result looks unchanged, but the amplitude of the Butterworth is reduced and matches the complex demodulation better when the semidiurnal tide has been removed.