%\documentclass{elsart}
% Use the option doublespacing or reviewcopy to obtain double line spacing
%\documentclass[doublespacing]{elsart}
\documentclass[doublespacing]{elsart}
% the natbib package allows both number and author-year (Harvard)
% style referencing;
\usepackage{natbib}
\usepackage{epsfig,amsmath,moreverb}
\journal{Computers and Geosciences}

%\oddsidemargin -0in %left margin on odd-numbered pages is 1in + 0in
%\topmargin -.5in %top margin is 1in -0.5in
\textwidth 6.375in %width of text
%\textheight 9in % size of page

\begin{document}

\begin{frontmatter}


\title{Classical Tidal Harmonic Analysis with Errors in MATLAB using T\_TIDE}

\author{Rich Pawlowicz\corauthref{cor1}}
\address{Dept. of Earth and Ocean Science, 
University of British Columbia, \\
6270 University Blvd., Vancouver, B.C., Canada, V6T 1Z4}
\ead{rich@ocgy.ubc.ca}
\corauth[cor1]{Corresponding author. Tel 604-822-1356, fax: 604-822-6091}

 
\author{Bob Beardsley, Steve Lentz}
\address{Woods Hole Oceanographic Institution, Woods Hole, MA, 02543 USA}

 

\begin{abstract} 
A standard part of any oceanic pressure gauge or current meter analysis is the
separation of tidal from non-tidal components of the signal. The tidal signal can
either be discarded, or its characteristics described in some fashion useful for
further analysis. Although tidal
signals can be removed by standard high or band-pass filtering techniques, their relatively
deterministic character and large amplitude make special techniques more
effective. In classical harmonic analysis, the tidal signal is modelled as the
sum of a finite set of sinusoids at specific frequencies related to
astronomical parameters. 
A set of programs has been written in MATLAB to a) perform 
classical harmonic analysis for periods of about one year or shorter, b) account 
for (some) unresolved constituents using nodal corrections, and c) compute 
confidence intervals for the analyzed components.
\end{abstract}

\begin{keyword}
% keywords here, in the form: keyword \sep keyword
tides \sep confidence interval \sep 


% PACS codes here, in the form: \PACS code \sep code
\PACS 92.10.Hm 

\end{keyword}

\end{frontmatter}

%%\section*{Introduction}

As the earth rotates on its axis, spatially varying gravitational forces from the moon and the 
sun act on the ocean, generating a forced elevation and current response primarily (but not solely) at
diurnal and semi-diurnal frequencies. Body forces act directly on deep oceanic
waters. Tidal effects in coastal regions are not directly forced by these astronomical
forces. Instead they arise as a side-effect of deep oceanic variability, 
propagating through shallower coastal waters as a wave.  
In a typical oceanic time series, tidal 
variability is often the largest signal.
Power spectra for such time series are often characterized by a broad hump with a low-frequency
maximum and a decline at higher frequencies. Superimposed are a number of sharp tidal peaks near
diurnal and semi-diurnal frequencies, and sometimes a 
broader peak associated with Coriolis or inertial effects. 
Dynamical analysis requires the separation of the tidal signal from sub or
super-tidal variations, or in some cases separation of tidal effects from inertial effects at
a nearby frequency. The tidal information is either discarded or kept for
further analysis.

Standard high/low/bandpass filtering techniques \citep[e.g.,][]{jackson:86} can be used but
in general these are inefficient because fairly narrow filters with a great deal of rejection
are needed. Also, although these are useful in analyzing non-tidal variability, they provide
no compression of the tidal information. Specialized techniques have been devised to take
advantage of the ``deterministic'' nature of tidal processes. In classical harmonic analysis,
the tidal forcing is modelled as a set of spectral lines, i.e., the sum of a finite set of 
sinusoids at specific frequencies. These frequencies are
are specified by 
various combinations of sums and differences 
of integer multiples of 6 fundamental frequencies arising from planetary motions 
\citep{godin:72}. These fundamental parameters  
represent effects of rotation of the earth (lunar day of 24.8 hours), the
orbit of the moon around the earth (lunar month of 27 days) and the earth 
around the sun (tropical year), and periodicities in the location of 
lunar perigee (8.85 years), lunar orbital tilt (18.6 years), 
and the location of perihelion ($\approx \! 21000$ years).
The set of 6 signed integers required to describe a particular
frequency are called the {\it Doodson numbers}. 
Many of the more important frequencies have names
such as ``$M_2$'', ``$K_1$'', etc. From astronomical considerations alone, 
an ``equilibrium'' response
can be predicted; this is the phase and amplitude that would be observed if the response of the
earth was fast enough that the surface deformation was effectively in equilibrium with the forcing at
all times. The real ocean is very definitely not in equilibrium with the tidal
forcing.  However, as tidal
amplitudes are small compared with total ocean depth, the dynamics are very nearly linear,
implying that the forced response contains only those frequencies present in the
forcing. A least squares fit can be
used to determine the relative phase and amplitude of each frequency in the response. This
phase/amplitude data thus provides a compression of the data in the complete time series,
which can then be compared with similar data at other locations to understand the characteristics
of tidal dynamics, or can be used to synthesize time series of tidal effects at other times
for predictive purposes. 

