Let us consider finite energy signals ∥f¯∥2=∫|f¯(x)|2dx∥f¯∥2=∫|f¯(x)|2dx of finite support, which is normalized to [0,1][0,1] or [0,1]2[0,1]2 for images. A sampling process implements a filtering of f¯(x)f¯(x) with a low-pass impulse response φ¯s(x)φ¯s(x) and a uniform sampling to output a discrete signal:

In two dimensions, n=(n1,n2)n=(n1,n2) and x=(x1,x2)x=(x1,x2). These filtered samples can also be written as inner products:

with φs(x)=φ¯s(-x)φs(x)=φ¯s(-x). Chapter 3 explains that φ_s is chosen, like in the classic Shannon–Whittaker sampling theorem, so that a family of functions {φs(x-ns)}1≤n≤N{φs(x-ns)}1≤n≤N is a basis of an appropriate approximation space U_N. The best linear approximation of f¯f¯ in U_N recovered from these samples is the orthogonal projection f¯Nf¯N of ff in U_N, and if the basis is orthonormal, then

A sampling theorem states that if f¯∈UNf¯∈UN then f¯=f¯Nf¯=f¯N so recovers f¯(x)f¯(x) from the measured samples. Most often, f¯f¯ does not belong to this approximation space. It is called aliasing in the context of Shannon–Whittaker sampling, where U_N is the space of functions having a frequency support restricted to the N lower frequencies. The approximation error ∥f¯-f¯N∥2∥f¯-f¯N∥2 must then be controlled.

The approximation error is computed by finding an orthogonal basis B={g¯m(x)}0≤m<+∞B={g¯m(x)}0≤m<+∞ of the whole analog signal space L2(R)[0,1]2L2(R)[0,1]2, with the first N vector {g¯m(x)}0≤m<N{g¯m(x)}0≤m<N that defines an orthogonal basis of U_N. Thus, the orthogonal projection on U_N can be rewritten as

Since f¯=∑m=0+∞〈f¯,g¯m〉g¯mf¯=∑m=0+∞〈f¯,g¯m〉g¯m, the approximation error is the energy of the removed inner products:

This error decreases quickly when N increases if the coefficient amplitudes |〈f¯,g¯m〉||〈f¯,g¯m〉| have a fast decay when the index m increases. The dimension N is adjusted to the desired approximation error.

Figure (a) shows a discrete image f[n]f[n] approximated with N=2562N=2562 pixels. Figure (c) displays a lower-resolution image fN/16fN/16 projected on a space UN/16UN/16 of dimension N/16N/16, generated by N/16N/16 large-scale wavelets. It is calculated by setting all the wavelet coefficients to zero at the first two smaller scales. The approximation error is ∥f-fN/16∥2/∥f∥2=14×10-3∥f-fN/16∥2/∥f∥2=14×10-3. Reducing the resolution introduces more blur and errors. A linear approximation space U_N corresponds to a uniform grid that approximates precisely uniform regular signals. Since images f¯f¯ are often not uniformly regular, it is necessary to measure it at a high-resolution N. This is why digital cameras have a resolution that increases as technology improves.

For M<NM<N, an approximation fMfM is computed by selecting the “best” M<NM<N vectors within BB. The orthogonal projection of ff on the space V_λ generated by M vectors {gm}m∈Λ{gm}m∈Λ in BB is

Since f=∑m∈γ〈f,gm〉gmf=∑m∈γ〈f,gm〉gm, the resulting error is

We write |λ||λ| the size of the set λ. The best M=|λ|M=|λ| term approximation, which minimizes ∥f-fλ∥2∥f-fλ∥2,  is thus obtained by selecting the M coefficients of largest amplitude. These coefficients are above a threshold T that depends on M:

This approximation is nonlinear because the approximation set λ_T changes with ff. The resulting approximation error is:

(Reference)(b) shows that the approximation support λ_T of an image in a wavelet orthonormal basis depends on the geometry of edges and textures. Keeping large wavelet coefficients is equivalent to constructing an adaptive approximation grid specified by the scale–space support λ_T. It increases the approximation resolution where the signal is irregular. The geometry of λ_T gives the spatial distribution of sharp image transitions and edges, and their propagation across scales. Chapter 6 proves that wavelet coefficients give important information about singularities and local Lipschitz regularity. This example illustrates how approximation support provides “geometric” information on ff, relative to a dictionary, that is a wavelet basis in this example.

(Reference)(d) gives the nonlinear wavelet approximation fMfM recovered from the M=N/16M=N/16 large-amplitude wavelet coefficients, with an error ∥f-fM∥2/∥f∥2=5×10-3∥f-fM∥2/∥f∥2=5×10-3. This error is nearly three times smaller than the linear approximation error obtained with the same number of wavelet coefficients, and the image quality is much better.

An analog signal can be recovered from the discrete nonlinear approxima-tion fMfM:

Since all projections are orthogonal, the overall approximation error on the original analog signal f¯(x)f¯(x) is the sum of the analog sampling error and the discrete nonlinear error:

