format: text

text:
# 1 Introduction

## 1 Introduction

Wavelets became a necessary mathematical tool in many investigations. They are used in those cases when the result of the analysis of a particular signal111The notion of a signal is used here for any ordered set of numerically recorded information about some processes, objects, functions etc. The signal can be a function of some coordinates, would it be the time, the space or any other (in general, $`n`$-dimensional) scale. should contain not only the list of its typical frequencies (scales) but also the knowledge of the definite local coordinates where these properties are important. Thus, analysis and processing of different classes of nonstationary (in time) or inhomogeneous (in space) signals is the main field of applications of wavelet analysis . The most general principle of the wavelet construction is to use dilations and translations. Commonly used wavelets form a complete orthonormal system of functions with a finite support constructed in such a way. That is why by changing a scale (dilations) they can distinct the local characteristics of a signal at various scales, and by translations they cover the whole region in which it is studied. Due to the completeness of the system, they also allow for the inverse transformation to be done. In the analysis of nonstationary signals, the locality property of wavelets leads to their substantial advantage over Fourier transform which provides us only with the knowledge of global frequencies (scales) of the object under investigation because the system of the basic functions used (sine, cosine or imaginary exponential functions) is defined on the infinite interval222Comparison of the wavelet transform with the so-called windowed (within the finite interval) Fourier transform will be briefly discussed below.. However, as we shall see, the more general definitions and, correspondingly, a variety of forms of wavelets are used which admit a wider class of functions to be considered. According to Y. Meyer , ”the wavelet bases are universally applicable: ”everything that comes to hand”, whether function or distribution, is the sum of a wavelet series and, contrary to what happens with Fourier series, the coefficients of the wavelet series translate the properties of the function or distribution simply, precisely and faithfully.”

The literature devoted to wavelets is very voluminous, and one can easily get a lot of references by sending the corresponding request to Internet web sites. Mathematical problems are treated in many monographs in detail (e.g., see ). The introductory courses on wavelets can be found in the books . The nice review paper adapted for beginners and practical users with demonstration of wavelet transform of some signals was published in this journal about four years ago  and attracted much attention. However continuous wavelets are mostly considered there whereas discrete wavelets are just briefly mentioned. This choice was dictated by the fact that the continuous wavelets admit somewhat more visual and picturesque presentation of the results of the analysis of a signal in terms of local maxima and sceleton graphs of wavelet coefficients with continuous variables.

At the same time, the main bulk of papers dealing with practical applications of wavelet analysis uses discrete wavelets which will be our main concern here. This preference of discrete wavelets is related to the fact that widely used continuous wavelets are not, strictly speaking, orthonormal because they are both infinitely differentiable and exponentially decreasing at infinity what violates orthonormalization property whereas there is no such problem for discrete wavelets. That is why the discrete wavelets allow often for more accurate transform and presentation of the analyzed signal, and, especially, for its inverse transform after the compression procedure. Moreover, they are better adapted to communication theory and practice. These comments do not imply that we insist on use of only discrete wavelets for signal analysis . On the contrary, the continuous wavelets may provide sometimes more transparent and analytical results in modelling the signal analysis than the discrete wavelets.

The choice of a specific wavelet , would it be discrete or continuous one, depends on the analyzed signal and on the problem to be solved. Some functions are best analyzed using one method or another, and it depends on the relative simplicity of the decomposition achieved. The researcher’s intuition helps resolve this ”miracle”. As an analogy, the example with number systems is often considered. It is a matter of tradition and convenience to choose the systems with the base 10, 2 or $`e`$. However, the Roman number system is completely excluded if one tries to use it for multiplication. At the same time, different problems can ask for more or less efforts both in their solution and in graphical presentation depending on the system chosen, and our intuition is important there.

The programs exploiting the wavelet transform are widely used now not only for scientific research but for commercial purposes as well. Some of them have been even described in books (e.g., see ). At the same time, the direct transition from pure mathematics to computer programming and applications is non-trivial and asks often for the individual approach to the problem under investigation and for a specific choice of wavelets used. Our main objective here is to describe in a suitable way the bridge that relates mathematical wavelet constructions to practical signal processing.

The discrete wavelets look as strangers to all those accustomed to analytical calculations because they can not be represented by analytical expressions (except of a simplest one) or by solutions of some differential equations, and instead are given numerically as solutions of definite functional equations containing rescaling and translations. Moreover, in practical calculations their direct form is not required even, and the numerical values of the coefficients of the functional equation are only used. Thus wavelets are defined by the iterative algorithm of the dilation and translation of a single function. This leads to a very important procedure called multiresolution analysis which gives rise to the multiscale local analysis of the signal and fast numerical algorithms. Each scale contains an independent non-overlapping set of information about the signal in form of wavelet coefficients, which are determined from an iterative procedure called fast wavelet transform. In combination, they provide its complete analysis and the diagnosis of the underlying processes.

After such an analysis was done, one can compress (if necessary) the resulting data by omitting some inessential part of the encoded information. This is done with the help of the so-called quantization procedure which commonly allocates different weights to various wavelet coefficients obtained. In particular, it helps erase some statistical fluctuations and, therefore, increase the role of dynamical features of a signal. It can however falsify the diagnostic if the compression is done inappropriately. Usually, the accurate compression gives rise to a substantial reduction of the required computer storage memory and transmission facilities, and, consequently, to a lower expenditure. The number of vanishing moments of wavelets is important at this stage. Unfortunately, the compression introduces unavoidable systematic errors. The mistakes one has made will consist of multiples of the deleted wavelet coefficients, and, therefore, the regularity properties of a signal play an essential role. Reconstruction after such compression schemes is then not perfect any more. These two objectives are clearly antagonistic. Nevertheless, when one tries to reconstruct the initial signal, the inverse transformation (synthesis) happens to be rather stable and reproduces its most important characteristics if proper methods are applied. The regularity properties of wavelets used become also crucial at the reconstruction stage. The distortions of the reconstucted signal due to quantization can be kept small, although significant compression ratios are attained. Since the part of the signal which is not reconstructed is often called noise, in essence, what we are doing is denoising signals. Namely at this stage the superiority of the discrete wavelets becomes especially clear.

Thus, the objectives of signal processing consist in the accurate analysis , effective coding, fast transmission and, finally, careful reconstruction (at the transmission destination point) of the initial signal. Sometimes the first stage of signal analysis and diagnosis is enough for the problem to be solved and anticipated goals to be achieved.

It has been proven that any function can be written as a superposition of admissible wavelets, and there exists a numerically stable algorithm to compute the coefficients for such an expansion. Moreover, these coefficients completely characterize the function, and it is possible to reconstruct it in a numerically stable way by knowing these coefficients. Because of their unique properties, wavelets were used in functional analysis in mathematics, in studies of (multi)fractal properties, singularities and local oscillations of functions, for solving some differential equations, for investigation of inhomogeneous processes involving widely different scales of interacting perturbations, for pattern recognition, for image and sound compression, for digital geometry processing, for solving many problems of physics, biology, medicine, technique etc (see the recently published books ). This list is by no means exhaustive.

One should however stress that, even though this method is very powerful, the goals of wavelet analysis are rather modest. It helps us describe and reveal some features, otherwise hidden in a signal, but it does not pretend to explain the underlying dynamics and physics origin although it may give some crucial hints to it. Wavelets present a new stage in optimization of this description providing the best known representation of a signal. With the help of wavelets, we merely see things a little more clearly. To understand the dynamics, standard approaches introduce models assumed to be driving the mechanisms generating the observations. To define the optimality of the algorithms of wavelet analysis , some (still debatable!) energy and entropy criteria have been developed. They are internal to the algorithm itself. However, the choice of the best algorithm is also tied to an objective goal of its practical use, i.e., to some external criteria. That is why in practical applications one should submit the performance of a ”theoretically optimal algorithm” to the judgements of experts and users to estimate its advantage over the previously developed ones.

Despite a very active research and impressive results, the versatility of wavelet analysis and existence of many different inversion formulas for wavelet transform (redundancy) imply that these studies are presumably not in their final form yet. We shall try to describe the situation in its status nascendi.

The main part of this paper (Sections 2–14) is devoted to description of general properties of wavelets and use of the wavelet transform in computer calculations. Some applications to different fields are briefly described in Section 15.

## 2 Wavelets for beginners

Each signal can be characterized by its averaged (over some intervals) values (trend) and by its variations around this trend. Let us call these variations as fluctuations independently of their nature, would it be dynamic, stochastic, psychological, physiological or any other origin. When processing a signal, one is interested in its fluctuations at various scales because therefrom one can learn about their origin. The goal of the wavelet analysis is to provide tools for such processing.

Actually, physicists dealing with experimental histograms analyze their data at different scales when averaging over different size intervals. This is a particular example of a simplified wavelet analysis treated in this Section. To be more definite, let us consider the situation when an experimentalist measures some function $`f(x)`$ within the interval $`0x1`$, and the best resolution obtained with the measuring device is limited by 1/16th of the whole interval. Thus the result consists of 16 numbers representing the mean values of $`f(x)`$ in each of these bins and can be plotted as a 16-bin histogram shown in the upper part of Fig. 1. It can be represented by the following formula

$$f(x)=\underset{k=0}{\overset{15}{}}s_{4,k}\phi _{4,k}(x),$$
(1)
where $`s_{4,k}=f(k/16)/4`$, and $`\phi _{4,k}`$ is defined as a step-like block of the unit norm (i.e. of a height 4 and a widths 1/16) different from zero only within the $`k`$-th bin. For an arbitrary $`j`$, one imposes the condition $`𝑑x|\phi _{j,k}|^2=1`$, where the integral is taken over the intervals of the lengths $`\mathrm{\Delta }x_j=1/2^j`$ and, therefore, $`\phi _{j,k}`$ have the following form $`\phi _{j,k}=2^{j/2}\phi (2^jxk)`$ with $`\phi `$ denoting a step-like function of the unit height in such an interval. The label 4 is related to the total number of such intervals in our example. At the next coarser level the average over the two neighboring bins is taken as is depicted in the histogram just below the initial one in Fig. 1. Up to the normalization factor, we will denote it as $`s_{3,k}`$ and the difference between the two levels shown to the right of this histogram as $`d_{3,k}`$. To be more explicit, let us write down the normalized sums and differences for an arbitrary level $`j`$ as

$$s_{j1,k}=\frac{1}{\sqrt{2}}[s_{j,2k}+s_{j,2k+1}];d_{j1,k}=\frac{1}{\sqrt{2}}[s_{j,2k}s_{j,2k+1}],$$
(2)
or for the backward transform (synthesis)

$$s_{j,2k}=\frac{1}{\sqrt{2}}(s_{j1,k}+d_{j1,k});s_{j,2k+1}=\frac{1}{\sqrt{2}}(s_{j1,k}d_{j1,k}).$$
(3)
Since, for the dyadic partition considered, this difference has opposite signs in the neighboring bins of the previous fine level, we introduce the function $`\psi `$ which is 1 and -1, correspondingly, in these bins and the normalized functions $`\psi _{j,k}=2^{j/2}\psi (2^jxk)`$. It allows us to represent the same function $`f(x)`$ as

$$f(x)=\underset{k=0}{\overset{7}{}}s_{3,k}\phi _{3,k}(x)+\underset{k=0}{\overset{7}{}}d_{3,k}\psi _{3,k}(x).$$
(4)
One proceeds further in the same manner to the sparser levels 2, 1 and 0 with averaging done over the interval lengths 1/4, 1/2 and 1, correspondingly. It is shown in the subsequent drawings in Fig. 1. The most sparse level with the mean value of $`f`$ over the whole interval denoted as $`s_{0,0}`$ provides

$`f(x)=s_{0,0}\phi _{0,0}(x)+d_{0,0}(x)\psi _{0,0}(x)+{\displaystyle \underset{k=0}{\overset{1}{}}}d_{1,k}\psi _{1,k}(x)`$
$`+{\displaystyle \underset{k=0}{\overset{3}{}}}d_{2,k}\psi _{2,k}(x)+{\displaystyle \underset{k=0}{\overset{7}{}}}d_{3,k}\psi _{3,k}(x).`$ (5)
The functions $`\phi _{0,0}(x)`$ and $`\psi _{0,0}(x)`$ are shown in Fig. 2. The functions $`\phi _{j,k}(x)`$ and $`\psi _{j,k}(x)`$ are the normalized by the conservation of the norm, dilated and translated versions of them. In the next Section we will give explicit formulae for them in a particular case of Haar scaling functions and wavelets. In practical signal processing, these functions (and their more sophisticated versions) are often called as the low and high-path filters, correspondingly, because they filter the large and small scale components of a signal. The subsequent terms in Eq. (5) show the fluctuations (differences $`d_{j,k}`$) at finer and finer levels with larger $`j`$. In all the cases (1)–(5) one needs exactly 16 coefficients to represent the function. In general, there are $`2^j`$ coefficients $`s_{j,k}`$ and $`2^{j_n}2^j`$ coefficients $`d_{j,k}`$, where $`j_n`$ denotes the finest resolution level (in the above example, $`j_n=4`$).

All the above representations of the function $`f(x)`$ (Eqs. (1)-(5)) are mathematically equivalent. However, the latter one representing the wavelet analyzed function directly reveals the fluctuation structure of the signal at different scales $`j`$ and various locations $`k`$ present in a set of coefficients $`d_{j,k}`$ whereas the original form (1) hides the fluctuation patterns in the background of a general trend. The final form (5) contains the overall average of the signal depicted by $`s_{0,0}`$ and all its fluctuations with their scales and positions well labelled by 15 normalized coefficients $`d_{j,k}`$ while the initial histogram shows only the normalized average values $`s_{j,k}`$ in the 16 bins studied. Moreover, in practical applications the latter wavelet representation is preferred because for rather smooth functions, strongly varying only at some discrete values of their arguments, many of the high-resolution $`d`$-coefficients in relations similar to Eq. (5) are close to zero and can be discarded. Bands of zeros (or close to zero values) indicate those regions where the function is smooth enough.

At first sight, this simplified example looks somewhat trivial. However, for more complicated functions and more data points with some elaborate forms of wavelets it leads to a detailed analysis of a signal and to possible strong compression with subsequent good quality restoration. This example provides also the illustration for the very important feature of the whole approach with successive coarser and coarser approximations to $`f`$ called the multiresolution analysis and discussed in more detail below.

## 3 Basic notions and Haar wavelets

To analyze any signal, one should, first of all, choose the corresponding basis, i.e., the set of functions to be considered as ”functional coordinates”. In most cases we will deal with signals represented by the square integrable functions defined on the real axis (or by the square summable sequences of complex numbers). They form the infinite-dimensional Hilbert space $`L^2(R)`$ ($`l^2(Z)`$). The scalar product of these functions is defined as

$$f,g=_{\mathrm{}}^{\mathrm{}}𝑑xf(x)\overline{g}(x),$$
(6)
where the bar stands for the complex conjugate. For finite limits in the sum it becomes a finite-dimensional Hilbert space. Hilbert spaces always have orthonormal bases, i.e., families of vectors (or functions) $`e_n`$ such that

$$e_n,e_m=\delta _{nm},$$
(7)
$$f^2=𝑑x|f(x)|^2=\underset{n}{}|f,e_n|^2.$$
(8)
In the Hilbert space, there exist some more general families of linear independent basis vectors called a Riesz basis which generalize the equation (8) to

$$\alpha f^2\underset{n}{}|f,e_n|^2\beta f^2$$
(9)
with $`\alpha >0,\beta <\mathrm{}`$. It is the unconditional basis where the order in which the basis vectors are considered does not matter. Any bounded operator with a bounded inverse maps an orthonormal basis into a Riesz basis.

Sometimes we will consider the spaces $`L^p(R)(1p<\mathrm{};p2)`$, where the norm is defined as

$$f_{L^p}=[𝑑x|f(x)|^p]^{1/p},$$
(10)
as well as some other Banach spaces333Those are linear spaces equipped with norm which generally is not derived from a scalar product and complete with respect to this norm ($`L^p(R)(1p<\mathrm{};p2)`$ is a particular example of it.)  (see Sections 11, 12).

The Fourier transform with its trigonometric basis is well suited for the analysis of stationary signals. Then the norm $`f`$ is often called energy. For nonstationary signals, e.g., the location of that moment when the frequency characteristics has abruptly been changed is crucial. Therefore the basis should have a compact support. The wavelets are just such functions which span the whole space by translation of the dilated versions of a definite function. That is why every signal can be decomposed in the wavelet series (or integral). Each frequency component is studied with a resolution matched to its scale. The above procedure of normalization of functions $`\phi _{j,k}`$ is directly connected with the requirement of conservation of the norm of a signal at its decompositions.

The choice of the analyzing wavelet is, however, not unique. One should choose it in accordance with the problem to be solved. The simplicity of operations (computing, in particular) and of representation (minimum parameters used) plays also a crucial role. Bad choice may even prevent from getting any answer as in the above example with Roman numbers. There are several methods of estimating how well the chosen function is applicable to the solution of a particular problem (see Section 6).

Let us try to construct functions satisfying the above criteria. ”An educated guess” would be to relate a function $`\phi (x)`$ to its dilated and translated version. The simplest linear relation with $`2M`$ coefficients is

$$\phi (x)=\sqrt{2}\underset{k=0}{\overset{2M1}{}}h_k\phi (2xk)$$
(11)
with the dyadic dilation 2 and integer translation $`k`$. At first sight, the chosen normalization of the coefficients $`h_k`$ with the factor $`\sqrt{2}`$ looks somewhat arbitrary. Actually, it is defined a’posteriori by the traditional form of fast algorithms for their calculation (see Eqs. (40) and (41) below) and normalization of functions $`\phi _{j,k}(x),\psi _{j,k}(x)`$. It is used in all the books cited above. However, sometimes (see , Chapter 7) it is replaced by $`c_k=\sqrt{2}h_k`$.

For discrete values of the dilation and translation parameters one gets discrete wavelets. The value of the dilation factor determines the size of cells in the lattice chosen. The integer $`M`$ defines the number of coefficients and the length of the wavelet support. They are interrelated because from the definition of $`h_k`$ for orthonormal bases

$$h_k=\sqrt{2}𝑑x\phi (x)\overline{\phi }(2xk)$$
(12)
it follows that only finitely many $`h_k`$ are nonzero if $`\phi `$ has a finite support. The normalization condition is chosen as

