Abstract

Orthonormal wavelet analysis (OWA) is a special form of wavelet analysis, especially suitable for analyzing spatial structures, such as atmospheric fields. For this purpose, OWA is much more efficient and accurate than the nonorthogonal wavelet transform (WT), which was introduced to the meteorological community recently and which is more suitable for time series analysis. Whereas the continuous WT is strictly correct only for infinite domains, OWA is derived from periodizing and discretizing the infinite-domain case and so is correct for periodic boundary conditions.

Unlike Fourier spectra, OWA is not shift invariant. Nor is it equivariant like the WT; that is, the OWA output does not shift as its input shifts. Two remedies are to combine all possible shifts, known as the overcomplete, nonorthogonal shift equivariant WT, or else to use a “best shift,” known as best shift wavelet analysis. Although shift invariant and orthonormal w.r.t. arbitrary inputs, the latter’s optimization generally depends nonlinearly on a reference structure.

OWA is a generalization of Fourier series, and the associated multiresolution analysis (MRA) generalizes Reynolds averaging. Like these, OWA/MRA on discrete and continuous domains satisfy analogous identities arithmetically exactly, unlike the WT. OWA/MRA is much more efficient than Fourier series for analyzing several examples, including nearly singular synthetic structures, and a 90-day Hovmöller diagram of observed geopotential height, containing five atmospheric blocking events. OWA’s variance compression is comparable to that of EOF analysis, but is much faster and less data-dependent.

Implementing the basic OWA operations is documented with WaveLab, a freely available package requiring MATLAB.

1. Introduction

Understanding atmospheric structure and dynamics requires both functional and statistical analysis. This paper is about relatively new approaches to these analyses and the situations in which they are preferable. Historically there have mainly been two approaches. One option is to use the observable field values at discrete locations. For example, Onsager (1949) introduced a statistical hydrodynamics based on the approximation of a continuous vorticity field by a set of discrete “point vortices” of constant circulation and variable location. Of course, another “spatial” representation is by values on a fixed grid, as in finite-difference numerical modeling. The alternative technique has been Fourier series, or similar orthogonal function series. Lorenz (1960) showed that analogs of the fields’ physical conservation laws were obeyed by truncations of Fourier series. Charney (1963) observed this to be the case also for general orthogonal function series and for point vortices, but not for fixed-grid discretization. He called the point-vortex and Fourier representations “dual” in the sense of the Heisenberg uncertainty principle, or wave particle duality. These terms usually refer to the reciprocal uncertainty of location and wavenumber information in quantum mechanics. Actually, an uncertainty principle applies to the statistical variances of spatial and wavenumber representations of any distribution (Chui 1992;Wickerhauser 1994; Prestin and Quak 1999). The implied trade-off of information can be a basis for choosing between the analysis methods that we will discuss.

The point-vortex technique is limited to relatively simple kinds of dynamics. Thus to maintain conservation laws, some kind of orthogonal function expansion is needed, in order to treat realistic dynamics. Because actual atmospheric fields may be spatially localized and can have large gradients, Fourier series and related global analyses are not ideally suited; some kind of compromise representation is called for.

a. Spatial versus wavenumber compromises

1) Before wavelets

The uncertainty principle suggests that there can be a mutual exchange of information content, or resolution, between the “extreme” cases of points or waves. In recent decades there have been various imperfect attempts to formulate compromise representations, which somehow combine localization in space and scale. Siggia (1977) introduced the use of “wave packets” to partially decouple position and scale in turbulent dynamics. Zimin (1981) introduced a disjoint wavenumber band representation, later identified by Zimin and Hussain (1995) as a nonorthogonal wavelet basis. Both these approaches required very complicated approximations. Besides these contributions from turbulence investigators, there was also an application to meteorology of bandpass filtering (BPF) designed for limited spatial domains, which had nonnegligible domain boundary effects (Bettge and Baumhefner 1980). These investigations, and those in other sciences, set the stage for the introduction of orthonormal wavelet analysis (OWA), as labeled, for example, by Mak (1995). There are already several reviews available; for example, Kumar and Foufoula-Georgiou (1997) review wavelet analysis for geophysical applications. Examples throughout broader areas of physics are collected by van den Berg (1999).

2) Relevant wavelet applications

OWA is especially relevant when dominant spatially “local” structures are present in a synoptic field, for example, during a blocking situation. [Following Hansen and Sutera (1984), we define blocking, in the 55°–80°N average, by slowly propagating ridges in the 50-kPa height, which exceed the zonal mean by at least 250 m for at least seven days.] The present paper is the first in a sequence that argues that in the case of compactly organized structures such as blocking, wavelet-based analysis techniques are more appropriate than Fourier analysis, to evaluate the relative importance of different-scale effects. [This outstanding issue has been investigated by numerous authors, e.g., Tsou and Smith (1990), using part of the present dataset.] Although wavelet analysis has been applied in the time domain to atmospheric boundary layer turbulence (e.g., Foufoula-Georgiou and Kumar 1994) and climatic time series (e.g., Bolton et al. 1995; Lau and Weng 1995; Mak 1995; Torrence and Compo 1998; Park and Mann 2000), and has been applied in the space domain to observed (Zubair et al. 1992; Zubair 1993) and numerically simulated turbulence (Meneveau 1991; Farge et al. 1992; Fournier and Stevens 1996), there has been little application to observed global synoptic meteorological data. The present author has endeavored to remedy this situation in Fournier (1996, 1999) and Fournier (1999, three manuscripts submitted to J. Atmos. Sci., hereafter F1, F2, and F3).

This paper is also designed as a complement to earlier tutorials by Lau and Weng (1995) and Torrence and Compo (1998), whose methods are suitable for analysis of time series, not spatial structures. The latter authors presented several formulas that should be identified as approximations to appropriate integrals, to lowest order in the time increment δt and/or the scale increment δj;they correctly emphasized the importance of choosing sufficiently small increments. [Shensa (1992) discusses exact methods for such integrals, in the course of his thorough discussion of the relationship of OWA to the continuous wavelet transform (CWT) and other related transforms.] This is adequate for sufficiently smooth signals but may lead to difficulties in reconstructing, filtering, decomposing the variance of, or otherwise analyzing less-smooth signals. By contrast, the methods described in the present paper are arithmetically exact, that is, accurate to computer precision, for the analysis of a uniform sampling of any spatial field. The analyses of the field’s samples arithmetically exactly satisfy equations analogous to those of the continuous field’s analyses.

b. A motivating example

The following further motivates the hypothesis that OWA should be advantageous for atmospheric problems. Figure 1a shows an image of the earth. Cloud patterns suggest the complicated dynamics and structure of the underlying fields, and exhibit a broad range of scales, anisotropy, and organized structures.

Fig. 1.

(a) Geostationary satellite earth image. (b) Box-averaged earth image. (c) Fourier-truncated earth image. (d) Wavelet-truncated earth image

Fig. 1.

(a) Geostationary satellite earth image. (b) Box-averaged earth image. (c) Fourier-truncated earth image. (d) Wavelet-truncated earth image

Often in image processing or geoscience one needs to use a relatively small number of values to represent such an earth image. Say we want to reduce the 218 pixel values to only 212 values of some kind. One method, corresponding to a purely spatial viewpoint, would be to simply average the pixel values in 23 × 23 squares, yielding the piecewise uniform image in Fig. 1b. Obviously, much important detail structure has been lost.

A method corresponding to a purely wavenumber viewpoint could be to truncate wavenumbers ≥25. The image reconstructed from the remaining (2 × 25)2 = 212 Fourier–sine–transform coefficients is shown in Fig. 1c. Although the image is much smoother and more realistic, there is still very little detail left.

The alternative method is to combine some of the spatial and wavenumber information, using coefficients of expansion in a wavelet basis, which is simultaneously localized in both space and wavenumber. The reconstruction from the 212 largest-magnitude coefficients is shown in Fig. 1d. Much of the structure lost by the previous two methods has been regained, such as the ITCZ and the vortices west of South America. This observation has been the main motivation to the investigations of Fournier (1995, 1996, 1999), F1, F2, F3, and Fournier and Stevens (1996) into how to apply OWA to atmospheric spatial structure and dynamics.

c. Outline