There are several drawbacks to classical harmonic analysis. The first is that, ignoring
the modulation of perihelion which is effectively constant over historical time, an 
$\approx18.6$ 
year time
series is required to resolve all of the listed frequencies (that is, the number of wavelengths of
each constituent in the record is at least 1 different from all other constituents). 
In practise,
record lengths are often one year or shorter.
In order to handle this issue
an assumption is made that the phase/amplitudes of response sinusoids with similar frequencies
(i.e., those whose first three Doodson numbers are identical)
are in the same proportion as those of the equilibrium response under the
reasonable premise that the ocean response should be similar at similar frequencies.
In such a cluster large equilibrium
peaks are surrounded by small subsidiary peaks in frequency-space which provide ``nodal
modulations'' (or more correctly, ``satellite modulations'') to 
the main peak. The appearance of the total signal will be a sinusoid whose phase and
amplitude varies slowly with time. These changes are slow enough to be considered effectively
constant for record lengths of up to one year.
At much shorter record lengths another problem arises. The frequency resolution further
degrades until even dissimilar constituents are unresolvable. The best solution is to
apply inference. This technique for finding the absolute phase/amplitude requires that the
relative differences in phase/amplitude between the two unresolved constituents is known
from other nearby data.
If this is not the case, it
is thought best to discard the smaller constituents and fit only to the largest in a given
frequency interval.

Another drawback of classical analysis is that it provides no easy way
to determine whether the resulting phase/amplitude of a given sinusoid is meaningful in a
deterministic way (i.e., it is truly a tidal line), or whether it results from fitting to a
component of the broad-spectrum
variability. In general a fit is likely to include elements of both and some kind of
confidence interval for the deterministic part would be useful.
In order to address this issue, the ``response'' method was invented \citep{response}.
Although this provides better results than classical harmonic analysis, it has
not found widespread use. 

Further problems with classical harmonic analysis arise in coastal regions where the tidal
response is in the form of a wave propagating onshore. In large estuaries,
the seasonal change in salinity and flow may change the dynamic response 
but as these changes can vary from year to year this is not really deterministic. 
Instead spectral peaks are broadened so that they are no
longer pure lines, but, depending on the situation, such variations may be
treated as lines in the analysis. Within smaller estuaries,
tidal height variations may be significant compared to water column depth and a variety
of nonlinear effects can occur. For example, flood periods shorten and intensify and ebbs lengthen. As 
long
as these effects are reasonably deterministic they may be handled by adding extra ``shallow
water'' constituents which occur at sum/difference frequencies of the major
constituents. More problematic in these regions are the effects of internal variability.
Tidal interactions with varying topography can produce large internal waves and bores whose
characteristics are highly sensitive to ambient stratification. In such cases the assumption
of ``line'' frequencies becomes questionable and other techniques such as wavelet analysis
have been suggested \citep{jay:99}. More comprehensive descriptions of analysis techniques,
their use, and their limitations is given in, e.g., \citet{foreman:95b} and 
\citet{godin:91}.

