Page updated: March 22, 2021
Author: Curtis Mobley
View PDF

Autocovariance Functions: Theory

The preceding pages have described the statistical properties of random sea surfaces using variance spectra. This page begins the exploration of an alternate path to the specification of surface roughness properties. Clearly, wind-dependent variance spectra are applicable only to surface waves that are generated by wind. Consider, however, the surface of a flowing river. The river’s surface can have ripples or waves generated by turbulence resulting from unstable shear flow induced by flow over a shallow bottom, or by eddies created as the water flows around rocks in the river. These water surfaces do not depend on the wind speed and can have different statistical properties, hence different optical properties, than wind-roughened water surfaces. Such surfaces can be described by their autocovariance functions.

Autocovariance functions can be converted to elevation variance spectra via the Wiener-Khinchin theorem. This page shows how that is done. Once a given autocovariance function has been converted to its equivalent elevation variance spectrum, the algorithms of the previous chapters are immediately applicable, even though the variance spectrum is not wind-dependent. Indeed, this conversion enables the Fourier transform methods of the previous chapters to be used to generate random realizations of any surface, not just water surfaces.

As is often the case, there is a large gap between textbook theory—usually developed for continuous variables or an infinite sample size of discrete values—and its implementation in a computer program for a finite sample size of discrete variables. In particular, careful attention must be paid to sampling of an autocovariance function in order to obtain the corresponding variance spectrum, or vice versa. I find it disappointing and frustrating (but not surprising) that numerical matters such as the effects of finite sample size, maximum lag size, and exactly how to sample spectra or autocovariances (in particular, the material on page Autocovariance Functions: Sampling) never seems to be mentioned in textbooks on digital signal processing or related subjects. It is left to the innocent student to spend a few weeks of unfunded time figuring out why various numerical results are not internally consistent or do not perfectly match the textbook theory.

This page begins with a discussion of autocovariances, and then the Wiener-Khinchin theorem is stated. The theorem is numerically illustrated on the next page first using the wind-dependent Pierson-Moskowitz elevation variance spectrum, for which certain values can be analytically calculated and used to check the numerical results. The modeling of water surfaces generated by shear-induced turbulence is then illustrated, again using analytical functions that allow for a rigorous check on the numerical results.


The autocovariance Czz() of z(r) is defined as

Czz() {z(r)z(r + )}, (1)

where denotes statistical expectation and is the spatial lag. This definition in terms of the expectation holds for both continuous and discrete variables. For the present discussion with z(r) being sea surface elevation, Czz() shows how strongly the sea surface elevation at one location is correlated to the elevation at a distance away. Czz() has units of m2, and Czz( = 0) is the variance of the surface elevation. The autocovariance is an even function of the lag: Czz() = Czz().

Consider an infinite sample of discrete zero-mean surface elevations spaced a distance Δx apart. The autocovariance is then computed by (e.g., Eq. 2.6.3 of Proakis and Manolakis (1996) with a minor change in their notation)

Czz() = r=+z(r)z(r + ),for = 0,±1,±2,....

Here is indexing the lag distance in units of the sample spacing Δx. For a finite sample of N discrete values, the same authors define the sample autocovariance by (their Eq. 2.6.11)

Czz() = r=0N||1z(r)z(r + ). (2)

As usual, there are competing definitions. For a finite sample of N discrete values, perhaps with a non-zero mean m, the IDL autocorrelation function (A_CORRELATE) uses

Czz(r) = 1 N r=0N||1[z(r)m][z(r+)m]for(N2) r N2, (3)


m = 1 N r=0Nz(r)

is the sample mean. Matlab computes the autocovariance via

Czz() = 1 N 1 r=0N||1[z(r) m][z(r + ) m].

Note that the lag must be less than the length of the sample. (That is, the sample locations are at xr = rΔx,r = 0,...,N 1 and the allowed lag distances are r = rΔx,r = 0,...,N 2.) Note also the factor of 1N in the IDL definition (3), which does not appear in Eq. (2), and which is a factor of 1(N 1) in the Matlab version.

Nor is there even any consensus on the terms “autocovariance” and “autocorrelation.” Some authors (and this page) define the nondimensional autocorrelation ρzz() as the autocovariance normalized by the variance, i.e.