An introduction to the basic elements of OWA follows a review of simple statistical ideas familiar to meteorologists, in section 2. Unlike most popular introductions, which were geared toward time series analysis, this includes a review of periodization, appropriate for analysis on latitude circles. Also briefly reviewed are extensions initiated by Beylkin (1992) and Shensa (1992) to compensate for the lack of shift invariance of the basic analysis. This important objection can be satisfactorily dealt with by either using all possible shifts simultaneously (with minor additional computation) or choosing a “best shift.” There are at this time few references to this subject available, and mostly they are aimed at an audience specializing in statistical mathematics or signal processing. The representation of blocking structures by best shift is investigated in F3 and summarized in section 2h. Section 3 contains a detailed survey of examples, directly comparing the accuracy of equally truncated Fourier and wavelet series, for analyzing various test signals. This tutorial is to familiarize the reader with OWA enough to understand its suitability for the physical problem treated in section 4, namely, multiresolution analysis (MRA) of Hovmöller diagrams of blocking.

To offer the reader “hands-on” practice, Table 1 summarizes the actual function syntaxes to implement the key algorithms. These functions are part of the free online MATLAB library, WaveLab (Buckheit et al. 1995).

Table 1.

WaveLab functions for computing several of the quantities described in the text. Affixes indicate multidimensional results or arguments, for example, (0)4ℓ=1 ≡ (0 0 0 0)

WaveLab functions for computing several of the quantities described in the text. Affixes indicate multidimensional results or arguments, for example, (0)4ℓ=1 ≡ (0 0 0 0)
WaveLab functions for computing several of the quantities described in the text. Affixes indicate multidimensional results or arguments, for example, (0)4ℓ=1 ≡ (0 0 0 0)

2. Average and detailed components of structures

a. Basic MRA concepts

Since the work of Reynolds, it has been a common practice in turbulence (large eddy simulation) and meteorological modeling to decompose a field f(x) into its average and deviation parts: f = 𝒲f + 𝒲f. Here 𝒲f is the result of some kind of averaging projection operation, 𝒲 is a smoothing operator, and the residual 𝒲f defines a detailing operator, 𝒲 ≡ 1 − 𝒲. The label ⊥ will appear in different ways here, to signify some kind of orthogonality. For functions and operators these orthogonalities are expressed by 〈𝒲f𝒲fx = 0 [where 〈fx10f(x) dx] and 𝒲𝒲 = 0, respectively. For example, low-pass filtering (LPF) (Blackmon 1976;Bettge and Baumhefner 1980; Saltzman 1990) gives a 𝒲 in the “ideal” LPF case, (⁠𝒲)2 = 𝒲, but the commonly used running average is not a projection: each application of running average further smoothes the field. For concreteness, henceforth 𝒲f and 𝒲f are taken to be 〈fx and ff − 〈fx, respectively.

This is generalized to MRA as follows. Begin with a trivial-looking “telescoping sum” decomposition:

 
formula

Note that this expansion is true regardless of the symbols’ definitions. If fJ stands for a “fine approximation” of f at some “finescale” j = J, then one can always derive it from a “gross approximation” f0 at some “coarse scale” j = 0 by adding a series of “details,”

 
formula

between successive approximations. Now let us combine these simple ideas with OWA.

We are interested in fs that depend on periodic coordinates such as

 
formula

so in the following sections we survey periodic OWA;see, for example, Perrier and Basdevant (1989) for details. (Henceforth “periodic” shall be subsumed into OWA, etc.) First we review Fourier analysis. For ease of discussion, all formulas and notation are summarized in Table 2, arranged in columns by analysis method and in rows by the ingredients of the analyses. Additional notations are 〈fλ ≡ (2π)−1 ππf(λ) and 〈anΣn=−∞an.

Table 2.

Fourier vs OWA formulas comparison. Column 1 describes the row contents. Columns 2–4 contain formulas for the Fourier and wavelet analyses described in the text

Fourier vs OWA formulas comparison. Column 1 describes the row contents. Columns 2–4 contain formulas for the Fourier and wavelet analyses described in the text
Fourier vs OWA formulas comparison. Column 1 describes the row contents. Columns 2–4 contain formulas for the Fourier and wavelet analyses described in the text

b. Fourier analysis

If f is 2π periodic and either absolutely integrable and has bounded ∂λf at λ, or else square integrable and continuous at λ, then f(λ) may be expanded in the Fourier basis. In practice, the limit ∞ in the Table 2, column 1 formulas is replaced by the number 2J of uniform function samples available; thus we get discrete analogs that are arithmetically exact.

c. Wavelet analysis

More generally than for Fourier series, any 2π periodic, continuous or absolutely integrable function f(λ) may be expanded in a periodic orthonormal wavelet basis. The wavelet coefficients j,k are defined using 〈 〉λ or, in Fourier space, 〈 〉n. Again, in practice ∞ is replaced by J − 1 so that the coefficients (including 〈fλ) total in number and are computed as described in section 2f. All the wavelet equations in Table 2 follow from section 9.3 of Daubechies (1992).

1) Different bases and the meaning of the labels

Another way to see how the OWA basis compromises between point and wave representations is given by Fig. 2a. The left panel depicts basis functions taking the value 1 at each of 2J longitudes,

 
λJ,1−Jπ,
(2.1)

and zero elsewhere. This basis underlies the standard uniform discretization of f(λ). It is entirely localized in λ. The middle panel depicts the sines and cosines with wavenumbers n from 0 to 2J−1, as used in the FFT. The Fourier basis is perfectly localized in wavenumber, sometimes equated with scale, as defined ∝n−1. In the right panel, obviously the basis functions are somewhat localized in space. The basis functions come in octaves or levels j = 0, . . . , J − 1, with 2j functions at each level. The functions at lower j are clearly less localized in λ, so j stands for resolution. It turns out that they are more localized in scale (properly defined as intervals of n), as we will see in Fig. 5b. In some sense, some scale information has been traded off to regain lost location information. Such trade-offs must always satisfy the uncertainty principle w.r.t. scale and location. We can visualize this with a unit box, Fig. 3a, symbolizing a 2D location–wavenumber plane to be covered by 2J = 8 basis elements. The space representation (no scale or wavenumber information), the wavenumber representation (no location information), and the wavelet representation (some of both kinds of information) are symbolized by the tilings in Figs. 3b–d, respectively. The uncertainty principle is expressed by each tile in every tiling having the same area, 2J = 1/8. As one tile dimension is reduced, to increase the resolution in that aspect, the other dimension must be increased, which reduces its resolution aspect.

Fig. 2.

(a) Three possible orthonormal bases plotted (vertically offset) as a function of λ. Left: The basis corresponding to gridpoint values. Center: The Fourier basis. Ordinate is wavenumber n. Right:The Daubechies-4 wavelet basis and W0,1 = 1. The ordinate indicates 2j + k, k = 1, . . . , 2j, j = 0, . . . , 4. (b) Different Daubechies-4 orthonormal bases possible for values (left to right) 0 to 4 of level of curtailing R (Table 2), as a function of λ.

Fig. 2.

(a) Three possible orthonormal bases plotted (vertically offset) as a function of λ. Left: The basis corresponding to gridpoint values. Center: The Fourier basis. Ordinate is wavenumber n. Right:The Daubechies-4 wavelet basis and W0,1 = 1. The ordinate indicates 2j + k, k = 1, . . . , 2j, j = 0, . . . , 4. (b) Different Daubechies-4 orthonormal bases possible for values (left to right) 0 to 4 of level of curtailing R (Table 2), as a function of λ.

Fig. 5.

