Page updated: May 23, 2021
Author: Curtis Mobley
View PDF

Surfaces to Spectra: 2D

The preceding pages have given a detailed look at one-dimensional surfaces. We now consider the more useful case of two-dimensional sea surfaces. The extension to two dimensions is mathematically straight forward.

Let z(x,t) = z(x,y,t) be the sea surface elevation in meters at point x = (x,y) at time t. The spatial extent of the sea surface is 0 x < Lx and 0 y < Ly. This surface is sampled on a rectangular grid of Nx by Ny points, where both Nx and Ny are powers of 2 for the FFT. The spatial sampling points are then

x(r) = [0, 1, 2,...,Nx 1] Lx Nx = rΔx,r = 0,...,Nx 1 y(s) = [0, 1, 2,...,Ny 1] Ly Ny = sΔy,s = 0,...,Ny 1.

This spatial sampling frequency gives Nx and Ny spatial frequencies kx and ky (in math frequency order)

kx(u) =[(Nx2 1),...,1, 0, 1,...,Nx2]2π Lx = uΔkx,u = (Nx2 1),...,Nx2 ky(v) =[(Ny2 1),...,1, 0, 1,...,Ny2]2π Ly = vΔky,v = (Ny2 1),...,Ny2 .

As before, the x-dimension 0 frequency is at array element u = Nx2 1 and the Nyquist frequency is at element u = Nx 1. Thus the positive and negative pairs of kx values are related by kx(u) = kx(Nx 2 u),u = 0,...,Nx 2, with the Nyquist frequency always being a special case because there is only a positive Nyquist frequency. Corresponding relations hold for the y direction.

Let k = (kx,ky) denote a spatial frequency vector, where kx and ky are frequencies in the x and y directions, respectively. For discrete values, we write kuv = (kx(u),ky(v)). The magnitude of k is k = kx 2 + ky 2. In our (x,y) coordinate system, let the wind blow in the + x direction. The x direction is then upwind, and the ± y directions are the cross-wind directions. With this choice, kx > 0 indicates frequencies of waves propagating more or less downwind, and kx < 0 for waves propagating against the wind. kx = 0 and ky0 indicates a wave propagating at exactly a cross-wind ± y direction. The angle of wave propagation relative to the downwind direction for a wave of frequency (kx(u),ky(v)) is given by

φ(kuv) = φ(u,v) = tan 1 ky(v) kx(u) . (1)

We can thus write the 2-D surface as z(x), its Fourier amplitude as (k), and the associated variance spectrum as Ψ(k). The discrete variance spectrum Ψ(kuv) = Ψ(u,v) gives the variance of the wave with wavelength 2πkuv propagating in direction φ(kuv) relative to the downwind direction.

Even though the mathematical transition from one to two dimensions causes no problems, it is again educational to take a careful look at a couple of contrived examples.

Example: A Random Sea Surface

For the first example, a sea surface area of size Lx × Ly = 10 × 10m is sampled using Nx = 16 and Ny = 8 points in the x and y directions. Sea surface elevations z(x,y) were created by drawing a 𝒩(0, 1) random number at each (x,y) value.

The upper-left panel of Fig. 1 shows a contour plot of this surface for a particular sequence of random numbers. The spatial periodicity of the Fourier representation is used to extend the contour plot to the full (Lx,Ly) range of the tile, which gives a good visual appearance. Thus the elevations at x = Lx are the same as those at x = 0, those at y = Ly equal those at y = 0, and the elevation at (Lx,Ly) duplicates the elevation at (0, 0). The Nx × Ny grid of sample points is shown by the solid silver dots. The points at x = Lx and y = Ly obtained by periodicity are shown by open silver circles at the right and top of the surface plot.


Figure 1: A two-dimensional random sea surface and its Fourier transform and two-sided variance spectrum. The surface elevations in the upper left panel were randomly drawn from a 𝒩(0, 1) distribution. The black arrows highlight a particular ±kuv frequency pair. The blue-blue and red-red symmetry of the real part of the Fourier amplitudes, and the red-blue symmetry of the imaginary part, shows the Hermitian nature of the amplitudes. The gray dots show the locations of the discrete values that were contoured to create the figures.

The Fourier amplitudes (kuv) = (kx(u),ky(v)) are obtained from the 2-D DFT of z(xrs) = z(x(r),y(s)):

(kuv) = 𝒟{z(xrs)}.

The usual warning on FFT frequency order applies here. The 2-D FFT gives back the amplitudes (kuv) with the array elements corresponding the the FFT frequency order of Eq. (15) of the Fourier Transforms page. Before plotting, the (kuv) array elements must be shifted into the math frequency order in both the kx and ky array directions using the appropriate shift function for the computer language used in the plotting. For the IDL routine used to generate Fig. 1, the 2-D circular shift is given by the IDL command

realzhatplot = SHIFT(REAL_PART(zhat),Nx2 1,Ny2 1)

where zhat is the complex 2-D array returned by the FFT routine, and realzhatplot is the 2-D array plotted in the upper right panel of Fig. 1.