$$_{\mathrm{}}^{\mathrm{}}𝑑x\phi (x)=1.$$
(13)
The function $`\phi (x)`$ obtained from the solution of this equation is called a scaling function444It is often called also a ”father wavelet ” but we will not use this term.. If the scaling function is known, one can form a ”mother wavelet ” (or a basic wavelet ) $`\psi (x)`$ according to

$$\psi (x)=\sqrt{2}\underset{k=0}{\overset{2M1}{}}g_k\phi (2xk),$$
(14)
where

$$g_k=(1)^kh_{2Mk1}.$$
(15)
The simplest example would be for $`M=1`$ with two non-zero coefficients $`h_k`$ equal to $`1/\sqrt{2}`$, i.e., the equation leading to the Haar scaling function $`\phi _H(x)`$:

$$\phi _H(x)=\phi _H(2x)+\phi _H(2x1).$$
(16)
One easily gets the solution of this functional equation

$$\phi _H(x)=\theta (x)\theta (1x),$$
(17)
where $`\theta (x)`$ is the Heaviside step-function equal to 1 at positive arguments and 0 at negative ones. The additional boundary condition is $`\phi _H(0)=1,\phi _H(1)=0`$. This condition is important for the simplicity of the whole procedure of computing the wavelet coefficients when two neighboring intervals are considered.

The ”mother wavelet ” is

$$\psi _H(x)=\theta (x)\theta (12x)\theta (2x1)\theta (1x).$$
(18)
with boundary values defined as $`\psi _H(0)=1,\psi _H(1/2)=1,\psi _H(1)=0`$. This is the Haar wavelet  known since 1910 and used in the functional analysis . Namely this example has been considered in the previous Section fot the histogram decomposition. Both the scaling function $`\phi _H(x)`$ and the ”mother wavelet ” $`\psi _H(x)`$ are shown in Fig. 2. It is the first one of a family of compactly supported orthonormal wavelets $`{}_{M}{}^{}\psi :\psi _H=_1\psi `$. It possesses the locality property since its support $`2M1=1`$ is compact.

The dilated and translated versions of the scaling function $`\phi `$ and the ”mother wavelet” $`\psi `$

$$\phi _{j,k}=2^{j/2}\phi (2^jxk),$$
(19)
$$\psi _{j,k}=2^{j/2}\psi (2^jxk)$$
(20)
form the orthonormal basis as can be (easily for Haar wavelets) checked555We return back to the general case and therefore omit the index $`H`$ because the same formula will be used for other wavelets.. The choice of $`2^j`$ with the integer valued $`j`$ as a scaling factor leads to the unique and selfconsistent procedure of computing the wavelet coefficients. In principle, there exists an algorithm of derivation of the compact support wavelets with an arbitrary rational number in place of 2. However, only for this factor, it has been shown that there exists the explicit algorithm with the regularity of a wavelet increasing linearly with its support. For example, for the factor 3 the regularity index only grows logarithmically. The factor 2 is probably distinguished here as well as in music where octaves play a crucial role. If the dilation factor is 2, then the Fourier transform of the ”mother” wavelet is essentially localized between $`\pi `$ and $`2\pi `$. However, for some practical applications, a sharper frequency localization is necessary, and it may be useful to have wavelet bases with a narrower bandwidth. The fractional dilation wavelet bases provide one of the solutions of this problem but there also exist other possibilities.

The Haar wavelet oscillates so that

$$_{\mathrm{}}^{\mathrm{}}𝑑x\psi (x)=0.$$
(21)
This condition is common for all the wavelets. It is called the oscillation or cancellation condition. From it, the origin of the name wavelet becomes clear. One can describe a ”wavelet” as a function that oscillates within some interval like a wave but is then localized by damping outside this interval. This is a necessary condition for wavelets to form an unconditional (stable) basis. We conclude that for special choices of coefficients $`h_k`$ one gets the specific forms of ”mother” wavelets and orthonormal bases.

One may decompose any function $`f`$ of $`L^2(R)`$ at any resolution level $`j_n`$ in a series

$$f=\underset{k}{}s_{j_n,k}\phi _{j_n,k}+\underset{jj_n,k}{}d_{j,k}\psi _{j,k}.$$
(22)
At the finest resolution level $`j_n=j_{max}`$ only $`s`$-coefficients are left, and one gets the scaling-function representation

$$f(x)=\underset{k}{}s_{j_{max},k}\phi _{j_{max},k}.$$
(23)
In the case of the Haar wavelets it corresponds to the initial experimental histogram with finest resolution. Since we will be interested in its analysis at varying resolutions, this form is used as an initial input only. The final representation of the same data (22) shows all the fluctuations in the signal. The wavelet coefficients $`s_{j,k}`$ and $`d_{j,k}`$ can be calculated as

$$s_{j,k}=𝑑xf(x)\phi _{j,k}(x),$$
(24)
$$d_{j,k}=𝑑xf(x)\psi _{j,k}(x).$$
(25)
However, in practice their values are determined from the fast wavelet transform described below.

In reference to the particular case of the Haar wavelet , considered above, these coefficients are often referred as sums ($`s`$) and differences ($`d`$), thus related to mean values and fluctuations.

For physicists familiar with experimental histograms, it generalizes the particular example discussed in the previous Section. The first sum in (22) with the scaling functions $`\phi _{j,k}`$ shows the average666The averaging is done with the weight functions $`\phi _{j,k}(x)`$. values of $`f`$ within the dyadic intervals $`[k2^j,(k+1)2^j)`$, and the second term contains all the fluctuations of the function $`f`$ in this interval. They come from ever smaller intervals which correspond to larger values of the scale parameter $`j`$. One would say that it ”brings into focus” the finer details of a signal. This touching in of details is regularly spaced, taking account of dimension - the details of dimension $`2^j`$ are placed at the points $`k2^j`$. At the lowest (most sparse) level $`j_0`$ the former sum consists of a single term with the overall weighted average $`f=s_{j_0,k_0}`$, where $`k_0`$ is the center of the histogram. The second sum in (22) shows fluctuations at all the scales. At the next, more refined level $`j_1>j_0`$, there are two terms in the first sum which show the average values of $`f`$ within half-intervals with their centers positioned at $`k_1,k_2`$. The number of terms in the second sum becomes less by one term which was previously responsible for the fluctuations at the half-interval scale. The total number of terms in the expansion stays unchanged. Here we just mention that according to (22) the number of terms in each sum depends on a definite resolution level. Changing the level index by 1, we move some terms to another sum, and each of these representations are ”true” representations of the histogram at different resolution levels.

Formally a similar procedure may be done the other way round by going to the sparser resolutions $`j<j_0`$. Even if we ”fill out” the whole support of $`f`$, we can still keep going with our averaging trick. Then the average value of $`f`$ diminishes, and one can neglect the first sum in (22) in the $`L^2`$ sense because its $`L^2`$-norm (8) tends to zero. In the histogram example, it decreases as $`|f|N^1`$ and $`|f|^2N^2`$ whereas the integration region is proportional to $`N`$, i.e., $`f^2N^10`$ for $`N\mathrm{}`$. That is why only the second term in (22) is often considered, and the result is often called as the wavelet expansion. It also works if $`f`$ is in the space $`L^p(R)`$ for $`1<p<\mathrm{}`$ but it can not be done if it belongs to $`L^1(R)`$ or $`L^{\mathrm{}}(R)`$. For example, if $`f=1`$ identically, all wavelet coefficients $`d_{j,k}`$ are zero, and only the first sum matters. For the histogram interpretation, the neglect of this sum would imply that one is not interested in its average value but only in its shape determined by fluctuations at different scales. Any function can be approximated to a precision $`2^{j/2}`$ (i.e., to an arbitrary high precision at $`j\mathrm{}`$) by a finite linear combination of Haar wavelets.

The Haar wavelets are also suitable for the studies of functions belonging to $`L^p`$ spaces, i.e., possessing higher moments.

Though the Haar wavelets provide a good tutorial example of an orthonormal basis, they suffer from several deficiences. One of them is the bad analytic behavior with the abrupt change at the interval bounds, i.e., its bad regularity properties. By this we mean that moments of the Haar wavelet are different from zero - only the integral (21) of the function itself is zero. The Haar wavelet does not have good time-frequency localization. Its Fourier transform decays like $`|\omega |^1`$ for $`\omega \mathrm{}`$.

It would be desirable to build up wavelets with better regularity. The advantage of them compared to the Haar system shows up in the smaller number of the wavelet coefficients which are large enough to account for and in their applicability to a wider set of functional spaces besides $`L^2`$. The former feature is related to the fact that the wavelet coefficients are significantly different from zero only near singularities of $`f`$ (strong fluctuations!). Therefore wavelet series of standard functions with isolated singularities are ”sparse” series in contrast to the Fourier series which are usually dense ones for rather regular functions. The latter feature allows us to get an access to local and global regularities of functions under investigation. The way to this program was opened by the multiresolution analysis .

## 4 Multiresolution analysis and Daubechies wavelets

The relation (22) shows that a general function $`f`$ can be approximated by a sequence of very simple functions $`\phi _{j,k},\psi _{j,k}`$. The above example has demonstrated that the Haar functions are local and cover the whole space $`L^2(R)`$ by using the translation $`k`$. They are orthogonal for different resolution scales $`j`$. The transition from $`j`$ to $`j+1`$ is equivalent to the replacement of $`x`$ by $`2x`$, i.e., to the rescaling which allows for the analysis to be done at various resolutions.

However the Haar wavelets are oversimplified and not regular enough. The goal is to find a general class of those functions which would satisfy the requirements of locality, regularity and oscillatory behavior. Note that in some particular cases the orthonormality property sometimes can be relaxed. They should be simple enough in the sense that they are of being sufficiently explicit and regular to be completely determined by their samples on the lattice defined by the factors $`2^j`$.

The general approach which respects these properties is known as the multiresolution approximation. A rigorous mathematical definition is given in Appendix 1. Here we just describe its main ingredients.

The multiresolution analysis consists of a sequence of successive approximation spaces $`V_j`$ which are scaled and invariant under integer translation versions of the central functional space $`V_0`$. To explain the meaning of these spaces in a simple example, we show in Fig. 3 what the projections of some function on the Haar spaces $`V_0`$, $`V_1`$ might look like. One easily recognizes the histogram representation of this function. The comparison of histograms at the two levels shows that the first sum in Eq. (22) provides the ”blurred image” or ”smoothed means” of $`f(x)`$ in each interval, while the second sum of this equation adds finer and finer details of smaller sizes. Thus the general distributions are decomposed into the series of correctly localized fluctuations having a characteristic form defined by the wavelet chosen.

The functions $`\phi _{j,k}`$ form an orthonormal basis in $`V_j`$. The orthogonal complement of $`V_j`$ in $`V_{j+1}`$ is called $`W_j`$. The subspaces $`W_j`$ form a mutually orthogonal set. The sequence of $`\psi _{j,k}`$ constitutes an orthonormal basis for $`W_j`$ at any definite $`j`$. The whole collection of $`\psi _{j,k}`$ and $`\phi _{j,k}`$ for all $`j`$ is an orthonormal basis for $`L^2(R)`$. This ensures us that we have constructed a multiresolution analysis approach, and the functions $`\psi _{j,k}`$ and $`\phi _{j,k}`$ constitute the small and large scale filters, correspondingly. The whole procedure of the multiresolution analysis can be demonstrated in graphs of Fig. 4.

In accordance with the above declared goal, one can define the notion of wavelets (see Appendix 1) so that the functions $`2^{j/2}\psi (2^jxk)`$ are the wavelets (generated by the ”mother” $`\psi `$), possessing the regularity, the localization and the oscillation properties.

At first sight, from our example with Haar wavelets, it looked as if one is allowed to choose coefficients $`h_k`$ at his own wish. This impression is, however, completely wrong. The general properties of scaling functions and wavelets define these coefficients in a unique way in the framework of the multiresolution analysis approach.

Let us show how the program of the multiresolution analysis works in practice when applied to the problem of finding out the coefficients of any filter $`h_k`$ and $`g_k`$. They can be directly obtained from the definition and properties of the discrete wavelets. These coefficients are defined by the relations (11) and (14)

$$\phi (x)=\sqrt{2}\underset{k}{}h_k\phi (2xk);\psi (x)=\sqrt{2}\underset{k}{}g_k\phi (2xk),$$
(26)
where $`_k|h_k|^2<\mathrm{}`$. The orthogonality of the scaling functions defined by the relation

$$𝑑x\phi (x)\phi (xm)=0$$
(27)
leads to the following equation for the coefficients:

$$\underset{k}{}h_kh_{k+2m}=\delta _{0m}.$$
(28)
The orthogonality of wavelets to the scaling functions

$$𝑑x\psi (x)\phi (xm)=0$$
(29)
gives the equation

$$\underset{k}{}h_kg_{k+2m}=0,$$
(30)
having a solution of the form

$$g_k=(1)^kh_{2M1k}.$$
(31)
Another condition of the orthogonality of wavelets to all polynomials up to the power $`(M1)`$, defining its regularity and oscillatory behavior

$$𝑑xx^n\psi (x)=0,n=0,\mathrm{},(M1),$$
(32)
provides the relation

$$\underset{k}{}k^ng_k=0,$$
(33)
giving rise to

$$\underset{k}{}(1)^kk^nh_k=0,$$
(34)
when the formula (31) is taken into account.

The normalization condition

$$𝑑x\phi (x)=1$$
(35)
can be rewritten as another equation for $`h_k`$:

$$\underset{k}{}h_k=\sqrt{2}.$$
(36)
Let us write down the equations (28),(34), (36) for $`M=2`$ explicitly:

$`h_0h_2+h_1h_3=0,`$
$`h_0h_1+h_2h_3=0,`$
$`h_1+2h_23h_3=0,`$
$`h_0+h_1+h_2+h_3=\sqrt{2}.`$
The solution of this system is

$$h_3=\frac{1}{4\sqrt{2}}(1\pm \sqrt{3}),h_2=\frac{1}{2\sqrt{2}}+h_3,h_1=\frac{1}{\sqrt{2}}h_3,h_0=\frac{1}{2\sqrt{2}}h_3,$$
(37)
that, in the case of the minus sign for $`h_3`$, corresponds to the well known filter

$$h_0=\frac{1}{4\sqrt{2}}(1+\sqrt{3}),h_1=\frac{1}{4\sqrt{2}}(3+\sqrt{3}),h_2=\frac{1}{4\sqrt{2}}(3\sqrt{3}),h_3=\frac{1}{4\sqrt{2}}(1\sqrt{3}).$$
(38)
These coefficients define the simplest $`D^4`$ (or $`{}_{2}{}^{}\psi `$) wavelet from the famous family of orthonormal Daubechies wavelets with finite support. It is shown in the upper part of Fig. 5 by the dotted line with the corresponding scaling function shown by the solid line. Some other higher rank wavelets are also shown there. It is clear from this Figure (especially, for $`D^4`$) that wavelets are smoother in some points than in others. The choice of the plus sign in the expression for $`h_3`$ would not change the general shapes of the scaling function and wavelet $`D^4`$. It results in their mirror symmetrical forms obtained by a simple reversal of the signs on the horizontal and vertical axes, correspondingly. However, for higher rank wavelets different choices of signs would correspond to different forms of the wavelet . After the signs are chosen, it is clear that compactly supported wavelets are unique, for a given multiresolution analysis up to a shift in the argument (translation) which is inherently there. The dilation factor must be rational within the framework of the multiresolution analysis . Let us note that $`{}_{2}{}^{}\phi `$ is Hölder continuous with the global exponent $`\alpha =0.55`$ (see Eq. (83) below) and has different local Hölder exponents on some fractal sets. Typically, wavelets are more regular in some points than in others.

For the filters of higher order in $`M`$, i.e., for higher rank Daubechies wavelets, the coefficients can be obtained in an analogous manner. It is however necessary to solve the equation of the $`M`$-th power in this case. Therefore, the numerical values of the coefficients can be found only approximately, but with any predefined accuracy. The wavelet support is equal to $`2M1`$. It is wider than for the Haar wavelets. However the regularity properties are better. The higher order wavelets are smoother compared to $`D^4`$ as seen in Fig. 5. The Daubechies wavelet with $`r`$ vanishing moments has $`\mu r`$ continuous derivatives where $`\mu 0.2`$ as was estimated numerically. This means that $`80\%`$ of zero moments are ”wasted”. As the regularity $`r`$ increases, so does the support in general. For sufficiently regular functions, Daubechies wavelet coefficients are much smaller ($`2^{Mj}`$ times) than the Haar wavelet coefficients, i.e., the signal can be compressed much better with Daubechies wavelets. Since they are more regular, the synthesis is more efficient also.

One can ask the question whether the regularity or the number of vanishing moments is more important. The answer depends on the application, and is not always clear. It seems that the number of vanishing moments is more important for stronger compression which increases for larger number of vanishing moments, while regularity can become crucial in inverse synthesis to smooth the errors due to the compression (omission of small coefficients).

In principle, by solving the functional equation (11) one can find the form of the scaling function and, from (14), the form of the corresponding ”mother wavelet ”. There is no closed-form analytic formula for the compactly supported $`\phi (x),\psi (x)`$ (except for the Haar case). Nevertheless one can compute their plots, if they are continuous, with arbitrarily high precision using a fast cascade algorithm with the wavelet decomposition of $`\phi (x)`$ which is a special case of a refinement scheme (for more details, see ). Instead of a refinement cascade one can compute $`\phi (2^jk)`$ directly from Eq. (11) starting from appropriate $`\phi (n)`$. However, in practical calculations the above coefficients $`h_k`$ are used only without referring to the shapes of wavelets.

Except for the Haar basis, all real orthonormal wavelet bases with compact support are asymmetric, i.e., they have neither symmetry nor antisymmetry axis (see Fig. 5). The deviation of a wavelet from symmetry is judged by how much the phase of the expresion $`m_0(\omega )=_kh_ke^{ik\omega }`$ deviates from a linear function. The ”least asymmetric” wavelets are constructed by minimizing this phase. Better symmetry for a wavelet necessarily implies better symmetry for the coefficients $`h_k`$ but the converse statement is not always true.

## 5 Fast wavelet transform and coiflets

