Harmonic Analysis: From Observations to Constituents

Harmonic Analysis: From Observations to Constituents#

Harmonic analysis is the inverse of tidal prediction: given a time series of observed sea levels, the goal is to recover the amplitudes and phase lags of the tidal constituents. This section presents the mathematical formulation implemented in PyFES.

The Observation Equation#

At discrete observation times \(t_j\) (\(j = 1, \ldots, M\)), the observed sea level can be modelled as a sum of \(N\) tidal constituents plus a residual:

\[h(t_j) = S_0 + \sum_{k=1}^{N} f_k(t_j) \left[ a_k \cos\!\left(V_k(t_j) + u_k(t_j)\right) + b_k \sin\!\left(V_k(t_j) + u_k(t_j)\right) \right] + \varepsilon_j\]

where:

The unknown parameters are \(S_0\), \(a_k\), and \(b_k\) for each constituent.

Least-Squares Formulation#

Defining the modulated phase of constituent \(k\) at time \(t_j\) as:

\[\phi_{k,j} = V_k(t_j) + u_k(t_j)\]

and collecting all unknowns into a vector \(\mathbf{x} = (S_0, a_1, b_1, a_2, b_2, \ldots, a_N, b_N)^T\) of dimension \(2N + 1\), the observation equation becomes the linear system:

\[\mathbf{h} = \mathbf{A} \mathbf{x} + \boldsymbol{\varepsilon}\]

where the design matrix \(\mathbf{A}\) has \(M\) rows (one per observation) and \(2N + 1\) columns:

\[\mathbf{A}_{j,:} = \left(1, \; f_1 \cos \phi_{1,j}, \; f_1 \sin \phi_{1,j}, \; f_2 \cos \phi_{2,j}, \; f_2 \sin \phi_{2,j}, \; \ldots, \; f_N \cos \phi_{N,j}, \; f_N \sin \phi_{N,j} \right)\]

The least-squares solution minimises \(\|\boldsymbol{\varepsilon}\|^2\) and is given by the normal equations:

\[\mathbf{A}^T \mathbf{A} \, \mathbf{x} = \mathbf{A}^T \mathbf{h}\]

This is the formulation implemented in PyFES’s harmonic analysis routines. The mean sea level \(S_0\) is estimated together with the constituents, so that it does not bias them, but only the constituents are returned.

Complex Amplitude Representation#

The harmonic constants are conventionally expressed as an amplitude \(H_k\) and phase lag \(G_k\):

\[H_k = \sqrt{a_k^2 + b_k^2}, \qquad G_k = \operatorname{atan2}(b_k, a_k)\]

PyFES stores these as complex numbers:

\[Z_k = a_k + i \, b_k = H_k \, e^{i G_k}\]

This complex representation is used throughout the tidal model data structures (e.g., CartesianComplex128, LGP1Complex128).

The Rayleigh Criterion#

Two constituents of angular speeds \(\omega_1\) and \(\omega_2\) can be independently resolved from a record of duration \(T\) only if:

\[T > \frac{2\pi}{|\omega_1 - \omega_2|} \quad \Leftrightarrow \quad T > \frac{1}{|f_1 - f_2|}\]

where \(f_i = \omega_i / 2\pi\) are the cyclic frequencies. This is the Rayleigh criterion.

In practice, a Rayleigh factor \(C_R \geq 1\) (typically \(C_R = 1\)) is used, requiring \(T > C_R / |f_1 - f_2|\). As the record length increases, more constituent pairs become separable. For example:

  • A 15-day record separates the species (diurnal from semidiurnal) and only just resolves \(M_2\) from \(S_2\), which requires \(1 / |f_{S_2} - f_{M_2}| \approx 14.77\) days.

  • A 30-day record separates most constituents within a species except closely spaced pairs like \(K_1\) and \(P_1\).

  • A year-long record resolves the vast majority of significant constituents.

The select_waves_for_analysis() method on WaveTableInterface automatically applies this criterion to select the constituents resolvable for a given record duration.

Role of Nodal Corrections in Analysis#

The nodal corrections \(f_k\) and \(u_k\) appear explicitly in the design matrix \(\mathbf{A}\). If they are omitted or computed with insufficient accuracy, the recovered amplitudes and phases will be biased. For records shorter than a year, the nodal corrections are essentially constant and can be evaluated at the midpoint of the record. For multi-year records, time-varying corrections should be applied at each observation time.

References#

  • Simon, B. (2013). Marées Océaniques et Côtières (943-MOC), Ch. VI, pp. 126–158.

  • Schureman, P. (1940). Manual of Harmonic Analysis and Prediction of Tides, SP 98, pp. 49–100.