Wave Spectrum in Wolfram Language

2025-01-30 • Last updated: 2025-03-21 • Posted in: Code, Wolfram Language

Warning/Update I published a post with an improvement and correction of the code below. The new post is Welch Periodogram in Wolfram Language. I’ll keep the old code below for historical reasons.

Surface waves in water bodies are usually described with the wave spectrum. This spectrum is erroneously referred to as a wave energy spectrum, but it is really a variance density spectrum. It describes how the variance of the water surface elevation is distributed over frequencies (and, optionally, directions). The plot has units of m²/Hz versus Hz, or equivalent.

Formally, water surface elevation is regarded as a stochastic process, and the wave spectrum is defined as the Fourier transform of its autocorrelation function. In practice, however, computing the autocorrelation function is expensive, so we define the spectrum with the Fourier transform of the surface elevation and apply the efficient Fast Fourier Transform (FFT) algorithm. This way, the wave spectrum is essentially a periodogram.

Since I acquired a Mathematica student license, I’m trying to do all of my Ph.D. analyses with the Wolfram Language. I was trying to compute the wave spectrum, but I found that the Fourier and PeriodogramArray functions don’t compute the wave spectrum directly. In Stanisław Massel’s Ocean Surface Waves (2017; DOI 10.1142/10666), the discrete Fourier transform, $X_{ℓ}$ , of a water surface elevation record, $ζ_{k}$ , of length $N$ sampled at a frequency $f_{s}$ is given by the expression below.

X_{ℓ} = \frac{1}{f_{s}} \sum_{k = 0}^{N - 1} ζ_{k} \exp (- 2 π i \frac{k ℓ}{N})

However, function Fourier in Wolfram Language has a different definition of the discrete Fourier transform $v_{ℓ}$ of a list $u_{k}$ of length $N$ :

v_{ℓ} = \frac{1}{N^{(1 - a) / 2}} \sum_{k = 1}^{N} u_{k} \exp (2 π i b \frac{(k - 1) (ℓ - 1)}{N})

By default, $(a, b) = (0, 1)$ , but these values can be changed to compute alternative definitions. Also note that this definition is independent of the sampling frequency.

I did some algebra and found that redefining the indices in Massel’s formulation results in Wolfram Language definition with $(a, b) = (1, - 1)$ and divided by $f_{s}$ :

X_{ℓ} = \frac{1}{f_{s}} \sum_{k = 1}^{N} ζ_{k} \exp (- 2 π i \frac{(k - 1) (ℓ - 1)}{N})

In addition, function PeriodogramArray returns the squared magnitude of the discrete Fourier transform, i.e., $| v_{ℓ} |^{2}$ , but the wave spectrum, $S_{ℓ}$ , is normalized as shown below.

S_{ℓ} = S (f_{ℓ}) = 2 \frac{f_{s}}{N} | X_{ℓ} |^{2}, ℓ = 0, 1, 2, \dots, \frac{N}{2}

The corresponding frequencies are:

f_{ℓ} = \frac{ℓ f_{s}}{N}

Note Massel’s equation 16.28 is actually for ordinary frequency, $f$ , even though it is written in terms of angular frequency, $ω$ . This was probably a typesetting error. If a spectrum in terms of $ω$ is desired, one must take $ω = 2 π f$ and $S (ω) = S (f) / (2 π)$ .

Wolfram Language code

In order to compute $| X_{ℓ} |^{2}$ with PeriodogramArray, the following expression must be evaluated for surface elevation list zeta and sampling frequency fs.

1/fs^2 PeriodogramArray[zeta, FourierParameters -> {1, -1}]

Then, if n is the length of zeta, we get the wave spectrum with the expression below.

2/(n fs) PeriodogramArray[zeta,
	FourierParameters -> {1, -1}][[1 ;; n/2 + 1]]

However, since n is already in Fourier definition, taking $(a, b) = (0, - 1)$ simplifies the expression.

2/fs PeriodogramArray[zeta,
	FourierParameters -> {0, -1}][[1 ;; n/2 + 1]]

Finally, an enhanced wave spectrum S can be computed as the mean of the spectra of partitions of length m overlapped by m/2 (Welch’s method). A Hann smoothing window can be applied to reduce spectral leakage.

S = 2/fs PeriodogramArray[zeta, m, m/2, HannWindow,
	FourierParameters -> {0, -1}][[1 ;; m/2 + 1]]

Frequencies are simply:

f = fs/m Range[0, m/2]

Example

Let’s apply this Wolfram Language code to some of my laboratory data. I have a record of water surface elevation from a wind wave experiment in a small acrylic tank. Data were sampled at 200 Hz for 2 minutes.

First, we detrend the data. We can fit a linear model to times, t, and trended elevations, zetaRaw. The linear model is subtracted from the raw data to get detrended elevations.

lm = LinearModelFit[Transpose[{t, zetaRaw]}], x, x];
zeta = zetaRaw - (lm /@ t)

Then, we compute frequencies, f, and spectrum, S, with the code from previous section and partition size of 2⁹ or 512 (according to documentation, Fourier is faster with powers of 2). A plot of the resulting spectrum is shown below.

Plot of a wave spectrum — Wave spectrum from laboratory data.

Finally, we can verify (approximately!) that the zeroth moment of the spectrum, $\int_{0}^{\infty} S (f) d f$ , equals the variance of the surface elevation. If we take the square root, we obtain a metric of wave height in meters.

{f, S} //
	Transpose //
	Interpolation //
	Integrate[#[x], {x, 0, 20}] & //
	Sqrt

(* 0.000560723 *)

StandardDeviation[zeta]

(* 0.00055394 *)