After calculation of the coefficients $`h_k`$ and $`g_k`$, i.e., the choice of a definite wavelet , one is able to perform the wavelet analysis of a signal $`f(x)`$ because the wavelet orthonormal basis ($`\psi _{j,k};\phi _{j,k}`$) has been defined. Any function $`fL^2(R)`$ is completely characterized by the coefficients of its decomposition in this basis and may be decomposed according to the formula (22). Let us make the sum limits in this formula more precise. The function $`f(x)`$ may be considered at any $`n`$-th resolution level $`j_n`$. Then the separation of its average values and fluctuations at this level looks like

$$f(x)=\underset{k=\mathrm{}}{\overset{\mathrm{}}{}}s_{j_n,k}\phi _{j_n,k}(x)+\underset{j=j_n}{\overset{\mathrm{}}{}}\underset{k=\mathrm{}}{\overset{\mathrm{}}{}}d_{j,k}\psi _{j,k}(x).$$
(39)
At the infinite interval, the first sum may be omitted as explained above, and one gets the pure wavelet expansion. As we stressed already, the coefficients $`s_{j,k}`$ and $`d_{j,k}`$ carry information about the content of the signal at various scales. They can be calculated directly using the formulas (24), (25). However this algorithm is inconvenient for numerical computations because it requires many ($`N^2`$) operations where $`N`$ denotes a number of the sampled values of the function. We will describe a faster algorithm. It is clear from Fig. 6, and the fast algorithm’ formulas are presented below.

In real situations with digitized signals, we have to deal with finite sets of points. Thus, there always exists the finest level of resolution where each interval contains only a single number. Correspondingly, the sums over $`k`$ will get the finite limits. It is convenient to reverse the level indexation assuming that the label of this fine scale is $`j=0`$. It is then easy to compute the wavelet coefficients for more sparse resolutions $`j1`$.

The multiresolution analysis naturally leads to an hierarchical and fast scheme for the computation of the wavelet coefficients of a given function. The functional equations (11), (14) and the formulas for the wavelet coefficients (24), (25) give rise, in case of Haar wavelets, to the relations (2), or for the backward transform (synthesis) to (3).

In general, one can get the iterative formulas of the fast wavelet transform

$$s_{j+1,k}=\underset{m}{}h_ms_{j,2k+m},$$
(40)
$$d_{j+1,k}=\underset{m}{}g_ms_{j,2k+m}$$
(41)
where

$$s_{0,k}=𝑑xf(x)\phi (xk).$$
(42)
These equations yield fast algorithms (the so-called pyramid algorithms) for computing the wavelet coefficients, asking now just for $`O(N)`$ operations to be done. Starting from $`s_{0,k}`$, one computes all other coefficients provided the coefficients $`h_m,g_m`$ are known. The explicit shape of the wavelet is not used in this case any more. The simple form of these equations is the only justification for introducing the factor $`\sqrt{2}`$ in the functional equation (11). In principle, the coefficients $`h_m,g_m`$ could be renormalized. However, in practice the Eqs. (40) and (41) are used much more often than others, and this normalization is kept intact. After choosing a particular wavelet for analysis ,i.e., choosing $`h_m,g_m`$, one uses only the Eqs. (40) and (41) for computing the wavelet coefficients, and additional factors in these equations would somewhat complicate the numerical processing.

The remaining problem lies in the initial data. If an explicit expression for $`f(x)`$ is available, the coefficients $`s_{0,k}`$ may be evaluated directly according to (42). But this is not so in the situation when only discrete values are available. To get good accuracy, one has to choose very small bins (dense lattice) which is often not accessible with finite steps of sampling. In such a case, the usually adopted simplest solution consists in directly using the values $`f(k)`$ of the available sample in place of the coefficients $`s_{0,k}`$ and start the fast wavelet transform using formulas (40), (41). This is a safe operation since the pyramid algorithm yields perfect reconstruction, and the coefficient $`s_{0,k}`$ essentially represents a local average of the signal provided by the scaling function.

In general, one can choose

$$s_{0,k}=\underset{m}{}c_mf(km).$$
(43)
The above supposition $`s_{0,k}=f(k)`$ corresponds to $`c_m=\delta _{0m}`$. This supposition may be almost rigorous for some specific choices of scaling functions named coiflets (by the name of R. Coifman whose ideas inspired I. Daubechies to build up these wavelets). It is possible to construct multiresolution analysis with the scaling function having vanishing moments, i.e., such that

$$𝑑xx^m\phi (x)=0;\mathrm{\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}0}<m<M.$$
(44)
To construct such wavelets (coiflets), one has to add to the equations for determining the coefficients $`h_k`$ a new condition

$$\underset{k}{}h_kk^m=0,\mathrm{\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}0}<m<M,$$
(45)
which follows from the requirement (44).

Coiflets are more symmetrical than Daubechies wavelets as is seen from Fig. 7 if compared to Fig. 5. The latter do not have the property (44). The price for this extra generalization is that coiflets are longer than Daubechies wavelets. If in the latter case the length of the support is $`2M1`$, for coiflets it is equal to $`3M1`$. The error in the estimation of $`s_{j,k}`$ decreases with the number of vanishing moments as $`O(2^{jM})`$.

There are other proposals to improve the first step of the iterative procedure promoted by using $`T`$-polynomials by W. Sweldens  or the so-called ”lazy” or interpolating wavelets by S. Goedecker and O. Ivanov . The latter is most convenient for simultaneous analysis at different resolution levels, in particular for non-equidistant lattices.

Inverse fast wavelet transform allows for the reconstruction of the function starting from its wavelet coefficients.

## 6 Choice of wavelets

Above, we demonstrated three examples of discrete orthonormal compactly supported wavelets. The regularity property, the number of vanishing moments and the number of wavelet coefficients exceeding some threshold value were considered as possible criteria for the choice of a particular wavelet not to say about computing facilities. Sometimes the so-called information cost functional used by statisticians is introduced, and one tries to minimize it and thus select the optimal basis. In particular, the entropy criterium for the probability distribution of the wavelet coefficients is also considered . The entropy of $`f`$ relative to the wavelet basis measures the number of significant terms in the decomposition (22). It is defined by $`\mathrm{exp}(_{j,k}|d_{j,k}|^2\mathrm{log}|d_{j,k}|^2)`$. If we have a collection of orthonormal bases, we will choose for the analysis of $`f`$ the particular basis that yields the minimum entropy.

The number of possible wavelets at our disposal is much larger than the above examples show. We will not discuss all of them just mentioning some and referring the reader to the cited books.

* First, let us mention splines which lead to wavelets with non-compact support but with the exponential decay at infinity and with some (limited) number of continuous derivatives. The special orthogonalization tricks should be used there. The splines are intrinsically associated with interpolation schemes for finding more precise initial values of $`s_{0,k}`$ relating them to some linear combinations of the sampled values of $`f(x)`$.
* To insure both the full symmetry and exact reconstruction, one has to use the so-called biorthogonal wavelets. It means that two dual wavelet bases $`\psi _{j,k}`$ and $`\psi ^{j,k}`$, associated with two different multiresolution ladders, are used. They can have very different regularity properties. A function $`f`$ may be represented in two forms, absolutely equivalent until the compression is done:

$$f=\underset{j,k}{}f,\psi _{j,k}\psi ^{j,k}$$
(46)
$$=\underset{j,k}{}f,\psi ^{j,k}\psi _{j,k},$$
(47)
where $`\psi `$ and its dual wavelet satisfy the biorthogonality requirement $`\psi _{j,k}|\psi ^{j^{},k^{}}=\delta _{j,k;j^{},k^{}}`$. In distinction to Daubechies wavelets, where regularity is tightly connected with the number of vanishing moments, biorthogonal wavelets have much more freedom. If one of them has the regularity of the order $`r`$, then its dual partner wavelet has automatically at least $`r`$ vanishing moments. If $`\psi ^{j,k}`$ is much more regular than $`\psi _{j,k}`$, then $`\psi _{j,k}`$ has many more vanishing moments than $`\psi ^{j,k}`$. It allows us to choose, e.g., very regular $`\psi ^{j,k}`$ and get many vanishing moments of $`\psi _{j,k}`$. The large number of vanishing moments of $`\psi _{j,k}`$ leads to better compressibility for reasonably smooth $`f`$. If the compression has been done, the formula (46) is much more useful than (47). The number of significant terms is much smaller, and, moreover, better regularity of $`\psi ^{j,k}`$ helps reconstruct $`f`$ more precisely. Biorthogonal bases are close to an orthonormal basis. Both wavelets can be made symmetric. Symmetric biorthogonal wavelets close to orthonormal basis are close to coiflets. The construction of biorthogonal wavelet bases is simpler than of orthonormal bases.
* The existence of two-scale relations is the main feature of the construction of wavelet packets. The general idea of the wavelet packets is to iterate further the splitting of the frequency band, still keeping the same pair of filters. The scaling function introduced above acquires the name $`w_0`$, and the packet is built up through the following iterations

$$w_{2n}(x)=\underset{k}{}h_kw_n(2xk),$$
(48)
$$w_{2n+1}(x)=\underset{k}{}g_kw_n(2xk).$$
(49)
The usual mother wavelet is represented by $`w_1`$. This family of wavelets forms an orthonormal basis in $`L^2(R)`$ which is called the fixed scale wavelet packet basis. Fig. 8 demonstrates the whole construction.
* One can abandon the orthonormality property and construct non-orthogonal wavelets called frames. An important special class of frames is given by the Riesz bases of $`L^2(R)`$. A Riesz basis is a frame, but the converse is not true in general. The frames satisfy the following requirement:

$$Af^2\underset{jJ}{}|f,\phi _j|^2Bf^2.$$
(50)
The constants $`A`$ and $`B`$ are called the frame bounds. For $`A=B`$ one calls them as tight frames. The case $`A=B=1`$ corresponds to orthonormal wavelets.
* When acting by (singular) operators, one sometimes gets infinities when usual wavelets are considered. Some suitable function $`b(x)`$ may be used to specify the extra conditions which will be necessary and sufficient for the result of the (singular integral) operator action to be continuous in $`L^2`$. In this case one chooses the so-called ”wavelets which are adapted to $`b`$”. Any function $`f`$ is again decomposed as

$$f(x)=\underset{\lambda }{}a(\lambda )\psi _\lambda ^{(b)}(x)$$
(51)
but the wavelet coefficients are now calculated according to

$$a_\lambda =𝑑xb(x)f(x)\psi _\lambda ^{(b)}(x).$$
(52)
They satisfy the normalization condition

$$𝑑xb(x)\psi _\lambda ^{(b)}(x)\psi _\lambda ^{}^{(b)}(x)=\delta _{\lambda ,\lambda ^{}}.$$
(53)
The cancellation condition now reads

$$𝑑xb(x)\psi _\lambda ^{(b)}(x)=0.$$
(54)
As we see, the cancellation is also adapted to the function $`b`$ (in general, to the ”complex measure” $`dxb(x)`$).
* Up to now we considered the wavelets with the dilation factor equal 2. It is most convenient for numerical calculations. However, it can be proved  that within the framework of a multiresolution analysis , the dilation factor must be rational, and no other special requirements are imposed. Therefore one can construct schemes with other integer or fractional dilation factors. Sometimes they may provide a sharper frequency localization. For wavelets with the dilation factor 2, their Fourier transform is essentially localized within one octave between $`\pi `$ and $`2\pi `$, whereas the fractional dilation wavelet bases may have a bandwidth narrower than one octave.
* Moreover, one can use the continuous wavelets as described at some length in Ref. . The wavelet is written in a form

$$\psi _{a,b}(x)=|a|^{1/2}\psi (\frac{xb}{a}).$$
(55)
The direct and inverse formulas of the wavelet transform look like

$$W_{a,b}=|a|^{1/2}𝑑xf(x)\psi (\frac{xb}{a}),$$
(56)
$$f(x)=C_\psi ^1\frac{dadb}{a^2}W_{a,b}\psi _{a,b}(x).$$
(57)
Here

$$C_\psi =\frac{d\omega }{|\omega |}|\psi (\omega )|^2=\frac{d\omega }{|\omega |}|𝑑xe^{ix\omega }\psi (x)|^2.$$
(58)
Herefrom one easily recognizes that the oscillation of wavelets required by Eq. (21) is a general property. The vanishing Fourier transform of a wavelet at $`\omega 0`$, which is just directly the condition (21), provides a finite value of $`C_\psi `$ in (58). One of the special and often used examples of continuous wavelets is given by the second derivative of the Gaussian function which is called the Mexican hat (MHAT) wavelet for its shape. Actually, it can be considered as a special frame as shown by Daubechies. The reconstruction procedure (synthesis) is complicated and can become unstable in this case. However it is widely applied for analysis of signals. The formula (56) is a kind of a convolution operation. That is why the general theory of the so-called Calderon-Zygmund operators  (see Appendix 2) is, in particular, applicable to problems of the wavelet decomposition.

## 7 Multidimensional wavelets

The multiresolution analysis can be performed in more than 1 dimensions. There are two ways  to generalize it to the two-dimensional case, for example, but we will consider the most often used construction given by tensor products. The tensor product method is a direct way to construct an $`r`$-regular multiresolution approximation which produces multidimensional wavelets of compact support. This enables us to analyze every space of functions or distributions in $`n`$ dimensions whose regularity is bounded by $`r`$.

The trivial way of constructing a two-dimensional orthonormal basis starting from a one-dimensional orthonormal wavelet basis $`\psi _{j,k}(x)=2^{j/2}\psi (2^jxk)`$ is simply to take the tensor product functions generated by two one-dimensional bases:

$$\mathrm{\Psi }_{j_1,k_1;j_2,k_2}(x_1,x_2)=\psi _{j_1,k_1}(x_1)\psi _{j_2,k_2}(x_2).$$
(59)
In this basis the two variables $`x_1`$ and $`x_2`$ are dilated independently.

More interesting for many applications is another construction, in which dilations of the resulting orthonormal wavelet basis control both variables simultaneously, and the two-dimensional wavelets are given by the following expression:

$$2^j\mathrm{\Psi }(2^jxk,2^jyl),j,k,lZ,$$
(60)
where $`\mathrm{\Psi }`$ is no longer a single function: on the contrary, it consists of three elementary wavelets. To get an orthonormal basis of $`W_0`$ one has to use in this case three families $`\phi (xk)\psi (yl),\psi (xk)\phi (yl),\psi (xk)\psi (yl)`$. Then the two-dimensional wavelets are $`2^j\phi (2^jxk)\psi (2^jyl),\mathrm{\hspace{0.33em}2}^j\psi (2^jxk)\phi (2^jyl),\mathrm{\hspace{0.33em}2}^j\psi (2^jxk)\psi (2^jyl)`$. In the two-dimensional plane, the analysis is done along the horizontal, vertical and diagonal strips with the same resolution in accordance with these three wavelets.

Fig. 9 shows how this construction looks like. The schematic representation of this procedure in the left-hand side of the Figure demonstrates how the corresponding wavelet coefficients are distributed for different resolution levels ($`j=1`$ and 2). The superscripts $`d,v`$ and $`h`$ correspond to diagonal, vertical and horizontal coefficients. In the right-hand side of the Figure, a set of geometrical objects is decomposed into two layers. One clearly sees how vertical, horizontal and diagonal edges are emphasized in the corresponding regions. One should notice also that the horizontal strip is resolved into two strips at a definite resolution level.

In the general $`n`$-dimensional case, there exist $`2^n1`$ functions which form an orthonormal basis and admit the multiresolution analysis of any function from $`L^2(R^n)`$. The normalization factor is equal to $`2^{nj/2}`$ in front of the function, as can be guessed already from the above two-dimensional case with this factor equal to $`2^j`$ in contrast to $`2^{j/2}`$ in one dimension.

There exists also a method to form the wavelet bases which are not reducible to tensor products of one-dimensional wavelets (see ). In dimension 1, every orthonormal basis arises from a multiresolution approximation. In dimensions greater than 1, it is possible to form an orthonormal basis such that there is no $`r`$-regular multiresolution approximation ($`r1`$) from which these wavelets can be obtained .

## 8 The Fourier and wavelet transforms

In many applications (especially, for non-stationary signals), one is interested in the frequency content of a signal locally in time, i.e., tries to learn which frequencies are important at some particular moment. As has been stressed already, the wavelet transform is superior to the Fourier transform, first of all, due to the locality property of wavelets. The Fourier transform uses sine, cosine or imaginary exponential functions as the main basis. It is spread over the entire real axis whereas wavelet basis is localized. It helps analyze the local properties of a signal using wavelets while the Fourier transform does not provide any information about the location where the scale (frequency) of a signal changes. Decomposition into wavelets allows singularities to be located by observing the places where the wavelet coefficients are (abnormally) large. Obviously, nothing of the kind happens for the Fourier transform. Once the wavelets have been constructed they perform incredibly well in situations where Fourier series and integrals involve subtle mathematics or heavy numerical calculations. But wavelet analysis cannot entirely replace Fourier analysis , indeed, the latter is often used in constructing the orthonormal bases of wavelets needed for analysis with wavelet series. Many theorems of wavelet analysis are proven with the help of the Fourier decomposition. The two kinds of analysis are thus complementary rather than always competing.

The Fourier spectrum $`f_\omega `$ of a one-dimensional signal $`f(t)`$ having finite energy (i.e., square-integrable) is given by

$$f_\omega =_{\mathrm{}}^{\mathrm{}}f(t)e^{i\omega t}𝑑t.$$
(61)
The inverse transform restores the signal

$$f(t)=\frac{1}{2\pi }_{\mathrm{}}^{\mathrm{}}f_\omega e^{i\omega t}𝑑\omega .$$
(62)
It is an unitary transformation

$$|f(t)|^2𝑑t=\frac{1}{2\pi }|f_\omega |^2𝑑\omega .$$
(63)
This is the so-called Parseval identity which states the conservation of energy between the time and the frequency domains. The formula (61) asks for information about the signal $`f(t)`$ from both past and future times. It is suited for a stationary signal when the frequency $`\omega `$ does not depend on time. Thus its time-frequency band is well located in frequency and practically unlimited in time, i.e., Eq. (61) gives a representation of the frequency content of $`f`$ but not of its time-localization properties. Moreover, the signal $`f(t)`$ should decrease fast enough at infinities the integral (61) to be meaningful.