In practice, N is imposed by the resolution of the signal-acquisition hardware, and M is typically adjusted so that ϵn(M,f)≥ϵl(N,f)ϵn(M,f)≥ϵl(N,f).

Sparse representations are obtained in a basis that takes advantage of some form of regularity of the input signals, creating many small-amplitude coefficients. Since wavelets have localized support, functions with isolated singularities produce few large-amplitude wavelet coefficients in the neighborhood of these singularities. Nonlinear wavelet approximation produces a small error over spaces of functions that do not have “too many” sharp transitions and singularities. Chapter 9 shows that functions having a bounded total variation norm are useful models for images with nonfractal (finite length) edges.

Edges often define regular geometric curves. Wavelets detect the location of edges but their square support cannot take advantage of their potential geometric regularity. More sparse representations are defined in dictionaries of curvelets or bandlets, which have elongated support in multiple directions, that can be adapted to this geometrical regularity. In such dictionaries, the approximation support λ_T is smaller but provides explicit information about edges' local geometrical properties such as their orientation. In this context, geometry does not just apply to multidimensional signals. Audio signals, such as musical recordings, also have a complex geometric regularity in time-frequency dictionaries.

To optimize the estimation operator D, one must take advantage of prior information available about signal ff. In a Bayes framework, ff is considered a realization of a random vector FF and the Bayes risk is the expected risk calculated with respect to the prior probability distribution π of the random signal model FF:

Optimizing D among all possible operators yields the minimum Bayes risk:

In the 1940s, Wald brought in a new perspective on statistics with a decision theory partly imported from the theory of games. This point of view uses deterministic models, where signals are elements of a set Θ, without specifying their probability distribution in this set. To control the risk for any f∈Θf∈Θ, we compute the maximum risk:

The minimax risk is the lower bound computed over all operators D:

In practice, the goal is to find an operator D that is simple to implement and yields a risk close to the minimax lower bound.

Figure 1

It is tempting to restrict calculations to linear operators D because of their simplicity. Optimal linear Wiener estimators are introduced in Chapter 11. Figure (a) is an image contaminated by Gaussian white noise. Figure (b) shows an optimized linear filtering estimation F˜=X☆h[n]F˜=X☆h[n], which is therefore diagonal in a Fourier basis BB. This convolution operator averages the noise but also blurs the image and keeps low-frequency noise by retaining the image's low frequencies.

If ff has a sparse representation in a dictionary, then projecting X on the vectors of this sparse support can considerably improve linear estimators. The difficulty is identifying the sparse support of ff from the noisy data X. Donoho and Johnstone (DonohoJ:94) proved that, in an orthonormal basis, a simple thresholding of noisy coefficients does the trick. Noisy signal coefficients in an orthonormal basisB={gm}m∈ΓB={gm}m∈Γ are

Thresholding these noisy coefficients yields an orthogonal projection estimator

The set Λ˜TΛ˜T is an estimate of an approximation support of ff. It is hopefully close to the optimal approximation support λT={m∈γ:|〈f,gm〉|≥T}λT={m∈γ:|〈f,gm〉|≥T}.

Figure 1(b) shows the estimated approximation set λ˜Tλ˜T of noisy-wavelet coefficients, |〈X,ψj,n|≥T|〈X,ψj,n|≥T, that can be compared to the optimal approximation support Λ_T shown in (Reference)(b). The estimation in Figure 1(d) from wavelet coefficients in λ˜Tλ˜T has considerably reduced the noise in regular regions while keeping the sharpness of edges by preserving large-wavelet coefficients. This estimation is improved with a translation-invariant procedure that averages this estimator over several translated wavelet bases. Thresholding wavelet coefficients implements an adaptive smoothing, which averages the data X with a kernel that depends on the estimated regularity of the original signal ff.

Donoho and Johnstone proved that for Gaussian white noise of variance σ², choosing T=σ2logeNT=σ2logeN yields a risk E{∥f-F˜∥2}E{∥f-F˜∥2} of the order of ∥f-fΛT∥2∥f-fΛT∥2, up to a logeNlogeN factor. This spectacular result shows that the estimated support λ˜Tλ˜T does nearly as well as the optimal unknown support λ_T. The resulting risk is small if the representation is sparse and precise.

The set λ˜Tλ˜T in Figure 1(b) “looks” different from the λ_T in (Reference)(b) because it has more isolated points. This indicates that some prior information on the geometry of λ_T could be used to improve the estimation. For audio noise-reduction, thresholding estimators are applied in sparse representations provided by time-frequency bases. Similar isolated time-frequency coefficients produce a highly annoying “musical noise.” Musical noise is removed with a block thresholding that regularizes the geometry of the estimated support λ˜Tλ˜T and avoids leaving isolated points. Block thresholding also improves wavelet estimators.

If W is a Gaussian noise and signals in Θ have a sparse representation in BB, then Chapter 11 proves that thresholding estimators can produce a nearly minimax risk. In particular, wavelet thresholding estimators have a nearly minimax risk for large classes of piecewise smooth signals, including bounded variation images.