(a) Wavelet basis functions Wj,k for j = 1, . . . , J − 1 = 6, as a function of λ (°) and k (ordinate). There are 2j basis functions at resolution j. The curves are offset in proportion to k to avoid overlapping, and only the nonvanishing parts are plotted. (b) Wavelet Fourier coefficients for j = 1, . . . , J − 1 = 6, as a function of wavenumber n (abscissa) and k (ordinate). Only 2j(maxn |'|)−1|| is plotted (black curve). The label “Mag.<” signifies maxn ||. The gray curves show (2π)−1 arg + n as a function of k, only for n such that || ≥ 0.1

Fig. 5.

(a) Wavelet basis functions Wj,k for j = 1, . . . , J − 1 = 6, as a function of λ (°) and k (ordinate). There are 2j basis functions at resolution j. The curves are offset in proportion to k to avoid overlapping, and only the nonvanishing parts are plotted. (b) Wavelet Fourier coefficients for j = 1, . . . , J − 1 = 6, as a function of wavenumber n (abscissa) and k (ordinate). Only 2j(maxn |'|)−1|| is plotted (black curve). The label “Mag.<” signifies maxn ||. The gray curves show (2π)−1 arg + n as a function of k, only for n such that || ≥ 0.1

Fig. 3.

(a) The location (abscissa λ)–wavenumber (ordinate n) plane; (b) space representation; (c) wavenumber representation; (d)–(f) OWA representations curtailed at levels R = 0, 1, 2, respectively

Fig. 3.

(a) The location (abscissa λ)–wavenumber (ordinate n) plane; (b) space representation; (c) wavenumber representation; (d)–(f) OWA representations curtailed at levels R = 0, 1, 2, respectively

2) Wavelets on the real line

The periodic OWA basis functions Wj,k are obtained from those on the real line, ψj,k, by periodization (see appendix A). The ψj,k are generated by 2j dilation and (k − 1) shift:

 
formula

There are many possible “mother wavelets” (Meyer 1993), W, which generate such a representation. Most Ws cannot be written in terms of explicit functions, but many can be defined by

 
formula

where the scaling function W (Daubechies 1992), also called “wavelet father” (Meyer 1993), is the solution of the dilation equation (e.g., Kaiser 1994; Strang and Nguyen 1995):

 
formula

Sometimes (2.3b) is called the “scaling identity” (e.g., Strichartz 1993), or the “two-scale equation” (e.g., Wickerhauser 1994), or (2.3a) and (2.3b) together are called the “two-scale relations” (e.g., Chui 1992). Ultimately, everything is determined by the two filter sequences and

 
formula

For any , (2.3c) (even with any odd integer replacing 3) implies for any integer ℓ⁠. This in turn implies as long as . Equations (2.2), (2.3a)–(2.3c) are derived from (1.3.4), (5.1.32), (5.1.15), (5.1.35), respectively, of Daubechies (1992). The significance of the cascade interpretation of (2.3a), (2.3b) to turbulence is discussed in appendix B.

There is one case for which W and W may be expressed explicitly. For = 2−1/2(−1, 1) the solutions to (2.3a), (2.3b) are

 
formula

(Daubechies 1992). These relations are schematically depicted in Figs. 4a,b. The orthonormal basis generated by these solutions via (A.1) and (2.2) was invented by Haar (1910).

Fig. 4.

Graphical representation of (a) the Haar solution of (2.3b); (b) the Haar solution of (2.3a); (c) as in (a), for Daubechies-4; (d) as in (b), for Daubechies-4

Fig. 4.

Graphical representation of (a) the Haar solution of (2.3b); (b) the Haar solution of (2.3a); (c) as in (a), for Daubechies-4; (d) as in (b), for Daubechies-4

Note that Haar’s solutions are compactly supported (CS, i.e., they vanish outside a closed, bounded set of t), discontinuous, and zero mean: . Daubechies (1988, 1992) generalized the solutions, showing that, with certain requirements on (reviewed by F3), solutions W exist that are CS and have smoothness and a number of vanishing moments that all increase roughly linearly with the support length 2N of . Figures 4c,d show the simplest (N = 2) Daubechies solutions, for

 
formula

In Fig. 4 the wavelet and scaling-function graphs are horizontally aligned to accurately show how features (such as peaks, dips, and cusps) combine to construct the larger structure. The size (inversion) of the th smaller scaling function indicates the magnitude (sign) of the coefficient or w−2N+2 that multiplies it.

The investigator should choose a particular wavelet type according to the problem to be addressed. For example, F1 uses the Daubechies-14 wavelet, with only 2N = 14 nonzero w. For the purpose of analyzing a partial differential equation solution, this wavelet was chosen to provide the most CS basis available that had continuous second derivatives. In the present section 4, Fournier (1999), and F2 and F3, the Daubechies-20 wavelet is used, with only 2N = 20 nonzero w. This was the smoothest CS wavelet available, selected on the assumption that atmospheric structures are at least twice differentiable on the scale observed.

3) Wavelet location and support

Figure 5a shows the Wj,ks (for λε = 0 in Table 2, and N = 10) in an arrangement that sets up Fig. 9. The Wj,ks are centered approximately at the pyramidal set of longitudes λj,k = λJ,j,k. Each λj,k is related to the grid (2.1) by the pyramidal index set

 
formula

and constructs a grid from −π to π with step 21−jπ. Pyramidal simply describes the constraint k = 1, . . . , 2j as exemplified in Table 3. The support of Wj,k has length equal to 2π min[1, 2j(2N − 1)].

Fig. 9.

(a) Hovmöller diagrams of Z(λ, 65°N, 70 kPa, t) (Dm), vs λ (°) and t (day). White bars indicate t intervals T and t average crest longitudes. The horizontal bar graph at right is ∝〈(∂λ/∂t)Zλt for each T. Dark (light) curves and lower-left (-right) to upper-right (-left) diagonal structures typify eastward (westward) progressing synoptic weather patterns. Vertical (horizontal) alignment indicates quasi steadiness (rapid motion). Lines parallel to the main diagonal correspond to the speed (90 day)−12πr[fy1,15]=Z cos 65° ≈ 8 km h−1. (b)–(g) Zj(λ, 65°N, 70 kPa, t), as in (a), arranged by j as in Fig. 5. (h)–(m) As in (b)–(g), except here each panel labeled “BPF j” corresponds to Fourier BPF (3.1). (n) Z(λ, 70 kPa, 44.5 day) and topography (Dm) at 65°N (ordinate) vs λ (°)

Fig. 9.

(a) Hovmöller diagrams of Z(λ, 65°N, 70 kPa, t) (Dm), vs λ (°) and t (day). White bars indicate t intervals T and t average crest longitudes. The horizontal bar graph at right is ∝〈(∂λ/∂t)Zλt for each T. Dark (light) curves and lower-left (-right) to upper-right (-left) diagonal structures typify eastward (westward) progressing synoptic weather patterns. Vertical (horizontal) alignment indicates quasi steadiness (rapid motion). Lines parallel to the main diagonal correspond to the speed (90 day)−12πr[fy1,15]=Z cos 65° ≈ 8 km h−1. (b)–(g) Zj(λ, 65°N, 70 kPa, t), as in (a), arranged by j as in Fig. 5. (h)–(m) As in (b)–(g), except here each panel labeled “BPF j” corresponds to Fourier BPF (3.1). (n) Z(λ, 70 kPa, 44.5 day) and topography (Dm) at 65°N (ordinate) vs λ (°)

Table 3.

Pyramidal indexes ℓj,k for J = 3. The horizontal spacing around ℓj,k indicates the support of Wj,k, while the 2D tableau arrangement illustrates the inheritance of scale/location information under the computations (2.5a)–(2.5e)

Pyramidal indexes ℓj,k for J = 3. The horizontal spacing around ℓj,k indicates the support of Wj,k, while the 2D tableau arrangement illustrates the inheritance of scale/location information under the computations (2.5a)–(2.5e)
Pyramidal indexes ℓj,k for J = 3. The horizontal spacing around ℓj,k indicates the support of Wj,k, while the 2D tableau arrangement illustrates the inheritance of scale/location information under the computations (2.5a)–(2.5e)

d. Generalization to arbitrary curtailing

The functions Wj,k defined by (A.1), replacing W by W in (2.2), satisfy the limited orthogonalities in Table 2. The WR,k are averaging functions, associated with LPF, and f0 = f0,1 = 〈fλ since W0,k(λ) = 1. This is called OWA “curtailed at level R” (Nason and Silverman 1995). Its physical interpretation is that the scaling coefficients fR,k now contain all information corresponding to a certain low-wavenumber interval, bounded ∝2R.

Figure 2b shows J = 5 possible curtailings. The leftmost panel copies the lower half (j < J − 1) of the right panel of Fig. 2a. Level j = J − 1 is unaffected by curtailing at levels R < J. Curtailing at level R = J corresponds to the identity transformation. From left to right, as R increases from 0 to J − 1 the first 2R basis functions shown in the leftmost panel are replaced by the corresponding scaling functions WR,k. The 2J − 2R wavelets above are unaffected. Each of these panels (taken together with the j = J − 1 half not shown) illustrates another possible orthonormal basis. For J = 3 the tilings corresponding to R = 1 and 2 are in Figs. 3e,f, respectively.

e. Fourier picture of OWA

OWA can also be visualized in the Fourier domain; the are shown in Fig. 5b. From the nth Fourier coefficient of (A.1) we see that the magnitude || is independent of location index k for a given resolution index j, so it is only plotted once, for k = 1 (black curve). Also, the phase arg is just a shift by 2π(1 − k)2jn from the continuous Fourier transform of W at frequency 2jn (gray curve for each n). Thus for n = 2j the phase shift is zero, while for n = 2j−1 every change of k by 2 sweeps a whole cycle.