The attempt to overcome these difficulties and improve time-localization while still using the same basis functions is done by the so-called windowed Fourier transform. The signal $`f(t)`$ is considered within some time interval (window) only. In practice, one has to restrict the search for the optimal window to the easily generated windows. A simple way of doing so amounts to multiplying $`f(t)`$ by a simplest compactly supported window function, e.g., by $`g_w=\theta (tt_i)\theta (t_ft)`$, where $`\theta `$ is the commonly used Heaviside step-function different from zero at positive values of its argument only, $`t_i,t_f`$ are the initial and final cut-offs of the signal (more complicated square-integrable window functions $`g`$, well concentrated in time, e.g., a Gaussian or canonical coherent states777In quantum mechanics, they are introduced for quantizing the classical harmonic oscillator. In signal analysis , they are known by the name of Gabor functions., can be used as well). Thus one gets

$$f_{\omega ,w}=_{\mathrm{}}^{\mathrm{}}f(t)g(t)e^{i\omega t}𝑑t.$$
(64)
In its discrete form it can be rewritten as

$$f_{\omega ,w}=f(t)g(tnt_0)e^{im\omega _0t}𝑑t,$$
(65)
where $`\omega _0,t_0>0`$ are fixed, and $`m,n`$ are some numbers which define the scale and location properties. The time localization of the transform is limited now but the actual window is fixed by different functions in time and frequency which do not depend on the resolution scale and have fixed widths. Moreover, the orthonormal basis of the windowed Fourier transform can be constructed only for the so-called Nyquist density (corresponding to the condition $`\omega _0t_0=2\pi `$; see Section 14), whereas there is no such restriction imposed on wavelets. At this critical value, frames are possible, but the corresponding functions are either poorly localized, or have poor regularity. This result is known as the Balian–Low phenomenon. In practical applications of the windowed Fourier transform, to achieve better localization one should choose $`\omega _0t_0<2\pi `$ thus destroying orthonormality.

The difference between the wavelet and windowed Fourier transforms lies in the shapes of the analyzing functions $`\psi `$ and $`g`$. All $`g`$, regardless of the value of $`\omega `$, have the same width. In contrast, the wavelets $`\psi `$ automatically provide the time (location) resolution window adapted to the problem studied, i.e., to its essential frequencies (scales). Namely, let $`t_0,\delta `$ and $`\omega _0,\delta _\omega `$ be the centers and the effective widths of the wavelet basic function $`\psi (t)`$ and its Fourier transform. Then for the wavelet family $`\psi _{j,k}(t)`$ (20) and, correspondingly, for wavelet coefficients, the center and the width of the window along the $`t`$-axis are given by $`2^j(t_0+k)`$ and $`2^j\delta `$. Along the $`\omega `$-axis they are equal to $`2^j\omega _0`$ and $`2^j\delta _\omega `$. Thus the ratios of widths to the center position along each axis do not depend on the scale. It means that the wavelet window resolves both the location and the frequency in fixed proportions to their central values. For the high-frequency component of the signal it leads to quite large frequency extension of the window whereas the time location interval is squeezed so that the Heisenberg uncertainty relation is not violated. That is why wavelet windows can be called as Heisenberg windows. Correspondingly, the low-frequency signals do not require small time intervals and admit wide window extension along the time axis. Thus wavelets localize well the low-frequency ”details” on the frequency axis and the high-frequency ones on the time axis. This ability of wavelets to find a perfect compromise between the time localization and the frequency localization by choosing automatically the widths of the windows along the time and frequency axes well adjusted to their centers location is crucial for their success in signal analysis . The wavelet transform cuts up the signal (functions, operators etc) into different frequency components, and then studies each component with a resolution matched to its scale providing a good tool for time-frquency (position-scale) localization. That is why wavelets can zoom in on singularities or transients (an extreme version of very short-lived high-frequency features!) in signals, whereas the windowed Fourier functions cannot. In terms of the traditional signal analysis , the filters associated with the windowed Fourier transform are constant bandwidth filters whereas the wavelets may be seen as constant relative bandwidth filters whose widths in both variables linearly depend on their positions.

In Fig. 10 we show the difference between these two approaches. It demonstrates the constant shape of the windowed Fourier transform region and the varying shape (with a constant area) of the wavelet transform region. The density of localization centers is homogeneous for the windowed Fourier transform whereas it changes for the wavelet transform so that at low frequencies the centers are far apart in time and become much denser for high frequencies.

From the mathematical point of view, it is important that orthonormal wavelets give good unconditional888For unconditional bases, the order in which the basis vectors are taken does not matter. All of the known constructions of unconditional bases of wavelets rely on the concept of the multiresolution analysis .  bases for other spaces than that of the square integrable functions, out-performing the Fourier basis functions in this respect. It is applied in the subsequent Sections to a characterization of such functions using only the absolute values of wavelet coefficients. In other words, by looking only at the absolute values of wavelet coefficients we can determine to which space this function belongs. This set of spaces is much wider than in the case of the Fourier transform by which only Sobolev spaces999 The function $`f`$ belongs to the Sobolev space $`W^s(R)`$ if its Fourier transform provides the finite integrals $`𝑑\omega (1+|\omega |^2)^s|f(\omega )|^2`$. can be completely characterized.

As we mentioned already, the wavelet analysis concentrates near the singularities of the function analyzed. The corresponding wavelet coefficients are negligible in the regions where the function is smooth. That is why the wavelet series with a plenty of non-zero coefficients represent really pathological functions, whereas ”normal” functions have ”sparse” or ”lacunary” wavelet series and easy to compress. On the other hand, Fourier series of the usual functions have a lot of non-zero coefficients, whereas ”lacunary” Fourier series represent pathological functions. Let us note at the end that, nevertheless, the Fourier transform is systematically used in the proof of many theorems in the theory of wavelets. It is not at all surprising because they constitute the stationary signals by themselves.

## 9 Wavelets and operators

The study of many operators acting on a space of functions or distributions becomes simple when suitable wavelets are used because these operators can be approximately diagonalized with respect to this basis. Orthonormal wavelet bases provide the unique example of a basis with non-trivial diagonal, or almost-diagonal, operators. The operator action on the wavelet series representing some function does not have uncontrollable sequences, i.e., wavelet decompositions are robust. One can describe precisely what happens to the initial series under the operator action and how it is transformed. In a certain sense, wavelets are stable under the operations of integration and differentiation. That is why wavelets, used as a basis set, allow us to solve differential equations characterized by widely different length scales found in many areas of physics and chemistry. Moreover, wavelets reappear as eigenfunctions of certain operators.

To deal with operators in the wavelet basis it is convenient, as usual, to use their matrix representation. For a given operator $`T`$ it is represented by the set of its matrix elements in the wavelet basis:

$$T_{j,k;j,^{}k^{}}=<\psi _{j,k}^{}|T|\psi _{j^{},k^{}}>.$$
(66)
For linear homogeneous operators, their matrix representation can be explicitly calculated .

It is extremely important that it is sufficient to first calculate the matrix elements at some ($`j`$-th) resolution level. All other matrix elements can be obtained from it using the standard recurrence relations. Let us derive the explicit matrix elements $`r_{j,l;j,l^{}}`$ of the homogeneous operator $`T`$ of the order $`\alpha `$:

$$r_{j,l;j,l^{}}=<\phi _{j,l}|T|\phi _{j,l^{}}>.$$
(67)
Using the recurrence relations between the scaling functions at the given and finer resolution levels, one gets the following equation relating the matrix elements at neighboring levels:

$$r_{j,l;j,l^{}}=<(\underset{k}{}h_k\phi _{j+1,2lk})|T|(\underset{k^{}}{}h_k^{}\phi _{j+1,2l^{}k^{}})>=\underset{k}{}\underset{k^{}}{}h_kh_k^{}r_{j+1,2lk;j+1,2l^{}k^{}}.$$
(68)
For the operator $`T`$ having the homogeneity index $`\alpha `$ one obtains

$$r_{j,l;j,l^{}}=2^\alpha \underset{k}{}\underset{k^{}}{}h_kh_k^{}r_{j,2lk;j,2l^{}k^{}}.$$
(69)
The solution of this equation defines the required coefficients up to the normalization constant which can be easily obtained from the results of the action by the operator $`T`$ on a polynomial of a definite rank. For non-integer values of $`\alpha `$, this is an infinite set of equations.

The explicit equation for the $`n`$-th order differentiation operator is

$`r_k^{(n)}=\phi (x)|{\displaystyle \frac{d^n}{dx^n}}|\phi (xk)=`$
$`{\displaystyle \underset{i,m}{}}h_ih_m\phi (2x+i)|{\displaystyle \frac{d^n}{dx^n}}|\phi (2x+mk)=`$
$`2^n{\displaystyle \underset{i,m}{}}h_ih_mr_{2kim}^{(n)}.`$ (70)
It leads to a finite system of linear equations for $`r_k`$ (the index $`n`$ is omitted):

$$2r_k=r_{2k}+\underset{m}{}a_{2m+1}(r_{2k2m+1}+(1)^nr_{2k+2m1}),$$
(71)
where both $`r_k`$ and $`a_m=_ih_ih_{i+m}`$ ($`a_0=1`$) are rational numbers in the case of Daubechies wavelets. The wavelet coefficients can be found from these equations up to a normalization constant. The normalization condition reads :

$$\underset{k}{}k^nr_k=n!.$$
(72)
For the support region of the length $`L`$, the coefficients $`r_k`$ differ from zero for $`L+2kL2`$, and the solution exists for $`Ln+1`$. These coefficients possess the following symmetry properties:

$$r_k=r_k$$
(73)
for even $`k`$, and

$$r_k=r_k$$
(74)
for odd values of $`k`$.

Here, as an example, we show in the Table the matrix elements of the first and second order differential operators in the Daubechies wavelet basis with four vanishing moments ($`D^8`$).

Table.

| k | $`h_k`$ | $`\phi (x)||\phi (xk)`$ | $`\phi (x)|^2|\phi (xk)`$ |
| --- | --- | --- | --- |
| -6 | 0 | 0.00000084 | 0.00001592 |
| -5 | 0 | -0.00017220 | -0.00163037 |
| -4 | 0 | -0.00222404 | -0.01057272 |
| -3 | 0 | 0.03358020 | 0.15097289 |
| -2 | -0.07576571 | -0.19199897 | -0.69786910 |
| -1 | -0.02963552 | 0.79300950 | 2.64207020 |
| 0 | 0.49761866 | 0 | -4.16597364 |
| 1 | 0.80373875 | -0.79300950 | 2.64207020 |
| 2 | 0.29785779 | 0.19199897 | -0.69786910 |
| 3 | -0.09921954 | -0.03358020 | 0.15097289 |
| 4 | -0.01260396 | 0.00222404 | -0.01057272 |
| 5 | 0.03222310 | 0.00017220 | -0.00163037 |
| 6 | 0 | -0.00000084 | 0.00001592 |

For a continuous linear operator $`T`$ represented by a singular integral

$$Tf(x)=K(x,y)f(y)𝑑y$$
(75)
with some definite conditions imposed on the kernel $`K`$ (see Appendix 2), there exists the important theorem (called $`T(1)`$ theorem) which states a necessary and sufficient condition for the extension of $`T`$ as a continuous linear operator on $`L^2(R^n)`$ (more details and the elegant wavelet proof of the theorem are given in Ref. ).

Let us note that a standard difficulty for any spectral method is a representation of the operator of multiplication by a function. As an example for physicists, we would mention the operator of the potential energy in the Schrödinger equation. However, it is well known that this operation is trivial and diagonal in the real space. Therefore for such operations one should deal with a real space, and after all operations done there return back. Such algorithm would need $`O(N)`$ operations only.

## 10 Nonstandard matrix multiplication

There are two possible ways to apply operators onto functions within wavelet theory. They are called the standard and nonstandard matrix forms.

For smooth enough functions most wavelet coefficients are rather small. For a wide class of operators, most of their matrix elements are also quite small. Let us consider the structure of the elements of the matrix representation of some operator $`T`$ that are large enough. The matrix elements satisfy the following relations

$$T_{j,k;j^{},k^{}}0,where|kk^{}|\mathrm{},$$
(76)
$$T_{j,k;j^{},k^{}}0,where|jj^{}|\mathrm{}.$$
(77)
The topology of the distribution of these matrix elements within the matrix can be rather complicated. The goal of the nonstandard form is to replace the latter equation (77) by another, more rigorous one:

$$T_{j,k;j^{},k^{}}=0,ifjj^{}.$$
(78)
It avoids taking the matrix elements between the different resolution levels. To deal with it, one should consider, instead of the wavelet space, the overfull space with the basis containing both wavelets and scaling functions at various resolution levels.

Let us consider the action of the operator $`T`$ on the function $`f`$ which transforms it into the function $`g`$:

$$g=Tf.$$
(79)
Both $`g`$ and $`f`$ have the wavelet representations with the wavelet coefficients $`(_fs_{j,k};_fd_{j,k})`$ and $`(_gs_{j,k};_gd_{j,k})`$. At the finest resolution level $`j_n`$ only $`s`$-coefficients differ from zero, and the transformation looks like

$${}_{g}{}^{}s_{j_n,k}^{}=\underset{k^{}}{}T_{SS}(j_n,k;j_n,k^{})_fs_{j_n,k^{}}.$$
(80)
At the next level one gets in both the standard and nonstandard approaches

$${}_{g}{}^{}s_{j_n1,k}^{}=\underset{k^{}}{}T_{SS}(j_n1,k;j_n1,k^{})_fs_{j_n1,k^{}}+\underset{k^{}}{}T_{SD}(j_n1,k;j_n1,k^{})_fd_{j_n1,k^{}},$$
(81)
$${}_{g}{}^{}d_{j_n1,k}^{}=\underset{k^{}}{}T_{DS}(j_n1,k;j_n1,k^{})_fs_{j_n1,k^{}}+\underset{k^{}}{}T_{DD}(j_n1,k;j_n1,k^{})_fd_{j_n1,k^{}},$$
(82)
where $`T_{SS}(j_n,k;j_n,k^{})=𝑑x\overline{\phi }_{j_n,k}(x)T\phi _{j_n,k^{}}(x)`$ and the replacement of subscripts $`SD`$ corresponds to the substitute $`\phi \psi `$ in the integrals. There is coupling between all resolution levels because all $`s`$-coefficients at this $`(j_n1)`$-th level should be decomposed by the fast wavelet transform into $`s`$ and $`d`$-coefficients at higher levels. Therefore even for the almost diagonal initial step the standard matrix acquires rather complicated form as demonstrated in Fig. 11 for 4-level operations (similar to those discussed above in the case of the Haar wavelets). It is thus inefficient for numerical purposes.

As we see in Fig. 11, at the final stage of the standard approach we have to deal with the wavelet representation corresponding to the formula (22) with only one $`s`$-coefficient left in the vectors which represents the overall weighted average of the functions (the $`SS`$-transition from $`f`$ to $`g`$ is given by the upper left box in the Figure). At the same time, in the process of approaching it from the scaling-function representation (23), (1) we had to deal with average values at intermediate levels decomposing them at each step into $`s`$ and $`d`$ parts of further levels. These intermediate $`s`$-coefficients have been omitted since they were replaced by $`s`$ and $`d`$-coefficients at the next levels. That is why the matrix of the standard approach looks so complicated.

To simplify the form of the matrix, it was proposed  to use the redundant set of the wavelet coefficients. Now, let us keep these average values in the form of corresponding $`s`$-coefficients of the intermediate levels both in initial and final vectors representing functions $`f`$ and $`g`$. Surely, we will deal with redundant vectors which are larger than necessary for the final answer. However, it should not bother us because we know the algorithm to reduce the redundant vector to a non-redundant form. At the same time, this trick simplifies both the form of the transformation matrix and the computation. This non-standard form is shown in Fig. 12. Different levels have been completely decoupled because there are no blocks in this matrix which would couple them. The block with $`SS`$-elements has been separated, and in place of it the zero matrix is inserted. The whole matrix becomes artificially enlarged. Correspondingly, the vectors, characterizing the functions $`f`$ and $`g`$, are also enlarged. Here all intermediate $`s`$-coefficients are kept for the function $`f`$ (compare the vectors in the right-hand sides of Figs. 11 and 12). Each $`S_{j+1}`$ is generated from $`S_j`$ and $`D_j`$. That is where the coupling of different levels is still present. In the transformation matrix all $`SS`$-elements are zero except for the lowest one $`S_0S_0`$. All other $`SD,DS,DD`$ matrices are almost diagonal due to the finite support of scaling functions and wavelets. The redundant form of the $`g`$-function vector of Fig. 12 may be reduced to its usual wavelet representation of Fig. 11 by splitting up any $`S_j`$ into $`S_{j1}`$ and $`D_{j1}`$ by a standard method. Then these $`S_{j1}`$ and $`D_{j1}`$ are added to the corresponding places in the vector. This iterative procedure allows, by going from $`S_{j1}`$ down to $`S_0`$, to get the usual wavelet series of the function $`g`$. We get rid of all $`s`$-coefficients apart from $`S_0`$. The computation becomes fast.

## 11 Regularity and differentiability

The analysis of any signal includes finding the regions of its regular and singular behavior. One of the main features of wavelet analysis is its capacity of doing a very precise analysis of the regularity properties of functions. When representing a signal by a wavelet transform, one would like to know if and at which conditions the corresponding series is convergent and, therefore, where the signal is described by the differentiable function or where the singularities appear. For certain singular functions, essential information is carried by a limited number of their wavelet coefficients. Such an information may be used in the design of numerical algorithms. We start with the traditional Hölder conditions of regularity and, in the next Section, proceed to their generalization following Ref. .

The Hölder definition of the pointwise regularity at a point $`x_0`$ of a real-valued function $`f`$ defined on $`R^n`$ declares that this function belongs to the space $`C^\alpha (x_0)`$ ($`\alpha >n`$) if there exists a polynomial $`P`$ of order at most $`\alpha `$ such that

$$|f(x)P(xx_0)|C|xx_0|^\alpha ,$$
(83)
where $`C`$ is a constant.