ρzz() Czz() Czz(0). (4)

However, Proakis and Manolakis (1996) call the autocovariance as used here the autocorrelation, and they call the autocorrelation of Eq. (4) the “normalized autocorrelation.” These sorts differences in the definitions and computations of autocovariances can can cause much grief when comparing the numerical outputs of different computer codes, or numerical outputs with textbook examples.

The Wiener-Khinchin Theorem

Now that the autocovariance has been defined, the Wiener-Khinchin theorem can be stated: The Fourier transform of the autocovariance equals the variance spectral density function. Symbolically,

ν{Czz()} = 𝒮2s(ν). (5)

Indeed, some texts define the spectral density as the Fourier transform of the autocovariance. The inverse is of course

ν1{𝒮 2s(ν)} = Czz(). (6)

Here 𝒮2s is a two-sided spectral density function as discussed on several previous pages (e.g., Surfaces to Spectra: 1D).

It is important to note (as emphasized by the ν subscript on the Fourier transform operator ) that the theorem is written with the ν version of the Fourier transform (Eq. 1 of the Fourier Transforms page), and the density function is written in terms of the spatial frequency ν, which has units of 1/meters. (In the time domain, the conjugate variables are time t in seconds and frequency f in 1/seconds = cycles/second = Hz.) The spectral density function for the angular spatial frequency k = 2πν (or angular temporal frequency ω = 2πf in radians per second in the time domain) can be obtained by noting that corresponding intervals of the spectral densities contain the same amount of variance:

𝒮2s(k)dk = 𝒮2s(ν)dν,

which gives

𝒮2s(k) = 1 2π𝒮2s(ν = k2π). (7)

Note that varies from to + and, likewise, ν or k ranges over all negative and positive values. The variance spectrum obtained from the Fourier transform of an autocovariance function is therefore a two-sided spectrum.

Comment: In light of Eq. (7), the theorem stated in terms of angular spatial frequency k appears to be

k{Czz()} = 2π𝒮2s(k), (8)

with the inverse

k1{𝒮 2s(k)} = 1 2πCzz(). (9)

I say “appears to be” because I’ve never actually seen the theorem written this way because the textbooks all seem to stick with x and ν (or t and f in the time domain). As Press et al. (1992) say in Numerical Recipes (p. 491), “There are fewer factors of 2π to remember if you use the (ν or) f convention, especially when we get to discretely sampled data....” In any case, Eqs. (8) and (9) are consistent with the k spectrum of Eq. (17) discussed below.

The theorem is usually proved in textbooks for continuous variables x and ν. However, in numerical application to a finite number of discrete samples, discrete variables xr and νu or ku are used, and proper attention must be paid to pesky factors of N, 2π, and bandwidth, and to the array ordering required by a particular FFT routine.

The continuous-variable Fourier transform in Eq. (5) gives a spectral density 𝒮2s(ν) with units of m2(1m). However, if the theorem is written for a DFT of discrete data Czz(r) = Czz(r),

𝔇ν{Czz(r)} = 𝒮2s(νu), (10)

then the resulting discrete spectrum 𝒮2s(νu) = 𝒮2s(u) has units of m2. Just as was discussed on the Fourier Transforms page, this discrete spectrum must be divided by the bandwidth Δν in order to obtain the density at ν = νu. That is,

𝒮2s(ν = νu) = 𝒮2s(u)Δν. (11)


Horoshenkov et al. (2013) (their Eq. 4) give an analytic formula for the autocorrelation function of surface waves generated by bottom-induced turbulence in shallow flowing water. In the notation of this page, this function is

ρzz() = exp 2 2σw2 cos 2π Lo. (12)

In their words, “σw relates to the spatial radius of correlation (correlation length)” and “Lo relates to the characteristic period in the surface wave pattern.” The average values for the physical conditions of the Horoshenkov et al. study are σw = 0.22m and Lo = 0.17m. [Note: Eq. (12) is Horoshenkov’s W(ρ) as shown in their abstract and in their conclusions, where it has a factor of 1/2 in the exponential. Their Eq (4) does not have the 1/2. This is probably a typo since Gaussians usually have the form exp[x2(2σ2)].]