Since both W and are localized in their respective spaces, j,k selects the contribution to f from location ∝2jk and scale ∝2j. The set defines an MRA of f. It is essentially a sophisticated kind of BPF technique, to use a term more familiar to climatologists (Blackmon 1976; Saltzman 1990). The are for BPF, while the are for LPF. The BPF bandwidth and band center are seen to increase with j as the spatial resolution increases (Fig. 5a), in accordance with the uncertainty principle (Meyer 1993; Wickerhauser 1994). Fourier and wavelet series are further compared in appendix C.

f. Practical OWA

In practice, 2J discrete uniform samples f(λJ,k) ≈ 2J/2fJ,k of f are given. Then discrete analogs of the 2R scaling coefficients fR,k+1, k = 0, . . . , 2R − 1, and the 2J − 2R wavelet coefficients j,k+1, k = 0, . . . , 2j − 1, j = R, . . . , J − 1, are computed in O(2J) convolution–decimation calculations:

 
formula

(Beylkin 1992; Nason and Silverman 1995, flipping the sign of εj to conform with WaveLab). (The decimation of ak is a2kε, ε = 0, 1.) For 2π periodic f and wℓ<1 = wℓ>2N = 0 (the periodic CS case) these become

 
formula

Equation (2.5) involves just , so Wj,k never needs to be evaluated. This efficient computation is due to the“self-similarity” (2.3a), (2.3b). Note that the corresponding FFT requires O(J2J) calculations and is therefore slower. The reconstruction is the inverse of this orthogonal transformation:

 
formula

In the periodic CS case this becomes

 
formula

A visual pattern corresponding to (2.5a)–(2.5c) corresponds to ascending a pyramid resembling Table 3; likewise (2.5d), (2.5e) correspond to descending. This discrete wavelet transform (DWT) satisfies the discrete analogs of all the previous equations, just by replacing the operation 〈fλ everywhere by 2J Σ2Jℓ=1f.

The binary sequence εj = 0, 1 appearing in (2.5) determines whether even- or odd-indexed samples are used at level j of convolution–decimation, as introduced independently by Shensa (1992) and (in more algorithmic detail) by Beylkin (1992). The reference longitude λε is related to in appendix D. Including all possible , the “overcomplete” set of corresponds to using all 2J possible periodized shifts of 2j-dilated W, instead of the orthonormal set of 2j shifts indexed by k = 1, . . . , 2j, but still only requires O(J2J) calculations (Beylkin 1992), the same cost as the FFT. Appropriately ordered, this larger set of analysis functions generates the nonorthogonal shift-equivariant WT (SEWT), sometimes called “undecimated WT” (Shensa 1992), “stationary WT” (Nason and Silverman 1995), “shift-invariant WT” (Lang et al. 1995), “translation-invariant WT” (Coifman and Donoho 1995; Del Marco and Weiss 1997) or “maximal-overlap DWT” (MODWT; Percival and Mofjeld 1997). The most physical description is equivariant (e.g., von Sachs et al. 1998), meaning a shift of the function produces the same shift of the SEWT [as opposed either to the relations (D.3), wherein the “detail bits” of the shift are truncated, or to (D.4)]. This term is distinguished from “shift invariant,” describing the Fourier spectrum |n|2 but not the Fourier transform itself, which merely diagonalizes the shift operator.

The unrestricted domain of location index k “fills in the gaps” left from the orthonormal basis. Figure 6 illustrates this overcomplete, or redundant, representation. The standard OWA (SOWA, defined by λε ≡ 0) basis functions are shown in heavy outline. The other 2J − 1 possible orthonormal sets are all contained in the redundant set.

Fig. 6.

Shift-equivariant wavelet basis based on the Haar basis (2.4) with R = 0. The J = 3 columns of panels indicate the resolution index j = 0, 1, 2 from left to right. The 2J = 8 rows indicate the position λJ,. The standard (λε ≡ 0) orthonormal set of basis functions is shown in heavy outline

Fig. 6.

Shift-equivariant wavelet basis based on the Haar basis (2.4) with R = 0. The J = 3 columns of panels indicate the resolution index j = 0, 1, 2 from left to right. The 2J = 8 rows indicate the position λJ,. The standard (λε ≡ 0) orthonormal set of basis functions is shown in heavy outline

Among the first SEWT applications were denoising (Coifman and Donoho 1995; Lang et al. 1995) and estimating noisy transients and time delays (Pesquet 1996). Further applications of SEWTs and their wavelet-packet generalizations (Pesquet 1996; Del Marco and Weiss 1997) include examining magnetic field shock waves (Walden and Contreras Cristan 1998), assessing formula errors with observed turbulent heat fluxes (Qi and Neumann 1997), synthesizing multifractal anisotropic 2D cloud profiles (Davis et al. 1999), and analyzing river water levels (B. Whitcher et al. 1999, manuscript submitted to Water Resour. Res.). In particular, Percival and Mofjeld (1997) demonstrated the power of their MODWT for detecting and meaningfully alligning MRA features with the original data, and efficiently decomposing data variance, and explained MODWT’s advantages over both short-time Fourier transform and complex demodulation.

With a minor abuse of notation, the coefficients may be denoted fj,2jJ, and j,2jJ where the fractional k index accurately reflects the location defined by λj,k→2jJ and allows = 1, . . . , 2J. Or, generalizing Shensa’s (1992) Eq. (3.7), use (D.1) to define

 
formula

which are easily shown to be shift equivariant, and by (D.4) reduce to tautologies of the form A = A as → 2Jjk. Although not orthonormal, the SEWT still satisfies a Parseval identity,

 
formula

which is essential for any kind of analysis of variance, for example, the wavelet energetics of F1, F2, and F3. The factors 2j compensate for an emphasis of the larger scales in this representation, in a sense, the “price” of spreading the large-scale representation over a finer location grid. In a scale analysis of covariance, omitting the factors 2j would cause larger-scale contributions to be emphasized, relative to the orthonormal analysis.

g. Best shift of the wavelet basis

A particular orthogonal representation, corresponding to a particular

 
formula

may be extracted from the overcomplete set by maximizing some measure of the representation’s efficiency, w.r.t. . Such a best-shift wavelet analysis (BSWA) was developed roughly contemporaneously by several investigators; see Liang and Parks (1996), Pesquet (1996), and their references. The λ+ε computation is shift equivariant but generally depends nonlinearly on one reference input. For the appropriate λ+ε, the BSWA is shift invariant and orthonormal w.r.t. arbitrary inputs.

One may find λ+ε in O(J2J) calculations by minimizing a cost functional such as () [the analysis entropy of Coifman and Wickerhauser (1992)] w.r.t. , where

 
formula

As shown by Wickerhauser (1994), is useful for at least two reasons. It is convex w.r.t. possible analyses (), that is,

 
formula

Hence any weighted averaging of analyses increases ℰ⁠, reducing information, so conversely reducing increases information, pointing to the most efficient representation. In the extreme cases, 0 ⩽ () ⩽ J, where () = 0 if and only if has at most one nonzero term (maximum representation efficiency), and () = J if and only if a is independent of (worst case: no information). Wickerhauser (1994) also shows that OWAs optimized in similar ways to this approach the optimal EOF analysis at substantially reduced computational cost. Liang and Parks (1996) classify several alternative cost functionals and evaluate their suitability for fast algorithms and effectiveness for compression. Del Marco and Weiss (1997) extended the method from base-2 to base-M in (2.3a), (2.3b), and to extra flexibility in minimizing by simultaneously combining multiple λεs.

h. Application of BSWA

Let us briefly review some recent applications of BSWA that are closely related to problems investigated by traditional techniques. BSWA is applied by F3 to observed 50-kPa geopotential height data. In four out of five cases of time- and latitude-averaged atmospheric blocking, only two optimally shifted wavelets contribute more variance (at least 89%) than four or more SOWA wavelets. For certain j = 1 type blocks, the observed height distribution remarkably resembles the wavelet W1,k(λ) itself. The study of F3 discusses physical reasons why this might be the case and other topics, for example, how to interpret physically the signal-processing motivations underlying wavelet design. Also, F3 derives a minimal Fourier coefficient model of blocking from W1,k, relates it to wavelet design considerations, and interprets the minimal model with regard to spectral dynamics models of Platzman (1962) and Wiin Christensen and Wiin-Nielsen (1996).

3. Examples of wavelet analysis

Now we shall compare examples of various N = 2 analyses of several test functions f accompanying the free WaveLab software (Buckheit et al. 1995). Each f from Table 4 was first centered and normalized by f − 〈fxf,f−1ff.

Table 4.