The supremum of all values of $`\alpha `$ such that (83) holds is called the Hölder exponent of $`f`$ at $`x_0`$ and denoted $`\alpha (x_0)`$101010In principle, the expression $`|xx_0|^\alpha `$ in the right-hand side of Eq. (83) can be replaced by a more general function satisfying definite conditions and called a modulus of continuity but we shall use only the above definition. In general, a modulus of continuity defines the largest deviation from that best polynomial approximation of a function $`f`$ which is characterized by the set of smallest deviations compared to other polynomial approximations..

Correspondingly, a point $`x_0`$ is called a strong $`\alpha `$-singularity $`(n<\alpha 1)`$ of $`f`$ if, at small intervals, the following inequality is valid

$$|f(x)f(y)|C|xy|^\alpha $$
(84)
for a relatively large set of $`x`$’s and $`y`$’s close to $`x_0`$.

Such definitions work quite well for a rather smooth function. However, in the case of a function drastically fluctuating from point to point they are difficult to handle at each point because, e.g., the derivative of $`f`$ may have no regularity at $`x_0`$, the condition (83) may become unstable under some operators such as Hilbert transform etc. The generalized definition of pointwise regularity will be given in the next Section by introducing two-microlocal spaces. Here, we concentrate on global properties of the function $`f`$.

The uniform regularity of a function $`f`$ at positive non-integer values of $`\alpha `$ consists in the requirement for the Eq. (83) to hold for all real $`n`$-dimensional $`x_0`$ with the uniform value of the constant $`C`$. At first sight, this condition looks rather trivial because for smooth functions the uniform regularity coincides with the pointwise regularity, which is everywhere the same. To understand that this is non-trivial, one could consider, e.g., the function $`f(x)=x\mathrm{sin}1/x`$ with pointwise exponents $`\alpha (0)=1,\alpha (x_0)=\mathrm{}`$ for any fixed $`x_0`$ different from 0, and uniform Hölder exponent $`\alpha =1/2`$ which appears for the set of points $`x_0|xx_0|^{1/2}|xx_0|`$.

The requirement of a definite uniform regularity of a $`n`$-dimensional function $`f`$ on the whole real axis at positive non-integer values of $`\alpha `$ can be stated in terms of wavelet coefficients as an inequality

$$|d_{j,k}|C2^{(\frac{n}{2}+\alpha )j}.$$
(85)
Thus from the scale behavior of wavelet coefficients one gets some characteristics of uniform regularity of the function. In particular, the linear dependence of the logarithms of wavelet coefficients on the scale $`j`$ would indicate on the scaling properties of a signal, i.e., on the fractal behavior whose parameters are determined from higher moments of wavelet coefficients (see below). For pointwise regularity determining the local characteristics, the similar condition is discussed in the next Section.

Now, let us formulate the conditions under which wavelet series converge at some points, or are differentiable. It has been proven that

* if $`f`$ is square integrable, the wavelet series of $`f`$ converges almost everywhere;
* the wavelet series of $`f`$ is convergent on a set of points equipotent to $`R`$ if

$$|d_{j,k}|=2^{j/2}\eta _j,$$
(86)
where $`\eta _j`$ tends to 0 when $`j`$ tends to $`+\mathrm{}`$;
* the function $`f`$ is almost everywhere differentiable and the derivative of the wavelet series of $`f`$ converges almost everywhere to $`f^{}`$, if the following condition is satisfied

$$|d_{j,k}|^22^{2j}\mathrm{};$$
(87)
* the function $`f`$ is differentiable on a set of points equipotent to $`R`$ and, at these points the derivative of its wavelet series converges to $`f^{}`$, if

$$|d_{j,k}|2^{3j/2}\eta _j,$$
(88)
where $`\eta _j`$ tends to 0 when $`j`$ tends to $`+\mathrm{}`$.

Limitations imposed on wavelet coefficients are generalized in the next Section to include the pointwise regularity condition when two-microlocal spaces are considered. These better estimates do not hold just uniformly, but hold locally (possibly, apart from a set of points of small Hausdorff dimension which can be determined by wavelet analysis ).

Note finally that one can also derive global regularity of $`f`$ from the decay in $`\omega `$ of the absolute value of its windowed Fourier transform, if the window function $`g`$ is chosen to be sufficiently smooth. In most cases, however, the value of the uniform Hölder exponent computed from the Fourier coefficients will not be optimal. Nothing can be said from Fourier analysis about the local regularity, in contrast to two-microlocal wavelet analysis discussed below.

## 12 Two-microlocal analysis

The goal of the two-microlocal analysis is to reveal the pointwise behavior of any function from the properties of its wavelet coefficients. The local regularity of the function is thus established.

The scalar product of an analyzed one-dimensional function and a wavelet maps this function into two-dimensional wavelet coefficients which reveal both scale and location properties of a signal. In the above definitions of the pointwise regularity of the function at $`x_0`$ determined by the Hölder exponent only local but not scaling properties of the signal are taken into account. The problem of determining the exact degree of the Hölder regularity of a function is easily solved if this regularity is everywhere the same, because in such a case it usually coincides with the uniform regularity. The determination of the pointwise Hölder regularity, however, becomes much harder if the function changes wildly from point to point, i.e., its scale (frequency) characteristics depend strongly on the location (time). In that case we have to deal with a very non-stationary signal. To describe the singularity of $`f(x)`$ at $`x_0`$ only by the Hölder exponent $`\alpha (x_0)1`$ one looks for the order of magnitude of the difference $`|f(x)f(x_0)|`$ when $`x`$ tends to $`x_0`$, without taking into account, e.g., the possible high-frequency oscillations of $`f(x)f(x_0)`$, i.e., its scale behavior.

The properties of the wavelet coefficients as functions of both scale and location111111Let us note that the word ”location” (and, correspondingly, ”scale”) does not necessarily imply a single dimension but could describe the position (and size or frequency) in any $`n`$-dimensional space. provide a unique way of describing the pointwise regularities. It is more general than the traditional Hölder approach because it allows us to investigate, characterize and easily distinguish some specific local behaviors such as approximate selfsimilarities and very strong oscillatory features like those of the indefinitely oscillating functions which are closely related to the so-called ”chirps” of the form $`x^\alpha \mathrm{sin}(1/x^\beta )`$ (reminding of bird, bat or dolphin sonar signals with very sharp oscillations which accelerate at some point $`x_0`$). Chirps are well known to everybody dealing with modern radar and sonar technology. In physics, they are known to the authors, e.g., in theoretical considerations of the dark matter radiation and of the gravitational waves. Sometimes, similar dependence can reveal itself even in the more traditional correlation analysis . The large value of a second derivative of the phase function of a frequency modulated signal is their typical feature. Such special behavior with frequency modulation laws hidden in a given signal can be revealed  with the help of the two-microlocal analysis . Actually, there is no universally accepted definition of a chirp. Sometimes, any sine of a non-linear polynomial function of time is also called a chirp.

The two-microlocal space $`C^{s,s^{}}(x_0)`$ of the real-valued $`n`$-dimensional functions $`f`$ (distributions) is defined by the following simple decay condition on their wavelet coefficients $`d_{j,k}`$

$$|d_{j,k}(x_0)|C2^{(\frac{n}{2}+s)j}(1+|2^jx_0k|)^s^{},$$
(89)
where $`s`$ and $`s^{}`$ are two real numbers. This is a very important extension of the Hölder conditions. The two-microlocal condition is a local counterpart of the uniform condition (85). It is very closely related to the pointwise regularity condition (83) because it expresses, e.g., the singular behavior of the function itself at the point $`x_0`$ in terms of the $`k`$-dependence of its wavelet coefficients at the same point. Such an estimate is stable under derivation and fractional integration.

For $`s^{}=0`$, the space thus defined is the global Hölder space $`C^s(R^n)`$.

If $`s^{}>0`$, these conditions are weaker at $`x_0`$ than far from $`x_0`$ so that they describe the behavior of a function which is irregular at $`x_0`$ in a smooth enviroment, i.e., the regularity of $`f`$ gets worse as $`x`$ tends to $`x_0`$. If $`s^{}`$ is positive, the inequality (89) implies that $`|d_{j,k}|C2^{(\frac{n}{2}+s)j}`$, so that $`f`$ belongs to $`C^s(R^n)`$. The Hölder norms do not give any information about this singularity because one has to consider them on domains that exclude $`x_0`$, and the Hölder regularity of $`f`$ deteriorates when $`x`$ gets close to $`x_0`$. That is why we need a uniform regularity assumption when $`s^{}`$ is positive.

If $`s^{}<0`$, the converse situation takes place when $`x_0`$ is a point of regularity in a nonsmooth enviroment.

Thus the inequality Eq. (89) clarifies the relationship between the pointwise behavior of a function and the size properties of its wavelet coefficients. If $`s>0`$ and $`s^{}>s`$, then $`C^{s,s^{}}(x_0)`$ is included in $`C^s(x_0)`$. One can prove that the elements of the space $`C^{s,s^{}}(x_0)`$, for $`s^{}>s`$, are functions whereas the elements of $`C^{s,s}(x_0)`$ are (in general) not functions but distributions (and, moreover, quite ”wild distributions” for which the local Hölder condition (83) cannot hold and should be generalized by multiplying the right-hand side by the factor $`\mathrm{log}1/|xx_0|`$). However, if $`f`$ belongs to the space $`C^s(x_0)`$ then it belongs to $`C^{s,s}(x_0)`$ too. In particular, if at $`x`$ close to $`x_0`$ the following inequality $`|f(x)|C|xx_0|^s`$ is valid for $`s0`$, imposing some limit on the nature of the singularity at this point, then the estimates for wavelet coefficients are given by

$$|d_{j,k}|C2^{(\frac{n}{2}+s)j}|k2^jx_0|^s,$$
(90)
if the support of $`\psi _{j,k}`$ is at least at a distance $`2^j`$ from $`x_0`$, and

$$|d_{j,k}|C2^{(\frac{n}{2}+s)j},$$
(91)
if the support of $`\psi _{j,k}`$ is at a distance less than $`2^j`$ from $`x_0`$, which demonstrate the above statement.

The two-microlocal spaces have some stability properties. In particular, the pseudodifferential operators of order 0 are continuous on these spaces (in contrast with the usual pointwise Hölder regularity condition which is not preserved under the action of these operators, for example, of the Hilbert transform). The position of the points of regularity of a function $`f`$ which belongs to the space $`C^{s,s^{}}(x_0)`$ is essentially preserved under the action of singular integral operators such as the Hilbert transform. This property leads to a pointwise regularity result for solutions of partial differential equations which is formulated  by the following

Theorem.
Let $`\mathrm{\Lambda }`$ be a partial differential operator of order $`m`$, with smooth coefficients and elliptic at $`x_0`$. If $`\mathrm{\Lambda }f=g`$ and $`g`$ belongs to $`C^{s,s^{}}(x_0)`$, then $`f`$ belongs to $`C^{s+m,s^{}}(x_0)`$.

Thus, if $`\mathrm{\Lambda }f=g`$, and if $`g`$ is a function that belongs to $`C^s(x_0)`$, then there exists a polynomial $`P`$ of degree less than $`s+m`$ such that, for $`|xx_0|1`$,

$`|f(x)P(x)|C|xx_0|^{s+m}`$.

All the above conditions look as if the wavelets were eigenvectors of the differential operators $`d^\alpha ,|\alpha |r`$, with corresponding eigenvalues $`2^{j|\alpha |}`$.

The more general functional spaces which admit the power dependence of wavelet coefficients on the scale $`j`$ defined by an extra parameter $`p`$ are considered in Ref. . In particular, those are the Sobolev spaces where the integrability is required up to some power of both the function itself and its derivatives up to the definite order. In this case, the additional factor $`j^{2/p}`$ appears in the right-hand side of the inequality (89), or, in other words, in addition to the linear term in the exponent one should consider another logarithmically dependent term. Thus there is no rigorous selfsimilarity (fractality) of the function any more. However, the wavelet analysis allows us to determine the fractal dimensions of the sets of points where the function $`f`$ is singular.

For continuous wavelets, the condition equivalent to Eq. (89) takes the form

$$|W(a,b)|Ca^s(1+\frac{|bx_0|}{a})^s^{}.$$
(92)
The definite upper limits on the behavior of $`|W(a,b)|`$ can also be obtained if $`f`$ is integrable in the neighborhood of the origin (see ). They are somewhat complicated, and we do not show them here.

It is instructive to examine the meaning of the conditions of the type (92). Let us consider the cone in the $`(b,a)`$ half-plane defined by the condition $`|bx_0|<a`$. Within this cone we have $`|W(a,b)|=O(a^s)`$ as $`a0`$. Outside the cone, the behavior is governed by the distance of $`b`$ to the point $`x_0`$. These two behaviors are generally different and have to be studied independently. However, it is shown in  that non-oscillating singularities may be characterized by the behavior of the absolute values of their wavelet coefficients within the cone. It is also shown in  that rapidly oscillating singularities, which we consider below, cannot be characterized by the behavior of their wavelet transform in the cone.

The two-microlocal methods proved especially fruitful in the analysis of oscillating trigonometric and logarithmic chirps.

A simplest definition of trigonometric chirps at the origin $`x=0`$ could simply be read

$$f(x)=x^\alpha g(x^\beta ),$$
(93)
where $`g`$ is a $`2\pi `$-periodic $`C^r`$ function with vanishing integral. This type of behavior leads to the following expansion of continuous wavelet coefficients at small positive $`b\delta `$ on curves $`a=\lambda b^{1+\beta }`$

$$W(\lambda b^{1+\beta },b)=b^\alpha m_\lambda (b^\beta ),$$
(94)
where $`m_\lambda `$ is a $`2\pi `$-periodic function with a definite norm. It shows that rapidly oscillating singularities cannot be characterized by the behavior of their wavelet transform in the cone as in the previous examples. In this case, the wavelet coefficients should be carefully analyzed not in the cone but on definite curves because they are concentrated along these ridges. For small enough scales, such a ridge lies outside the influence cone of the singularity. A typical example is given by the function $`f(x)=\mathrm{sin}1/x`$ whose instantaneous frequency tends to infinity as $`x0`$. The wavelet coefficient modulus is maximum on a curve of equation $`b=Ca^2`$ for some constant $`C`$ depending only on the wavelet and this curve is not a cone. Therefore it is not sufficient to study the behavior of wavelet coefficients inside the cone to characterize the rapidly oscillating singularities. Let us note, that in case of noisy signals, the contribution of the deterministic part of the signal may be expected to be much larger than that of a noise just near the ridges of the wavelet transform because the wavelet transform of the noise is spread in the whole $`(j,k)`$-plane. It can be thus used for denoising the signal.

Logarithmic chirps have an approximate scaling invariance near the point $`x_0`$. A function $`f`$ has a logarithmic chirp of order ($`\alpha ,\lambda `$) and of regularity $`\gamma 0`$ at the origin $`x=0`$ if for $`x>0`$ there exists a $`\mathrm{log}\lambda `$-periodic function $`G(\mathrm{log}x)`$ in $`C^\gamma (R)`$ such that

$$f(x)P(x)=x^\alpha G(\mathrm{log}x).$$
(95)
The continuous wavelet coefficients of it satisfy the condition

$$W(a,b)=a^\alpha H(\mathrm{log}a,b/a),$$
(96)
where $`H`$ is $`\mathrm{log}\lambda `$-periodic in the first argument and its behavior in $`b/a`$ is restricted by some decay conditions. In other words, if the following scaling property is observed

$$C(\lambda a,\lambda b)=\lambda ^\alpha C(a,b)$$
(97)
for $`\lambda <1`$ and small enough $`a`$ and $`b`$, then $`f`$ has a logarithmic chirp of order ($`\alpha ,\lambda `$). More general and rigorous statements and proofs can be found in Ref. . Thus we conclude that the type and behavior of a chirp can be determined from the behavior of its wavelet coefficients.

As an example of the application of the above statements, let us mention the continuous periodic function $`\sigma (x)`$ represented by the Riemann series

$$\sigma (x)=\underset{1}{\overset{\mathrm{}}{}}\frac{1}{n^2}\mathrm{sin}\pi n^2x,$$
(98)
which belongs to the space $`C^{1/2}`$. It is best analyzed by a continuous wavelet transform using the specific complex wavelet $`\psi (x)=(x+i)^2`$ proposed by Lusin in 1930s in the functional analysis. Then the wavelet transform of the Riemann series is given by the well known Jacobi Theta function. The behavior of wavelet coefficients determines the chirp type and its parameters. With the help of the two-microlocal analysis it was proven that this function has trigonometric chirps at some rational values of $`x=x_0`$ which are ratios of two odd numbers (also, its first derivative exists only at these points) and logarithmic chirps at quadratic irrational numbers (near these points its wavelet transform possesses the scaling invariance properties).

For practical purposes, however, it may happen that very large values of $`j`$ are needed to determine Hölder exponents reliably from the above conditions. For the global Hölder exponent, no assumptions about regularity of wavelets is required whereas determination of the local Hölder exponents requires a more detailed approach.

A nice illustration of application of wavelet analysis for studies of local regularity properties of the function $`f(x)`$ is given in  for the following function

$$f(x)=2e^{|x|}\theta (x1)+e^{|x|}\theta (x+1)\theta (1x)+e^x[(x1)^2+1]\theta (x1),$$
(99)
which is infinitely differentiable everywhere except at $`x=1,0,1`$ where, respectively, $`f,f^{},f^{\prime \prime }`$ are discontinuous. One computes for each of the three points the maxima of wavelet coefficients $`A_j=max_k|W_{j,k}|`$ at various resolution levels $`3<j<10`$ and plots $`\mathrm{log}A_j/\mathrm{log}2`$ versus $`j`$. The linear dependences on $`j`$ are found at these points. The slopes at $`x=1,0,1`$ are, correspondingly, -0.505 ; -1.495 ; -2.462, leading with pretty good accuracy of 1.5$`\%`$ to rather precise estimates of the Hölder exponents 0, 1, 2.

Even better accuracy in determining the local regularity of a function can be achieved with the help of redundant wavelet families (frames) where the translational non-invariance is much less pronounced, the wavelet regularity plays no role and only the number of its vanishing moments is important.