Here we describe T\_TIDE, a package of routines that can be used to
perform classical harmonic analysis with nodal corrections,
inference, and a variety of user-specified options. There are several novel
features. First, although the harmonic analysis algorithm with nodal corrections, etc., 
itself is not
original (other than conversion to complex algebra), it 
is implemented in MATLAB, an analysis package
widely used by oceanographers. This allows for easy use within the framework of a complete
analysis involving plotting of raw data, scatter plots, and so forth. Second, the code is
written directly in matrix terms and thus relatively easy to understand and modify
if required. Finally, in order to differentiate between true deterministic
(line) frequencies and broad-spectrum variability, confidence intervals for the estimated
tidal parameters are computed using one of several user-selectable algorithms. The
package is made up of a number of files each of which contains one or more functions.
User-callable functions generally have a ``{\it t\_}'' prefix to prevent namespace
collisions.

The paper is composed of four parts. In the first the form of the equilibrium
potential is described. In the second part the mathematical basis of the technique
for making phase/amplitude estimates is described. In the third part the generation of confidence
intervals is outlined. Finally an illustrative example is discussed. 

\section{Tidal Potential}

The effect of gravitational force vectors from the sun and moon, $\bf F$, can be written as the gradient of
a scalar potential $V$, ${\bf F}= -\nabla V$. The magnitude of this potential at the earth's
surface at any time obviously is dependent on the relative positions of the earth, moon,
and sun. In Doodson's development (\citealp{doodson}, described in \citealp{godin:72}) 
the potential is 
written as a function of
lunar time $\tau$ (defined to begin at ``lunar midnight'') 
and other astronomical variables (which are also functions of time):
\begin{eqnarray}
 s & = & \text{mean longitude of moon,} \nonumber \\
 h & = & \text{mean longitude of sun,} \nonumber \\
 p & = & \text{longitude of perigee,} \nonumber \\
 N' & = & \text{negative of the longitude of the ascending node,} \nonumber \\
 p' & = & \text{longitude of perihelion,}\nonumber 