Errors for analyses of synthetic fs (column 1) in terms of Fourier (columns 2, 5, 8), SOWA (columns 3, 6, 9), and BSWA (columns 4, 7, 10) truncated series. Column triplets 2–4, 5–7, and 8–10 show mean absolute, root-mean-square, and maximum errors, respectively. The font indicates least (italic) and greatest (bold) of each triplet. Thresholds and numbers of retained coefficients are given in Fig. 7, and the parameters of rows e and f are given in Table 5

Errors for analyses of synthetic fs (column 1) in terms of Fourier (columns 2, 5, 8), SOWA (columns 3, 6, 9), and BSWA (columns 4, 7, 10) truncated series. Column triplets 2–4, 5–7, and 8–10 show mean absolute, root-mean-square, and maximum errors, respectively. The font indicates least (italic) and greatest (bold) of each triplet. Thresholds and numbers of retained coefficients are given in Fig. 7, and the parameters of rows e and f are given in Table 5
Errors for analyses of synthetic fs (column 1) in terms of Fourier (columns 2, 5, 8), SOWA (columns 3, 6, 9), and BSWA (columns 4, 7, 10) truncated series. Column triplets 2–4, 5–7, and 8–10 show mean absolute, root-mean-square, and maximum errors, respectively. The font indicates least (italic) and greatest (bold) of each triplet. Thresholds and numbers of retained coefficients are given in Fig. 7, and the parameters of rows e and f are given in Table 5

a. Fourier analysis versus SOWA and BSWA

Except where noted, both the BSWA and Fourier series were truncated at the number of SOWA coefficients exceeding the magnitude shown, to illustrate the relative merits of the methods. At the left of Figs. 7e,f,k,z, the wavelets are vertically sorted by the variance percentage (VP) each contributes to the reconstruction, cumulatively labeled for each (j, k). After summing some (usually small) number of terms, for example, 3(8) for Fig. 7e (f), the eventual inferiority of the Fourier components (not graphed) is indicated by the right-hand labels’ (n) lower scores in the VP measure. The VP measure also indicates analysis of variance efficiency.

Fig. 7.

Comparison of truncated wavelet vs Fourier analyses of (a)–(e) Ramp, Table 4, row a, with best-shift wavelet; (f) as in (e), with standard wavelet; (g)–(k) as in (a)–(e), for the singular function, Table 4, row c; (l)–(p) the bumps function, Table 4, row e; (q)–(u) the block function, Table 4, row f; (v)–(z) the discontinuous sine function, Table 4, row i

Fig. 7.

Comparison of truncated wavelet vs Fourier analyses of (a)–(e) Ramp, Table 4, row a, with best-shift wavelet; (f) as in (e), with standard wavelet; (g)–(k) as in (a)–(e), for the singular function, Table 4, row c; (l)–(p) the bumps function, Table 4, row e; (q)–(u) the block function, Table 4, row f; (v)–(z) the discontinuous sine function, Table 4, row i

Figure 7a shows the “Ramp” function, such as might model a breaking wave. [In fact, t−1f(x) solves the nonlinear advection equation Ut = −UUx; see F1.] The truncated Fourier series exhibits the Gibbs phenomenon (Fig. 7b), while the truncated BSWA yields a much better approximation (Fig. 7c). For most x, the BSWA errors are much smaller than the Fourier series errors (Fig. 7d). Note in Fig. 7e how the smaller-scale (higher j) wavelets appear only at the discontinuity, and how the cusps of the larger-scale (lower j) wavelets align and sum to a smooth line. Figures 7e,f show the improvement of BSWA over SOWA by the VP measure. Note that VP compression does not necessarily follow visually perceived accuracy; pointwise SOWA errors (not shown) are less than BSWA’s shown in Fig. 7d, as indicated in row a of Table 4.

How well do the methods reconstruct a point singularity? Table 4, row b, shows the truncation errors for the Kronecker (smoothed Dirac) function. Graphically this is similar to the less singular f analyzed in row c and Figs. 7g–k. With only 23 coefficients the Fourier series has barely started to build up this singularity (Fig. 7h) and performs terribly w.r.t. VP (Fig. 7k). The OWA errors are significantly smaller (Fig. 7j). The Fourier series is better for Table 4, row c, than row b, reflecting the weaker localization in x, but in neither case does it perform nearly as well as OWA.

Further delocalizing in x, consider Table 4, row d. Again the Fourier series improves but is still inferior to OWA, especially for describing the cusp (not shown). In this case the OWA actually yields two more coefficients larger than the threshold (0.01), relative to the Fourier analysis.

Functions a–d might represent various kinds of vortices, examples of localized phenomena occurring in geophysical fluid dynamics. This suggests that wavelets are useful for representing vortices, as has been demonstrated by Farge (1992), Farge et al. (1992, 1999), and Schneider et al. (1997). Wavelets are also advantageous for representing topography-like structures such as shown in Figs. 7l–p. In order to capture the very localized “bumps,” Fourier analysis must include high wavenumbers, which contaminate the level “plains” (Fig. 7m) and increase the overall truncation error (Fig. 7o). The small-scale wavelet components only occur where needed in Fig. 7p. Similar behavior appears for the piecewise uniform function shown in Figs. 7q–u.

We now consider functions in which Fourier analysis has some advantage. Clearly for the purely harmonic cases, Table 4, rows g and h, the Fourier analysis does a somewhat better job. In these two analyses there were many more significant OWA coefficients compared to Fourier coefficients, so for the sake of comparison, truncation was made on the magnitude of the latter (as not in the other examples).

Suppose we add a little discontinuity to a sine wave. Figures 7w,x show the reconstructions using 35 coefficients for Fourier or BSWA analysis, determined from the BSWA coefficient magnitudes. Figure 7z shows better VP compression by the Fourier coefficients. As seen in Table 4, row i, Fourier analysis also has slightly smaller mean absolute error, because the wavelets in this case are visibly not smooth enough and so introduce some roughness (Fig. 7x). However, OWA reproduces the discontinuities more sharply and so has smaller rms and max errors. Note: in this case BSWA does not improve over SOWA, since the two discontinuities are separated in x by 0.42, which is nondyadic and so not amenable to a single best shift. We speculate that the multiple-λε approach of Del Marco and Weiss (1997) may help here.

We conclude this survey with two cases related to harmonic functions, yet in a sense more amenable to OWA. Table 4, row j, shows the truncation errors for spatially varying wavenumber, nD(x) ≈ 21[x(20x + 1)]−1 (away from x = 0, 1). For the low-wavenumber part Fourier analysis does better, but at high wavenumbers the truncation greatly increases the Fourier series error, and OWA does better. OWA also does better with the linearly varying wavenumber “chirp” function, Table 4, row k. OWA outperforms Fourier analysis in this case at both the extreme large and small wavenumbers.

b. MRAs and SEWTs versus ideal BPF

Now we revisit some of the wavelet coefficients just discussed, combined scale-by-scale for λε = 0 (SOWA case) to form MRAs, or sorted by each λε to form SEWTs. Figure 8 shows three analyses of the functions in Table 4, rows e and j. The sum of the fj (or the fbpj below, but not the SEWT) regains f, with the indicated error ɛ. These cases illustrate how MRA and SEWT isolate local features at the appropriate scale. Figures 8d–f exhibit nD(x): observe how the “activated” wavelets proceed to smaller scales as x decreases. Given that n ∝ 2j for a given wavelet level (see Fig. 5b), it is satisfying to observe that the maximum amplitude follows very closely to the indicated (gray) curve Ds ≡ {(x, j)|j − log2nD(x) = s = (small constant)}. (The nD = 2j curve D0 becomes less least squares fit in going from MRA through SEWT to BPF.) Similarly for Table 4, row k, the wavenumber varies as nc(x) ≡ 2J−2x, and the maximum amplitude mostly follows the curve (not shown) j − log2nc(x) = (small constant).

Fig. 8.

Test MRA analyses fj [(a), (d)], SEWTs j,2jJℓ=x [(b), (e)], and dyadic BPFs (3.1) fbpj [(c), (f)] vs x (abscissa). The analyzed functions are from (a)–(c) row e, (d)–(f) row j of Table 4. Each curve is offset by j (ordinate) and rescaled to avoid overlap

Fig. 8.

Test MRA analyses fj [(a), (d)], SEWTs j,2jJℓ=x [(b), (e)], and dyadic BPFs (3.1) fbpj [(c), (f)] vs x (abscissa). The analyzed functions are from (a)–(c) row e, (d)–(f) row j of Table 4. Each curve is offset by j (ordinate) and rescaled to avoid overlap