If orthonormal wavelets are used, then their regularity can become essential. Typically, they are continuous, have a non-integer uniform Hölder exponent and different local exponents at different points. In fact, there exists a whole hierarchy of (fractal) sets in which these wavelets have different Hölder exponents. There is a direct relation between the regularity of the function $`\psi `$ and the number of its vanishing moments. The more moments are vanishing, the more smooth is the function $`\psi `$. For example, for widely used Daubechies wavelets with a compact support the degree of regularity increases linearly with the maximum number of vanishing moments $`N`$ for higher-tap wavelets as $`\mu N`$ with $`\mu =0.2`$ and, correspondingly, with the support width, as was already mentioned in Section 4. These properties justify the name ”mathematical microscope”, which is sometimes bestowed on the wavelet transform. Let us note that the choice of the $`h_k`$ that leads to maximal regularity is, however, different from the choice with maximal number of vanishing moments for $`\psi `$. How well the projections of a multiresolution approximation $`V_j`$ converge to a given function $`f`$ depends on the regularity of the functions in $`V_0`$.

There exists a special class of wavelets called vaguelets. They possess the regularity properties characterized by the following restrictions:

$$|\psi _{j,k}(x)|C2^{nj/2}(1+|2^jxk|)^{n\alpha },$$
(100)
$$|\psi _{j,k}(x^{})\psi _{j,k}(x)|C2^{j(n/2+\beta )}|x^{}x|^\beta $$
(101)
with $`0<\beta <\alpha <1`$ and a constant $`C`$. Surely, the cancellation condition (21) must be also satisfied. The norm of any function $`f`$ is then limited by its wavelet coefficients as

$$fC^{}(\underset{j,k}{}|d_{j,k}|^2)^{1/2}$$
(102)
with a constant $`C^{}`$. Any continuous linear operator $`T`$ on $`L^2`$ which satisfies the condition $`T(1)=0`$ transforms an orthonormal wavelet basis into vaguelets.

In practice, the regularity of a wavelet can become especially important during the synthesis when after omission of small wavelet coefficients it is better to deal with a rather smooth $`\psi `$ to diminish possible mistakes at the restoration stage. On the contrary, for analysis it seems to be more important to have wavelets with many vanishing moments to ”throw out” the smooth polynomial trends and reveal potential singularities. Also, the large number of vanishing moments leads to better compression but can enlarge mistakes at inverse procedure of reconstruction because of worsened regularity of wavelets. The use of the biorthogonal wavelets helps a lot at this stage because among two dual wavelets one has many vanishing moments whereas another possesses good regularity properties. By choosing the appropriate inversion formula one can minimize possible mistakes.

## 13 Wavelets and fractals

Some signals (objects) possess the self-similar (fractal) properties (see, e.g., ). It means that by changing the scale one observes at a new scale the features similar to those previously noticed at other scales. This property leads to power-like dependences. The formal definition of a (mono)fractal Hausdorff dimension $`D_F`$ of a geometrical object is given by a condition

$$0<lim_{ϵ0}N(ϵ)ϵ^{D_F}<\mathrm{},$$
(103)
which states that $`D_F`$ is the only value for which the product of the minimal number of the (covering this object) hypercubes $`N(ϵ)`$ with a linear size $`l=ϵ`$ and of the factor $`ϵ^{D_F}`$ stays constant for $`ϵ`$ tending to zero. In the common school geometry of homogeneous objects, it coincides with the topological dimension. The probability $`p_i(ϵ)`$ to belong to a hypercube $`N_i(ϵ)`$ is proportional to $`ϵ^{D_F}`$, and the sum of moments is given by

$$\underset{i}{}p_i^q(ϵ)ϵ^{qD_F}.$$
(104)
The fractal dimension is directly related to the Hölder exponents.

Moreover, for more general objects called multifractals, the ”fractal exponents” $`\alpha (x_0)`$ (see Eq. (83)) vary from point to point. The Hausdorff dimension of the set of points $`x_0`$, where $`\alpha (x_0)=\alpha _0`$, is a function of $`d(\alpha _0)`$ whose graph determines multifractal properties of a signal, let it be Brownian motion, a fully developed turbulence or a purely mathematical construction of the Riemann series. Thus the weights of various fractal dimensions inside a multifractal differ for different multifractals, and the value $`D_F`$ is now replaced by $`D_{q+1}`$ which depends on $`q`$. It is called the Renyi (or generalized) dimension . Usually, it is a decreasing function of $`q`$. The fractal (Hausdorff), information and correlation dimensions are, correspondingly, obtained from the Renyi dimension at $`q=1;0;1`$. The difference between the topological and Renyi dimensions is called the anomalous dimension (or codimension). Fractals and multifractals are common among the purely mathematical constructions (the Cantor set, the Serpinsky carpet etc) and in the nature (clouds, lungs etc).

The pointwise Hölder exponents are now determined using wavelet analysis . As we have seen, all wavelets of a given family $`\psi _{j,k}(x)`$ are similar to the basic wavelet $`\psi (x)`$ and derived from it by dilations and translations. Since wavelet analysis just consists in studying the signal at various scales by calculating the scalar product of the analyzing wavelet and the signal explored, it is well suited to reveal the fractal peculiarities. In terms of wavelet coefficients it implies that their higher moments behave in a power-like manner with the scale changing. The wavelet coefficients are less sensitive to noise because they measure, at different scales, the average fluctuations of the signal.

Namely, let us consider the sum $`Z_q`$ of the $`q`$-th moments of the coefficients of the wavelet transform at various scales $`j`$

$$Z_q(j)=\underset{k}{}|d_{j,k}|^q,$$
(105)
where the sum is over the maxima of $`|d_{j,k}|`$. Then it was shown  that for a fractal signal this sum should behave as

$$Z_q(j)2^{j[\tau (q)+\frac{q}{2}]},$$
(106)
i.e.,

$$\mathrm{log}Z_q(j)j[\tau (q)+\frac{q}{2}].$$
(107)
Thus the necessary condition for a signal to possess fractal properties is the linear dependence of $`\mathrm{log}Z_q(j)`$ on the level number $`j`$. If this requirement is fulfilled the dependence of $`\tau `$ on $`q`$ shows whether the signal is monofractal or multifractal. Monofractal signals are characterized by a single dimension and, therefore, by a linear dependence of $`\tau `$ on $`q`$, whereas multifractal ones are described by a set of such dimensions, i.e., by non-linear functions $`\tau (q)`$. Monofractal signals are homogeneous, in the sense that they have the same scaling properties throughout the entire signal. Multifractal signals, on the other hand, can be decomposed into many subsets characterized by different local dimensions, quantified by a weight function. The wavelet transform removes lowest polynomial trends that could cause the traditional box-counting techniques to fail in quantifying the local scaling of the signal. The function $`\tau (q)`$ can be considered as a scale-independent measure of the fractal signal. It can be further related to the Renyi dimensions, Hurst and Hölder (at $`q=1`$ as is clear from the examples in the previous Sections) exponents (for more detail, see Refs. ). The range of validity of the multifractal formalism for functions can be elucidated  with the help of the two-microlocal methods generalized to the higher moments of wavelet coefficients. Thus, wavelet analysis goes very far beyond the limits of the traditional analysis which uses the language of correlation functions (see, e.g., ) in approaching much deeper correlation levels.

Let us note that $`Z_2(j)`$ is just the dispersion (variance) of wavelet coefficients whose average is equal to zero. For positive values of $`q`$, $`Z_q(j)`$ reflects the scaling of the large fluctuations and strong singularities, whereas for negative $`q`$ it reflects the scaling of the small fluctuations and weak singularities, thus revealing different aspects of underlying dynamics.

## 14 Discretization and stability

In signal analysis , real-life applications produce only sequences of numbers due to the discretization of continuous time signals. This procedure is called the sampling of analog signals. Below we consider its implications. At first sight, it seems that in this case the notions of singularities and Hölder exponents are meaningless. Nevertheless, one can say that the behavior of wavelet coefficients across scales provides a good way of describing the regularity of functions whose samples coincide with the observations at a given resolution.

First, we treat the doubly infinite sequence $`f^{(d)}=\{f_n\}=\{f(n\mathrm{\Delta }t)\}`$ for $`\mathrm{}<n<\mathrm{}`$ obtained by sampling a continuous (analog) signal at the regularly spaced values $`t_n=n\mathrm{\Delta }t`$. If within each $`n`$-th interval $`\mathrm{\Delta }t`$ the function $`f`$ can be replaced by the constant value $`f(n\mathrm{\Delta }t)`$ then for small enough $`\mathrm{\Delta }t`$ one gets

$$f_\omega =_{\mathrm{}}^{\mathrm{}}f(t)e^{i\omega t}𝑑t\mathrm{\Delta }t\underset{n=\mathrm{}}{\overset{\mathrm{}}{}}f(n\mathrm{\Delta }t)e^{in\omega \mathrm{\Delta }t}.$$
(108)
The inversion formula reads

$$f_n=\frac{1}{2\pi }_{\mathrm{}}^{\mathrm{}}f_\omega e^{in\omega \mathrm{\Delta }t}𝑑\omega .$$
(109)
The function $`f_\omega `$ is periodic with the period $`2\pi /\mathrm{\Delta }t`$. It means that one can consider it only within the frequency interval $`[\pi /\mathrm{\Delta }t,\pi /\mathrm{\Delta }t)`$. Moreover, for real signals even the interval $`[0,\pi /\mathrm{\Delta }t)`$ is sufficient. For time-limited signals the summation in Eq. (108) is done from $`n_{min}=0`$ to $`n_{max}=N1`$ for $`N`$ sampled values of a signal. For fast Fourier transform algorithms, it asks for $`O(N\mathrm{log}N)`$ computations.

To get the relation between Fourier transforms of the continuous time signal $`f`$ and of the sequence of samples $`f^{(d)}`$ we rewrite the formula (109)

$`f_n={\displaystyle \frac{1}{2\pi }}{\displaystyle \underset{k=\mathrm{}}{\overset{\mathrm{}}{}}}{\displaystyle _{(2k1)\pi /\mathrm{\Delta }t}^{(2k+1)\pi /\mathrm{\Delta }t}}f_\omega e^{in\omega \mathrm{\Delta }t}𝑑\omega =`$
$`{\displaystyle \frac{1}{2\pi }}{\displaystyle _{\pi /\mathrm{\Delta }t}^{\pi /\mathrm{\Delta }t}}{\displaystyle \underset{k=\mathrm{}}{\overset{\mathrm{}}{}}}f_{\omega +2k\pi /\mathrm{\Delta }t}e^{in\omega \mathrm{\Delta }t}d\omega ={\displaystyle \frac{1}{2\pi }}{\displaystyle _{\pi /\mathrm{\Delta }t}^{\pi /\mathrm{\Delta }t}}f_\omega ^{(d)}e^{in\omega \mathrm{\Delta }t}𝑑\omega ,`$ (110)
wherefrom one gets

$$f_\omega ^{(d)}=\underset{k=\mathrm{}}{\overset{\mathrm{}}{}}f_{\omega +2k\pi /\mathrm{\Delta }t}.$$
(111)
Thus the Fourier transform of the discrete sample contains, in general, the contributions from the Fourier transform of the continuous signal not only at the same frequency $`\omega `$ but also at countably infinite set of frequencies $`\omega +2k\pi /\mathrm{\Delta }t`$. These frequencies are called aliases of the frequency $`\omega `$. For ”undersampled” sets, the aliasing phenomenon appears, i.e., the admixture of high frequency components to lower frequencies. The frequency $`\omega _{(N)}=\pi /\mathrm{\Delta }t`$ is called the Nyquist frequency. To improve the convergence of the series, the ”oversampling” is used, i.e., $`f`$ is sampled at a rate exceeding the Nyquist rate.

The band-limited function $`f`$ can be recovered from its samples $`f_n`$’s whenever the sampling frequency $`\omega _s(\mathrm{\Delta }t)^1`$ is not smaller than the band-limit frequency $`\omega _f`$, i.e., whenever $`\omega _s\omega _f`$. This statement is known under the name of the sampling (or Shannon-Kotelnikov) theorem.

In this case, there is no aliasing since only one aliased frequency is in the band limits $`[\omega _f,\omega _f]`$. If $`\omega _s<\omega _f`$, then the function $`f`$ cannot be recovered without additional assumptions.

The rigorous proof of the theorem exists, but it is quite clear already from intuitive arguments that one cannot restore the high-frequency content of a signal by sampling it with lower frequency. Therefore, to get knowledge of high-frequency components one should sample a signal with frequency exceeding all frequencies important for a given physical process. Only then the restoration of a signal is stable. One can reduce the sampling frequency $`\omega _s`$ all the way down to $`\omega _f`$ without losing any information. It allows to subsample the signal by keeping only smaller sets of data, i.e., to get a shorter sampled signal.

In case of wavelet analysis , to have a numerically stable reconstruction algorithm for $`f`$ from $`d_{j,k}`$ one should be sure that $`\psi _{j,k}`$ constitute a frame. For better convergence, one needs frames close to tight frames, i.e., those satisfying the condition $`(\frac{B}{A}1)1`$. The orthonormal wavelet bases have good time-frequency localization. In principle, if $`\psi `$ itself is well localized both in time and in frequency, then the frame generated by $`\psi `$ will share that property as well.

Discrete wavelets are quite well suited for the description of functions by their mean values at the equally spaced points. However in practical real-life applications, apart from this projection of a function, one has to deal with a finite interval of its definition and with a finite number of resolution levels. As was mentioned in Section 5, one usually rescales the ”units” of levels by assuming that the label of the finest available scale is $`j=0`$, and the coarser scales have positive labels $`1jJ`$. The fine-level coefficients are defined by the sampled values of $`f(x)`$ as $`s_{0,k}=f(k)`$, i.e., instead of the function $`f(x)`$ one considers its projection $`P_0f`$. The higher-level coefficients are found from the iterative relations (40), (41), i.e., with the help of fast wavelet transform without direct calculation of integrals $`f(x)\psi _{j,k}(x)𝑑x`$ and $`f(x)\phi _{j,k}(x)𝑑x`$. Therefore the approximate representation of the function $`f(x)`$ which corresponds to the redefined versions of Eqs. (5) and (22) can be written as

$$f(x)P_0f=\underset{j=1}{\overset{J}{}}\underset{k}{}d_{j,k}\psi _{j,k}+\underset{k}{}s_{J,k}\phi _{J,k},$$
(112)
where the sum over $`k`$ is limited by the interval in which $`f(x)`$ is defined. Moreover, to save the computing time, one can use not a complete set of wavelet coefficients $`d_{j,k}`$ but only a part of them omitting small coefficients not exceeding some threshold value $`ϵ`$. This standard estimation method is called estimation by coefficient thresholding. If the sum in (112) is taken only over such coefficients $`|d_{j,k}|>ϵ`$ and the number of the omitted coefficients is equal to $`n_o`$, then the function $`f_ϵ`$ which approximates $`f(x)`$ in this situation will differ from $`f(x)`$ by the norm as

$$f(x)f_ϵ(x)<ϵn_o^{1/2}.$$
(113)
Thus for a function, which is quite smooth in the main domain of its definition and changes drastically only in very small regions, the threshold value $`ϵ`$ can be extremely small. Then it admits large number of wavelet coefficients to be omitted with rather low errors in final approximations. Instead of this so-called hard-shrinkage procedure, one can use a soft-shrinkage thresholding  which shifts positive and negative coefficients to their common origin after the omission procedure, i.e., replaces non-omitted $`d_{j,k}`$ in the formula (112) by

$$d_{j,k}^{(ϵ)}=sign(d_{j,k})(|d_{j,k}|ϵ).$$
(114)
It has been proven that such an approach leads to optimal min-max estimators.

One can find the coefficients of the expansion (112) with use of the fast wavelet transform since the coefficients $`s_{0,k}`$ are fixed as the discrete values of $`f(x)`$. In iterative schemes, the error will however accumulate and their precision will not be sufficiently high. Much better accuracy can be achieved if the interpolation wavelets are used. In this case, the values of the function on the homogeneous grid $`f(k)`$ are treated as $`s`$-coefficients for the interpolation basis, and initial values of $`s_{0,k}`$ are formed by some linear combinations of them with coefficients which depend on the shapes of wavelets considered.

The more elaborate procedure of weighting different wavelet coefficients before their omission called quantization is often used instead of this simplified procedure of the coefficient thresholding (see also Section 15.4). According to experts estimates, different weights (importance) are ascribed to different coefficients from the very beginning. These weights depend on the final goals of compression. Only then the thresholding is applied. Many orthogonal wavelet bases are infinitely more robust with respect to this quantization procedure than the Fourier trigonometric basis. Inspite of this, the lack of symmetry for Daubechies wavelets with compact support does not always satisfy the experts because some visible defects appear. This problem can be cured when one uses symmetric biorthogonal wavelets having compact support.

## 15 Some applications

In this section we describe mostly those applications of wavelets which are closest to our personal interests (the brief summary is given in Ref. ). Even among them we have to choose those where, in our opinion, the use of wavelets was crucial for obtaining principally new information unavailable with other methods (see Web site www.awavelet.com).

### 15.1 Physics

Physics applications of wavelets are so numerous that it is impossible to describe even the most important of them here in detail (see, e.g., ). They are used both in purely theoretical studies in functional calculus, renormalization in gauge theories, conformal field theories, nonlinear chaoticity, and in more practical fields like quasicrystals, metheorology, acoustics, seismology, nonlinear dynamics of accelerators, turbulence, structure of surfaces, cosmic ray jets, solar wind, galactic structure, cosmological density fluctuations, dark matter, gravitational waves etc. This list can easily be made longer. We discuss here two problems related to use of wavelets for solving the differential equations with an aim to get the electronic structure of an extremely complicated system, and for pattern recognition in multiparticle production processes in high energy collisions.

#### 15.1.1 Solid state and molecules

The exact solution of a many-body problem is impossible, and one is applying some approximate methods for solid state problems, e.g., the density functional theory . However, electronic spectra of complex atomic systems are so complicated that, even within this approach, it is impossible, in practice, to decipher them by commonly used methods. For example, it would require to calculate about $`2^{100}`$ Fourier coefficients to represent the effective potential of the uranium atom (and even more, in case of uranium dimers). This is clearly an unrealistic task. Application of the wavelet analysis methods makes it possible to resolve this problem . The potential of the uranium dimer is extremely singular. It varies by more than 10 orders of magnitude. Its reconstruction has now become a realistic problem. Quite high precision has been achieved with the help of wavelets as seen in Fig. 13.

