April 15, 2020
Author: Curtis Mobley
The representation of sea surfaces as sums of sinusoids suggests further analysis using Fourier transforms, which will be a fundamental mathematical tool for our subsequent study of sea surfaces. There are many texts on Fourier transforms. Bracewell (1986) is a standard reference, and excellent sets of lecture notes and videos of lectures can be found on the web. The results needed for the subsequent discussions will therefore be stated without proof.
Given a real function of a continuous variable , the Fourier transform of is deﬁned as
The inverse Fourier transform is given by
It can be shown that if we insert the integral of Eq. (1) into Eq. (2), we recover . (This is known as Fourier’s integral theorem, the proof of which is not trivial.) Equations (1) and (2) are termed a Fourier transform pair.
Understanding the units of a Fourier transform is important. In our case, is sea surface elevation, and both and have units of meters. Equation (1) shows that thus has units of . The variance of also has units of , which gives the ﬁrst hint at a profound connection between the variance of a physical quantity and its Fourier transform. The units of in the Fourier transform can be rewritten as m/(1/m), which is units of divided by units of spatial frequency . Indeed, the transform can be interpreted as showing ”how much of there is per unit frequency interval.” The inverse transform then has units of ( over frequency) times frequency, which returns the original units of .
A Fourier transform is a spectral density function. The integral of a spectral density function over a given frequency interval gives the variance in the physical quantity contributed by the frequencies in the integration interval. Density functions are rather peculiar mathematical creatures compared to point functions, which simply give the value of something at a given value of the independent variable (e.g. the temperature at location and time ). The blackbody radiation function is another example of a spectral density function. The blackbody function shows how much energy is emitted (at a given temperature) per unit frequency interval of the emitted electromagnetic radiation. The blackbody function is discussed on the page A Common Misconception. If you are not familiar with the distinction between point and density functions, especially regarding how to change variables in density functions, you should take a look at the page on blackbody spectra before continuing with the present discussion.
The Fourier transform deﬁnitions above with the in the exponents are those of the ”Stanford school” of Bracewell (1986) and Goodman (1996). You will see others in the literature. For example, if we use as the frequency variable, then Eqs. (1) and (2) become
This reappearance of the in the second equation is required so that the inverse transform of the transform gets you back to where you started. In the version, some people put the in front of the other integral, and some put a in front of each integral. Some authors, e.g. Press et al. (1992), put the in Eq. (1) and the in Eq. (2). The choice of which sign to use on the and where to put the is almost a religion—most people stay with what they ﬁrst learned, are convinced of the superiority of their deﬁnition, and are willing to die rather than change. Fortunately, it doesn’t matter which deﬁnitions you use, so long as you are consistent in how the transform pair is deﬁned so that you get back to where you started if you inverse transform a transform, or vice versa.
In our work, is the sea surface elevation, which is a real number. However, even though is real, (or ) is complex. Expanding the complex exponential in Eq. (1) as the sum of a cosine and a sine, it is easy to see that . Such functions are called Hermitian. A real function has a Hermitian Fourier transform. Conversely, if a function is Hermitian, it has a real inverse Fourier transform. The Hermitian property will be an important constraint on the next pages when we wish to generate random realizations of a sea surface by computing the inverse Fourier transform of a complex function: we will have to conjure up a Hermitian function so that we end up with a real sea surface.
Now suppose that we have sampled the sea surface at a set of evenly spaced points , in a region of size ; . We want to describe this sampled sea surface as a sum of sinusoids. In general, these values can be represented as a sum of a constant term, cosine terms, and sine terms (recall that there is no two-point sine term):
with . Note that the sum runs from the fundamental frequency () through the Nyquist frequency (), with only a cosine term for the two-point wave. This sum is equivalent to
which also contains a total of independent real and imaginary parts of the coeﬃcients. (Recall from Eq. (5) of the previous page that , so these coeﬃcients are not independent for and pairs.) These equations are the discrete-variable forms of Eqs. (6) and (7) of the previous page.
Comment on notation: It is common to use the letters for dummy summation indices. However, we’ve already used for and for angular wavenumber, so the preceding equations would be hopelessly confusing if we reused and for summation indices. We will therefore use and for indices on spatial variables, e.g., , and and for indices on frequency variables, e.g., or . and also will be used as needed for dummy indices.
We now have a ﬁnite number of discrete samples of the sea surface, so we need a discrete form of the Fourier transform. The discrete Fourier transform (DFT) of is deﬁned as
for . Recalling that and gives . It is also common to write and , in which case the previous equation becomes
The corresponding inverse discrete Fourier transform is given by
It was emphasized above that the continuous function deﬁned by the continuous Fourier transform (1) is a density function with units of per spatial frequency interval, e.g., m/(1/m) if is sea surface height in meters. However, the discrete function deﬁned by the discrete Fourier transform (7) has the same units as . The discrete Fourier transform is a point function that shows how much of is contained in a ﬁnite frequency interval centered at frequency . Discrete Fourier transforms convert point functions to point functions .
The and notations will be used to distinguish continuous vs. discrete Fourier transforms. As just seen, the continuous and discrete transforms are diﬀerent mathematical constructs with diﬀerent units and interpretations; they must not be confused. Likewise, if necessary, a subscript can be appended to show the frequency variable, e.g., as in Eq. (1) or as in Eq. (3). As always, there are competing deﬁnitions. Equations (7) and (8) are used in Bracewell and the IDL computer language. Numerical Recipes interchanges the and . Matlab puts the factor on the inverse transform. Also, Matlab does not support array indices of 0, so the summation indices are shifted from 0 to to 1 to , with a corresponding shift from to in the exponentials. The devil is in the details, and details like where to put the factor and diﬀerences in array indexing in diﬀerent computer languages can cause great misery when it comes time to actually write computer programs or to compare results computed by diﬀerent canned subroutines.
The sums in the last two equations require computing complex exponentials (i.e., sines and cosines), multiplying by the corresponding values of or and adding up the results. These equations can be evaluated for any value of . The number of computations required to do this is of order . The computation time thus increases very rapidly for large .
However, a classic paper by Cooley and Tukey (1965) showed how these sums can be computed in order computations, if is a power of 2. Their technique is now called the Fast Fourier Transform or FFT. The diﬀerence in computer time becomes enormous for large . For example, if , then , and the diﬀerence in computation time is a factor of . Thus in the case of , a roughly 6 minute computer run becomes a 1 second run. The development (or, perhaps, reinvention, since the basic idea goes all the way back to Gauss) of the FFT was a major advance in numerical analysis, which enables the computations on the following pages to be performed extremely eﬃciently. Subroutines for computing FFTs and inverse FFTs are available in all computer languages commonly used in science (Fortran, C, IDL, Matlab, etc). Fortunately we do not need to concern ourselves here with the details of how the FFT algorithm actually works, any more than we need to worry about how a canned subroutine actually computes the cosine of a number. If you are interested in how the FFT works, a web search will yield many detailed explanations. It is important to remember that the FFT is not another type of transform; the FFT is a numerically eﬃcient way to evaluate the DFT if the number of data values is a power of two.
The one-dimensional (1-D) equations seen above are easily extended to two or more dimensions. For two dimensions , we can sample a region of size by meters over points in the direction and points in the direction, with and both powers of 2 so we can use FFTs. Equations (5) and (6) then become
The corresponding 2-D DFT pair is
It will be important below to keep notational track of continuous vs. discrete versions of various functions. For any variable , will denote a continuous function of , will denote the continuous function evaluated at the discrete value , and will denote a discrete point function. Keep in mind that the density function and the point function have diﬀerent units. In the following pages and will denote the continuous and discrete versions, respectively, of wave elevation variance spectra.
The diﬀerences in units between continuous and discrete Fourier amplitudes sometimes makes it tricky to make the transition between discrete and continuous versions of the same quantity. In particular, it will be necessary to explicitly include the frequency intervals in some of the later calculations that involve both continuous and discrete variables. For example, if we have a continuous density function and we need to convert to a corresponding discrete function, we must multiple the continuous function by the ﬁnite frequency interval, e.g.
Conversely, if we have discrete amplitudes and we wish to estimate the continuous spectral density , then we must divide by the frequency interval:
(If you are an optical oceanographer familiar with the scattering phase function, you can ﬁnd an analogous situation in the estimation of the scattering phase function from measurements of scattered light. The scattering phase function is a measure of how much light energy is scattered from an incident direction into a particular direction, per unit solid angle; it therefore has units of 1/steradian. If you measure the scattered light using an instrument with a ﬁnite solid angle , then you get the total amount of energy scattered into the solid angle . To estimate the phase function from this measurement, you must divide the measured value by the solid angle of the instrument; this gets you back to units of 1/steradian.)
There is a peculiarity to most (perhaps all) FFT subroutines. The discrete FFT of Eq. (7) returns an array of complex numbers , which are associated with discrete spatial frequencies. What is peculiar is the order in which the elements of the array are returned by an FFT subroutine.
Let represent the fundamental frequency. If wavenumber is the frequency variable, then . If angular spatial frequency is the frequency variable, then ; for temporal angular frequency , . In any case, the discrete frequencies associated with the discrete Fourier amplitudes are equally spaced at intervals of and are in the negative-to-positive order
I’ll call this the ”math frequency order” because this is the natural order of arranging values in mathematics. This frequency order is convenient for plotting all of the amplitudes, as will be seen in the following pages. Plots showing both negative and positive frequencies are called ”two-sided,” and examples will be seen below.
However, FFT routines generally return their amplitudes corresponding to frequencies in the order of 0 ﬁrst, then the positive frequencies from the smallest to the Nyquist frequency, then the negative frequencies in reverse order:
I’ll call this the ”FFT frequency order.” Given the Hermitian symmetry of the amplitudes about the 0 frequency, the FFT order is convenient for plotting amplitudes only for the positive frequencies, with the negative-positive symmetry of understood. Such plots are called ”one-sided” or ”folded.”
Either frequency order can be obtained from the other by a repeated circular shift, which moves an array element oﬀ of one end of an array and copies it to the other end of the array, shifting all elements right or left by one position in the process. The detail to watch is that diﬀerent computer languages deﬁne a circular shift in diﬀerent ways. For example, the IDL routine SHIFT (and the Matlab routine CIRCSHIFT) moves the array elements to the right for a ”positive” shift (a negative shift moves elements to the left), whereas the Fortran 95 CSHIFT routine moves the array elements to the left for a positive shift (a negative shift moves elements to the right). Thus
- In IDL:
- math order = SHIFT(FFT order, N/2-1)
- FFT order = SHIFT(math order, -(N/2-1))
- In Fortran95:
- math order = CSHIFT(FFT order, -(N/2-1))
- FFT order = CSHIFT(math order, N/2-1)
Some FFT routines allow the user to set a ﬂag specifying which frequency order is to be returned. In any case, sorting out the frequency order of the amplitudes returned by a particular FFT routine, and ﬁguring out how to shift between math and FFT frequency orders in a particular computer language, can drive you to tears.
Note ﬁnally that in Eq. (7) you are simply providing an array of values and getting back the same number of values. What values correspond to the values is irrelevant. That is, might correspond to the spatial range from to , or from to , or to any other range. It is only the number of samples and the corresponding values that matters. In other words, the amplitudes coming out of the Fourier transform are independent of the origin of the spatial coordinate system. The frequencies depend on both the number of points and the physical size of the sampled region via the presence of (or time interval ) in the fundamental frequency . Thus the size of the region sampled and the number of samples, along with the sample values themselves, fully deﬁne the associated Fourier transforms.
The appearance of negative frequencies in Fourier transforms may seem somewhat mysterious. Frequency, after all, is a physical quantity that is inherently a positive number, e.g., the number of wavelengths in a given distance. However, one way to interpret mathematically negative frequencies is that they are simply the mathematical price we pay for the convenience of using complex numbers. Consider, for example, the representation of the cosine as a sum of complex exponentials for the frequency:
We can interpret the complex representation of the real cosine as having one term with a positive frequency and one term with a negative frequency . A similar equation holds for the complex representation of . The Fourier transform of a real function always contains both negative and positive frequencies, which arise from the complex exponentials in the deﬁnition of the transform.
Additional comments about negative frequencies will be made on later pages, where it will be seen that positive frequencies can be associated with waves propagating in the direction, and negative frequencies correspond to waves propagating in the opposite, , direction.
The physical and spectral variables of a continuous Fourier transform pair satisfy
The extension to the two-dimensional case is straightforward:
The discrete forms of Parseval’s relations provide important checks on numerical calculations. For example, it is easy to misplace factors of N, which appear in diﬀerent places depending on the exact form used for the deﬁnition of the discrete transforms.
You can take a class in Fourier transforms and prove many more pretty theorems about their properties. However, we have now assembled the mathematical tools needed to describe wind-blown sea surfaces, and so we can now get back to oceanography.