For comparison with traditional techniques, attempt the same effect with the ideal BPF achieving roughly the same dyadic partitioning in n: Compare fj with

 
formula

( being the integer part of t), using nonorthogonal wavelet

 
formula

whose n′th Fourier coefficient − 1), 0 (otherwise). This filter is wavenumber orthogonal (〈ψnψnλ = δn,n) and zero mean (〈ψnλ = 0), and decomposes the zonal deviation f = Σj=0fbpj and zonal variance f2λ = Σj=0 (fbpj)2. The results of this BPF of the cases from Table 4, rows e and j, are shown in Figs. 8c,f. In every case some of the features revealed by SOWA–MRA and SEWT still appear, but with much poorer spatial resolution, especially at large scales. There are contaminating oscillations where the graphs should be flat. The fbpj are smoother, but the fj may be made smoother by using smoother wavelets. In terms of analysis of variance (label to right of jth curve indicating contribution from jth component), SOWA–MRA achieves the most compression for both functions (Figs. 8a,d). SEWT has some variance compression advantage over BPF for the variable wavenumber function (Figs. 8e,f) but, as noted before, not for the function containing discontinuities with nondyadic separation (Figs. 8b,c).

4. Time-dependent analyses of geopotential height

a. MRA

In these sections five blocking periods are analyzed, of 7–13-day duration (Hansen and Sutera 1984; Fournier 1996, 1999; F2; F3). After F2, ϕ = 65°N was judged to be a representative latitude to discriminate blocking from nonblocking, since it cuts through the block structure. It is standard to depict blocking event sequences with a Hovmöller (1949) diagram. This is shown in Fig. 9a for Z(λ, t) at 70 kPa, with t = 0 ≐ 0000 UTC 1 December 1978. The light regions of this figure clearly depict four blocking events in the time intervals T = [0, 6.5], [19, 26.5], [44, 56.5], [77, 83] days near the Atlantic and one in T = [28.5, 40] days near the Pacific. The third Atlantic blocking event is very strong and clearly moves in a retrograde (westward) manner. Apart from these features it is difficult to extract other insight from the complicated patterns. The time and zonal average phase velocity 〈(∂λ/∂t)Zλt does not consistantly show the slowing commonly associated with blocking. There is range of scales displayed, since there are both regions with many close contour curves and broad, more uniform areas.

A standard way to proceed is to filter the Hovmöller diagram in Fourier space for each t (Wiin-Nielsen and Chen 1993). For comparison, take the N = 10 MRA

 
formula

shown in Figs. 9b–g. The sum of these six panels reproduces Fig. 9a. The j = 1 by itself clearly shows the same blocking ridge pattern as Fig. 9a. This suggests that information encoded in just the two numbers 1,1 and 1,2 may instantaneously describe the presence of blocking, at least for these data. Although j = 2 evidently partly contributes to the block, generally the intermediate scales 2j = 1/4, 1/8, 1/16 capture patterns associated with synoptic cyclone waves, from Atlantic size (j = 2) down to Hudson Bay size (j = 4). Generally speaking, there are more prograde (eastward tending, dark curves) structures evident during nonblocking. For example, Z1 for |λ| ≲ 90 during the first nonblocking interval, [7, 18.5] days, and Z2 for 0 ≲ λ ≲ 100° (over Europe and Russia) during the fourth nonblocking interval, [57, 76.5] days; Z3 at these times and also during the third nonblocking interval, [40.5, 43.5] days. The observable tendencies are usually reflected by 〈(∂λ/∂t)Zjλt.

It should be noted that there are exceptions to these observations. That is, occasionally Zj may appear steady during nonblocking or prograde during blocking. Likewise, strong local anomalies prevent 〈(∂λ/∂t)Zjλt from being a foolproof indicator of blocking action. This is partly why the very definition of blocking is still a subject of much debate (Liu 1994), and why part of the sequence of studies by the present author is aimed at characterizing blocking using MRA techniques.

b. Comparison with Fourier BPF

The apparent low-amplitude (a few meters), small-scale (j > 4) steady structures over the eastern United States and Europe in Figs. 9f,g seem unphysical. These are probably artifacts but not due to this analysis; almost identical patterns are seen in high-j Fourier BPF (3.1) of these data, shown in Figs. 9h–m. This figure also shows how blocks, belonging to the large-scale (low j) band, are more poorly represented by simple BPF than by MRA. For example, although Zbp1 exhibits greater amplitude during some blocking events, it cannot indicate in which hemisphere the block lies. Figure 9n shows Z(λ, 44.5 day) and the surface topography at 65°N. This time was chosen from the high-j panels of previous figures as exhibiting relatively high-amplitude small-scale structure. It is clear that the high-j steady structures must be associated with orography, especially mountain ranges. The perturbations are probably smaller than the original measurement accuracy and so probably indicate mountain-generated artifacts from the National Meteorological Center (now the National Centers for Environmental Prediction) analysis procedure (K. E. Trenberth 1997, personal communication).

c. Wavelet coefficient phase space

To investigate whether the presence of blocking might be described by just the two coefficients 1,1 and 1,2, their time series are presented in Fig. 10a. For the most part, the curves tend to stay at large magnitudes, apart from each other during blocking, and to approach each other at small magnitudes during nonblocking. As might be expected, the Eastern and Western Hemisphere wavelet modes exchange signs during the Pacific blocking relative to the Atlantic blocking.

Fig. 10.

(a) Time series of 1,1 (black curve) and 1,2 (gray curve, Dm), at (ϕ, p) = (65°N, 70 kPa). Abscissa is t (day). White, gray, and light gray backgrounds indicate nonblocking, Atlantic blocking, and Pacific blocking, respectively. (b) Phase space of 1,1 (ordinate) vs 1,2 (abscissa, Dm). Dark gray, gray, and white strips of width 2 Dm for nonblocking, Atlantic blocking, and Pacific blocking, respectively. Light gray background indicates quadrants I and III (like-signed 1,k). Every third day an arrow to the right of the trajectory shows the 12-h average direction of d(1,1, 1,2)/dt. Later segments overlie earlier ones at crossings

Fig. 10.

(a) Time series of 1,1 (black curve) and 1,2 (gray curve, Dm), at (ϕ, p) = (65°N, 70 kPa). Abscissa is t (day). White, gray, and light gray backgrounds indicate nonblocking, Atlantic blocking, and Pacific blocking, respectively. (b) Phase space of 1,1 (ordinate) vs 1,2 (abscissa, Dm). Dark gray, gray, and white strips of width 2 Dm for nonblocking, Atlantic blocking, and Pacific blocking, respectively. Light gray background indicates quadrants I and III (like-signed 1,k). Every third day an arrow to the right of the trajectory shows the 12-h average direction of d(1,1, 1,2)/dt. Later segments overlie earlier ones at crossings

The evolution of 1,k(t), k = 1, 2, can also be described in a dynamical phase space; see Fig. 10b. Ordering the quadrants in the standard way indicated (Shapiro 1977), it may be seen that Atlantic blocking is roughly characterized by clockwise orbits, 1,1 > 0 leading 1,2 < 0, mostly in quadrant II, while Pacific blocking is characterized by counterclockwise orbits, 1,1 < 0 lagging 1,2, in quadrants III and IV. Nonblocking points tend to cluster nearer to the origin in all quadrants. The empty regions enclosed by the orbits suggest centers of instability, although either longer time series or projection onto higher-dimensional phase space, or both, would be required to further support this. If such diagnostic qualities were robust for more data, they could be useful in predicting blocking, an outstanding problem in extended-range weather prediction.

5. Conclusions

Orthonormal wavelet analysis (OWA) yields the scale separation, or bandpass filtering, of Fourier analysis, while retaining spatial information. The combined accuracy of retained information is limited by the uncertainty principle. The examples in this paper, whether from images, explicit functions, or observed geopotential height, illustrate the greater efficiency of OWA over Fourier analysis, in many cases.

The author offers this tutorial and examples, including the guide (Table 1) to the use of freely available software, with three purposes in mind:

  1. to supplement the well-established material reviewing wavelet analysis of time series (e.g., Lau and Weng 1995; Torrence and Compo 1998) with techniques that are suitable for analysis of spatial structures;

  2. to help a broader part of the meteorological community to go beyond the limitations of pure location or wavenumber representations; and

  3. to illustrate, using atmospheric data, the advantages and complexities that await the user of these OWA techniques, in regard to physical interpretation.

Fig. 5.

(Continued)

Fig. 5.

(Continued)

Fig. 7.

(Continued)

Fig. 7.

(Continued)

Fig. 7.