The solution of the density functional theory equations by the wavelet methods was obtained by several groups . The method described in Refs  has been applied to a large variety of different materials. It possesses good convergence properties as is shown in Figs. 14 and 15. They demonstrate how fast is the decrease of energy (or wave function) recipies with respect to the number of iterations (Fig. 14) or of the absolute values of the wavelet coefficients (Fig. 15) with respect to their index (when the coefficients are ordered according to their magnitude). Such a fast decrease allows us to consider much more complicated systems than it was possible before. This method has been successfully tested for solid hydrogen crystals at high pressure, manganate and hydrogen clusters , 3d-metals and their clusters etc. The density matrix can be effectively represented with the help of Daubechies wavelets .

#### 15.1.2 Multiparticle production processes

First attempts to use wavelet analysis in multiparticle production go back to P. Carruthers  who used wavelets for diagonalisation of covariance matrices of some simplified cascade models. The proposals of correlation studies in high multiplicity events with the help of wavelets were promoted , and also, in particular, for special correlations typical for the disoriented chiral condensate . It was recognized  that wavelet analysis can be used for pattern recognition in individual high multiplicity events observed in experiment.

High energy collisions of elementary particles result in production of many new particles in a single event. Each newly created particle is depicted kinematically by its momentum vector, i.e., by a dot in the three-dimensional phase space. Different patterns formed by these dots in the phase space would correspond to different dynamics. To understand this dynamics is a main goal of all studies done at accelerators and in cosmic rays. Especially intriguing is a problem of the quark-gluon plasma, the state of matter with deconfined quarks and gluons which could exist during an extremely short time intervals. One hopes to create it in collisions of high energy nuclei. Nowadays, the data about Pb-Pb collisions are available where, in a single event, more than 1000 charged particles are produced. We are waiting for RHIC accelerator in Brookhaven and LHC in CERN to provide events with up to 20000 new particles created. However we do not know yet which patterns will be drawn by the nature in individual events. Therefore the problem of phase space pattern recognition in an event-by-event analysis becomes meaningful.

It is believed that the detailed characterization of each collision event could reveal the rare new phenomena, and it will be statistically reliable due to a large number of particles produced in a single event.

When individual events are imaged visually, the human eye has a tendency to observe different kinds of intricate patterns with dense clusters and rarefied voids. Combined with search for maxima (spikes) in pseudorapidity 121212The pseudorapidity is defined as $`\eta =\mathrm{log}\mathrm{tan}(\theta /2)`$, where $`\theta `$ is a polar angle of particle emission. (polar angle) distributions, it has lead to indications on the so-called ring-like events in cosmic ray studies  and at accelerators . However, the observed effects are often dominated by statistical fluctuations. The method of factorial moments was proposed  to remove the statistical background but it is hard to apply in event-by-event approach. The wavelet analysis avoids smooth polynomial trends and underlines the fluctuation patterns. By choosing the strongest fluctuations, one hopes to get those which exceed the statistical component. The wavelet transform of the pseudorapidity spectra of cosmic ray high multiplicity events of JACEE collaboration was done in Ref. .

In Ref.  wavelets were first applied to analyze patterns formed in the phase space of the accelerator data on individual high multiplicity events of Pb-Pb interaction at energy 158 GeV per nucleon. With emulsion technique used in experiment the angles of emission of particles are measured only, and the two-dimensional (polar+azimuthal amgles) phase space is considered therefore. The experimental statistics is rather low but acceptance is high and homogeneous that is important for the proper pattern recognition. To simplify the analysis , the two-dimensional target diagram representing the polar and azimuthal angles of created charged particles was split into 24 azimuthal angle sectors of $`\pi /12`$ and in each of them particles were projected onto the polar angle $`\theta `$ axis. Thus one-dimensional functions of the polar angle (pseudorapidity) distribution of these particles in 24 sectors were obtained. Then the wavelet coefficients were calculated in all of these sectors and connected together afterwards (the continuous MHAT wavelet was used). The resulting pattern showed that many particles are concentrated close to some value of the polar angle, i.e., reveal the ring-like structure in the target diagram. The interest to such patterns is related to the fact that they can result from the so-called gluon Cherenkov radiation  or, more generally, from the gluon bremsstrahlung at a finite length within quark-gluon medium (plasma, in particular). More elaborate two-dimensional analysis with Daubechies ($`D^8`$) wavelets was done recently  and confirmed these conclusions with jet regions tending to lie on some ring-like formations. It is seen, e.g. in Fig. 16 where dark regions correspond to large wavelet coefficients of the large-scale particle fluctuations in two of the events analyzed. The resolution levels $`6j10`$ were left only after the event analysis was done to store the long-range correlations in the events and get rid of short-range ones and background noise. Then the inverse restoration was done to get the event images with these dynamic correlations left only, and this is what is seen in Fig. 16. It directly demonstrates that large-scale correlations chosen have a ring-like (ridge) pattern. With larger statistics, one will be able to say if they correspond to theoretical expectations. However preliminary results favor positive conclusions . It is due to the two-dimensional wavelet analysis that for the first time the fluctuation structure of an event is shown in a way similar to the target diagram representation of events on the two-dimensional plot.

Previously, some attempts  to consider such events with different methods of treating the traditional projection and correlation measures revealed just that such substructures lead to spikes in the polar angle (pseudorapidity) distribution and are somewhat jetty. Various Monte Carlo simulations of the process were compared to the data and failed to describe this jettiness in its full strength. More careful analysis  of large statistics data on hadron-hadron interactions (unfortunately, however, for rather low multiplicity) with dense groups of particles separated showed some ”anomaly” in the angular distribution of these groups awaited from the theoretical side. Further analysis using the results of wavelet transform is needed when many high multiplicity events become available. The more detailed review of this topic is given in Ref. .

### 15.2 Aviation (engines)

Multiresolution wavelet analysis has proved to be an extremely useful mathematical method for analyzing complicated physical signals at various scales and definite locations. It is tempting to begin with analysis of signals depending on a single variable131313The trick of reducing the two-dimensional function to a set of one-dimensional ones was described in the previous subsection devoted to particle production.. The time variation of the pressure in an aircraft compressor is one of such signals. The aim of the analysis of this signal is motivated by the desire to find the precursors of a very dangerous effect (stall+surge) in engines leading to their destruction.

An axial multistage compressor is susceptible to the formation of a rotating stall which may be precipitated by a distorted inlet flow induced by abrupt manoeuvres of an aircraft (helicopter) or flight turbulence. The aerodynamic instability behind the rotating blades is at the origin of this effect. The instability regions called stall cells separate from the blades and rotate with a speed about 60$`\%`$ of the rotor speed. Thus they are crossed by the blades approximately one time every 1.6 rotor revolutions. It limits abruptly the flow and thus the pressure delivery to the downstream plenum (the combustion chamber) where the fuel is burnt up and gas jet is formed. The high prestall pressure build-up existing in the combustion chamber tends to push the flow backwards through the compressor. If the flow is reversed one calls it ”deep surge”. In many cases it is attested by flames emerging from the engine inlet. It can have serious consequences for the engine life and operation, not to say about an aircraft and its passengers if it happens at in-flight conditions. That is why the search for precursors of these dangerous phenomena is very important. Attempts to predict the development of a rotating stall and a subsequent surge with a velocity measuring probe such as a hot wire anemometer  provide a precursor or warning of only about 10 msec what is not enough for performing any operations which would preclude surge development.

Multiresolution wavelet analysis of pressure variations in a gas turbine compressor reveals much earlier precursors of stall and surge processes . Signals from 8 pressure sensors positioned at various places with the compressor were recorded and digitized at 3 different modes of its operation (76$`\%`$, 81$`\%`$, 100$`\%`$ of nominal rotation speed) in stationary conditions with interval 1 msec during 5–6 sec before the stall and the surge developed. An instability was induced by a slow injection of extra air into the compressor inlet. After a few minutes, this led to a fully blown instability in the compressor. The interval of last 5–6 sec was wavelet-analyzed. Since the signal fluctuates in time, so too does the sequence of wavelet coefficients at any given scale. A natural measure for this variability is the standard deviation (dispersion) of wavelet coefficients as a function of scale. A scale of $`j=5`$ showed the remarkable drop about 40$`\%`$ of the dispersion (variance) of the wavelet coefficients for more than 1 sec prior the malfunction develops. The dispersion is calculated according to the standard expression

$$\sigma (j,M)=\sqrt{\frac{1}{M1}\underset{k=0}{\overset{M1}{}}[d_{j,k}d_{j,k}]^2},$$
(115)
where $`M`$ is the number of wavelet coefficients at the scale $`j`$ within some time interval which was chosen equal to 1 sec.

In Fig. 17 the time variation of the pressure in the compressor at one of the indicated above modes of engine operation is shown by a very irregular line. Large fluctuations of the pressure at the right-hand side denote the surge onset. It is hard to guess about any warning of this drastic instability from the shape of this curve. The wavelet coefficients of this function were calculated with Daubechies 8-tap ($`D^8`$) wavelets. The dashed curve in the same Figure shows the behavior of the dispersion of wavelet coefficients as a function of time at the scale level $`j=5`$ which happens to be most sensitive to the surge onset141414It shows that the correlations among successive values of the sampled signal are likely to occur on intervals of the order of 30 msec. Let us note that similar correlation lengths are typical for a sampled speech signal as well.. Its precursor is seen as the maximum and the subsequent drop of the standard deviation by about 40$`\%`$ which appears about 1 - 2 sec prior the malfunction denoted by the large increase of both the pressure and dispersion at the right-hand side. Initial part of the dispersion plot is empty because the prescribed interval for the compilation of the representative distribution of wavelet coefficients is chosen equal to 1 sec. The randomly reordered (shuffled) sample of the same values of the pressure within the pre-surge time interval does not show such a drop of the dispersion as demonstrated by the dot-dashed line in the Figure. It proves the dynamic origin of the effect. Further analyzed characteristics of this process (fractality, high-moments behavior) are discussed in . No fractal properties of the signal have been found. The scale $`j=5`$ falls off the linear dependence of the partition function on $`j`$ necessary for fractal behavior. Let us note that the multifractality of the signal may fail because of accumulation of points where a chirp-like behavior happens. If it is the origin of this effect, it should provide a guide to the description of the dynamics of the analyzed processes. The time intervals before and after appearence of the precursor were analyzed separately. Higher moments of wavelet coefficients as functions of their rank $`q`$ behave in a different way at pre- and post-precursor time. It indicates on different dynamics in these two regions.

This method of wavelet analysis can be applied to any engines, motors, turbines, pumps, compressors etc. It provides significant improvement in the diagnostics of the operating regimes of existing engines, which is important for preventing their failure and, consequently, for lowering associated economic losses. Two patents on diagnostics and automatic regulation of engines dated 19.03.1999 have been obtained by the authors of Ref. .

There are other problems in aviation which could be approached with wavelet analysis , e.g., combustion instabilities, analysis of ions in the jet from the combustion chamber or more general problem of metal aging and cracks etc.

### 15.3 Medicine and biology

Applications of the wavelet analysis to medicine and biology (see, e.g., ) follow also the same procedures as discussed above . They are either deciphering information hidden in one-dimensional functions (analysis of heartbeat intervals, ECG, EEG, DNA sequences etc) or pattern recognition (shapes of biological objects, blood cell classification etc).

#### 15.3.1 Hearts

Intriguing results were obtained  when the multiresolution wavelet analysis was applied to the sequence of time intervals between human heartbeats. It has been claimed that the clinically significant measure of the presence of heart failure is found just from analysis of heartbeat intervals alone without knowledge of the full electrocardiogram plot while the previous approaches provided statistically significant measures only. Series of about 70000 interbeat intervals were collected for each of 27 patients treated. They were plotted versus interval number. Signal fluctuations were wavelet transformed and dispersions of wavelet coefficients at different scales were calculated as in the case of aviation engines considered above but with averaging over the whole time interval because it is not necessary now to study the evolution during this time interval. Thus each patient was characterized by a single number (dispersion) at a given scale. It happened that the sets of numbers for healthy and heart failure patients did not overlap at the scale151515By coincidence it is the same as for aviation engines. $`j=5`$ as seen in Fig. 18. This is considered as the clinically significant measure in distinction to the statistically significant measure when these sets overlap partly. The dispersions for healthy patients are larger probably corresponding to a higher flexibility of healthy hearts. The fluctuations are less anti-correlated for heart failure than for healthy dynamics.

This analysis was questioned by another group  stating that the above results on the scale-dependent separation between healthy and pathologic behavior depend on the database analyzed and on the analyzing wavelet . They propose new scale-independent measures – the exponents characterizing the scaling of the partition function of the wavelet coefficients of the heartbeat records. These exponents reveal the multifractal behavior for healthy human heartbeat and loss of multifractality for a life-threatening condition, congestive heart failure. This distinction is ascribed to nonlinear features of the healthy heartbeat dynamics. The authors claim that the multifractal approach to wavelet analysis robustly discriminates healthy subjects from heart-failure subjects. In our opinion, further studies are needed.

Earlier, the Fourier analysis of the heart rate variability disclosed the frequency spectrum of the process. Three frequency regions play the main role. The high frequency peak is located about 0.2 Hz, and the low frequency peak is about 0.1 Hz. The superlow frequency component reminds $`1/f`$-noise with its amplitude strongly increasing at low frequencies. There are some attempts  to develop theoretical models of such a process with the help of wavelet transform.

The wavelet analysis of one-dimensional signals is useful also for deciphering the information contained in ECG and EEG. The functions to be analyzed there are more complicated than those in above studies. Some very promising results have been obtained already. In particular, it was shown that anomalous effects in ECG appear mostly at rather large scales (low frequencies) whereas the normal structures are inclined to comparatively small scales (high frequencies). It corresponds to the above results on heartbeat intervals. Such analysis of ECG was reported, e.g., in Ref. . We shall not describe it in detail here because these studies are still in their initial stage only. The time-frequency analysis of EEG can be found, e.g., in Ref. . It can locate the source of the epileptic activity and its propagation in the brain. The general description of EEG methodics see, e.g., in Refs. .

The wavelet analysis has been used also for diagnostics of the embrion state during the pregnancy period . The matching pursuit method was developed to fit the properties of the signal as well as possible.

#### 15.3.2 DNA sequences

The wavelet transform with continuous wavelets has been used for fractal analysis of DNA sequences  in attempt to reveal the nature and the origin of long range correlations in these sequences. It is still of debate whether such correlations are different for protein-coding (exonic) and noncoding (intronic, intergenic) nucleotide sequences. To graphically portray these sequences in the form of one-dimensional functions, the so-called ”DNA walk” analysis was applied with the conversion of the 4-letter text of DNA into a binary set . Applying wavelet transform to 121 DNA sequences (with 47 coding and 74 noncoding regions) selected in the human genome, the authors of Refs.  have been able to show that there really exists the difference between the two subsequences. It is demonstrated by the presence of long range correlations in noncoding sequences while the coding sequences look like uncorrelated random walks. To show this, the averaged partition (generating) functions of wavelet coefficients (see Eq.(105)) over these two statistical samples were calculated. Both of them scaled with the exponent predicted for homogeneous Brownian motions, i.e., $`\tau (q)=qH1`$. The main difference between the noncoding and coding sequences is the value of $`H`$, namely, $`H_{nc}=0.59\pm 0.02`$ and $`H_c=0.51\pm 0.02`$ which distinguish uncorrelated and correlated subsamples. Moreover, it has been shown that $`\tau (q)`$ spectra extracted from both sets are surprisingly in remarkable agreement with the theoretical prediction for Gaussian processes if the probability density of wavelet coefficients at fixed scales are plotted versus these coefficients scaled with their r.m.s. value at the same scale. The results of wavelet analysis clearly show that the purine (A, G) versus pyrimidine (C, T) content161616A, C, G, T denote the usual four letter alphabet of any DNA text. of DNA is likely to be relevant to the long range correlation properties observed in DNA sequences.

#### 15.3.3 Blood cells

Another problem which can be solved with the help of wavelet analysis is the recognition of different shapes of biological objects. By itself, it has very wide range of applicability. Here, we consider the blood cell classification scheme according to the automatic wavelet analysis developed by us, and a particular illustration is given for erythrocytes classification. The automatic search, stability in determining the cells shapes and high speed of processing can be achieved with the computer wavelet analysis . The main idea of the method relies on the fact that at a definite resolution scale wavelet analysis reveals clearly the contours of the blood cells what allows to classify them. In Fig. 19 we demonstrate how it became possible to improve the image resolution by such a method. The low quality image of blood cells has been transformed in its wavelet image with clearly seen contours of the cells.

The shapes of erythrocytes differ for various types of them. Depending on their shape, they can play either positive or negative role for a human being. Therefore their differentiation is very important. After microscope analysis of a dry blood smear and registration of blood cells in a computer, the wavelet analysis of the set of registered blood cells has been done.

However the extreme random irregularities at some points of a cell contour can prevent from performing the analysis . Therefore a special smearing procedure was invented to avoid such points without loss of crucial peculiarities typical for a definite type of cells. After that, the set of blood cells with somewhat smeared contours is ready for wavelet analysis . It consists in the wavelet correlation analysis which shows different behavior of the correlation measures depending on the particular cell shapes. Using the expert classification of cells in a definite sample, the wavelet characteristics of correlations typical for a particular class were found and then inserted in the computer program. This software was used for classification of blood cells obtained from other patients, and results were cross-checked by experts. Their positive conclusions are encouraging. The whole procedure is now done fast and automatically without human intervention. In Fig. 20 the blood cells of various kinds are shown.

### 15.4 Data compression

Data compression is needed if one wants, e.g., to store the data spending as low memory capacity as possible or to transfer it at a low cost using smaller packages. The example of FBI using wavelet analysis for pattern recognition and saving in that way a lot of money on computer storage of fingerprints is well known. This is done by omitting small wavelet coefficients after the direct wavelet transform was applied. Surely, to restore the information one should be confident in stability and good quality of the inverse transformation. Therefore, both analysis and synthesis procedures are necessary for the data compression and its subsequent reconstruction. Above, we have shown how successful is wavelet analysis in solving many problems. Due to the completeness of the wavelet system it is well suited for the proper inverse transform (synthesis) as well (see, e.g., Ref. ).