\end{eqnarray}
where all terms are in units of cycles. These variables can be evaluated for a 
given Julian date
using function {\it t\_astron} which implements formulas in \citet{ephemerides}.
Their effects are combined with the Doodson number set for a particular
constituent $\{i', j', k', l', m', n' \}$
into the astronomical argument $V_a = i' \tau + j' s + k' h + l' p + m' N' + n' p'$. Sets with a common $i'$ are
called a species (thus the slow, diurnal, and semidiurnal species for $i'=0,1,2$ respectively), 
and sets with common $i'j'k'$ are called a subgroup. The constituent 
frequency $\sigma$ is defined as $\sigma = 2\pi d V_a / dt $.
The tidal potential is then written in the form
\begin{equation}
V = \sum_{i' = 0}^3 \left[ G_{i'} (\theta) \sum_{j' k' l' m' n'}  A'_{i'j'k'l'm'n'} \cos(2 \pi V_a)
 + G'_{i'} (\theta) \sum_{j' k' l' m' n'}  B'_{i'j'k'l'm'n'} \sin(2 \pi V_a) \right]. \label{equilib}
\end{equation}
For a given Doodson number set either $A'$ or $B'$ is nonzero, but not both. 
These constants are 
tabulated and stored in data structures that can be loaded using {\it t\_getconsts}. 
The geodetic 
functions $G_{i'}$ and $G'_{i'}$ vary with species $i'$ and latitude $\theta$, and also 
depend on
such constants as the radius of the earth and the masses and separations of the earth, 
moon, and
sun. The equilibrium amplitude for a constituent is defined as $G A'/g$ or $G' B'/g$, where
$g$ is the gravitational acceleration, and can
be generated for a particular latitude using {\it t\_equilib}.


\section{Phase/Amplitude Estimates}

The algorithm used here for making phase and amplitude estimates is 
based on algorithms and FORTRAN code described by \citet{godin:72}, 
\citet{foreman:77},  and
\citet{foreman:78}. 
However, unlike those authors we use complex algebra directly rather than
deal with sine and cosine fits separately. This has the advantage of
unifying the treatment for scalar (e.g., pressure) and vector (e.g.,
horizontal currents) time series which are represented as complex numbers $u+iv$.

Consider a time series of observations $y(t), t=t_1, t_2, ... t_M$ arranged in a vector, 
where the observation
times are regularly spaced at an interval $\Delta t$ (default 1 hour) and $M$ is an odd number
(an endpoint is discarded if required).
The time axis is defined such that the origin (or central time) is at $t_{(M+1)/2}$.
Some missing observations can be handled
by using a ``missing data'' marker in the input vector (by MATLAB convention this is NaN, 
the IEEE arithmetic representation for Not-a-Number).
This regular interval restriction does not arise from the least-squares
fit itself but rather from the automated constituent-selection algorithm and is also a
requirement when spectra are estimated in one of the confidence interval algorithms.
This time series may be  composed of either real or complex numbers. The time series is
passed to the analysis program {\it t\_tide} along with a variety of (mostly optional) parameters.
 
The tidal response is modelled as
\begin{equation}
x(t) = b_0 + b_1 t + \sum_{k=1...N} \!\! a_k e^{i \sigma_k t } + a_{-k} e^{-i \sigma_k t} \label{complex}
\end{equation}
where $N$ constituents (each with unique Doodson number sets) are used. 
Each constituent has a frequency $\sigma_k$ which is known from the development of the potential,
and a complex amplitude $a_k$ which is not known, although if $y(t)$ is a real time series
$a_k = \overline{a_{-k}}$. A possible offset and (optional) linear drift is handled by the first two terms.
The traditional approach uses real sinusoids:
\begin{equation}
x(t) = b_0 + b_1 t + \sum_{k=1...N} \!\! A_k \cos (\sigma_k t) + B_k \sin (\sigma_k t) \label{real}
\end{equation}
and is related to (\ref{complex}) by $A_k = a_k + a_{-k}$ and $B_k = i (a_k - a_{-k})$.
The real representation is more convenient for the linear error analysis described later. 

Constituents can be chosen from a list of 45 astronomical and 101 shallow-water
constituents. Data structures containing names and other information about these constituents are
loaded using {\it t\_getconsts}. There are several alternatives for
selecting constituents. For general use, an automated selection
algorithm \citep[following][]{foreman:77} is in place, which works as follows. A basis 
of all astronomical and 24 of the most
important shallow-water constituents are gathered together. 
All constituents are listed in order of
predefined importance. Less important constituents whose frequencies are less than
a Rayleigh resolution limit $\alpha (N \Delta t)^{-1}$ (with default $\alpha = 1$)
apart from more important constituents in frequency are 
discarded. Additional shallow-water constituents can be specified if required.
If the relative phase/amplitude of two constituents that are otherwise unresolvable is
known from other sources, then an inference procedure can be carried out.
Alternatively, constituent lists can be explicitly specified.

The least squares fit are the coefficients minimizing
\begin{equation}
E = \sum_m | x(t_m) - y(t_m) |^2 = || T a - y ||^2
\end{equation}
where $y = [ y(t_1), y(t_2), ... y(t_M)]'$, $a=[b_0,b_1,a_1,a_{-1},a_2,a_{-2},...a_{-N}]'$, and
T is an $M\times 2N+2$ matrix of linear and sinusoidal basis functions 
evaluated at observation times. The solution is found using the 
Matlab ``\verb+\+'' matrix division operator.

Once the fit has been performed, various corrections are applied. These are generated in {\it t\_vuf}.
First, the phase of the constituent response is usually reported as ``Greenwich phase'' $g_k$, that is, phase 
referenced to the phase of
the equilibrium response at $0^{\circ}$ longitude (the Greenwich meridian). This
can be interpreted
as reporting the phase of the response at the time when the equilibrium forcing is at its largest 
positive value at $0^{\circ}$ longitude.
It is simplest to find the fitted phase at the central time of the record ($t=0$); the equilibrium phase $v_k$ 
is then just $V_a$ for the given constituent computed at the Julian date corresponding to this
central time, with possible 
adjustments of $1/4, 1/2,$ or $3/4$ cycle depending on whether $A$ or $B$ is nonzero, and their signs.

Second, if a latitude is specified, then nodal or satellite corrections are computed as follows. Consider
a main peak of index $k$ with satellites with indices $kj$. The effect of the different satellites will 
be to slowly modulate the phase/amplitude of the main peak over various periods, all more than 8 years. 
Our fitted
response $\hat{a}_k$ over some period can then be written as a modification of the ``true'' response of
the main constituent $a_k$, 
in which the amplitude is changed by a factor $f_k$ and the phase by an angle $u_k$ due to the
presence of the satellites. $f_k$ and $u_k$ are called 
the nodal correction amplitude and phase respectively.
That is,
\begin{equation}
\hat{a}_k e^{i \sigma_k t} = f_k a_k e^{i \sigma_k t + i u_k} = a_k e^{i \sigma_k t } + \sum_j a_{kj} e^{i \sigma_{kj} t} 
\end{equation}
Cancelling common terms, we have
\begin{eqnarray}
f_k e^{i u_k }  &=&  1 + \sum_j \frac{a_{kj}}{a_k} e^{i (\sigma_{kj} - \sigma_k) t }  \nonumber \\
f_k e^{i u_k }  &\approx&  1 + \sum_j \frac{a_{kj}}{a_k}   
\end{eqnarray}
The final approximation will hold as long as $(\sigma_{kj} - \sigma_k) t$  remains 
``small'' (i.e., $N\delta t << 8$ years).
In general the true phases and amplitudes of the satellites are not known. However, since their
frequencies are very similar to that of the main peak it is standard to assume the ratio of true amplitudes
is the same as the ratio of amplitudes in the equilibrium response,
and the difference in true phases will be equal to the difference in equilibrium phases. The nodal
corrections are thus computed from the equilibrium response (\ref{equilib}). 
A latitudinal 
dependence arises from the geodetic functions. $G'_1$ is zero at the
equator and a crude limiting is used to prevent some diurnal corrections from getting overly large.
The validity of using the latitude-dependent equilibrium response to predict an aspect of the dynamic behavior
in one part of an ocean basin in such a simple way is not clear. If the record length is longer
than one year the comparison of successive one year analyses with and without nodal correction
can be used to test this validity of this process.



The products of the analysis above are a pair of complex values $\{ a_k, a_{-k} \}$,
possibly corrected for nodal
modulation,  for
each constituent $k$.  These are converted into standard parameters:
\begin{eqnarray}
U_{k} & = &  |a_k| + | a_{-k}|, \label{ell1} \\
u_{k} & = &  |a_k| - | a_{-k}|, \\
\theta_k & = & \frac{ang( a_k ) + ang( a_{-k} )}{2}, \\
g_k      & = & v_k  - \frac{ang( a_k ) - ang( a_{-k}) }{2}.  \label{ell2}
\end{eqnarray}
For horizontal currents, these parameters describe the features of an ellipse traced out by
the tip of the velocity vector:
the length of the semi-major and semi-minor axis of the ellipse ($U_{k}, u_{k}$ 
respectively), the inclination
of the northern semi-major axis counter-clockwise from due east $\theta_k$, and the Greenwich phase 
$g_k$. If $u_k>0$ ($<0$) then the ellipse is traced in a counter-clockwise (clockwise) direction.
For scalar time series, the parameter $U_k$ is merely called the amplitude, and 
$u_k, \theta_k \equiv 0$.
Note that the restriction of the definition of inclination to use the northern axis 
means that analyses of constituents
whose ellipses are aligned in an east/west direction may have inclinations that fluctuate between
close to $0^{\circ}$ and nearly $180^{\circ}$ due to noise. These apparently large jumps are an
artifact of the restriction but do not represent similarly large changes in the physical
behavior. 

Once ellipse parameters are found, these can be used for further analysis. They can also be used to
generate predictions at other times using {\it t\_predic}. Nodal corrections in {\it t\_predic}
are computed at the time series mid-point so that it is an exact inverse of {\it t\_tide}.

\section{Confidence Intervals}

One drawback of classical harmonic analysis is that the degree to which a given constituent
represents true tidal energy as opposed to the energy of a broad-band nontidal process is
not determined. This is useful information for two reasons: first, it allows one to make better
estimates of the tidal behavior, and second, it can allow one to quantitatively compare
different analyses. There are two steps to producing confidence intervals. First, we must
form an estimate of the characteristics of non-tidal or residual ``noise'' affecting 
the $a_k$ (or $A_k, B_k$). 
Second, we must convert these estimates into confidence intervals for the standard parameters
through a nonlinear mapping. We discuss the case of real time series first.

\subsection{Residual noise (real)}

After the harmonic analysis for an N-point real time 
series $y(t)$ is performed, we examine the structure of the residual
series. In the simplest case, the residuals are statistically  Gaussian and uncorrelated in time.
If this is the case then the total residual power $P_T = \sigma_x^2 = P/\Delta t$, where
$P$ is the two-sided spectral density. The amplitude of the fit to sine and cosine terms ($A$ and $B$ respectively) 
will be contaminated by errors arising from
unresolved noise components within a frequency interval of $\Delta f = (N \Delta t)^{-1}$
around the line. Thus $\sigma^2_{A} = \sigma^2_{B} = P\Delta f =  \sigma_x^2 /N$. It is unlikely that
a geophysical series will be spectrally flat, and a more sophisticated approach 
used in {\it t\_tide} is to find
a local value of $P$ by making a spectral estimate from the residual time
series (i.e., after the removal of all fitted constituents) and averaging the power over
frequency bins in a window around the frequency of any constituent, neglecting bins
in which fitted constituents reside. Here we chose a sequence of windows of width
0.4 cpd centered on $1,2,3,...$ cpd.

\subsection{Conversion to standard parameters (real)}

A conversion from errors in the cos/sine amplitudes to errors
in standard parameters (amplitude and phase) can be done  using
a linearized analysis. Consider a constituent $k$. Let
$\xi = F(A_k,B_k) $ be a nonlinear function of these
parameters, either the amplitude or the Greenwich phase. 
Then if $\{A_k,B_k\}$ are independent random variables, we
can find a linearized estimate of the standard error of $\xi$ in terms of
the standard errors of the sinusoid amplitudes:
\begin{equation}
\sigma_\xi^2 = \left( \frac{\partial F}{\partial A_k} \right)^2 \sigma_A^2 +
\left( \frac{\partial F}{\partial B_k} \right)^2 \sigma_B^2, 
\end{equation}
where the partial derivatives can be derived exactly (but tediously) 
from (\ref{ell1}-\ref{ell2}). 

Alternatively the nonlinear mapping can be handled directly using a ``parametric bootstrap'' 
\citep{bootstrap}.
In this case the residual variance estimates are used by the code
to simulate a number of realizations or replications of the 
analysis by taking the estimates of the
sinusoid amplitudes and adding Gaussian noise with the appropriate variance to them. 
All of these realizations are then converted
nonlinearly to standard parameters using (\ref{ell1}-\ref{ell2}) and an estimate of the 
standard error computed from this replicate dataset directly.

Once a standard error is found, 95\% confidence intervals can be estimated
using standard techniques. Alternatively, a signal-to-noise power ratio (SNR) can be computed
based on the ratio of amplitude to amplitude error.
Simulations performed in {\it t\_synth} (and described in the text file {\it t\_errors})
in which the variability of analyses carried out
on a fixed data set with different noise realization are compared with estimated confidence
intervals show that 
the linear procedure appears to be adequate for real time series (e.g. 
tidal height), as long as the SNR $>10$, and is probably 
not bad for SNR as low as 2 or 3. The nonlinear procedure gives similar results to the
linearized procedure at high SNR, and is more accurate at low SNR.

\subsection{Residual noise (complex)}

A complex residual time series $u+iv$ can be modelled as bivariate white noise and variances
$\sigma^2_u, \sigma^2_v$ and covariance $\sigma_{uv}$ computed. 
If we assume further
that the noise in both components is not correlated 
($\sigma_{uv} \approx 0$), then a coloured bivariate noise model
can be used and variances assigned to real and imaginary parts of constituent amplitude
separately on the basis of local spectral
densities as described above. If it is suspected that $\sigma_{uv} \ne 0$ then it is recommended that
the time series be rotated into a coordinate system in which this is true (e.g., into
along/across channel axes or into the principal axes).


\subsection{Conversion to standard parameters (complex)}

The linearized
analysis now involves functions of 4 variables, since both $A_k = A_r + iA_i$ and $B_k = B_r+iB_i$ have
real and imaginary parts, and analytic expressions for partial derivatives 
\begin{equation}
\sigma_\xi^2 = \left( \frac{\partial F}{\partial A_r} \right)^2 \sigma_{A_r}^2 +
\left( \frac{\partial F}{\partial A_i} \right)^2 \sigma_{A_i}^2 +
\left( \frac{\partial F}{\partial B_r} \right)^2 \sigma_{B_r}^2 +
\left( \frac{\partial F}{\partial B_i} \right)^2 \sigma_{B_i}^2 
\end{equation}
become large.
Some analytical simplification is possible by
assuming that all 4 variables are independent.

The bootstrap approach can also be applied to the complex coefficients $a_k$. One minor
complication that arises is that unless the noise is circular ($\sigma^2_u = \sigma^2_v, \sigma_{uv}=0$),
the errors in $a_k$ and $a_{-k}$ are correlated with each other. The bootstrap process requires
the generation of correlated noise replicates.




\section{Example}

The analysis of an example dataset provided in \citet{foreman:77} is shown in Figure \ref{f1}
and Table \ref{t1}. This example is included in datafile {\it t\_example.mat}. 
The example dataset consists of 66 days of hourly 
elevations with a 3 day gap. A tidal variation is visible superimposed on subtidal
variability. The time series can be loaded and analyzed using the demonstration script
{\it t\_demo}. Code within this script illustrates how the programs are called.
In this example,
the automated
constituent selection algorithm is used and it selects 35 constituents. 
In addition, one shallower water constituent
(M10) is manually added and two other constituents analyzed via inference.
The $P_1$ constituent is inferred from $K_1$, and $K_2$ is inferred from $S_2$.
Nodal corrections are performed. A linear trend is not included in the analysis.
The coloured bootstrap analysis is used to determine significance and confidence
intervals. 
Table \ref{t1} gives the output of the program. In the first column the name of
the constituent is given. Significant constituents (those with SNR in the last column
$>1$) are marked with a ``*''. The SNR is the ratio of amplitude (3rd column) to the 
error in amplitude (4th column).
Frequencies (first column) are listed in cph and Greenwich phase/phase 
error (5th and 6th columns) in degrees.
11 constituents
were judged to be
significant (only 6 would be significant at an SNR cutoff of 2). Figure \ref{f1}b
shows the residual series and an elevation series synthesized from the significant
constituents (note that it spans the 3-day gap). Results of the analysis are shown
in a spectrum in Figure \ref{f1}c. Most of the significant constituents are in the
diurnal and semidiurnal bands ($\approx0.04$ cph and $\approx0.08$ cph respectively) 
although several
higher-frequency constituents also appear to be significant. In spite of the large
amount of energy in the fortnightly band ($\approx 0.002$ cph) the fitted constituents
there are apparently not significant.
The
analyzed phases of significant constituents are shown in Figure \ref{f1}c. 
The significant constituents generally have reasonably
small phase errors.

\section{Summary}

Separation of tidal and nontidal energy  is an important task in
any analysis of oceanic time series. Here we discuss the 
theoretical foundation and implementation details
of a MATLAB package for classical harmonic analysis. The package can also compute
confidence intervals for the tidal parameters using one of three different sets of
assumptions about the structure of residual noise. An example is provided to show typical
results. The code is available at \verb+http://www.ocgy.ubc.ca/~rich+.

\section*{acknowledgements}
Elements of this MATLAB toolbox for tidal harmonic analysis were developed
over the last few years as part of the authors' research into tidal
phemonena.  We want to first acknowledge Mike Foreman's major contribution
in providing the FORTRAN codes and documentation to the community in
1977-78 that served as our starting point.  Rich Signell first derived the
formulae for the linearized error analysis used here, and XXX also
provided useful input. Support for this effort was provided by Natural Sciences 
and Engineering 
Research Council of Canada under grant OGPO194270 (RP)
and the National Science Foundation and Office of Naval Research (SL and
RB).  


%\bibliographystyle{/users/rich/doc/local/natsci}
\bibliographystyle{/users/rich/doc/local/elsart-harv}
\bibliography{/users/rich/doc/bib/oceanog,/users/rich/doc/bib/math,/users/rich/doc/bib/tide}


%\psfig{file=tidepic.eps,width=5in}
\begin{table}
\begin{center}
\psfig{file=PAPEROUT.ps,width=6in}
\caption{Analysis Results for Tuktoyuktuk dataset}
\label{t1}
\end{center}
\end{table}

\begin{figure}
\begin{center}
\psfig{file=tidepic.eps,width=5in}
\caption{Tuktoyuktuk analysis example. (a) Raw time series. 
(b) Upper curve is residual time series after removal of tidal signal.
Lower curve is a synthesized tidal series using significant constituents.
(c) Amplitude of all analyzed components with 95\% significance level. Note
frequency dependence. Significant constituents are marked with a solid circle.
(d) The phase of significant constituents with a 95\% confidence interval. 
}
\label{f1}
\end{center}
\end{figure}



\end{document}