(Continued)

Fig. 7.

(Continued)

Fig. 7.

(Continued)

Fig. 7.

(Continued)

Fig. 7.

(Continued)

Fig. 7.

(Continued)

Fig. 9.

(Continued)

Fig. 9.

(Continued)

Fig. 9.

(Continued)

Fig. 9.

(Continued)

Fig. 9.

(Continued)

Fig. 9.

(Continued)

Table 5.

Parameters of Table 4, rows e and f. Width, amplitude, and location are parameterized by ;utr, ;utc, and ;utx, respectively

Parameters of Table 4, rows e and f. Width, amplitude, and location are parameterized by ;utr, ;utc, and ;utx, respectively
Parameters of Table 4, rows e and f. Width, amplitude, and location are parameterized by ;utr, ;utc, and ;utx, respectively

Acknowledgments

I thank A. B. Davis, A. R. Lupo, and an anonymous reviewer for much useful input; B. Saltzman, R. B. Smith, K. R. Sreenivasan, R. R. Coifman, A. R. Hansen, P. D. M. Parker, J. C. van den Berg, and J. D. Berner for comments on earlier versions; and the National Center for Atmospheric Research for providing the observational data. This material is based upon work supported by the National Science Foundation under Grant 9420011 and through NCAR Project 36211014.

REFERENCES