The approach to the solution of the problem strongly depends on the actual requirements imposed on the final outcome. There are at least three of them. If one wants to keep the quality of the restored image (film) practically as good as the initial one, the compression should not be very strong. It is required, e.g., if the experts should not distinguish the compressed and uncompressed copies of a movie when they are shown on two screens in parallel. Another requirement would be important if one wants to compress an image as strongly as possible leaving it still recognizable. It is required, e.g., if one needs to transfer the information in a line with limited capacity. Finally, one can require the whole procedure of analysis and synthesis to be done as fast as possible. It is necessary, e.g., if the information must be obtained immediately but at lower cost. These three situations ask for different choice of wavelets to optimize analysis and synthesis. In all cases wavelet analysis has an advantage over the coding methods which use the windowed Fourier transform but the quantitative estimate of this advantage varies with the problem solved.

Let us remind that any image in the computer should be digitized and saved as the bitmap or, in other words, as a matrix, each element of which describes the color of the point in the original image. The number of elements of the matrix (image points) depends on the resolution chosen in the digital procedure. It is a bitmap that is used for the subsequent reproduction of the image on a screen, a printer etc. However, it is not desirable to store it in such a form because it would ask for the huge computer capacity. That is why, at present, the numerous coding algorithms (compression) of a bitmap are developed, whose effectiveness depends on the image characteristics171717It is evident that to store the ”Black square”, one does not have to deal with a matrix of all black dots but must store just the three numbers showing the width, the height and the color. This is the simplest example of the image compression.. All these algorithms belong to the two categories - they are either the coding with a loss of information or without any loss of it (in that case the original bitmap can be completely recovered by the decoding procedure). In a more general applications, the last algorithm is often called as the data arxivation.

As an example, we consider the compression algorithm with a loss of information. ”The loss of information” means in this case that the restored image is not completely identical to the original one but this difference is practically indistinguishable by a human eye181818The ”compression quality” is usually characterized by a parameter which varies from 0 to 100. Here 100 implies the minimal compression (the best quality), and the restored image is practically indistinguishable from the initial one, while 0 means the maximum compression which still allows to distinguish some details of the original image in the recovered one..

At present, most of the computer stored images (in particular, those used in Internet) with the continuous tone191919that is the images contained many slightly different colors. Ordinarily to store such images in a computer 16 millions colors per pixel are used are coded with the help of the algorithm JPEG (for the detailed description, see ). Main stages of this algorithm are as follows. One splits the image in the matrices with 8 by 8 dots. For each matrix, the discrete cosine-transform is performed. The obtained frequency matrices are subjected to the so-called ”quantization” procedure when the most crucial for the visual recognition elements are chosen according to a special ”weight” table compiled beforehand by specialists after their collective decision was achieved. This is the only stage where the ”loss” of information happens. Then the transformed matrix with the chosen frequencies (scales) becomes compact and it is coded by the so-called entropy method (also called arithmetic or Huffmann method).

Applied above algorithm differs from the described one by use of the wavelets instead the windowed cosine-transform and by the transform of the whole image instead the 8 by 8 matrix only. Fig. 21 demonstrates the original image and its two final images restored after the similar compression according to JPEG and wavelet algorithms. It is easily seen that the quality of the wavelet image is noticeably higher than for JPEG at practically the same size of the coded files. The requirement of the same quality for both algorithms leads to the file sizes 1.5 - 2 times smaller for the wavelet algorithm that could be crucial for the transmission of the image, especially if the transmission line capacity is limited.

### 15.5 Microscope focusing

Surely, the problem of microscope focusing is tightly related to pattern recognition. One should resolve well focused image from that with diffused contours. It is a comparatively easy task for wavelets because in the former case the image gradients at the contour are quite high while in the latest one they become rather vague. Therefore the wavelet coefficients are larger when microscope is well focused to the object and drastically decrease with defocusing. At a definite resolution level corresponding to the contour scale the defocusing effect is strongest. It is demonstrated in Fig. 22. The peak of the wavelet coefficients shows the most focused image. After doing the wavelet analysis of an image, the computer sends a command to shift the microscope so that larger values of the wavelet coefficients and thus better focusing are attained. At other levels it is somewhat less pronounced and, moreover, it is asymmetric depending on whether the microscope is positioned above the focus location or below it. This asymmetry has been used for the automatic microscope focusing with a well defined direction of its movement toward the focus location.

## 16 Conclusions

The beauty of the mathematical construction of the wavelet transformation and its utility in practical applications attract nowadays researchers from both pure and applied science. Moreover, commercial outcome of this research becomes quite important. We have outlined a minor part of activity in this field. However we hope that the general trends in the development of this subject became comprehended and appreciated.

Unique mathematical properties of wavelets made them a very powerful tool in analysis and subsequent synthesis of any signal. The orthogonality property allows to get an independent information from different scales. Normalization assures that at different steps of this transformation the value of information is not changed and mixed up. The property of localization helps get knowledge about those regions in which definite scales (frequencies) play most important role. Finally, the completeness of the system of wavelet functions formed by dilations and translations results in validity of the inverse transformation.

The multiresolution analysis leads to fast wavelet transform and together with the procedure of the nonstandard matrix multiplication to rather effective computing. Analytic properties of functions, local and global Hölder indices, multifractal dimensions etc are subject to wavelet analysis . The natural accomodation of differential operators in this enviroment opens a way to effective solution of differential equations.

All these properties enable us to analyze through the wavelet transform complex signals at different scales and locations, to solve equations describing extremely complicated nonlinear systems involving interactions at many scales, to study very singular functions etc. The wavelet transform is easily generalized to any dimension, and multidimensional objects can be analyzed as well. Therefore, wavelets are indispensable for pattern recognition.

Thus the wavelet applications in various fields are numerous and give nowadays very fruitful outcome. In this review paper we managed to describe only some of these applications leaving aside the main bulk of them. The potentialities of wavelets are still not used at their full strength.

However one should not cherish vain hopes that this machinery works automatically in all situations by using its internal logic and does not require any intuition. According to Meyer , ”no ’universal algorithm’ is appropriate for the extreme diversity of the situations encountered”. Actually it needs a lot of experience in choosing the proper wavelets, in suitable formulation of the problem under investigation, in considering most important scales and characteristics describing the analyzed signal, in the proper choice of the algorithms (i.e., the methodology) used, in studying the intervening singularities, in avoiding possible instabilities etc. By this remark we would not like to prevent newcomers from entering the field but, quite to the contrary, to attract those who are not afraid of hard but exciting research and experience.

Appendix 1

The general approach which respects all properties required from wavelets is known as the multiresolution approximation and defined mathematically in the following way:

DEFINITION

A multiresolution approximation of $`L^2(R^n)`$ is, by definition, an increasing sequence $`V_j,jZ`$, of closed linear subspaces of $`L^2(R^n)`$ with the following properties:

1. $`_{\mathrm{}}^{\mathrm{}}V_j=\{0\},_{\mathrm{}}^{\mathrm{}}V_j`$ is dense in $`L^2(R^n)`$;
2. for all $`fL^2(R^n)`$ and all $`jZ^n`$

$$f(x)V_jf(2x)V_{j+1};$$
3. for all $`fL^2(R^n)`$ and all $`kZ^n`$

$$f(x)V_0f(xk)V_0;$$
4. there exists a function, $`g(x)V_0`$, such that the sequence $`g(xk),kZ^n`$, is an orthonormal (or Riesz) basis of the space $`V_0`$. It is clear, that scaling versions of the function $`g(x)`$ form bases for all $`V_j`$.

Thus the multiresolution analysis consists of a sequence of successive approximation spaces $`V_j`$ which are scaled and invariant under integer translation versions of the central space $`V_0`$. There exists an orthonormal or, more generally, Riesz basis in this space. The Haar multiresolution analysis can be written in these terms as

$$V_j=\{fL^2(R);forallkZ:f|_{[k2^j,(k+1)2^j)}=const\}.$$
(116)
In Fig. 3 we showed what the projections of some $`f`$ on the Haar spaces $`V_0`$, $`V_1`$ might look like. The general distributions are decomposed into series of correctly localized fluctuations of a characteristic form defined by the wavelet chosen.

The functions $`\phi _{j,k}`$ form an orthonormal basis of $`V_j`$. The orthogonal complement of $`V_j`$ in $`V_{j+1}`$ is called $`W_j`$. The subspaces $`W_j`$ form a mutually orthogonal set. The sequence of $`\psi _{j,k}`$ constitutes an orthonormal basis for $`W_j`$ at any definite $`j`$. The whole collection of $`\psi _{j,k}`$ and $`\phi _{j,k}`$ for all $`j`$ is an orthonormal basis for $`L^2(R)`$. This ensures us that we have constructed a multiresolution analysis approach, and the functions $`\psi _{j,k}`$ and $`\phi _{j,k}`$ constitute the small and large scale filters, correspondingly. The whole procedure of the multiresolution analysis has been demonstrated in graphs of Fig. 4.

In accordance with the above formulated goal, one can define the notion of wavelets in the following way:

DEFINITION

A function $`\psi (x)`$ of a real variable is called a (basic) wavelet of class $`m`$ if the following four properties hold:

1. if $`m=0`$, $`\psi (x)`$ and $`\phi (x)`$ belong to $`L^{\mathrm{}}(R)`$; if $`m1,\psi (x),\phi (x)`$ and all their derivatives up to order $`m`$ belong to $`L^{\mathrm{}}(R)`$;
2. $`\psi (x),\phi (x)`$ and all their derivatives up to order $`m`$ decrease rapidly as $`x\pm \mathrm{}`$;
3. $`_{\mathrm{}}^{\mathrm{}}𝑑xx^n\psi (x)=0`$ for $`0nm`$, and $`_{\mathrm{}}^{\mathrm{}}𝑑x\phi (x)=1`$;
4. the collection of functions $`2^{j/2}\psi (2^jxk),\mathrm{\hspace{0.33em}2}^{j/2}\phi (2^jxk),j,kZ`$, is an orthonormal basis of $`L^2(R)`$.

Then the equation (22) is valid. If $`\psi `$ and $`\phi `$ both have compact support, then it gives a decomposition of any distribution of order less than $`m`$. Moreover, the order of the distribution $`f`$ (the nature of its singularities) can be calculated exactly and directly from the size of its wavelet coefficients as has been shown in Sections 11, 12. The functions $`2^{j/2}\psi (2^jxk)`$ are the wavelets (generated by the ”mother” $`\psi `$), and the conditions 1., 2., 3. express, respectively, the regularity, the localization and the oscillatory character. One sees that the property 3. is satisfied for the Haar wavelets for $`m=0`$ only, i.e., its regularity is $`r=0`$. In general, for each integer $`r1`$, there exists a multiresolution approximation $`V_j`$ of $`L^2(R)`$ which is $`r`$-regular and such that the associated real-valued functions $`\phi `$ and $`\psi `$ have compact support. As the regularity increases, so do the supports of $`\phi `$ and $`\psi `$. The wavelet $`2^{j/2}\psi (2^jxk)`$ is ”essentially concentrated” on the dyadic interval $`I=[k2^j,(k+1)2^j)`$ for $`j,kZ`$. Its Fourier transform is supported by $`2^j(2\pi /3)|\omega |2^j(8\pi /3)`$. In fact, it has a one octave frequency range.

Appendix 2

In the wavelet analysis of operator expressions the integral Calderon-Zygmund operators are often used. There are several definitions of them (see, e.g., the monograph ). We give here their definition used by Daubechies .

DEFINITION

A Calderon-Zygmund operator $`T`$ on $`R`$ is an integral operator

$$Tf(x)=𝑑yK(x,y)f(y)$$
(117)
for which the integral kernel satisfies

$$|K(x,y)|\frac{C}{|xy|},$$
(118)
$$|\frac{}{x}K(x,y)|+|\frac{}{y}K(x,y)|\frac{C}{|xy|^2},$$
(119)
and which defines a bounded operator on $`L^2(R)`$.

With such a definition the special care should be taken at $`x=y`$.

Appendix 3

In the literature devoted to the signal processing, the so-called Littlewood-Paley decomposition is often used. It is closely related to the wavelet transform. Therefore we give here the dictionary which relates the Littlewood-Paley coefficients denoted as $`\mathrm{\Delta }_j(f)(x)`$ to both discrete $`(d_{j,k})`$ and continuous $`(W(a,b))`$ wavelet coefficients.

$$\mathrm{\Delta }_j(f)(x)=2^{nj/2}d_{j,2^jx}=W(2^j,x),$$
(120)
or

$$d_{j,k}=2^{nj/2}W(2^j,2^jk)=2^{nj/2}\mathrm{\Delta }_j(f)(2^jk).$$
(121)
Acknowledgements

We are grateful to all our colleagues with whom we worked on wavelet problems. Our special thanks are to A. Leonidov who read this paper and made important suggestions. O.I. is especially grateful to S. Goedecker with whom many methods were developed and discussed. V.N. is also grateful to INTAS for the support (grant 97-103).

Figure Captions

Fig. 1. The histogram and its wavelet decomposition.
The initial histogram is shown in the upper part of the Figure. It corresponds to the level $`j=4`$ with 16 bins (Eq. (1)). The intervals are labelled on the abscissa axis at their left-hand sides. The next level $`j=3`$ is shown below. The mean values over two neighboring intervals of the previous level are shown at the left-hand side. They correspond to eight terms in the first sum in Eq. (4). At the right-hand side, the wavelet coefficients $`d_{3,k}`$ are shown. Other graphs for the levels $`j=2,1,0`$ are obtained in a similar way.

Fig. 2. The Haar scaling function $`\phi (x)\phi _{0,0}(x)`$ and ”mother” wavelet $`\psi (x)\psi _{0,0}(x)`$.

Fig. 3. The analyzed function (a) and its Haar projections onto two subsequent spaces V.

Fig. 4. The graphical representation of the multiresolution analysis with decomposition of $`V_{j+1}`$ space onto its subspace $`V_j`$ and the orthogonal complement $`W_j`$ iterated to the lower levels.

Fig. 5. Daubechies scaling functions (solid lines) and wavelets (dotted lines) for $`M=2,4`$.

Fig. 6. The fast wavelet transform algorithm.

Fig. 7. Coiflets (dotted lines) and their scaling functions (solid lines) for $`M=4`$.

Fig. 8. The wavelet packet construction.

Fig. 9. The example of the wavelet analysis of a two-dimensional plot.
One sees that either horizontal or vertical details of the plot are more clearly resolved in the corresponding coefficients. Also, the small or large size (correlation length) details are better resolved depending on the level chosen.

Fig. 10. The lattices of time-frequency localization for the wavelet transform (left) and windowed Fourier transform (right).

Fig. 11. The matrix representation of the standard approach to the wavelet analysis .
The parts containing non-zero wavelet coefficients are shaded.

Fig. 12. The nonstandard matrix multiplication in the wavelet analysis .

Fig. 13. The Coulomb potential V (the solid line) as a function of the distance and the precision of its calculation $`\mathrm{\Delta }V`$ (the dashed line) for the Uranium dimer in arbitrary units (a.u.).
It is drawn for the direction from one atom to another one in the left-hand side and for the opposite direction in the right-hand side. The complete symmetry is seen with a larger distances shown in the right-hand side.

Fig. 14. Fast decrease of energy recipies with respect to the number of iterations demonstrates good convergence of the solution of the density functional theory equations by the wavelet methods.

Fig. 15. The same as in Fig. 14 for the absolute values of the wavelet coefficients.

Fig. 16. The restored images of long-range correlations in experimental target diagrams.
They show the typical ring-like structure in some events of central Pb-Pb interactions.

Fig. 17. The signal of the pressure sensor (the solid line) and the dispersion of its wavelet coefficients (the dashed line).
The time variation of the pressure in the engine compressor (the irregular solid line) has been wavelet analyzed. The dispersion of the wavelet coefficients (the dashed line) shows the maximum and the remarkable drop prior the drastic increase of the pressure providing the precursor of this malfunction. The shuffled set of the data does not show such an effect for the dispersion of the wavelet coefficients (the upper curve) pointing to its dynamic origin.

Fig. 18. The sets of the values of the dispersions of the wavelet coefficients for the heartbeat intervals of the healthy (white circles) and heart failure (black circles) patients are shown.
They do not overlap at $`j=4`$.

Fig. 19. a) The photo of blood cell, obtained from the microscope.
b) The same photo after the wavelet analysis done.
The blurred image of a blood cell (left) becomes clearly visible one (right) after the wavelet transform.

Fig. 20. The classification of the erythrocytes cells.

Fig. 21. a) The original photo (the file size is 461760 bytes).
b) The photo reconstructed after compression according to the JPEG-algorithm.(the file size is 3511 bytes).
c) The photo reconstructed after compression according to the wavelet algorithm (the file size is 3519 bytes).
Better quality of the wavelet transform is clearly seen when comparing the original image (left) and two images restored after the similar compression by the windowed Fourier transform (middle) and the wavelet transform (right).

Fig. 22. The focusing line.
The values of the wavelet coefficients calculated for different positions of the microscope are shown on the vertical axis. The corresponding positions are numbered on the horizontal axis. The best focusing is obtained at the maximum of the curve (positions 6 - 10). Large wavelet coefficients correspond to the better focused image. The microscope moves to the focus position being driven by computer commands for increase of the wavelet coefficients of the analyzed objects.

|  |  |  |  | j |
| --- | --- | --- | --- | --- |
|  |  |  |  | 4 |
|  |  |  |  | 3 |
|  |  |  |  | 2 |
|  |  |  |  | 1 |
|  |  |  |  | 0 |

Fig. 1

Fig.2

|  |
| --- |
|  |
|  |

Fig. 3

Fig.4

|  |  | M=2 |
| --- | --- | --- |
|  |  | M=4 |

Fig.5

Fig.6

|  |  |
| --- | --- |
|  |  |

Fig.7

Fig.8

Fig.9

|  |  |
| --- | --- |

Fig.10

Fig.11

Fig.12

Fig.13

Fig.14

Fig.15

|  |  |
| --- | --- |

Fig.16

Fig.17

Fig.18

|  |  |
| --- | --- |

Fig.19

Fig.20

| See JPEG-files: | | |
| --- | --- | --- |
| fig21o.jpg | fig21j.jpg | fig21w.jpg |

Fig.21

Fig.22