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\left(x,t\right)=z\left(x,y,t\right)$ be the sea surface elevation in meters at point $x=\left(x,y\right)$ at time $t$. The spatial extent of the sea surface is $0\le x<L_x$ and $0\le y<L_y$. This surface is sampled on a rectangular grid of $N_x$ by $N_y$ points, where both $N_x$ and $N_y$ are powers of 2 for the FFT. The spatial sampling points are then

This spatial sampling frequency gives $N_x$ and $N_y$ spatial frequencies $k_x$ and $k_y$ (in math frequency order)

As before, the x-dimension 0 frequency is at array element $u=N_x∕2-1$ and the Nyquist frequency is at element $u=N_x-1$. Thus the positive and negative pairs of $k_x$ values are related by $k_x\left(u\right)=-k_x\left(N_x-2-u\right),u=0,...,N_x-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=\left(k_x,k_y\right)$ denote a spatial frequency vector, where $k_x$ and $k_y$ are frequencies in the $x$ and $y$ directions, respectively. For discrete values, we write $k_{uv}=\left(k_x\left(u\right),k_y\left(v\right)\right)$. The magnitude of $k$ is $k=\sqrt{k_x^2+k_y^2}$. In our $\left(x,y\right)$ coordinate system, let the wind blow in the $+x$ direction. The $-x$ direction is then upwind, and the $\pm y$ directions are the cross-wind directions. With this choice, $k_x>0$ indicates frequencies of waves propagating more or less downwind, and $k_x<0$ for waves propagating against the wind. $k_x=0$ and $k_y\ne 0$ indicates a wave propagating at exactly a cross-wind $\pm y$ direction. The angle of wave propagation relative to the downwind direction for a wave of frequency $\left(k_x\left(u\right),k_y\left(v\right)\right)$ is given by

We can thus write the 2-D surface as $z\left(x\right)$, its Fourier amplitude as $ẑ\left(k\right)$, and the associated variance spectrum as $\Psi \left(k\right)$. The discrete variance spectrum $\Psi \left(k_{uv}\right)=\Psi \left(u,v\right)$ gives the variance of the wave with wavelength $2\pi ∕k_{uv}$ propagating in direction $\phi \left(k_{uv}\right)$ 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 $L_x\times L_y=10\times 10m$ is sampled using $N_x=16$ and $N_y=8$ points in the x and y directions. Sea surface elevations $z\left(x,y\right)$ were created by drawing a $\mathcal{𝒩}\left(0,1\right)$ random number at each $\left(x,y\right)$ 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 $\left(L_x,L_y\right)$ range of the tile, which gives a good visual appearance. Thus the elevations at $x=L_x$ are the same as those at $x=0$, those at $y=L_y$ equal those at $y=0$, and the elevation at $\left(L_x,L_y\right)$ duplicates the elevation at $\left(0,0\right)$. The $N_x\times N_y$ grid of sample points is shown by the solid silver dots. The points at $x=L_x$ and $y=L_y$ obtained by periodicity are shown by open silver circles at the right and top of the surface plot.

The Fourier amplitudes $ẑ\left(k_{uv}\right)=ẑ\left(k_x\left(u\right),k_y\left(v\right)\right)$ are obtained from the 2-D DFT of $z\left(x_{rs}\right)=z\left(x\left(r\right),y\left(s\right)\right)$:

The usual warning on FFT frequency order applies here. The 2-D FFT gives back the amplitudes $ẑ\left(k_{uv}\right)$ with the array elements corresponding the the FFT frequency order of Eq. (15) of the Fourier Transforms page. Before plotting, the $ẑ\left(k_{uv}\right)$ array elements must be shifted into the math frequency order in both the $k_x$ and $k_y$ 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

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 $ẑ\left(k_{uv}\right)$, and the lower left panel plots the imaginary part. The Nyquist frequency $k_x^{Ny}=2\pi ∕\left(2\Delta x\right)=5.03rad∕m$ lies along the right side of the contour plot. The Nyquist frequency $k_y^{Ny}=2\pi ∕\left(2\Delta y\right)=2.51rad∕m$ 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 $\pm k_{uv}$ values is indicated by the black arrows; the $\left(0,0\right)$ frequency is shown by a black dot. Note that in the plot of the real part, $Re\left\{ẑ\left(-k_{uv}\right)\right\}=Re\left\{ẑ\left(+k_{uv}\right)\right\}$, whereas in the plot of the imaginary part, $Im\left\{ẑ\left(-k_{uv}\right)\right\}=-Im\left\{ẑ\left(+k_{uv}\right)\right\}$. The contouring is rather low quality for so few points, but it is easy to see in the digital output that when the $k_x$ is the Nyquist frequency (the points along the right column of points in the plot), the symmetries are given by $Re\left\{ẑ\left(k_x^{Ny},-ky\right)\right\}=Re\left\{ẑ\left(k_x^{Ny},+ky\right)\right\}$ and $Im\left\{ẑ\left(k_x^{Ny},-ky\right)\right\}=-Im\left\{ẑ\left(k_x^{Ny},+ky\right)\right\}$. A corresponding relation holds for $\pm k_x$ when $k_y=k_y^{Ny}$ (the points along the top row of the plot). The point at the upper right of the plot corresponds to both $k_x$ and $k_y$ 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: ${ẑ}^{\ast }\left(-k_{uv}\right)=ẑ\left(k_{uv}\right)$. The discrete 2-D variance spectrum $\Psi \left(k_{uv}\right)=\Psi \left(u,v\right)=ẑ\left(u,v\right){ẑ}^{\ast }\left(u,v\right)$ is contoured in the lower right panel of the figure. Note the $\Psi \left(-k_{uv}\right)=\Psi \left(k_{uv}\right)$ symmetry.

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

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:

where as usual $r=0,...,N_x-1;s=0,...,N_y-1;$ and where

$A_1=1.0m$ is the amplitude of the first wave

$N_{x1}=2$ is the number of wave lengths in the x direction in $\left[0,L_x\right]$ for the first wave

$k_{x1}=2.0\pi N_{x1}∕L_x=1.257rad∕m^2$ is the $k_x$ frequency of the first wave

$N_{y1}=1$ is the number of wave lengths in the y direction in $\left[0,L_y\right]$ for the first wave

$k_{y1}=2.0\pi N_{y1}∕L_y=0.628rad∕m^2$ is the $k_y$ frequency of the first wave

${\varphi }_1=0$ is the phase of the first wave; 0 gives a cosine wave

$A_2=0.5m$ is the amplitude of the second wave

$N_{x2}=4$ is the number of wave lengths in the x direction in $\left[0,L_x\right]$ for the second wave

$k_{x2}=2.0\pi N_{x2}∕L_x=2.513rad∕m^2$ is the $k_x$ frequency of the second wave

$N_{y2}=3$ is the number of wave lengths in the y direction in $\left[0,L_y\right]$ for the second wave

$k_{y2}=-2.0\pi N_{y2}∕L_y=-1.885rad∕m^2$ is the $k_y$ frequency of the first wave

${\varphi }_2=\pi ∕2$ is the phase of the second wave; $\pi ∕2$ gives a sine wave

The wavelength of the first wave is ${\Lambda }_1=2\pi ∕\sqrt{k_{x1}^2+k_{y1}^2}=4.47m$, and that of the second wave is ${\Lambda }_2=2.00m$. The direction of propagation of the first wave relative to the $+x$ axis is

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

The upper left panel of Fig. 2 shows this surface elevation pattern when Eq. (2) is sampled with $N_x=N_y=16$. The dominant red-blue pattern shows the first wave oriented with the direction of propagation along either the $+k_1$ direction at ${\phi }_1=26.57deg$ or the $-k_1$ direction at ${\phi }_1=26.57+180=206.57deg$. We cannot of course determine the actual direction $+k_1$ or $-k_1$ 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.

The choice above of ${\varphi }_1=0$ in Eq. (2) makes the dominate wave a cosine in our coordinate system, and ${\varphi }_2=\pi ∕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 $\pm k_1$ in the plot of the real part, and one pair at $\pm k_2$ 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 ${\Psi }_{2s}\left(k\right)$ 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. ${\Psi }_{2s}\left(-k\right)={\Psi }_{2s}\left(+k\right)$. Note also that you can look at a 2-D variance spectrum and see how much energy is propagating in a given $\pm 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

This value agrees exactly with the corresponding sums of the surface elevations, ${\sum }_r{\sum }_s{\left[z\left(r,s\right)\right]}^2$, and variance values, $N_x{N}_y{\sum }_u{\sum }_v\Psi \left(u,v\right)$.