REFERENCES
Arneodo, A., S. Manneville, J. F. Muzy, and S. G. Roux, 1999: Revealing a lognormal cascading process in turbulent velocity statistics with wavelet analysis. Philos. Trans. Roy. Soc. London, Series A, 357, 2415–2438
.
Bettge, T., and D. Baumhefner, 1980: A method to decompose the spatial characteristics of meteorological variables within a limited domain. Mon. Wea. Rev., 108, 843–854
.
Beylkin, G., 1992: On the representation of operators in bases of compactly supported wavelets. SIAM J. Numer. Anal., 29, 1716–1740
.
Blackmon, M., 1976: A climatological spectral study of the 500 mb geopotential height of the Northern Hemisphere. J. Atmos. Sci., 33, 1607–1623
.
Bolton, E. W., K. A. Maasch, and J. M. Lilly, 1995: A wavelet analysis of Plio-Pleistocene climate indicators: A new view of periodicity evolution. Geophys. Res. Lett., 22, 2753–2756
.
Buckheit, J., S. Chen, D. Donoho, I. Johnstone, and J. Scargle, cited 1995: About WaveLab. Stanford University Tech. Report. [Available online at http://www-stat.stanford.edu/∼wavelab/.]
.
Charney, J., 1963: Numerical experiments in atmospheric hydrodynamics. Experimental Arithmetic, High Speed Computing and Mathematics Vol. 15, Proc. Symp. Appl. Math. Chicago, IL, and Atlantic City, NJ, Amer. Math. Soc., 289–309
.
Chui, C. K., 1992: An Introduction to Wavelets. Vol. 1, Wavelet Analysis and its Applications, Academic Press, 266 pp
.
Coifman, R., and M. V. Wickerhauser, 1992: Entropy-based algorithms for best basis selection. IEEE Trans. Inform. Theory, 32, 712–718
.
——, and D. Donoho, 1995: Translation-invariant denoising. Wavelets and Statistics, A. Antoniadis and G. Oppenheim, Eds., Lecture Notes in Statistics, Vol. 103, Springer-Verlag, 125–150
.
Daubechies, I., 1988: Orthonormal bases of compactly supported wavelets. Commun. Pure Appl. Math., 41, 909–996
.
——, 1992: Ten Lectures on Wavelets. SIAM, 357 pp
.
Davis, A. B., A. Marshak, and E. E. Clothiaux, 1999: Anisotropic multiresolution analysis in 2D: Application to long-range correlations in cloud millimeter-radar fields. Proc. Wavelet Applications VI, Orlando, FL, SPIE, 194–207
.
Del Marco, S., and J. Weiss, 1997: Improved transient signal detection using a wavepacket-based detector with an extended translation-invariant wavelet transform. IEEE Trans. Sig. Proc., 45, 841–850
.
Farge, M., 1992: Wavelet transforms and their applications to turbulence. Ann. Rev. Fluid Mech., 24, 395–457
.
——, E. Goirand, Y. Meyer, F. Pascal, and M. Wickerhauser, 1992: Improved predictability of two-dimensional turbulent flows using wavelet packet compression. Fluid Dyn. Res., 10, 229–250
.
——, N. K.-R. Kevlahan, V. Perrier, and K. Schneider, 1999: Turbulence analysis, modelling and computing using wavelets. Wavelets in Physics, J. C. van den Berg, Ed., Cambridge, 117–200
.
Federbush, P., 1993: Navier and Stokes meet the wavelet. Commun. Math. Phys., 155, 219–248
.
Foufoula-Georgiou, E., and P. Kumar, Eds., 1994: Wavelets in Geophysics. Vol. 4, Wavelet Analysis and its Applications, Academic Press, 373 pp
.
Fournier, A., 1995: Wavelet representation of lower atmospheric long nonlinear wave dynamics, governed by the Benjamin-Davis-Ono-Burgers equation. Proc. Wavelet Applications II, Orlando, FL, SPIE, 672–681
.
——, 1996: Wavelet analysis of observed geopotential and wind: Blocking and local energy coupling across scales. Proc. Wavelet Applications in Signal and Image Processing IV, Denver, CO, SPIE, 570–581
.
——, 1999: Transfers and fluxes of wind kinetic energy between orthogonal wavelet components during atmospheric blocking. Wavelets in Physics, J. C. van den Berg, Ed., Cambridge, 263–298
.
——, and D. E. Stevens, 1996: Wavelet multiresolution analysis of numerically simulated 3d radiative convection. Proc. Wavelet Applications III, Orlando, FL, SPIE, 642–653
.
Frisch, U., 1995: Turbulence: The Legacy of A. N. Kolmogorov. Cambridge University Press, 296 pp
.
Haar, A., 1910: Zur Theorie der orthogonalen Funktionen-Systeme (On the theory of orthogonal function systems). Math. Ann., 69, 331–371
.
Hansen, A. R., and A. Sutera, 1984: A comparison of the spectral energy and enstrophy budgets of blocking versus nonblocking periods. Tellus, 36A, 52–63
.
Hovmöller, E., 1949: The trough-and-ridge diagram. Tellus, 1, 62–66
.
Kaiser, G., 1994: A Friendly Guide to Wavelets. Birkhauser, 300 pp
.
Kumar, P., and E. Foufoula-Georgiou, 1997: Wavelet analysis for geophysical applications. Rev. Geophys., 35, 385–412
.
Lang, M., H. Guo., J. E. Odegard, C. S. Burrus, and R. O. Wells, 1995: Nonlinear processing of a shift invariant DWT for noise reduction. Proc. Wavelet Applications II, Orlando, FL, SPIE, 640–651
.
Lau, K., and H. Weng, 1995: Climate signal detection using wavelet transform: How to make a time series sing. Bull. Amer. Meteor. Soc., 76, 2391–2402
.
Liang, J., and T. W. Parks, 1996: A translation-invariant wavelet representation algorithm with applications. IEEE Trans. Sig. Proc., 44, 225–232
.
Liu, Q., 1994: On the definition and persistence of blocking. Tellus, 46A, 286–298
.
Lorenz, E. N., 1960: Maximum simplification of the dynamic equations. Tellus, 12, 243–254
.
Mak, M., 1995: Orthogonal wavelet analysis: Interannual variability in the sea surface temperature. Bull. Amer. Meteor. Soc., 76, 2179–2186
.
Meneveau, C., 1991: Analysis of turbulence in the orthonormal wavelet representation. J. Fluid Mech., 232, 469–520
.
Meyer, Y., 1993: Wavelets: Algorithms and Applications. SIAM, 133 pp
.
Nason, G., and B. Silverman, 1995: The stationary wavelet transform and some statistical applications. Wavelets and Statistics, A. Antoniadis and G. Oppenheim, Eds., Lecture Notes in Statistics, Vol. 103, Springer-Verlag, 281–300
.
Nelkin, M., 1994: Universality and scaling in fully developed turbulence. Adv. Phys., 43, 143–181
.
Onsager, L., 1949: Statistical hydrodynamics. Suppl. Serie IX Nuovo Cimento, 6, 279–287
.
Park, J., and M. Mann, 2000: Interannual temperature events and shifts in global temperature: A “multiwavelet” correlation approach. Earth Interactions, 4, 4-001. [Available online at http://EarthInteractions.org.]
.
Percival, D. B., and H. O. Mofjeld, 1997: Analysis of subtidal coastal sea level fluctuations using wavelets. J. Amer. Stat. Assoc., 92, 868–880
.
Perrier, V., and C. Basdevant, 1989: Periodical wavelet analysis, a tool for inhomogeneous field investigation. Theory and algorithms. Rech. Aérosp., 3, 53–67
.
Pesquet, J.-C., 1996: Time-invariant orthonormal wavelet representations. IEEE Trans. Sig. Proc., 44, 1964–1970
.
Platzman, G. W., 1962: The analytical dynamics of the spectral vorticity equation. J. Atmos. Sci., 19, 313–327
.
Prestin, J., and E. Quak, 1999: Optimal functions for a periodic uncertainty principle and multiresolution analysis. Proc. Edinburgh Math. Soc., 42, 225–242
.
Qi, Y., and H. Neumann, 1997: Wavelet analysis on errors of the bulk aerodynamic flux formula over canopy for GCMs. Mon. Wea. Rev., 125, 2238–2246
.
Richardson, L. F., 1922: Weather Prediction by Numerical Process. Cambridge, 236 pp
.
Saltzman, B., 1990: Three basic problems of paleoclimatic modeling:A personal perspective and review. Climate Dyn., 5, 67–78
.
Schneider, K., N. K.-R. Kevlahan, and M. Farge, 1997: Comparison of an adaptive wavelet method and nonlinearly filtered pseudo-spectral methods for the two-dimensional Navier-Stokes equations. Theor. Comput. Fluid Dyn., 9, 191–206
.
Shapiro, M. S., Ed., 1977: Mathematics Encyclopedia. Doubleday, 289 pp
.
Shensa, M. J., 1992: The discrete wavelet transform: Wedding the à trous and Mallat algorithms. IEEE Trans. Sig. Proc., 40, 2464–2482
.
Siggia, E. D., 1977: Origin of intermittency in fully developed turbulence. Phys. Rev. A, 15, 1730–1750
.
Sreenivasan, K. R., 1991: Fractals and multifractals in fluid turbulence. Annu. Rev. Fluid Mech., 23, 539–600
.
Strang, G., and T. Nguyen, 1995: Wavelets and Filter Banks. Wellesley-Cambridge, 490 pp
.
Strichartz, R. S., 1993: How to make wavelets. Amer. Math. Mon., 100, 539–556
.
Torrence, C., and G. Compo, 1998: A practical guide to wavelet analysis. Bull. Amer. Meteor. Soc., 79, 61–78
.
Tsou, C.-H., and P. J. Smith, 1990: The role of synoptic/planetary-scale interactions during the development of a blocking anticyclone. Tellus, 42A, 174–193
.
van den Berg, J. C., Ed., 1999: Wavelets in Physics. Cambridge University Press, 440 pp
.
von Sachs, R., G. P. Nason, and G. Kroisandt, 1998: Spectral representation and estimation for locally stationary wavelet processes. CRM Proc. Lecture Notes, 18, 381–397
.
Walden, A. T., and A. Contreras Cristan, 1998: The phase-corrected undecimated discrete wavelet packet transform and its application to interpreting the timing of events. Proc. Roy. Soc. London Ser. A, 454, 2243–2266
.
Wickerhauser, M. V., 1994: Adapted Wavelet Analysis from Theory to Software. A. K. Peters, 504 pp
.
Wiin Christensen, C., and A. Wiin-Nielsen, 1996: Blocking as a wave-wave interaction. Tellus, 48A, 254–271
.
Wiin-Nielsen, A., and T.-C. Chen, 1993: Spectral energetics of blocking. Fundamentals of Atmospheric Energetics, Oxford, 162–176
.
Zimin, V., 1981: Hierarchic model of turbulence. Izv. Atmos. Oceanic Phys., 17, 941–946
.
——, and F. Hussain, 1995: Wavelet based model for small-scale turbulence. Phys. Fluids, 7, 2925–2927
.
Zubair, L. M., 1993: Studies in turbulence using wavelet transforms for data compression and scale separation. Ph.D. dissertation, Yale University, 221 pp
.
——, K. Sreenivasan, and M. V. Wickerhauser, 1992: Characterization and compression of turbulent signals and images using wavelet-packets. Studies in Turbulence, T. B. Gatski, S. Sarkar, and C. G. Speziale, Eds., Springer-Verlag, 489–513
.

APPENDIX A

Periodization

The periodization (e.g., Wickerhauser 1994) of the orthonormal wavelet basis functions ψj,k on the real line is given in Table 2, where the periodization operator ( )per acts on a function A defined on the real line as

 
formula

and λε is a reference origin longitude that may be chosen to optimally align the OWA to the analyzed function (cf. section 2g).

APPENDIX B

Significance of OWA for Turbulence Modelers

All of the graphical depictions of (2.3a), (2.3b) in section 2c(2) illustrate another property of OWA that has attracted researchers in turbulence: a fractal-like self-similarity. Notice in Fig. 4c how the scaling function reproduces itself as a linear combination of smaller copies of itself (2.3b). A related linear combination generates the mother wavelet in Fig. 4d (2.3a). Turbulence cascade modelers have invoked less explicit (i.e., only statistical) self-similarity in order to obtain the observed scaling in turbulence, based in part on Richardson’s (1922) cascade hypothesis. For instance, the ratio of turbulent kinetic energy dissipation on two scales (wavenumbers 2jn, 2j−1n) is modeled by a probability distribution that, if independent of j, yields (statistically) self-similar cascades (Sreenivasan 1991). The continuous Fourier transform of (2.3b) reveals a ratio of values at wavenumbers 2jn′, 2j−1n′ equal to 2−1/2wF∗(2jπn′)〉n (discussed further by F3). Although this ratio is deterministic for a given wavelet, the fact that wavelets are generated by such a process has encouraged their application to turbulence studies, such as the CWT studies of Arneodo et al. (1999). Wavelets are explicitly self-similar by (2.2) and generated by a scaling process (2.3a), akin to a cascade, and have been conclusively shown to yield theoretical or numerical advantage when applied to the Navier–Stokes equations that govern turbulence (e.g., Federbush 1993; Schneider et al. 1997). Wavelets’ efficient representation of intermittent, fractal, or multifractal structures such as those observed in turbulent signals is another indication of their relevance to high Reynolds number flows. See Nelkin (1994) and Frisch (1995) for a review of turbulence theory and empirical results. Farge (1992) and Farge et al. (1999) review the application of wavelet analysis and synthesis to turbulence.

APPENDIX C

Formal Comparison of Fourier and Wavelet Expansions

It is interesting to compare the completely expanded Fourier and SOWA expressions for f★(x), that is,

 
formula

The top form is an expansion in amplitudes |n| times scalings by a factor n of special shifts by −argn of the mother wave cos(2πx) (Chui 1992); the bottom form is an expansion in amplitudes 2j/2j,k times scalings by a factor 2j of general shifts by −2j + k − 1 of the mother wavelet W(x). Thus in both cases scale is introduced as a multiplier of the argument x, and clearly the linear index n affords more scale resolution than the exponential index j. However, because of the global structure of cos, the individual shifts (phases) argn do not localize structures in x, but must be taken all together to do so. Because W is CS, each individual location index k does localize structures in λ, although it can be seen in Fig. 5a that for low j the periodization (A.1) slightly diffuses this localization.

APPENDIX D

Relation of the εj to the Choice of λε

The εj are related to the choice of λε and the shift equivariance (limited to a given resolution j) of the basis functions Wj,k as follows. Let

 
εfλfλλε
(D.1)

be a shift by λε, where is used to represent λε on the grid (2.1) as

 
formula

The limiting cases correspond to = (0, . . . , 0) and = (1, . . . , 1), respectively. It has been shown (e.g., Nason and Silverman 1995) that the coefficients j,k and fj,k are shifted in k by ε with the “detail bits” εj′≥j set to zero, that is, by

 
formula

from the coefficients and computed with = (0, . . . , 0). (Examples are provided by F1.) A corollary generalizes Shensa’s (1992) Eq. (3.6): for any ,

 
formula

precisely expressing the shift variance of OWA. Ultimately this implies that there are only (JR)2J different j,k, and 2JsgnR different fR,k, for all possible , that is, all possible λε.

Footnotes

* Current affiliation: Department of Meteorology, College of Computer, Mathematical, and Physical Sciences, University of Maryland, College Park, Maryland.

Corresponding author address: Aimé Fournier, National Center for Atmospheric Research, Boulder, CO 80307-3000.