The upper right panel of Fig. 1 plots the real part of (kuv), and the lower left panel plots the imaginary part. The Nyquist frequency kxNy = 2π(2Δx) = 5.03radm lies along the right side of the contour plot. The Nyquist frequency kyNy = 2π(2Δy) = 2.51radm lies along the top of the contour plot. The white space at the left and bottom highlights that there are no negative Nyquist frequencies. In each amplitude plot a particular pair of ±kuv values is indicated by the black arrows; the (0, 0) frequency is shown by a black dot. Note that in the plot of the real part, Re{(kuv)} = Re{(+kuv)}, whereas in the plot of the imaginary part, Im{(kuv)} = Im{(+kuv)}. The contouring is rather low quality for so few points, but it is easy to see in the digital output that when the kx is the Nyquist frequency (the points along the right column of points in the plot), the symmetries are given by Re{(kxNy,ky)} = Re{(k xNy, +ky)} and Im{(kxNy,ky)} = Im{(k xNy, +ky)}. A corresponding relation holds for ± kx when ky = kyNy (the points along the top row of the plot). The point at the upper right of the plot corresponds to both kx and ky being at their respective Nyquist frequencies. As always, the array elements at the Nyquist frequencies must be treated as special cases when writing computer programs. These symmetries show that the 2-D amplitudes are Hermitian: (k uv) = (kuv). The discrete 2-D variance spectrum Ψ(kuv) = Ψ(u,v) = (u,v)(u,v) is contoured in the lower right panel of the figure. Note the Ψ(kuv) = Ψ(kuv) symmetry.

For this example, the standard check on the 2-D discrete Parseval’s relation, Eq. (16) on page Fourier Transforms, gives

r=0Nx1 s=0Ny1|z(r,s)|2 = N xNy u=0Nx1 v=0Ny1|(u,v)|2 = NxNy u=0Nx1 v=0Ny1Ψ(u,v) = 130.90m2.

In all of these plots, it should be remembered that the discrete values are known only at the locations of the silver dots. The contouring routine simply interpolates between these points to create a visually appealing figure. The continuous color in the plots does not imply that the values are continuous and known in between the discrete points.

Example: A Sea Surface of Crossing Sinusoids

For a second example, let us define a surface from a pair of crossing sinusoids as follows:

z(xr,ys) = A1 cos(kx1xr + ky1ys + ϕ1) + A2 cos(kx2xr + ky2ys + ϕ2), (2)

where as usual r = 0,...,Nx 1; s = 0,...,Ny 1; and where

A1 = 1.0m is the amplitude of the first wave
Nx1 = 2 is the number of wave lengths in the x direction in [0,Lx] for the first wave
kx1 = 2.0πNx1Lx = 1.257radm2 is the kx frequency of the first wave
Ny1 = 1 is the number of wave lengths in the y direction in [0,Ly] for the first wave
ky1 = 2.0πNy1Ly = 0.628radm2 is the ky frequency of the first wave
ϕ1 = 0 is the phase of the first wave; 0 gives a cosine wave
A2 = 0.5m is the amplitude of the second wave
Nx2 = 4 is the number of wave lengths in the x direction in [0,Lx] for the second wave
kx2 = 2.0πNx2Lx = 2.513radm2 is the kx frequency of the second wave
Ny2 = 3 is the number of wave lengths in the y direction in [0,Ly] for the second wave
ky2 = 2.0πNy2Ly = 1.885radm2 is the ky frequency of the first wave
ϕ2 = π2 is the phase of the second wave; π2 gives a sine wave

The wavelength of the first wave is Λ1 = 2πkx1 2 + ky1 2 = 4.47m, and that of the second wave is Λ2 = 2.00m. The direction of propagation of the first wave relative to the + x axis is

φ1 = tan 1 ky1 kx1 = 26.57deg,

and the direction of propagation of the second wave relative to the + x axis is

φ2 = tan 1 ky2 kx2 = 36.87deg.

The upper left panel of Fig. 2 shows this surface elevation pattern when Eq. (2) is sampled with Nx = Ny = 16. The dominant red-blue pattern shows the first wave oriented with the direction of propagation along either the + k1 direction at φ1 = 26.57deg or the k1 direction at φ1 = 26.57 + 180 = 206.57deg. We cannot of course determine the actual direction + k1 or k1 of propagation from sea surface elevations at a single time. The dominant wave pattern is modulated by the second wave, which has one half the amplitude of the first wave.


Figure 2: A two-dimensional sea surface composed of two crossing sinusoids, and the resulting Fourier amplitudes and variance.

The choice above of ϕ1 = 0 in Eq. (2) makes the dominate wave a cosine in our coordinate system, and ϕ2 = π2 makes the second wave a sine. As we saw in the 1-D examples, variance associated with cosine waves appears in the real part of the Fourier amplitudes, and the variance in sine waves appears in the imaginary parts. We see this again here for the first (cosine) and second (sine) waves of the surface. Since each wave pattern has only one frequency, there is only one pair of points at ±k1 in the plot of the real part, and one pair at ±k2 in the plot of the imaginary part. The amplitudes at all other frequencies are zero. The symmetries of these points again show the Hermitian nature of the amplitudes.

The lower right panel shows the two-sided variance function Ψ2s(k) for this surface. Note that the first wave has four times the variance of the second wave because the amplitude of the first wave is twice that of the second wave. Ψ2s(k) = Ψ2s(+k). Note also that you can look at a 2-D variance spectrum and see how much energy is propagating in a given ±k direction.

For this simple example involving just four frequencies it is easy (from the digital output) to hand check that the right-hand side of the 2D discrete Parseval’s relation, Eq. (18) of the Fourier Transforms page, is

NxNy u v|(u,v)|2 = 1616 (0.5)2 + (0.5)2 + (0.25)2 + (0.25)2 = 160m2.

This value agrees exactly with the corresponding sums of the surface elevations, r s[z(r,s)]2, and variance values, NxNy u vΨ(u,v).

Comments for Surfaces to Spectra: 2D:

Loading Conversation