This autocorrelation function provides a nice example of how to use the Wiener-Khinchin theorem to obtain the corresponding variance spectrum. Equation (12), when converted to an autocovariance via a factor of the variance, z2 = C zz(0), has the form

Czz() = Czz(0) exp 2 2σw2 cos 2π Lo = Czz(0) exp a22 cos(q o), (13)

where a = 1(2σw) and qo = 2πLo. This function has an easily computed analytical Fourier transform.

The continuous-variable Fourier transform of this Czz() is given by Eq. (1) of the Fourier Transforms page:

𝒮2s(ν) = ν{Czz()} =C zz()ei2πνd. (14)

Here and ν are continuous variables; 𝒮2s(ν) has units of m3, which is interpreted as m2(1m) as explained before. Note that this variance spectral density is two-sided, i.e., < ν < . Expanding the complex exponential via ei𝜃 = cos 𝜃 i sin 𝜃 gives

𝒮2s(ν) = C zz(0) exp a22 cos(q o) cos(2πν)d iC zz(0) exp a22 cos(q o) sin(2πν)d.

The imaginary term is zero because the integrand is an odd function of . Using the identity

cos x cos y = 1 2 cos(x + y) + cos(x y)


𝒮2s(ν) = 2Czz(0)0 exp a22 1 2 cos[(2πν + qo)] + cos[(2πν qo)] d.

The integral

0 exp a22 cos(b)d = π exp[b2(4a2)] 2a

gives the Fourier transform of the Czz() of Eq. (13):

𝒮2s(ν) = π 2σwCzz(0) exp 1 2(2πσw)2(ν + 1L o)2 + exp 1 2(2πσw)2(ν 1L o)2 . (15)

This variance spectrum has maxima at ν = ±1Lo, where the value is very close to π 2 σwCzz(0). Figure 1 plots this Czz() (Eq. 13) and 𝒮2s(ν) (Eq. 15) for typical values of σw = 0.22m, Lo = 0.17m, and Czz(0) = 2.5 × 107m2. Note that the sub peaks of the autocovariance lie at integer multiples of ± Lo, and that the peaks of the spectrum are at ± 1Lo.


Figure 1: The Horoshenkov autocovariance Czz() and elevation variance spectral density 𝒮2s(ν) for typical parameter values taken from their Table 3.

By definition, the integral over all frequencies of an elevation variance spectral density gives the total elevation variance z2:

z2 =𝒮 2s(ν)dν.

This can be computed analytically for the spectrum of Eq. (15). The 𝒮2s(ν) spectrum of Eq. (15) is the sum of two identical Gaussians centered at different ν values. Consider the one centered at ν = 1Lo, which involves the exponential with the ν 1Lo term. The total variance is then twice the integral of this Gaussian:

z2 = 2π 2σwCzz(0) exp 1 2(2πσw)2(ν 1L o)2 dν.

Letting x = ν 1Lo gives

z2 = 2π 2σwCzz(0) exp c2x2 dx,

where c2 = 1 2(2πσw)2. The integral

0 exp c2x2 dx = π 2c

then gives the final result:

z2 = 4π 2σwCzz(0) π 2 1 22πσw = Czz(0). (16)

Thus starting with a variance of Czz(0) in the autocovariance of Eq. (13), obtaining the variance spectral density from the Fourier ν-transform of the autocovarience, and then integrating the variance spectral density over ν thus gives back the variance as originally specified in the autocovariance function.

However, if the above process is naively carried through starting (as in Eq. 14) with the k-transform of Eq. (3) on the Fourier Transforms page, the end result (as in Eq. 16) is 2πCzz(0). This extra factor of 2π is rectified by the 12π factor seen in Eq. (7). Thus the k version of the Horeshenkov spectral density is

𝒮2s(k) = 1 2ππ 2σwCzz(0) exp 1 2σw2(k + q o)2 + exp 1 2σw2(k q o)2 . (17)

Integration of this 𝒮2s(k) over all k then results in z2 = C zz(0), as required. The inverse k transform of the spectral density (17) gives Czz()(2π), as expected from Eq. (9).

Comments for Autocovariance Functions: Theory:

Loading Conversation