.
For 1f31 < rr /2, we have cos f3 > 0, and the function in brackets decays exponentially to zero as R increases, thus offsetting the polynomial increase of the factor R a  I . This guarantees that, for any a, the integral along the large arc converges to 0 as R 00. A useful comment is in order here: For f3 = rr /2, when the first summand in brackets is equal to 1, and for a < 1, the integral over the large arc also converges to 0 because of the factor Ra  I . So, in the limits e 0 and R 00, equality (4) takes the form
f
00
o
f
00
taIeptdt = eipa
0
t a  I exp(peiPt)dt.
Chapter 4. Asymptotics ofFourier transforms
118
Multiplying both sides by pa, introducing on the lefthand side a new variable of integration t' = pt, and noticing that the transformed integral on the lefthand side coincides with the Gamma function (3), we get an equality
f
00
r(a)
=
pa eifJa
t a  1 exp(  peifJt)dt,
o
which, after a change of variables, u=peifJ=h+iw,
becomes an equality r(a) = u a
w=psin{J,
h=pcos{J,
1
00
t a 1 exp(ut)dt,
(6)
(7)
valid for a > 0, and h = Re u > O. For a fractional power function u a , its principal branch must be selected. For 0 < a < 1, we can put h = 0 ({J = 7r /2) (see the comment following the analysis of the integral over the large arc) to obtain the formula
f
00
r(a) = (iw)a
t a  1 exp(iwt)dt.
(8)
o
Formulas (7) and (8) have an easy interpretation in terms of Fourier transforms. The integral on the righthand side of (8), up to a factor 1/27r, coincides with the Fourier transform of function g(t; a)
so, its Fourier transform
= x(t)t a 1 ,
0 < a < 1,
_ r(a) g(w, a) = 27r(iw)a
(9)
(10)
In particular, for a = 1/2, we have a remarkable formula
1
x(t)
  t+
./i
(1
+ i)./27rw
•
(11)
Remark 1. Note, that it is easy to obtain the power law g(w) '" l/wa in (10) from dimensional analysis, although this type of argument will not yield the precise value r(a) /27r of the numerical coefficient.
4.4. Gamma function and Fourier transform of power functions
119
Remark 2. The validity of the same power law across the whole frequency range of the Fourier transform reflects two, fundamentally different, behaviors of the original time function. Its validity for w+oo indicates the presence of a discontinuity of the second kind of the original function at time t = 0, while its validity for w + 0, is a consequence of moderate ('" t a  1) decay of the tail of the original function as t + 00.
Remark 3. Equalities (78) were derived with the help of an integration contour in the upper halfplane (see Fig. 4.4.1). Hence (see (6)), frequency w = pcosfJ appearing in these formulas is positive. However, it is easy to show that the above proof remains in force if the contour is reflected into the lower halfplane. Thus, the formulas (711) remain valid for all frequencies 00 < w < 00. In particular,
x(t)
..fi.........
{ (1  i)/.J8rrlWT, (1 + i)/.J8rrlWT,
for w > 0; for w < o.
(12)
Remark 4. Time reversal in the original function results in the frequency changing sign in the Fourier transform (see (3.1.3c)). Hence, it follows from (12) that
x(t) .Jjtf
+
......... { (1 i)1 ,J81l' Iwl, (1  i) I ,J81l' Iwl,
for w > 0; for w < o.
Combining this relationship with (12), we find the Fourier transform of a symmetric
in time function 1/.Jjtf, displaying a discontinuity of the second kind: 1 1   ......... === .Jjtf ,J21l' Iwi·
(13)
Now, let us return to the discussion of the formula (7). With the help of notation introduced in (6), we can rewrite (7) in the form
f
00
r(a) = (h
+ iw)a
t a  1 exp (ht  iwt) dt.
(14)
o
In other words, function (15)
has the Fourier transform
_
r(a)
g(w; a, h) = 21l'(h
+ iw)a
(16)
Chapter 4. Asymptotics of Fourier transfonns
120
Its principal asymptotics at infinity is described by the relation r(a)
_
g(w; a, h) '" 27r(iw)IX'
(Iwl
+ (0),
(17)
which has a form identical to (10), but is valid for any a > O. For fractional a this asymptotic formula is a consequence of the original function's (see (15» discontinuities of the second kind, or of similar discontinuities of its derivatives. For the integer values of a the formula agrees with the asymptotic behavior of Fourier transforms of functions with explicit (in the function itself), or hidden (in the derivatives) discontinuities of the first kind which were discussed in Section 4.3. Despite asymmetry of the function (15) (it is identically equal to 0 for t < 0) one can still consider its even and odd components. For an arbitrary function g(t), these two components are given by 1
geven(t) = 2[g(t)
+ g(t)],
Clearly, geven (t)
1
godd(t) = 2[g(t)  g(t)].
(18)
+ godd(t) = g(t).
It follows from properties (3.1.3c) and (3.1.4) of the Fourier transform of real functions that the Fourier transforms of even and odd components of g correspond, respectively, to the real and the imaginary parts of the Fourier transform of the original function g. More formally, geven(w)
= Re g(w),
godd(W)
= i 1m g(w).
(19)
Separating the real and imaginary parts of the Fourier transform (16), we obtain that (20) The argument K of a complex number h + iw introduced above depends on the dimensionless frequency y via the formula K
= arctany,
y =
w/ h.
(21)
For y + ±oo and K + ±7r /2, the Fourier transforms of even and odd parts of the original function (15) have the following principal asymptotics: _
geven(w) '"
r(a)cos(a7r/2) 2 IX ' 7rW
(22a)
4.4. Gamma function and Fourier transform ofpower functions
_
godd
() w
. r(a) sin(a1l' /2).
"V
1,,_
kJ,Wa
()
Slgn W ,
(Iwl
121
00).
(22b)
Remark 5. In engineering and physical applications, these asymptotic formulas are much more useful and important than the exact formulas (20). Formulas (20) give Fourier transforms of a narrow class of gauge functions, whereas formulas (22) describe asymptotics of Fourier transforms of a much broader class of functions which, at some arbitrary instants of time tk, have local singularities (t  tk)al. "V
Remark 6. For odd values of a, the cosine in asymptotic formula (22a) becomes 0, and for even a the sine in (22b) vanishes. This means that the asymptotics of the corresponding functions is of order smaller than l/wa • This phenomenon is similar to the one already encountered for functions C and S in (4.3.45). Its essence can be explained with the help of two functions: X (t)t 2 and X (t)t 3. The former, extended to an odd function becomes sign (t) t 2 , which has discontinuities of the second derivative, whereas its even extension is an infinitely differentiable function t 2 • The latter becomes infinitely differentiable under the odd extension but has a discontinuity in the third derivative under the even extension. The above functions have Fourier transforms only in the generalized, djstributional sense. However, the general principle stating that Fourier transforms of infinitely differentiable functions decay, for w 00, faster than any power of w, extends to them as well. Consequently, the generalized Fourier transforms of an even function t 2 and the odd function t 3 do not have power asymptotics. In other words, all the coefficients in their asymptotic expansions in powers of l/w are equal to O. This will become clear when we recall the generalized Fourier transforms (3.3.5) of the above two functions. 00, of Fourier In the case of fractional a, the principal asymptotics, as Iwl transforms of both even and odd parts of function (9) are the same and of order l/wa . Either extension of function (9) to the negative halfline does not remove its characteristic absence of smoothness at the origin. Here, the reader may feel that the asymptotics of the Fourier transform of function (9) established rigorously above for any a, contradicts geometric common sense which seems to be telling us that the graphs of even and odd components of (9) are qualitatively different (see Fig. 4.4.3 (a) and (b». The graph of the even component of (23) with a characteristic cusp at t = 0, gives us an impression of a function that is much less smooth than that of the odd component, although the latter also has a vertical slope at t = O. Nevertheless, in spite of our geometric intuition, both of them have the same Fourier transform asymptotics 1/w3/ 2 • The Fourier transform of function (23) can be calculated explicitly with the help of formula (20), but an alternative derivation provided below is more direct.
122
Chapter 4. Asymptotics of Fourier transforms
g(t)
g(t)
even  a
oddb
+t
4t
FIGURE 4.4.3 Graphs ofeven (a) and odd (b) components of function g(t) = X (t).../ie ht , corresponding to a = 3/2. Since
r(3/2) = (1/2)! = ..;;r/2, it follows from (16) that
_(
g w) =
/1 1 V;J;3 4(1 + iy)../1 + iy'
y = w/h.
Separating the real and imaginary parts is easy if we use algebraic identities instead of the general relation (22) which contains trigonometric functions. Squaring both sides of equality ";1 + iy = x + iy, and solving the resulting algebraic equations with respect to x and y, we obtain the main branch of the radical
Further, standard but tedious calculations give that g(w) is equal to
x [(1 _
1+
y2
Jl + y2
) _ iy
(1 + _::1==)] 1+
Jl + y2
.
123
4.5. Generalized Fourier transfonns ofpower functions
4.5
Generalized Fourier transforms of power functions
In the previous section we demonstrated (4.4.1516) that
0 the above formula generates the whole family of generalized In the limit h Fourier transforms, many of them not encountered thus far. We shall study them in this section. Recall (Section 3.2), that g(w) is said to be a generalized Fourier transform of function g(t) if, for any test function t/J(t) E S, and its Fourier transform (also in S),
1
;T(
=
1
g(t)f(J(t)dt.
Note that the above formula corresponds to formula (3.2.6) with 'l' = 0, but is sufficient to uniquely determine distribution g(w). By definition, the Fourier transform of function g from (4.4.15) is equal to
1 00
g(w; a, h) =
talehtiwt dt.
(1)
o Multiplication of both sides of (1) by a function entire waxis, gives that
21f
1
1 00
g(w; a,
=
ta1e ht
E
S, and integration over the
[I
dt.
o
The integral in the brackets is equal to function t/J(t) E S which is absolutely integrable and rapidly (that is, faster than any power) decreases to 0 at infinity. For that reason, the integral on the righthand side of the equality
2T(
1
g(w; a,
1 00
=
ta1ehtt/J(t)dt
o
exists in the classical sense, for any h 2: 0 and a > O. In particular, for h = 0 we
Chapter 4. Asymptotics of Fourier transforms
124
arrive at a symbolic equality
21l'
f
f
00
=
g(w;
taI O. Separating their real and imaginary parts we get that
The regular limit, as h 0, of functions on the righthand side does not exist. However, this is not a serious obstacle as we are interested in the weak convergence, that is, in the result of integration of these functions against an arbitrary test function (iJ(w). In that sense, the real part, which is the familiar Lorentz function, converges to 8(w)/2. The imaginary part, when integrated against an arbitrary smooth and absolutely integrable function (iJ(w), gives that · 11m
h+o
fh
wiP(w) d W= 2
+w2
pvf
rp(w)d w,

w
where PV f stands for the principal value of the integral. At this point, we will not dwell on the notion of the principal value of the integral as it is going to be discussed in depth (together with its physical applications) in Chapter 6. We shall only show that the above principal value integral induces a new distribution which, symbolically, will be denoted
pv.!..w
In this notation, distribution g(w; 1) is given by equality 
g(w; 1)
1 11' IpV ;, = 2°(w) + 2 11: IW
(6)
which is often written in the form
1 1 . 0 = PV;IW+
IW
+ 11:8(w).
(7)
For a = 1, our general original function X(t)t a  1 degenerates to the Heaviside function X (t), so that (6) becomes the generalized Fourier transform of Heaviside function. In this case, its real part is equal to the Fourier transform of the even component of the Heaviside function, which is simply equal to 1/2, and the imaginary part is the Fourier transform of the odd component, which is sign(t)/2.
126
Chapter 4. Asymptotics ofFourier transforms
Utilizing the recurrence formulas (4), for any integer m, we can express the generalized Fourier transforms of functions X(t)t m  1 via the Fourier transform of Heaviside function
For m = 2, this formula gives the Fourier transform
_ g(cu;3)
=
I" . l d2 1 (cu) +zPV. 2 2 21l' dcu cu
of function x(t)t 2 discussed in Remark 4.4.6. Its real part is the generalized Fourier transform of infinitely differentiable function t 2 • It is identically equal to 0 for arbitrary cu "# 0, and , as was discussed before, does not have a power asymptotics for cu + 00.
The case a = 0; Fourier transform offunction X(t)/ta physical approach. The case of the Fourier transform of function X (t) / t, corresponding to the value a = 0, has to be considered separately from other generalized Fourier transforms of power functions. Copying formally the approach that was so successful in analyzing of generalized Fourier transforms for a > 0, we will try to determine the Fourier transform of X(t)/t as the weak limit (for h + 0) of the integral
f
00
g(cu; 0, h) =
exp(ht  icut)dt.
(8)
o
Just a passing glance at (8) permits an observation that, for any cu and any h, the integral on the righthand side is infinite, in view of the nonintegrable singularity at the lower limit. Nevertheless, neither mathematicians nor physicists throw up their hands in despair in such a situation. Mathematicians introduce new distributions that assign well defined values to integrals (8). Physicists quote additional physical arguments, which also give a finite answer. The ideas of mathematicians and physicists, although different in details and method of argumentation, are similar in essence. Various ways of computing integrals of type (8) can be grouped under a unifying umbrella of renormalization techniques. Most often, renormalization techniques are applied in quantum electrodynamics where the computation of physical quantities by the perturbation method leads to divergent integrals. Without reference to the physical processes that are described by function X (t) / t, we will use only very general renormalization ideas. A reasonably accurate measurement of the Fourier transform of a real physical process at frequency cu requires much longer time than the period T = 21l' / cu of the corresponding oscillation. For
127
4.5. Generalized Fourier transforms ofpower functions
that reason, the Fourier transform at frequency w = 0 is not an observable quantity as it would require for its measurement an infinitely long observation interval. Hence, it is natural to exclude it from considerations, reading its value off other values of the Fourier transform. For w = 0, the Fourier transform (8) 00
8(0; 0, h) =
/
exp( ht)dt
o
is infinite. Nevertheless, subtracting it from (8), we arrive at a renormalized Fourier transform 00
l(w; h)
/
=
1] dt,
[e iwt 
(9)
o which assumes finite values and correctly reflects the dependence of the Fourier transform 8 on frequency. Passing in (9) to a new variable of integration or = wt, 00
1/1 .  1] dor, /(w; h) = 21r ;eILT: [ e1T: o
where f,L = h / w is a dimensionless parameter. Differentiating both sides of the above equality with respect to that parameter, we obtain that _
d/ = d f,L
00
21T
/ [eT:(IL+i) _ eT:IL] dor = o
Utilizing the fact that if f,L +
00
21T
IL
S
[.!. __+
1_.] .
f,L
f,L
I
then 1 + 0, we can compute 1 from
 1/00[1.1]
/=
21T
S+l
_ _
ds=Re/+iIm/,
(10)
where (Ha)
1m /
=  41[1
1  21r arctg (f,L)
]
sign (f,L).
(Hb)
Chapter 4. Asymptotics of Fourier transfonns
128
The absolute value under the logarithm and function sign (It) in the imaginary part make these expression valid for negative It as well. Now, to find a generalized renormalized Fourier transform of function X (t)/t, it suffices to let h + 0 in expressions (10) and (11). At the beginning, we will do that for the imaginary part. Observing that in this case It + 0+ for w > 0, and It + 0 for w < 0, we get that 1m

41 sign (w).
f =
The real part of the renormalized Fourier transform X(t) / t requires a more thoughtful treatment. Recall that It = h/w, and observe that h in the first equality in (11) can not be taken to converge to 0 as its right hand side then becomes infinite. For that reason, we will impose a restriction w » h. Then It « 1, and we can utilize a simpler approximate expression Re

1
1
21r
1r
f = In(llti) = 2 lnOwi) + C,
1
C = In(h). 21r
The constant C above will be called the calibrating constant, as its value should be selected on the basis of comparison with results of the measurements and the choice of a frequency scale. Combining the last two formulas we arrive at a remarkable relation x(t)

t
1
t+ In(lwi)
21r
+C 
1 . isIgn (w).
4
(12)
The above "physical" approach to calculating the Fourier integral proved successful even in the case which, taken at its face value, diverged for any frequency.
The case a = 0; Fourier transform of 1/ ta mathematical approach. We shall now show how one can deal with this situation from the mathematical viewpoint. For simplicity, we shall restrict ourselves to a calculation of the Fourier transform of an even function l/ltl. In this situation, one defines a new distribution 1
T = PVit!, directly as a functional on test functions. Its continuity will be guaranteed by an exclusion of singularities in the corresponding integral. In the case under consideration, this functional is defined by the equality
f
1
T[4>] =
1
4>(t)  4>(0) dt +
Itl
f
Itl>1
4>(t) dt.
It I
(13)
129
4.5. Generalized Fourier transforms ofpower functions
Such regularization of divergent integrals is justified, from the viewpoint of mathematicians, by the fact that for any test function t/J(t) E S with t/J(O) = 0, the value of the function T[t/J] (13) coincides with the original integral t/J(t)/Itl dt. Let T(w) be the Fourier transform of our distribution. By Parseval formula (3.2.6) (for 1" = 0), it has to satisfy equality
J
T[4J] =
21l'
/
PV.!:..t/J(t)dt It I
=
[/1
21l'
t/J(t)  t/J(O) dt It I
1
+/ Itl>1
t/J(t) dt] . It I
If we transform the integrals in the brackets by expressing test function
1
in terms of its own Fourier transform, and change the order of integration, we get that
[i COS("';) 
f
=
Passing to the new variable of integration 1"
T[4J] =
/
1dt
008:"") dtJ dw.
= wt in the inner integrals gives
iP(w) [Cin(w)
where
+
+ Ci(w)] dw,
00
Ci(z) =  /
co:(S) ds
z is a special function called the integral cosine, and
· )= Cm(z
l
Z
o
1  cos(s)d s s
is another related special function. Neither of them, separately, can be expressed in terms of elementary functions, but their sum, up to a constant, is equal to the logarithmic function. Indeed, if we write function Cin of a real argument w in the form w
Cin(w) = /
o
1
1 ;oss ds = / 1 ;ass ds 0
w
+/
1
1 ;oss ds,
Chapter 4. Asymptotics of Fourier transforms
130
and the last integral above in the form
f
w
=f
f
w
1 coss S
1
ds
w
1 ;ds 
1
= In(lwl) 
1
then we see that Cin(w) where constant
f
1
y=
coss . sds  Cl(W),
1
= 1n(lwi) + y 
1 coss d
s
o
f
00
coss sds
f
00
s
1
Ci(w),
coss  d s. s
One can prove that constant y coincides with the Euler constant (4.1.5). Hence Cin(w)
+ Ci(w) = 1n(lwl) + y,
which gives the following addition to our tables of generalized Fourier transforms:
4.6 Discontinuities of the second kind Let us return to one of the main topics of this chapter: analysis of the asymptotics of Fourier transforms of functions with discontinuities of the second kind through a study of the gauge functions (4.4.9) and (4.4.15). Consider the integral
f
00
f(t)eiwtdt,
o where, for t > 0, f(t) is a sufficiently smooth function, and for t singularity of the type f(t) '" Ata  1 , a >0.
0, it has a
To make the situation more concrete assume that f(t) == 0 for t < 0, so that tQe above integral is equal, up to a factor of 1j2rr, to the Fourier transform of function f(t). In the previous two sections, a detailed analysis of gauge functions (4.4.9) and (4.4.15) gave us some hints that the Fourier transform j(w) of function f(t)
4.6. Discontinuities of the second kind
131
has, for w + 00, principal asymptotics of the order l/wa. This fact has yet to be proved rigorously, and explicit necessary conditions on function f(t) have to be spelled out. We will adopt a "patchingup" method which relies on removing singularities from the original function by superposing on it a "patchingup" function with the same singularity. More exactly, we will consider an auxiliary function
v(t)
= f(t) 
Ag(t; a, h).
(1)
Since g(t; a, h) '" t a  1 (t + 0), function v(t) is of order smaller than f(t), that is, v(t) = o{t a  1 } (t + 0). Hence, it is natural to expect asymptotics of its Fourier transform to be of order smaller than the expected principal asymptotics of the Fourier transform j(w), that is, v(w) = o{1/wa } (w + (0). If this is indeed the case, then the principal asymptotics of function
f(t) = Ag(t; a, h)
+ v(t)
coincides with the main asymptotics (4.4.17) of the gauge (or patchingup) function g(t; a, h) multiplied by A, and the Fourier transform of function f(t) has asymptotics
 w  A r(a)
f( ) 
21f(iw)a
+ 0 { wa
(w + (0).
(2)
The rigorous proof is based on the following modification of the more general RiemannLebesgue Lemma:
Let function v (t) be continuous for t > 0 and absolutely integrable on the infinite interval (0, (0). Furthermore, assume the same is true for all the derivatives of v(t) of order up to n = La + IJ a > O. Additionally, assume that, for t + 0, function v(t), and all its derivatives up to order n, satisfy the asymptotic relations v(m) = o{t a  m 1 } Then the integral
f
(t + 0).
(3)
00
F(w) =
v(t)e iwt dt
(4)
o
has the asymptotics (5)
Chapter 4. Asymptotics of Fourier transforms
132
First, let us sketch a proof of this modification. Asymptotic relations (3) are equivalent to the following statement: For any e > 0 we can find a K > 0 such that, for t < K, m = 1,2, ... ,n. (6) We shall apply these inequalities to the integral (4). Select a> > 11K, and split the interval of integration into subintervals (0, 1/a» and (1/a>, (0). For the integral over the first subinterval, . I I1o1/fJJ v(t)e'fJJtdt
11/fJJ 0
Iv(t)ldt
e
11/fJJ a 1 e 1 0 t  dt = ; wa'
which means that this piece of integral (4) is, as 1a>1
00,
(7)
of order smaller than
1/a>a.
The second piece of the integral (4) over interval (1/a>, (0) will be transformed by repeated integration by parts. The first integration by parts gives that
1
. e i ( 1 ) v(t)e1fJJtdt = .v l/fJJ la> a> 00
1 + ;
1
00
la> l/fJJ
. v'(t)e1fJJtdt.
In view of (6), the first summand on the righthand side is ela>a and, as a result, it is of order smaller than 1/a>a for 1a>1 00. Repeating the integration by parts another n  1 times, and making a similar observation to the effect that the nonintegral terms are o{1/a>a}, we are lead in the end to an asymptotic equality
1
00
l/fJJ
. dt = v(t)e 1fJJt
(
1 ;la>
)n 1 v(n) (t)e1fJJt dt + 00
.
l/fJJ
0 { 1
wa
}.
Let us split the integral on the righthand side
In view of (6), the contribution of the first summand
1 ( ;la>
)n 1
1//( v(n)(t)eifJJtdt
l/fJJ
I
1 ron
_IE
1//( t a  n 1dt
l/fJJ
So (00
J1/fJJ
v (t)e ifJJt dt =
(00
la>
A//(
V(n) (t)eifJJt dt
<  IE . 1.
+0
n
CJt
a>a
{2..} . wa
4.6. Discontinuities of the second kind
133
It follows from the RiemannLebesgue Lemma of Section 4.2 that the remaining integral, which contains a continuous and absolutely integrable function v(n)(t), converges to 0 as Iwl 00. This, together with (7), proves the validity of asymptotic relation (5). • Now, the proof of (2) immediately follows from the above modification of the RiemannLebesgue Lemma, since the auxiliary function v(t) in (1) satisfies conditions of the lemma. Example 1. Consider the asymptoties of Fourier transform of function I(t)
= exp(Itl a ),
cx >
o.
(8)
In contrast to functions considered above, for t = 0, it assumes a nonzero finite value 1(0) = 1. As a result, the auxiliary integral
has the principal asymptoties of the order O(l/w) which is absent in the actual asymptoties of the full Fourier transform of function (8). To exclude this asymptoties and to find the asymptotic behavior of the full Fourier transform, first consider the derivative of function (8) I'(t)
= _cxltl a  1 exp(Itla) sign (t).
(9)
Let us form an auxiliary function equal to 0 for t < 0 and, for t > 0, given by equality v(t) = I' (t) + cxg(t; cx, h). It is easy to see that it satisfies all the requirements of the lemma proved above.
Hence, the principal asymptoties, for Iwl 00, of the Fourier transform of the "onesided" derivative I'(t)x(t) is described by formula (2) with A = cx. The actual derivative (9) of the original function (8) is an odd function of t. Consequently, the asymptotics of its Fourier transform is twice the imaginary part of the asymptoties (2), that is . r(cx) sin(1l'cx/2) . 2ICX 2 Sign (w).
1l'lwl a
Finally, the asymptotics of Fourier transform of the original function (8) can be found by dividing the above expression by iw, thus obtaining
Chapter 4. Asymptotics of Fourier transform
134
r(a + 1) sin(1l"a/2) j (w) '" ,,
1l"lwl a +1
(10)
.
• 4.7 Exercises 1. Find the main power asymptotics (as x + O)offunction/(x) = (sin2x2sinx)/x. 2. Investigate asymptotic behavior, as x + 0, and as x + 00, of function / (x) = (x  tanhx)/x 2 • 3. Investigate asymptotic behavior, as x + 0, of function / (x) = (1  x cot x) Ix.
°
4. Utilizing answer to the above exercise find asymptotic behavior, as a + 00, of the root of the transcendental equation a  x cot x, < x < 1l" 12.
=
5. Investigate asymptotic behavior, as N + 00, of the expression /(N) ani N), where {an} is a bounded sequence.
= n:=1 (1 +
6. Find the asymptotics (as N + 00) of expression n:=1 (1 + where = tiN, and ql(1') is a function integrable in the interval l' e (0, t). Provide an upper estimate of the remainder term in the obtained asymptotic formula. 7. Determine the character of convergence to 0, as w + 00, of the following two integrals:
S(w) = rooSinwtSin(
10
atp4 )dt,C(W)= 1+ t
10roo
1 + pt
8. Find the principal asymptotics, for w + 00, of the Fourier transform of function exp( altI 3 ), and evaluate the infinitesimal order of the remainder term.
Jet) =
9. Assume that
Jet)
is an even, infinitely differentiable function on the interval t e
°
(1', 1'), which is identically equal to zero outside this interval. Furthermore, suppose that there exists limit limH .. o /(t)/(1'  t)n = A, where A > and n is positive integer. Find the principal asymptotics of the Fourier transform of this function for w + 00.
10. Find the principal asymptotics, for w + 00, of the integral
J(w) =
1 1
1
In
(2+t2) coswtdt. 1 +2t 2
11. Find the principal asymptotics, for w + 00, of the Fourier transform of function 4.7.1.
Jet), with the graph shown in Fig.
135
4.7. Exercises
fit)
FIGURE 4.7.1 12. What is the principal asymptotics (as (J) /( )
t =e
_r2 {
....
I,
l(t/(l+t»tI,
00) of the Fourier transform of function
t < 0; t>O,
(0 <
P<
1).
Its graph, reminiscent of the profile of an ocean wave (or a sand dune), appears in Fig. 4.7.2.
fit)
FIGURE 4.7.2 13. Find the principal asymptotics, for (J) .... 00, of the Fourier transform of the semi circle function /(t), equal to zero outside the interval t E (1,1), and equal to inside that interval.
ChapterS Stationary Phase and Related Methods
In this chapter we will use methods developed in Chapter 4 to provide a general scheme for finding asymptotics. The remarkable Kelvin's method of stationary phase will be employed as well.
5.1
Finding asymptotics: a general scheme
Consider the integral 1= I(x) =
foB I(t)exp[ixp(t)]dt,
(1)
where 1 (t) is a continuously differentiable function on the interval [0, B) and such that 1(0) = 1 =F O. The function p(t) appearing in the exponent will be assumed twice differentiable on (0, B) and monotonically increasing with p'(t) > O. To apply to (1) the standard methods of asymptotic analysis described in Chapter 4, we will change the variable of integration to s
= p(t) 
p(O).
Denote the monotonically increasing inverse of the above function by t = q(s)
(q(Q) = B).
After this change of variables, the integral (1) assumes familiar form of the Fourier integral 1= exp[ixp(O)]
10fQ F(s)q'(s)e .
© Springer Nature Switzerland AG 2018 A. I. Saichev and W. Woyczynski, Distributions in the Physical and Engineering Sciences, Volume 1, Applied and Numerical Harmonic Analysis, https://doi.org/10.1007/9783319979588_5
IXS
ds,
Chapter 5. Stationary phase and related methods
138
where F (s) = f (q (s)) is continuously differentiable on the interval [0, Q] function such that F (0) = f =1= O. Let us investigate the asymptotic behavior of I as x 00. First, observe that the factor in front of the integral has modulus 1 and has no effect on asymptotics. So, from now on, we will omit it (putting p(O) = 0) and Q consider 1= F(s)q'(s)e ixs ds. (2)
fo
The function q' (s) is also continuously differentiable on the interval (0, Q). If, in addition, it had a finite limit for s 0, then a further study of the integral I would repeat the previously developed asymptotic analysis of Fourier images of functions with discontinuities of the first kind. More interesting, and physically more important, is the case of functions p(t) which have the asymptotics (3)
0)
(t
for a certain ex > 0, which we will consider in some detail. In this situation q'(s) '" GsfJ l
where
(s
(4)
0),
1
f3 =. ex
The above asymptotics (s 0) of the integrand and experience gained in the previous chapter suggest the following asymptotics for the integral (2): 1 )fJ = f I'" r(f3)fG ( :IX
1 (1) (1.
r ex
ex
)l/a
P,X
(Ixl
00).
(5)
On the other hand, if F(Q) = f(B) =1= 0, then the behavior of the integrand close to the upper limit gives the asymptotics (Ix I
I'" f(B)q'(Q)fix
00),
which, for ex > 1, is of order smaller than (5). Thus, the main term of the asymptotic expansion of integral (2) is given by formula (5). Reinserting the factor omitted earlier, and returning to the notation of the original integral (1), we can finally write 00, that, for any ex > 1, as Ix I 1 ) iA(Xl f(t)exp[ixp(t)]dt",r ( ;;+1
(
1 ) l/a Pix f(A)exp[ixp(A)]. (6)
139
5.1. Finding asymptotics: a general scheme
The upper limit in (6) was deliberately set to be infinite to avoid distraction caused by smaller order asymptotics generated by a finite upper limit. The lower limit was kept arbitrary for the sake of generality. Observe that the assumption that function p(t) be strictly increasing is not necessary for the above result. A strict monotonicity suffices as the two cases can be transformed into each other by a nonessential replacement of i into i. Let us indicate a number of consequences of formula (6) that are important in applications: (1) If function p(t) is symmetric in the neighborhood of A, then the asymptotics
p(t) '" Pit  Ala,
(7)
for an a > 0, implies the doubled asymptotics (6) for the integral with infinite limits:
f
f(t)exp[ixp(t»)dt '"
2r
1)
f(A)exp[ixp(A»). (8)
(2) If p(t) is antisymmetric in the neighborhood of A, then only the real part is preserved in the asymptotics and
f
f(t) exp[ixp(t»)dt '"
+
r
1)
Re
I/a f(A) exp[ixp(A»).
(9) (3) In terms of distribution theory the equality (8) means that, for any a > 1, the family of functions of variable t
exp[ixp(t)It 
Ala] /
2r
1)
x eR,
I/a ,
weakly converges to the Dirac delta  A) as x 00. The particular case a = 2 corresponds to the familiar function (1.3.5). (4) If f(t) is constant and p(t) = P(t  A)a, then the asymptotic relation (6) becomes an exact equality. Specifying f = 1, A = 0, and x = 1, and separating the real and imaginary parts we arrive at the following standard integral formulas valid for a > 1 and P > 0:
10tx) cos(Pta ) dt =
r
(1 + ) (1 )l/a 1
10rOO sin(Pta)dt = r ( 1 + 1)
P
cos
( 1 ) I/a P sin
(2a) ,
(lOa)
(2a).
(lOb)
7r
7r
Chapter 5. Stationary phase and related methods
140
5.2 Stationary phase method Assumptions (5.1.3) and (5.1.7) which secured the above asymptotics may seem artificially chosen, just to make mathematics rigorous. Actually, many of them, and in particular the case a = 2, emerge perfectly naturally in the physical phenomena. Consider the integral
LB f(t)exp[ixp(t)]dt.
(1)
where p (t) is an arbitrary function twice differentiable on the interval of integration. Function f(t) will be assumed continuously differentiable. It turns out that in this fairly general situation the asymptotics of (1) corresponds to the special case a = 2. We shall begin the asymptotic (x 00) analysis of (1) by finding stationary points of p(t) where p'(t) = O.
Denote the roots of this equation by rm. m = 1, ...• N < 00. Assume that all of them correspond to simple extrema of function p(t) with p" (rm) =f. O. In their neighborhood, p(t) has a parabolic behavior
(2) Consequently, the integral (1) has automatically the symmetric asymptotics (5.1.7) with
a = 2,
P
= p"(rm )/2.
Let us partition the interval of integration into disjoint intervals, each containing just one of the stationary points. In our case a = 2 > 1, and the asymptotics contributed by the boundary points of the intervals are of order smaller than the asymptotics generated by the stationary points (5.1.8). For this reason, in final formulas only the latter appear. Summing contributions of all the stationary points we arrive at the asymptotic (x 00) formula
L B
N
f(t)eXp[iXP(t)]dt,...,'?; f(rm)
21r
. "( ) exp[ixp(rm)]. IXP rm
(3)
If some stationary points coincide with the endpoints of the interval of integration then the corresponding summands will appear with coefficient 1/2.
5.3. Fresnel approximation
5.3
141
Fresnel approximation
Let us take a look at the asymptotics (5.2.3) from a slightly different viewpoint and focus our attention on the case where there is only one stationary point T and the limits of integration are infinite. Then (5.2.3) reduces to the asymptotic equality
f
f(t) exp[ixp(t)] dt "" f(T)
.
2rr "
IXP (T)
exp[ixp(T)],
(x + (0).
(1)
We shall attempt to find the "hidden springs" of the stationary phase method by analyzing this example in some depth. The totally rigorous mathematical derivation of (1) seems bland and incomplete to physicists and engineers if it is given without that extra insight that comes from a perhaps imprecise but revealing heuristic arguments. Actually, mathematicians also often gain a deeper understanding of their subject by accumulating a store of sometimes imprecise analogies acquired in "reallife" experiences and physical "thought" experiments. The stationary phase method can also be elucidated by such "reallife" arguments: the fast oscillation (2) exp[ixp(t)] has current (timedependent) frequency w = xp' (t) and period T = 2rr /Iwl which decays like 1/x. If f(t) has a characteristic scale a then, in the domain of integration where T « a, the adjacent crests and troughs of the integrated process compensate, the better the bigger x, and only close to the stationary point t = T, where W(T) = 0, does that compensation becomes less effective. As a result, a small and shrinking with the growth of x neighborhood of point T gives the main contribution to the integral. In the neighborhood of the simple stationary point T, function p(t) is well approximated by the parabola p(t) = p(T)
+ r(t 
r = p"(T).
T)2/2,
Outside that small neighborhood, function (2) can be replaced by exp [ ip(T)  ixr
(t  T)2]
2
'
without changing the value of the integral significantly, since both functions oscillate quickly and give a small contribution to the integral. For this reason, the original integral (1) can be replaced by an asymptotically equivalent expression
. f
exp[lp(T)]
[.
f(t) exp lxr (tT)2] 2 dt.
(3)
142
Chapter 5. Stationary phase and related methods
Furthermore, observe that for x 00, in the neighborhood of the stationary point essential for the integral, function /(t) practically coincides with constant /(1:). The latter can be taken outside the integral sign and the remaining integral can be calculated with the help of the standard formula
Remark 1. Physicists often stop short of the asymptotic relation (1) and operate with the integral (3). Analogously with optics, where such integrals appear in the socalled "Fresnel approximation", the approximation of the integral on the lefthand side of (1) by the integral (3) will be also called the Fresnel approximation. Later on, discussing optics applications, we shall show that the asymptotic formulas (5.2.3) and (1) correspond to the crude geometric optics approximation. Remark 2. The Fresnel approximations of integrals of type (3) are closely related to the Fresnel sine and cosine integral special functions (4a)
(4b)
often encountered in wave problems. They are both odd functions of z, with limits C(oo) = S(oo) = 1/2.
(5)
The graphs of Fresnel integral functions are shown in Fig. 5.3.1.
5.4 Accuracy of the stationary phase approximation To enhance the practical value of the mainasymptotics formulas (5.2.3) and (5.3.1) for integrals of rapidly oscillating functions we need estimates of the remainder terms. Their magnitude strongly depends on functions /(t) and p(t). Moreover, no universal method of finding precise estimates is known. For that reason we will analyze the accuracy of asymptotics (5.3.1) in just two generic cases, hoping that detailed analysis of a few concrete examples will better illuminate the essence of the problem than plowing through a laborious general argument.
5.4. Accuracy of the stationary phase approximation
143
0.6 0.5 0.4 0.3 0.2
0.1 3
4
5
FIGURE 5.3.1 Graphs of the Fresnel integral functions C(z) and S(z).
Example 1. First, let us find the magnitude of error created by replacing the integral (5.3.3) by the righthand side of (5.3.1) in the case when /(t) is a continuously differentiable (smooth) Gaussian function. Then (5.3.3) becomes the standard integral
f [
2
ex  t2 i x r t ] dt= p 20 2 2
Iff.J1 
ixr
1 (ija 2xr)
.
(1)
For simplicity's sake put 'l' = 0; consideration of the more complex general case contributes little to the understanding of the essence of the situation. First factor on the righthand side of (1) corresponds to the main asymptotics (5.3.1), and the second describes the deviation from it. Thus, the relative error is of the order (obtained via formula (5.3.1»
1 I" "2012xr, I1 .J1  (ija 2xr)
(x + 00).
(2)
For a general smooth function / (t) the accuracy offormulas like (5.3.1) can be estimated by replacing a with the characteristic scale of function / (t )an admittedly nonrigorous but heuristically useful approach. •
144
Chapter 5. Stationary phase and related methods
Example 2. Consider the jump function /(t) = x(a Itl). Then the calculation is reduced to evaluation of the integral
1=
i:
exp [ixr
2/f.(C(z) 
dt =
is(z)),
(3)
where Z=
.
7r
To investigate the asymptotic behavior of the Fresnel integrals that appear in (3) we shall write the first of them in the form C(z)
1
= 2
where c(Z)
= C(oo) 
C(z)
=
c(z),
100
cos(7rt 2 /2)dt.
Changing to a new variable of integration y = t 2 , we get that c(z) =
1
2
[00 v'Y1 cos(7rY/2) dy, Z2
and integration by parts gives c(z) = 
7rZ
sin(7r Z2 /2)
+
27r
[00 y",1",Y sin(7rY/2) dy. Z2
Another integration by parts shows that the remaining integralis 0(1/z) as z + 00. Therefore, C(z) satisfies the asymptotic relation
+ o(l/z),
(z + 00).
1 = 21  7rZ cos(7rz 2 /2) + o(l/z),
(z + 00).
C(z) =
2
+
7rZ
sin(7rz 2/2)
An analogous relation is valid for S(z): S(z)
5.5. Method of steepest descent
145
Substituting these asymptotics into (3) we get that
1=
f{f .7r lxr
(
. 2 /2] + o(l/z) ) , 1 + .j2[ exp[l'lrXZ 7rZ
(Z * (0).
Now it is clear thatthe relative error of replacing I by the righthand side of (5.3.1) is
Note that it is of the order'" 1/,.fX and not'" 1 / x, the latter being the case for smooth functions f(t) (see (2». •
5.5
Method of steepest descent
Formula (5.2.3) contains the main ingredient of the stationary phase method. It is related to the steepest descent method (or Laplace's method) which is applicable to purely real integrals
i
B
f(t) exp[ xp(t)] dt.
(1)
Without repeating considerations that led us to (5.2.3) we shall give the final formula for the main asymptotics of integral (1). Let p(t) be a sufficiently smooth function which has only simple minima at points'l"m located inside the interval (A, B). It turns out that just rewriting (5.2.3) without the imaginary unit i gives the correct asymptotics
l
B
A
f(t)exp[xp(t)]dt'"
L f('l"m) m=l N
27r
" ) exp[xp('l"m)], xp ('l"m
(Z * (0).
(2)
Despite its superficial similarity, the above formula differs from formula (5.2.3) in an essential way. First of all, the summation in it is not over all extrema but just over all minima of the exponent. Secondly, and this is the main point here, for different values of the minima, the exponential factors in the sum have different magnitudes, and these differences increase with the growth of x. For that reason the main asymptotics of integral (2) contains only one term corresponding to the absolute minimum of function p(t) in the interval of integration
146
l
Chapter 5. Stationary phase and related methods B
A
f(t) exp[xp(t)]dt
'V
 , ,  exp[xp(r)], xp (r)
f(r)
(x + 00),
(2)
where p(t) ::: p(r), for all t E [A, B].
5.6 Exercises 1. The stationary phase method is often useful in problems of wave propagation. In particular, in Chapter 9 we will encounter the integral G(p) =  1
1
00
21ro
exp[ikpcosh(t)]dt.
which describes complex amplitude of a cylindrical wave. Analyze this integral using the stationary phase method. 2. The real Anger function J,Az) and Weber function E,Az) (here, we assume that w and z are real) are uniquely determined by the equation D",(z)
= J",(z) 
iE",(z)
=.!.. 11'
Jor exp[i(erxP 
z sin 1/1)] dl/l
(1)
and are often encountered in problems of mathematical physics. Find the main asymptotics of Anger and Weber functions for a fixed z and w 00. 3. Find the main asymptotics of Anger and Weber functions (see Exercise 2) for fixed w and
z
00.
4. Find the Fresnel approximation of function D",(z) for z » 1. 5. In Exercises 3 and 4 we have explored the asymptotic behavior of Anger and Weber functions along the w and z axes ofthe (w. z)plane. What happens in the rest ofthe plane? More precisely, study the asymptotic behavior of Anger and Weber functions along the raysz = pwforO < p < 1. 6. Study the asymptotics of (1) for p = 1. w 00. Hint: As p 1 the expression on the righthand side of the relation (5.5) in the Answers and Solutions chapter diverges to infinity. This fact indicates that for p = 1 the asymptotics is of a different order.
=
(1/11') 7. Complete investigation of the integral D",(pw) 1/1  p sin 1/1, by checking its asymptotic behavior for w
p(1/1) =
J; e;"'P(q,)dl/l. where 00
if p > 1.
Remarks on Exercises 27: The integral in Exercise 7 clearly has two different types of asymptotic behavior in the (z. w)plane: one in the octant 0 < z < w (p < 1) and another in the octant 0 < w < z (p > 1). Moreover, its asymptotic behavior on the boundary z = w (p = 1) of these octants is qualitatively different from its asymptotic behavior in
147
5.6. Exercises
either of them. The total picture can be summarized as follows: the integral in Exercise 7 obeys the asymptotic power law of order a (p) where
I,
a(p)
= { 1/3, 1/2,
if 0 :::: p < 1; if p = 1; if 1 < p.
(2)
The function a (p) which determines asymptotics of the integral in Exercise 7 has jumps as we move from one (z, w)region to another. An infinitesimal change of parameter p can cause a major change in the asymptotic behavior of that integral. Such phenomena are called phase transitions and they correspond to the physical phase transitions like melting, evaporation or crystallization, where small changes in temperature (or other physical parameters) can cause large and sudden changes in the physical properties of matter. At first sight, the phase transition for the integral in Exercise 7 in the vicinity of the critical point p = 1 could be puzzling. The next exercise provides and additional insight into why it occurs.
1
8. Consider the integral
1 21r
00
/(w)=
It follows from (4.3.3) that, for any
0
t"
I(w) '"
eit»tdt
,
t"
o.
(3)
> 0,
1 . !' 21fIW", t"
(4)
(w + 00)
i.e., (3) obeys the asymptotic power decay law of the order a = 1. On the other hand, it follows from (4.6.2) that, for t" = 0, T(w)=
1
(5)
2",1fiw
i.e., (3) obeys the asymptotic power decay law of the order a transition" in the asymptotic behavior of (3) as t" + o.
= 1/2.
Study the "phase
9. Utilize the method of steepest descent sketched in Section 5.5 to derive the Stirling's approximate formula for the factorial:
n! = n . (n  1) ..... 2 . 1 '" ../21fn n" e"
(n + 00).
10. Consider the Riemann equation
av at
+ v av
ax
= 0,
v(x, t
= 0) = vo(x),
(6)
where Vo (x) is an infinitely differentiable and absolutely integrable function whose derivative attains its minimum at a certain point z, i.e.,
inf
oo Rand
'PV
f
f
R
ds =
'PV
ds.
R
Representing function qJ(x) in the form qJ(X) = [qJ(x)  qJ(O)]
+ qJ(O),
and noticing that the principal value of the integral containing the constant term qJ(O) is zero, we get that
'PV
f
f
R
ds =
qJ(s)
qJ(O) ds.
(2)
R Now the integral on the righthand side can be understood in the ordinary sense since its singularity at 0 has been removed. In other words, the integrand has been regularized. This immediately follows from the fact that any test function qJ(s) e V satisfies the Lipschitz condition IqJ(a)  qJ(b) I < Kia  bl.
(3)
If function qJ(s) does not have compact support, then one can still use equality
'PV
f
f
R
qJ(s) ds =
s
lim
Roo
qJ(s)  qJ(O) ds,
s
R
instead of (2) to define the principal value of integrals with infinite limits. Finally, let us point out another obvious but useful representation of the principalvalue integral (1):
pvf
qJ(s) ds =
s
(00
10
qJ(s)  qJ(s) ds.
s
6.1. Principal value distribution
151
Standard operations on "wellbehaved" convergent integrals, such as the change of variables, differentiation with respect to a parameter, etc., can also be used in analysis of principalvalue integrals. The following example illustrates the situation.
Example 1. Consider the singular integral I(x, b) = 'PV
f
exp(b2s2) ds, sx
(4)
and observe that a change of variables 'l' = bs transforms it into the integral
I(x, b) = I(bx), where
I(x) = 'Pvf exp('l'2) d'l'. 'l'X
Introducing variable of integration y = 'l'  x, we obtain that
I(x) = exp(x 2)J(x), where
= 'PV
J(x)
f
exp(y:  2xy) dy.
Differentiation of the above expression with respect to parameter x gives
J' (x)
= 2
f
exp( _y2  2yx)dy
= 2,Jii exp(x 2).
Taking into account the fact that J (0) = 0, one can get that
I(x) = 2,JiiD(x), where
D(x) = exp(x 2) foX exp(l)dy
is the socalled Dawson integral. Thus, finally, we arrive at the formula
'PV
f
exp(b2s2) sx
ds = 2,JiiD(bx).
Let us remark that for Ixl + 0, we have that D(x) '" x, and that for Ixl + Dawson integral has the asymptotics D(x) '" 1/2x.
00
the
152
Chapter 6. Singular integrals and fractal calculus
D(x)
1
2
3
5 X
4
FIGURE 6.1.1 Graph of the Dawson integral.
6.2 Principal value of Cauchy integral Another approach to evaluation of integral (6.0.1) depends on its interpretation as the limit rp(s) ds = lim rp(s) ds, (1)
f
s x
where
y+o
f
S 
Z
z = x +iy
is a complex parameter. Moving into the complex plane vicinity of the real axis removes the singularity. Let us find the above limit by separating explicitly the real and the imaginary parts of the expression
1 s z
iy
s x
= (s  x)2
+ y2 +
(s  x)2
+ y2·
Substituting this sum into (1) and noticing that the integral of the real part converges to the principal value as y + 0, we obtain that lim
y+O
f
rp(s) ds s Z
= pVf
rp(s) ds s Z
+i lim
y+O
f
Y2 (s  x)
+y
2rp(s)ds.
(2)
153
6.3. A study of monochromatic wave
Notice that the factor in front of function q>(s) in the last integral coincides, up to rr, with the familiar Lorentz curve (1.3.3)
1
Y
; (x  s)2
+ y2 '
which weakly converges to 8(x  s) as y 0+. Hence, the evaluation of the integral (6.0.1) using the limit procedure (1) leads to an identity
f
q>(s)   d s = PV sx
f
q>(s) .   d s ± l7rq>(x). sx
The plus sign corresponds to the limit y 0+ with y's restricted to the upper halfplane and the minusto the limit y 0 with y's restricted to the lower halfplane. Thus equality (2) determines two distributions: 1 1    .  = PV
and
+ irr8(x 
s)
(3)
1 1    .  = PV  irr8(x  s).
(4)
s
X
s x
,0
+,0
S 
x
sx
Although assigning complex values to realvalued integrals may seem strange at the first glance, these formulas often give the correct physical answer. The point is that their imaginary parts reflect the causality principle which was not spelled out explicitly when the original physical problem was posed but which, as we will see later on, plays an important role. Obvious physical arguments then permit us to indicate which of the formulas (34) exactly corresponds to the physical problem under consideration.
6.3 A study of monochromatic wave Formulas of Section 2 give, for example, the correct physical answer in the problem of radiation by a monochromatic wave source with complex amplitude. For simplicity, we will discuss only the ID case. Then, wave radiation is described by the nonhomogeneous wave equation
Chapter 6. Singular integrals and fractal calculus
154
If the wave source is monochromatic, that is D(x, t)
= w(x)coswt = Rew(x)eiwt ,
then the radiated wave is also monochromatic and can be written in the form E(x,t) = Reu(x)i wt ,
where the complex amplitude u (x) of the propagating wave satisfies the Helmholtz equation (1)
Here, k = wlc is the socalled wavenumber. The equation can be solved with the help of Green's functionan approach that will be discussed later on. Here, we will solve it by passing to the frequency domain and considering the Fourier transform 1 U(K) = 2Jr u(x)e1Kxdx.
1 .
The inverse Fourier transform is then given by u(x) =
1
u(x)eiKXdx.
Taking the Fourier transforms of both sides of equation (1) we get U(K) _ W(K) _ W(K)  k 2  K2 2k
[_1+___1_] K
k
K k
'
where W(K) is the Fourier transform ofthe source function w(x). We will assume that function W(K) is sufficiently smooth and rapidly decaying to 0 for IKI 00. The sought complex amplitude of the propagating wave can now be found by the inverse Fourier transform:
u () x = 1
2k
[I
W(K)  eiKXd K K+k

1
W(K)  eiKXd K] . Kk
Changing the variable of integration K in the first integral to K and observing that W(K) = W*(K), where the asterisk denotes the complex conjugate, we obtain
that
11
u(x) = k
Re[w(K)eiKX ] dK. Kk
(2)
155
6.3. A study of monochromatic wave
First of all, let us calculate the principal value
1 [ PV PV[U(X)] = Re k
f
W(K). _ _ e'KxdK ] Kk
of that integral by splitting it into the sum of two components pVf W(K) eiKXdK Kk =eikx[pV f
f
(3)
and studying the asymptotic behavior of each of them separately as Notice that, according to (3.3.7), sin[(K  k)x] Kk
weakly as
Ixl
00.
00.
• k)' () JrO(K SIgn X ,
Consequently,
lim PV
x+oo
Ixl
f
W(K)
sin[(K  k)x] . dK = Jrw(k) sIgn (x). K  k
Now consider the first integral on the righthand side of (3). Adding and subtracting number w(k) from function W(K) we arrive at the equality pVf W(K) COS[(K  k)x]dK Kk
=
where
f
1{I(K, k) COS[(K  k)x]dK
 f
+ w(k)PV
k)x] dK, Kk
COS[(K 
W(K)  w(k) k 1{I(K, ) = k' K
By previous assumptions, 1{1 (K, k) is a continuous function of K. The second integral on the righthand side of the above equality vanishes because the integrand is odd. The first integral converges uniformly for Ix I > 0, and by the RiemannLebesgue Lemma, its value converges to 0 as x 00. Thus we arrive at the asymptotic formula
(Ixl
00),
Chapter 6. Singular integrals and fractal calculus
156
which permits us to drop the corresponding, converging to zero, term to get that 'PV[U(x)]""" IIm[w(k)e ikx sign (x)],
(Ixl
+ (0).
(4)
The above expression contradicts the radiation condition and is physically not acceptable. 1 To save the situation we will tum to formulas of Section 2 which give that 1T ·k (5) u(x) = 'PV[u(x)] =f i"kRe[w(k)e' X]. Formulas (4) and (5) imply that if we select the plus sign in the formulas of Section 2 then we shall arrive at the asymptotic formula u(x)
{
i ?!.w*(k)e ikx
.! (k) i k x ' '7Cw e ,
x > 0,
0 x < ,
which does satisfy the radiation condition. Here, the physically acceptable answer corresponds to the distribution
11. Kk+iO Kk
   = 'PV  l1T8(K  k).
(6)
Let us consider the physical arguments in favor of such a choice. There are no purely monochromatic wave sources in nature emitting radiation for infinitely long times. All realworld phenomena begin and end in finite time. The fact that the source was turned on sometime in the past can be taken into account by assuming that its time dependence reflects asymptotically negligible intensity of the source at time 00. That is exactly what the replacement of k = w/c in the preceding formulas by w  1. Y = k  I·0, k =(7) C
C
accomplishes and what justifies utilization of the distribution (6). Our choice was based on the causality principle which asserts that it is impossible to receive the wave before the wave source is turned on. The same result can be obtained if a principle of infinitesimal relaxation is applied. According to this principle, any medium (even the vacuum) damps waves. This means that, for example, the propagating to the right monochromatic wave exp(iwt  ikx) is attenuated with the growth of x. Again we are led to the conclusion that the real k in (2) should be replaced by the formula (7). INotice that the radiation condition demands that far from the source (i.e., for only exist waves running away from the source, like exp(i(wt  klxl».
Ixl
(0) there
157
6.4. The Cauchy fonnula
6.4 The Cauchy formula The results of Section 2 are closely related to the Cauchy formula from the theory of functions of a complex variable. The formula asserts that for any function f (z), analytic in a simply connected domain D in the complex plane C and continuous on its closure iJ including the boundary contour C, f( z) = _1_ "'_.
f
f(s)d s
,.
 z
C
,
where z is a point in the interior of D and contour C is oriented counterclockwise. If z is an arbitrary point of the complex plane then the above integral defines a new function F(z) =
2m
f
C
f(s)d s .
(1)
S z
For z inside the contour of integration F(z) = f(z),
zED.
If z ¢ iJ then the integrand is analytic everywhere in D and, by the Cauchy's Theorem, F(z) == O. For a boundary point z = so E C the integral (1) is singular and F(z) will be understood in the principalvalue sense. In the present context, this will mean that F(so)
. 1 = PV[F(so)] = r.O hm "'_. I
f
Cr
f(S)d
S
s,
so
(2)
where Cr is a curve obtained from contour C by removing its part contained in a disc of radius r + 0 and centered at so (see Fig. 6.4.1). Observe that, for any function f continuous on contour, function F(z) is well defined everywhere. Formula (2) generalizes the concept of the principal value of a singular integral on the real axis and we will study it for a function f(z) analytic inside contour C and continous on it. To that end substitute an identity f(s) = [f(s)  f(so)]
+ f(so),
158
Chapter 6. Singular integrals and fractal calculus
FIGURE 6.4.1 A schematic illustration of countour Cr.
into (2) and split the integral into two parts:
PV[F(so)] =
2m
f c
[f(s)  f(so))d s
S
so
+ lim
r ...... O
2m
s  so
We deliberately replaced Cr by C in the first integral since for f(z) satisfying the Lipschitz condition (6.1.3) the integral is no longer singular and one can integrate over the whole closed contour C. Since its integrand is analytic inside the contour and continous on the contour, the first integral vanishes by Cauchy's theorem. Thus
PV[F(so)] = lim
r ...... O
2m
f
c,
S
so
(3)
The last integral can be evaluated assuming that contour C is smooth in the vicinity of point so, which is called the regular point ofthe contour. Let us add and subtract from (3) the integral over portion Cr of the located in D circle of radius r with center at so. Since the integral over the closed contour Cr + Cr is equal to zero, equality (3) can be rewritten in the form
PV[F(so)) =  lim
r ...... O
2m
f c,
S
so
6.4. The Cauchy formula
159
The latter integral is easy to evaluate:
f
.!!L = 
f
o
. rel'fJ
= i1r.
11:
Cr
Hence, we obtain that
= and the principal value of the analytic function is equal to
'Pvf c
=

(4)
that is, the value of function f at the singular point of the integrand multiplied by i 1r. Thus, the behavior of the Cauchy integral in the neighborhood of a regular point of contour C can be summarized as follows:
f
.{
Z
= 2m
c
Z
0,
=
(5)
Z
where the plus sign corresponds to the limit value of the integral while is approached from the inside of contour C, and the minus sign corresponds to the approach from outside.
Example 1. Let us consider a contour C consisting of the interval [ R, R] on the real axis and the semicircle CR of radius R with the center at point Z = located in the upper halfplane. If f(z) is analytic in the upper halfplane and such that the integral over C R uniformly converges to zero as R _ 00, then the Cauchy formula is transformed into
°
f sz f(s)
.
y > 0,
ds = 2mf(z),
(6)
and equality (4) assumes the form 'PV
f sx  d s = l1rf(x). f(s)
.
(7)
Recall that, by the well known Jordan Lemma in the complex functions theory, function f(z) = exp(iAZ), (8)
Chapter 6. Singular integrals and fractal calculus
160
satisfies the conditions mentioned earlier so that the above formula implies that
PV
f
eiAs
.
  d s = 7rie'AX.
In particular, it follows that for x = 0
iPV
f
(9)
sx
eiAS
;ds
=
f
sin{}l.s) sds
= 7r.
(10)
Notice that the Cauchy formula (6) interpreted in the spirit ofthe distribution theory defines a new distribution All T(s  z) =  .   , 27r1 S  Z
(11)
which is called the analytic representation of the Dirac delta. Its functional action assigns to a function !(s), analytic in the upper halfplane and rapidly decaying at infinity, its value at the point z. Applied to any usual "wellbehaved" function of the real variable s, it defines a new function of complex variable z which is analytic everywhere with the possible exception of the real axis. Crossing the real axis at the point z = x the functional T(s  z)[fl has a jump of size !(x); this follows from formulas of Section 2. We will extract this jump by introducing a new distribution A
A
8(s  z) = T(s  z)
+T
A*
1 (s  z) = 7r
y (s  x)
2
+ y 2'
which is harmonic for y =f:. 0 and which converges to the usual Dirac delta 8 (s x) as y 0+. As we will see in Volume 2, the corresponding functional 8(s  z)[! (s)] solves the Dirichlet problem for the 2D Laplace's equation in the upper halfplane y>Q
•
6.S The Hilbert transform The integral Hilbert transform 1/I(t)
1 = PV 7r
f
rp(s) ds s t
(1)
161
6.5. The Hilbert transform
of function q;(t) is also defined in terms of the principal value of the integral involved. We shall find the inversion formula for the Hilbert transform by applying the Fourier transform to both sides of definition (1). The lefthand side is transformed into 1/I(w) = 1 1/I(t)e lwt dt, 211"
f
.
and the righthand side, after a change of the integration order and other simple manipulations, assumes the form
1
2
211"
Since PV
f
f
.[ f
q;(s)e IWS
eiw(st)
.
s t
dt
=i
PV
f
eiw(st)
s t
sin(wt)
dt t
]
dt ds.
= i1l" sign (w),
we get that t(w) = iip(w) sign (w).
(2)
Now let g(t) be the Hilbert transform of function 1/I(t). Then, according to the above formula, its Fourier transform is g(w) = it(w) sign (w) = ip(w)
so that g(t)
= q;(t).
Thus, the inverse Hilbert transform has the form
1 q;(t) = PV 11"
f
1/I(s) ds.
s
t
(3)
One of the most important applications of the Hilbert transform is related to the causality principle. We shall illustrate it in the example of an absolutely integrable function h(t) describing a response of a linear physical system to the Dirac delta impulse 8(t). The Fourier transform h(w) appears in many physical applications. Let us write it with the real frequency w replaced by a complex number)... = w+ia:
h()"') = 1
211"
1
00
0
h(t)e l')"tdt.
(4)
The zero lower limit takes into account the causality principle which requires that the system's response cannot appear before the action of the impulse: h (t < 0) == O. It is obvious from (4) that function h()"') is analytic in the lower halfplane and
Chapter 6. Singular integrals and fractal calculus
162
continuous on the real axis a = relation

o.
This, in tum, means that h(w) satisfies the
i7rh(w) =
'PV
f
h(K) dK, KW
similar to (6.4.7). The complex function h(w) can be represented as a sum of its real and imaginary parts: h(w) = fP(w) + iy,(w). Substituting them into the last equality and comparing separately the real and the imaginary parts, we discover that fP and y, are related to each other by the Hilbert transform. The connecting formulas (1) and (3) applied to the real and imaginary parts of the Fourier transform of response function are called in physics the dispersion relations and are widely used in the theory of wave propagation in dispersing media.
6.6 Analytic signals Another physical application of the Hilbert transform is related to the notion of analytic signal, which appears in various areas of physical sciences, from electrical engineering to quantum optics. It can be introduced in the following way. To each real process g(t), which will be assumed to be absolutely integrable, we will assign a complex signal with the real part equal to g(t) and the imaginary part TJ(t) defined by the condition that is an analytic function of the complex variable z = t + ifJ in the upper halfplane. The analyticity can be achieved by making the Fourier transform of the signal g(t) vanish for negative frequencies, that is, by replacing (w) by 2 X (w ). The last expression can be also written in the form
=
+
=
sign (w).
(1)
The first component on the righthand side is the Fourier transform of the original signal and the second is the Fourier transform of the imaginary part of the analytic signal Formula (6.5.2), relating the Fourier transforms of a function and its Hilbert transform, implies that the imaginary part of the analytic signal of real variable is expressed through its real part by
1 TJ(t) = 'PV 7r
f
g(s) ds. st
(2)
The concept of analytic signal helps to solve the crucial electrical engineering problem of definition of amplitude and phase of the narrow band signal whose Fourier transform is concentrated in the small neighborhood of the carrying fre
6.7. Fourier transform of Heaviside function
163
quency Q. Such a signal is usually represented in the form g(t)
= A(t) cos[Qt + 1jF(t)],
(3)
where A(t) and 1jF(t) are slowly varying within the period T = 2Jr IQ. Thestandard engineering problem of finding A(t) and 1jF(t) from the known form of Ht) is not well posed mathematically, since it reduces to solving one equation for two unknowns A and ((i. That fundamental difficulty makes, for example, comparing accuracy for phase measurements by different phase detectors questionable. From the theoretical viewpoint the best prescription is to uniquely define the amplitude and the phase via the concept of the analytic signal, the imaginary part providing the missing second equation
+ 1jF(t)].
11(t) = A(t) sin[Qt
(4)
6.7 Fourier transform of Heaviside function Having armed ourselves with the notion of the principal value of singular integral, we are now in a position to explore Fourier transforms of the Heaviside function and related distributions. This topic was only briefly mentioned in Chapter 4. We shall begin by rewriting the Hilbert transform (6.5.1) in the language of convolutions: 1jF(t)
= ((i(t) *
Comparing the Fourier transform (6.5.2) of this function with the formula f(t) 1+ 2Jr /(w)ip(w) (3.2.8), after simple transformations we obtain that
*
((i(t)
PV
1+
sign (w).
Inverting this expression with the help of formula (3.2.5) we get that sign (t)
1+
.!..Jr PV
lW
.
(1)
Now, we are ready to recover the Fourier transform of the Heaviside function. To do that we will represent the latter in the form X (t) =
1 1. ( ) 2 + 2 SIgn t .
Chapter 6. Singular integrals and fractal calculus
164
The Fourier transform of the first summand is equal to (w ) /2 so that for the Fourier transform of the Heaviside function we have
1 ( :1) + 'PV 21l lW
1
X(t)
2
.
It is convenient to write this relation with the help of the distributions (6.2.3):
X(t)
1 1 .. .. 2m W  10

(2)
Once the Fourier transform of the Heaviside function has been calculated we can evaluate Fourier transforms of a wide class of functions (3)
representable by integrals with variable upper limit of absolutely integrable functions f(t). Indeed, if we represent F(t) as a sum F(t) = F(oo)X(t)
where
+ G(t),
G(t) = F(t)  F(oo)X(t),
then, according to (3), the Fourier transform of F(t) can be expressed in terms of the Fourier transform of G(t) by

F(w)
1 1 = F(00)2 .   . + G(w). 111 w  10
(4)
Example 1. Consider an absolutely integrable function f (t) = X(t )ePI, P > O. Then G(t) = x (t)e pl / p. Its Fourier transform exists in the classical sense. Thus, equation (4) implies that
F(w)
1 [1  1 =] 21lip w  iO w  ip .
(5)
•
Integration of (3) leads, in turn, to a new function which linearly increases as t 00. We shall learn how to find the Fourier transform of such functions by first evaluating the Fourier transform of the absolute value function It I.
165
6.7. Fourier transform of Heaviside function
Example 2. Let us write the absolute value function in the product form It I = t sign (t). The Fourier transforms of each of the two factors are already known to us. Recall that t t+ i8'(w) and that the Fourier transform of sign (t) is given by formula (6.7.1). In this fashion, with the help of formula (3.2.10) according to which the Fourier transform of a product is equal to convolution of the Fourier transforms of factors, we get that It I
. = tSlgn(t)
., (w) 18
d(l)  .
1 (:1 ) = 'PV1 * t'V 7r IW 7r dw
W
(6)
Let us identify the new distribution arising on the righthand side through its action as a functional on an arbitrary test function rp E S :
(1) dw
d I rp(w) dw t'V ;
= t'V Irp'(W) ;;;dw.
Recall that the principal value of the above integral is, by definition, equal to
t'VI rp'(w) dw = lim W
[/£ rp'(w) dw +1rp'(w) dW] . 00
W
00
W
£
For the first integral, integrating by parts expression in the brackets, we obtain that e
e
rp'(w)d  I wI W
00
00
[rp(w)  rp(s)]d 2 w. w
A similar transformation of the second integral yields eventually that
The latter limit defines a new distribution
which functionally acts via the formula
t'V
(:2) [rp(w)]
166
Chapter 6. Singular integrals and fractal calculus
= lim [/ [4>(w)  4>( e)] dw + eM w2 e
00
1 00
e
[4>(w)  4>(e)] dW] . w2
Obviously, the last equality can be written in the form
[4>(w)] = 'PV
'PV
1
[4> (W);:z 4>(O)]dw,
(7)
or in an equivalent regularized form
'PV
)
[4>(w)]
=
1 00
4>(w)
+ 4>(:::)  24>(0) dw.
(8)
o Thus, we derived a new distributiontheoretic formula
which yields the Fourier transform of function It I:
It I t+
1f
w2
(9)
• 6.8 Fractal integration In the last three sections of this chapter we develop another class of important singular integrals which arise when one tries to extend the notion of ntuple integrals and of nth order derivatives of classical calculus to noninteger (or fractional) n. We begin with the concept of fractal or fractional integration. It is natural to introduce it as a generalization of the Cauchy formula
= (n _1 1)!
j'
00 (t
 s)nlg(s) ds,
(1)
167
6.S. Fractal integration
which expresses the result of ntuple integration of function g (t) of a single variable via the single integration operator. Before we move on to fractal integrals, let us take a closer look at the Cauchy formula (1). It is valid for absolutely integrable functions g(t) which decay for t + 00 sufficiently rapidly to guarantee the existence of the integral on the righthand side of (1). Assuming that the integrand g(t) vanishes for t < 0, the Cauchy formula can be rewritten in the form
=
1
(n I)!
lot 0
(ts)nlg(s)ds.
(2)
Remark 1. The Cauchy formula can be viewed as an illustration of the general Riesz Theorem about representation of any linear continuous (in a certain precise sense) operator L transforming function g of one real variable into another function L[g] as an integral operator L[g](t) =
f
h(t, s)g(s) ds
(3)
with an appropriate kernel h(t, s). In our case, the ntuple integration linear operator in (1) has a representation via the single integral operator with kernel h(t, s) = (t  s)nl/(n  I)!. Let us check the validity of the Cauchy formula (1) by observing that the ntuple integral in (1) is a solution of the differential equation
dn
x(t) = g(t) dt n
(4)
satisfying the causality principle. Such a solution, in view of (2.2.4), can be written as convolution x(t)
= f kn{t 
where kn{t) =
s)g(s) ds,
x (t)y(t),
and y(t) is the solution of the corresponding homogeneous equation
dn
x(t) dtn
=0
(5)
Chapter 6. Singular integrals and fractal calculus
168
with the initial conditions
y(O)
= y' (0) = ... = In2) (0) = 0,
y(nl)
= 1.
Solving the above initialvalue problem we get
kn(t) =
1 (n I)!
t
nl
X(t),
the kernel that appears in the Cauchy formula (1). This is a good point to introduce fractal integrals. Replacing integer n in kernel kn by an arbitrary positive real number a and the factorial (n  I)! by the gamma function r(a) (see (4.4.3» we arrive at the generalized kernel
(6) So it is natural to call the convolution operator
(la g)(t) = ka(t)
* g(t),
(7)
the fractal integration operator of order a. In the case when function g(t) == 0 for t < 0 the convolution (7) reduces to the integral
(la g)(t) = 1
lot (t  s)al g(s) ds
r(a) 0
(8)
where the upper integration limit reflects the causality property of the operator of fractal integration. Let us establish some of the important properties of the fractal integration operator assuming, for simplicity, that function g(t) in (8) is bounded and continuous.
Existence of fractal integrals. For a 1 the integrand in (8) is bounded and continuous, and the integral exists in the Riemann sense. For 0 < a < 1, the kernel (t  s)al is singular but the singularity is integrable and the integral is an absolutely convergent improper integral. Zeroorder integration. As a 0+ the operators fa tend to the identity operator. Indeed, in view of the recurrent formula r(a + 1) = ar(a) for the gamma function and the fact that r(l) = 1 we have the asymptotics r(a) ...... 1ja, (a 0+). Hence, lim (lag)(t) =
a ..... 0+
lim !g(S)X(ts)a(tS)aldS.
a ..... O+
6.8. Fractal integration
169
Example 1 in Section 1.9 shows that the function X(t  s)a(t  s)a1 weakly converges to the Dirac deha8(t  s  0) as a + 0+. 2 Consequently,
or, equivalently, in the distributional language, ko(t)
= 8(t).
(9)
Iteration o/fractal integrals. As in the case of usual ntuple integrals, repeated application of fractal integrals is subject to the rule (10)
To see this it suffices to check that, for any a, fJ > 0, ka
* kfJ = ka+fJ·
(11)
Indeed, the lefthand side of (11), in view of (6), equals
so that, passing to the new dimensionless variable of integration r = s / t, ka
* kfJ =
where B(a, fJ) =
x(t)t a+fJ  1
r(a)r(fJ) B(a, fJ),
10 1 r a  1 (1 
(12)
r)fJ 1dr
is the beta function. It can be expressed in terms of the gamma function by B(
a,
fJ)
= r(a)r(fJ) r(a
+ fJ)·
Substituting it into (12) we obtain equality (11). 2This notation emphasizes that the support of this Dirac delta lies inside the interval (0, t) so that 8(t  s  0) ds 1 and not 1/2 as in (2.9.6). A similar situation was encountered in formula (2.9.4) for ex + +00
=
170
Chapter 6. Singular integrals and fractal calculus
Fractal integrals as continuous operators. The following two inequalities show that integration of fractal order has some continuity properties as a linear operation. These properties will find an application in our construction of Brownian motion in Chapter 14 of Volume 2, and for those purposes it will suffice to assume that o < ex < 1. First, observe that, for ex > I Ct is a continuous operator from L2[0, 1] into LP[O, 1] for each p < 00. 3 Indeed, by the Schwartz Inequality,
!,
{1(10t
:SJo
f2(s)ds
:sc(ex,p) (
)P/2 . (t 10
)P/2dt
(13)
10(1 f2(s)ds )P/2
where c(ex, p) is a constant depending only on ex and p. Additionally, by a similar argument, but using the Holder Inequality with 1/ p + l/q = 1, IIIP flloo = sup
It
0:9:::1 10
t
{1 )l/q (10 :s ( 10 Ikp(s)lqds
kp(t  s)f(s)dsl
If(s)jPds
)1/P = c({3, p)lIf11
p'
(14)
so that for any {3 > 1/ p, the operator IP is continuous from V[O, 1] into the space of continuous functions qo, 1].
6.9
Fractal differentiation
The operator D Ct of fractal or fractional differentiation is defined as the inverse of the operator I Ct of the fractal integration, that is, via the operator equation (1)
3Recall the LP[o. 1] denotes the Lebesgue space of functions
f
on the interval [0, 1] which have
pth powers integrable, Le., for which the norm IIfllp := (Jo1If(sW ds)l/p <
00.
171
6.9. Fractal differentiation
where Id denotes the identity operator. Similarly to the fractal integration operator, the fractal differentiation operator has an integral representation
(2) where ra (t) is the convolution kernel which will be identified next. To do that notice that the operator equation (1) is equivalent to the convolution algebra equation (3)
where ka is the fractal integration kernel (6.8.6). Denote by y the solution of the equation a + y = n, where n = a 1 is the smallest integer greater than or equal to a. In other words,
r
n 1 < a
n,
(4)
y =n a,
Applying the fractal integration operator /Y to both sides of (3) we get
The expression on the lefthand side represents the usual ntuple integral of the fractal differentiation kernel r a. Thus, if we differentiate it n times we arrive at the explicit formula y =n a> O. (5) for the fractal differentiation kernel. The corresponding fractal differentiation operator is then given by the convolution (D a g)(t) =
* g(t).
(6)
For t > 0, the nth derivative of the kernel ky(t) appearing in (5) exists in the classical sense and
1
ra(t) = _ _ t a  1 = La(t), r(a)
t > O.
(7)
So, the fractal differentiation kernel is equal to the fractal integration kernel with the opposite index a. Here, the gamma function r (a) of negative noninteger variable is defined via the above mentioned recurrence property as follows:
r(a)
= r(n 
a)/(n  a 1)· ... · (a).
172
Chapter 6. Singular integrals and fractal calculus
Therefore, the fractal differentiation operator can be treated as the fractal integration operator of the negative order:
(8) which adds attractive symmetry to the fractal calculus. Consequently, for a function g(t) which vanishes on the negative halfaxis, equation (6) can be rewritten in the following symbolic integral form:
(9) For a 0, the above improper integral diverges in view of the nonintegrable singularity of the integrand in the vicinity of the upper limit of integration. Therefore, its values have to be taken as the values of the corresponding regularized integral which can be found treating equality (6) as the convolution of distributions. In view of properties of the distributional convolution, the operation of ntuple integration can be shifted from the first convolution factor to the second so that
(D a g)(t) = ky (t)
* g(n)(t),
(10)
with converging integral on the righthand side. In particular, for integer a = n, taking (6.8.9) into account, we obtain (as expected) that
Example 1. Let us now consider the special case of a function g(t) which vanishes identically for t < 0 and is of the form g(t) = X (t)q,(t),
(11)
where q,(t) is an arbitrary infinitely differentiable function. Differentiating (11) n times and taking into account the multiplier probing property (1.5.3) of the Dirac delta, we obtain
= X (t)q,(n)(t) + L
nl
g(n)(t)
8(m) (t)q,(nml) (0).
m=O
Substituting the above formula in (10), and remembering that y = n  a and t
> 0,
6.9. Fractal differentiation
173
we finally get that
it
1
(D a g)(t) =
f(n  a) 0
(t  s)nalt/J(n) (s) ds
+ L t/J(nml) (O)kn m a (t), nl
t >
o.
(12)
m=O
In particular, D a kfJ = kfJa.
Also Da(X(s)sfJ log lsI) =
f(fJ + 1) X(s)sfJ a [ log lsi f(fJ a + 1)
+ C] ,
where the constant (see the literature in the Bibliography for its derivation and other formulas of the fractal calculus) C = (log f)' (fJ
+ 1) 
(log f)' (fJ  a
+ 1).
Furthermore, for 0 < a < 1, we have a
(D g)(t)
1
= f(1 _ a)
f' t/J'(s) I_a Jo (t _ s)a ds + t/J(O) f(1 _ a) t ,
t > O. (13)
Observe that, in contrast to (9), the singular integral on the righthand sides of (12) and (13) converges absolutely, and that the regularizations (10), (12) define a new distributionthe principal of function Ta(t)X(t) (7):
'PV x(t)_I_ t  a 1 , f(a)
Its convolution (D a g)(t) = 1 'PV f(a)
it 0
a ::: O.
(14)
(t  s)al g(s) ds
with any function of the form (11) has a distributional interpretation via the righthand side of formula (12). • The following properties of the operation of differentiation of fractal order highlight its peculiarities.
174
Chapter 6. Singular integrals and fractal calculus
Nonloeal character. Values D n get) of the usual derivatives of integer orders depend only on values of function get) in the immediate and arbitrarily small (infinitesimal) vicinity of the point t. By contrast, the fractal (noninteger) derivatives are nonloeal operators since the value Da get) depends on the values of get") for all t" < t. In particular, this fact explains why a function's discontinuity at a certain point (t = 0 for function (11» generates slowly decaying "tails" in its fractal derivatives (the last sum on the righthand side of formula (12». Causality. Fractal derivatives enjoy the causality property: If function get) is identically equal to zero for t < to then so does its fractal derivative. Scale irwarianee. Like usual derivatives, fractal derivatives are scale invariant. This means that differentiation of the compressed (K > 1) function gK (t) = g(Kt) just requires mUltiplication of the compressed derivative by the compression factor:
(15)
Fourier transform. Under Fourier transformation, fractal derivatives behave just like the ordinary derivatives. In the distributional sense (16)
This formula follows directly from (10) and from the results of Section 4.4. Remark 1. Fractal Laplaeians. The above definitions of fractal differential operators for functions of one variable can be extended to fractal partial differential operators for functions of several variables (see the literature at the end for further details). For example the fractal Laplacian can be defined through the Fourier transform approach as follows: For any
The following integration by parts formula is then obtained via the Parseval equality
It is also clear that the fundamental solution (Green's function) G a of the equation
l1au = 8 has the Fourier transform
175
6.10. Fractal relaxation
The explicit inversion depends on the dimension d of the space. If a = 2, d or 0 a 2, d 2, or 0 a 1, d = 1, then
3,
6.10 Fractal relaxation Recently, more and more frequently, physicists find applications for the fractal calculus which permit construction of generalized mathematical models of such phenomena as relaxation, diffusion and wave propagation. In this section we will illustrate these possibilities in the example of relaxation processes. Informally, one says that a physical system has the relaxation property if within finite time. (called the relaxation time) of cessation of external perturbation the system "forgets" the perturbation and returns to its original state. Such systems are encountered in the wide spectrum of applied problems from physics and electrical engineering to biology and economics. In the simplest case, the mathematical model of a linear relaxing system is the firstorder ordinary differential equation X'
+ px =
g(t),
(1)
where function g(t) describes the external perturbation of the system and p = 1/. > 0 is called the relaxation frequency. The response x(t) of the system (1) to the external perturbation which satisfies the causality principle is given by formula (2.2.4): x(t)
= H(t) * g(t).
(2)
The fundamental solution H (t) entering (2) satisfies equation (2.2.5) which, in our particular case, is of the form H'
+ pH =
8(t),
H(t) = 0 for t < O.
(3)
Its well known solution is H(t) = x(t)exp(pt).
(4)
From the physicist's perspective the main feature ofthe relaxation model (1) is the presence of a unique (for this system) characteristic relaxation time • = 1/ p. The model itself is just a special case of the whole family of justifiable models
176
Chapter 6. Singular integrals and fractal calculus
described by fractal differential equations Dax + pax = g(t),
O0,
(1)
where functions I(t) and g(t) have finite support Fourier images, identically equal to zero for Iwl O. Then, the corresponding analytic signal is equal to
= [/(t)  ig(t)]e mt .
(2)
Prove it. Remark 1. Recall that the imaginary part of the analytic signal (6.6.2) is equal to minus the Hilbert transform (6.5.1) of This implies the following corollary to the above result: If (j/ 1 cos Ot + g sin Ot then its Hilbert transform is y, g cos Ot  1 sin Ot. Remark 2. Signals with finitesupport Fourier image seldom appear in electrical engineering applications. However, for narrowband signals the replacement of the actual analytic signal by the expression (6.9.3) gives a rather good approximation. For example, the signal from Exercise 6 has an unbounded Fourier image but is narrowband if O. » 1. In this case, it is easy to see that the approximate expression
=
=
a
(t)
1 = __ e,..... ..t t 2 +.2
is very close to 10. Let
= t sin Ot, 0 > O. Find
11. Use the concept of analytic signal to find the instantaneous amplitude, phase and frequency of 12. Let signal. 13. Let
= X(t) sin Ot. Find the imaginary part 1/(t) ofthe corresponding analytic = X(t)t a 
1,
0 < a < 1. Find
14. Find the analytic signal 15. Let 16. Let
(t) which corresponds to the even function
(t)
= It lal .
= 8' (t). Find the imaginary part 1/(t) ofthe corresponding analytic signal.
Ii be the Hilbert transform operator liq,(t)
= !..1r 'PV
f
q,(s) ds.
st
182
Chapter 6. Singular integrals and fractal calculus
Find the set of functions 4J(t) on which H acts as just the shift operator, i.e., H4J(t) = 4J(t + T), for a certain T. Fractal calculus 17. Extend the Cauchy formula to the case of the ntuple integral x(t) =
dt1a1(t1)
£:
dt2a2(t2) ...
i:
1
dtnan(tn)g(tn),
where a1 (t), a2 (t), ... , an (t), are known functions such that the above ntuple integral converges absolutely for any integrable function g(t). 18. Let A C Rn. Express the multiple integral 1=
f·'?· i
a(x1,x2, ... ,xn)g(b(X1,X2, ... ,xn»dx1 ... dxn
via a single integral of function g (u).
Chapter 7 Uncertainty Principle and Wavelet Transforms
The method of wavelet transforms, which provides a decomposition of functions in terms of a fixed family of functions of constant shape but varying scales and locations, recently acquired broad significance in the analysis of signals and of experimental data from various physical phenomena. It is clear that the potential of this method has not yet been fully tapped. Nevertheless, its value for the whole spectrum of problems in many areas of science and engineering, including the study of electromagnetic and turbulent hydrodynamic fields, image reconstruction algorithms, prediction of earthquakes and tsunami waves, and statistical analysis of economic data, is by now quite obvious. Although the systematic ideas of wavelet transforms have been developed only since the early 80s, to get the proper intuitions about sources of their effectiveness it is necessary to become familiar with a few more traditional ideas, tools and methods. One of those is the celebrated uncertainty principle for the Fourier transforms which will be given special attention in this chapter. A close relative of the wavelet transformthe windowed Fourier transform, will also be studied in this context. We begin though with a brief sketch of the notion of the functional Hilbert space which provide a convenient framework for our analysis.
7.1 Functional Hilbert spaces The extension of classic 3D Euclidean geometry concepts such as the space, vector, composition, multiplication of vectors by scalars, inner product of vectors, angle, orthogonality and parallelness, to a broad class of mathematical objects, was one of the success stories of twentieth century mathematics. As a result, a multitude offunctional spaces were introduced, studied and added to our permanent arsenal. The linear topological spaces of distributions briefly described in Section 1.9 are one such example. In this section, we discuss another class of functional spaces called Hilbert spaces. © Springer Nature Switzerland AG 2018 A. I. Saichev and W. Woyczynski, Distributions in the Physical and Engineering Sciences, Volume 1, Applied and Numerical Harmonic Analysis, https://doi.org/10.1007/9783319979588_7
184
Chapter 7. Uncertainty principle and wavelet transforms
FIGURE 7.1.1 Composition of 2D vectors. At first, recall the basic notions of the usual 3D geometry where each point of the space is identified, in a fixed Cartesian coordinate system, with a vector a anchored at the origin and with a tip at a given point. Each such vector is uniquely described by its coordinates (al, a2, a3)an ordered triple of real numbers. If b = (bI. b2, b3) is another 3D vector then the inner or scalar product of these two vectors is defined by the equality (1)
Since, alternatively, (a, b) =
lIallllbll cos a,
where, by Pythagoras' theorem,
lIall 2 =
(2)
(a, a)
is the square of the vector's norm (length, magnitude) and a is the angle between the two vectors, the inner product of two vectors clearly depends on their mutual orientation. In particular, aJ.b
if and only if
(a,b) = O.
The geometric composition of vector a with vector b (see, Fig. 7.1.1) corre
7.1. Functional Hilbert spaces
185
sponds to the algebraic operation
(3) of vector addition. Vector c = (Ct. c2. C3) is called the sum of vectors a and b if its coordinates are sums of corresponding coordinates of the summand vectors: Cn = an + bn • '! = 1.2.3. Such addition operation is obviously commutative, that is
a+b= b+a.
(4)
Besides addition, one introduces the operation of multiplication of a vector by a scalar: (5)
which geometrically represents vector contraction for Irl < 1, vector dilation if Ir I > 1 and vector reflection in case of r = 1. The following three properties of the inner product and the norm are fundamental for the geometric properties of the Euclidean space: (i) The norm is homogeneous, that is, for any scalar r and vector a,
IIrall
=
Iriliali.
(6)
(ii) The norm and the addition operation are related via the triangle inequality, that is, for any two vectors a and b,
lIa + bll
lIall + IIbll·
(7)
(iii) The norms and the inner product are related by the Schwartz inequality
I(a. b)1
lIalillbli.
(8)
The first step in the generalization of the geometry of 3D Euclidean spaces to abstract functional Hilbert spaces are two observations: (a) The concept of the inner product (and related geometry of the space) can be immediately extended to ddimensional Euclidean spaces by defining for any a = (at •...• ad) and b = (bt. ...• bd)
(b) If one wants to operate with vectors with complex (rather than real) coordinates and preserve the positivity of the norm the only adjustment in the definition of the inner product is as follows:
(9)
186
Chapter 7. Uncertainty principle and wavelet transforms
where, as usual, the asterisk denotes the complex conjugation. Then, the square
lIall 2 =
d
d
= n=l
defines a positive norm
lIall
L Ian 12 ::: 0
(10)
n=l
= J(a, a).
(11)
In the complex case, the symmetry of the inner product is replaced by the Hermitian property (a, b) = (b, a)*. (12) The linearity with respect to the second variable in the inner product is preserved, though, as for any complex number z, (a, zb) = z(a, b).
(13)
Example 1. Inner product space ofpolynomials. Let us consider the set all polynomials
of
of degree at most N  1 with complex coefficients. Each of these polynomials is uniquely determined by its coefficients for different powers of x, that is by the (complex) vector a = (ao, ... , aNI). The sum of such polynomials is again a polynomial of the above type and the same is true for a product of a complex number and such polynomial. Moreover, the summation and multiplication by scalars in the family of such polynomials corresponds to the analogous operations on the of all coefficient vectors and identifies the vector space structure of the family polynomials of degree at most N  1 as that of an Ndimensional complex vector space. The natural inner product leads to the notion of the "distance" I = lIa  bll between polynomials with coefficient vectors a and b. •
Example 2. Inner product space of complex exponentials. Consider the set of all infinite sums of complex exponentials 1 E(x) =  
L
00
...tiiin=oo
ane inx ,
x E [1l',1l']
(14)
with complex coefficients. This set forms a natural vector space under termwise addition and multiplication by scalars and can be identified with the (infinitedimensional) vector space of coefficient vectors (sequences) a = (... , aI, ao, aI,
7.1. Functional Hilbert spaces
187
a2, •. .). However, if we wanted to introduce the inner product in such space
associated with the norm
liE II =
we immediately run into the question of convergence of the above series and to proceed we have to assume additionally that
L
00
n=oo
lan l2
<
(15)
00.
This is the first fundamental difference with the finite dimensional spaces. The of sums E satisfying condition (15) remains closed under operations of subset termwise addition and multiplication by scalars since IIzE II = Izlil E II and
liE + FII
IIEII + IWII,
•
for any sums E, F and scalar z.
Attempts to generalize the above examples immediately lead us to the idea of the functional Hilbert space L 2 (R) of complexvalued functions I (x) defined on the entire real axis R and such that (16)
Now, we can introduce the inner product in L2(R) via the formula (f, g) =
I
(17)
I*(x)g(x) dx
and the related norm
11111 = J(f, f) =
(
l'/(X),2dX )
In view of the classical integral Schwartz inequality
l(f, g)1 =
II
f*(X)g(X)dXI
1/2
(18)
Chapter 7. Uncertainty principle and wavelet transfonns
188
:: f (
1/2 (
If(x)1 2dx )
f
Ig(x)1 2dx
)1/2
= IIfllllgll,
(19)
condition (16) assures that the inner product (17) is well defined (i.e., that f*g is integrable). The Schwartz inequality also immediately leads to (see Exercises) the triangle inequality
IIf + gil ::: IIfil + IIgll,
f, g
E L 2 (R),
(20)
which, incidentally implies that the Hilbert space L2(R) is closed under the usual pointwise addition of functions. It is also closed under multiplication by scalars since, by (18), The above inner product (17) is Hermitian, that is (f, g) = (g, 1)*
and homogeneous in the second variable, as (f, Zg) = z(f, g)
for any complex constant z and f, g E L2(R). Of course, the norm in L2(R) is nonnegative, that is for any
f
E L2(R)
IIfil and if f(x)
== 0 then IIfli
=
o.
A mathematical aside: functions with vanishing norm and the Lebesgue integral. 1 It is quite clear that beside the function equal to zero identically there are other functions f(x) :F 0 such O. For example, any that II f II O. In other words, II!II 0 does not necessarily imply f (x) function different from zero at a finite or countable number of points would have norm zero. This creates somewhat unpleasant situation of having two different functions f(x) :F g(x) for which o. The satisfactory resolution of this problem is not possible within the the distance II f framework of the Riemann integral which we implicitly used throughout the preceding chapters (and which is sufficient for our other purposes). It requires introduction of the more general Lebesgue integral (hence letter L in the notation of the functional Hilbert space) which permits integration of a much broader class of functions than the Riemann integral (see the bibliographical notes at the end). For example the Dirichlet function
=
=
=
g" =
D(x)
= {0,
1,
x
if x
IS
rational;
1This material may be skipped by the first time reader.
189
7.1. Functional Hilbert spaces
is not integrable on [0,1] in the Riemann sense since the upper approximating sums (always equal to 1) and the lower approximating sums (always equal to 0) do not converge to the same number. However, it is integrable in the Lebesgue sense and its Lebesgue integral is equal to O. Interpreting the integrals in (1618) as Lebesgue integrals, it is customary to formally define the functional Hilbert space L 2 (R) not just as the space of squareintegrable functions but as the space of equivalence classes of squareintegrable functions where two functions I, g are understood as equivalent (written I = g) if and only if III  gil = O. It is easy to see that, if we define the Lebesgue measure IA I of the set A eRas the Lebesgue integral of its indicator function I A (x) (equal to 1 on A and to 0 off A), then two functions belong to the same equivalence class if and only if they differ on a set of Lebesgue measure zero (or in measure theory jargon, are equal almost everywhere). Now, the norms and inner products are the same for all functions in a given equivalence class, so practically one always does computations on concrete functions, but if elements I of the Hilbert space L2(R) are meant as equivalence classes of functions equal almost everywhere then the desired strict positivity of the norm is achieved, that is
11/11
if, and only if
=0
1=0.
The norm IIfII in the functional Hilbert space L 2 (R), and the related distance II I  g II of (equivalence classes of) functions I, g permit us to introduce the notion of the limit of functions in L2(R). Namely, we say that
In if lim
n+oo
I
(21)
II/n  III = o.
This notion permits us to study approximation problems in the functional Hilbert space. Having introduced in the functional Hilbert space the algebraic structure (addition and multiplication by scalars), the compatible inner product structure and the related metric (norm) structure one could sensibly introduce in L2(R) and study the geometric concepts such as the angles between functions, their orthogonality, etc. This very fruitful approach is the essence of the branch of mathematics called functional analysis. Remark 2. Completeness of the functional Hilbert space. The finite dimensional inner product spaces introduced at the beginning of this sections enjoyed the important property of being complete, that is, the Cauchy criterion of convergence remained valid for them. The same criterion happens to hold true for the space L 2 (R). In other words, the functional Hilbert space is a complete inner product space, that is for any sequence of functions In E L 2 (R), n = 1, 2, ... , such that
lim
n,m+OO
II/nImil=O
190
there exists a function
Chapter 7. Uncertainty principle and wavelet transforms
I
E
L2(R) such that
lim
n+oo
II/n  III
= O.
Remark 3. Other Hilbert and normed functional spaces. The above discussion applies, with obvious adjustments, to Hilbert spaces of functions over other subsets of the real axis such as L2([a, b]), L 2([0,00», etc., as well as to their natural analogues L 2 (Rd ) defined for functions of several variables. On the other hand, it is often necessary to consider functional spaces where the introduction of the inner product structure is impossible and one has to be satisfied only with the norm structure. Examples of nonHilbertian functional normed spaces, such as LP(R), 1::: p < 00, p =f:. 2, which consists of all (equivalence classes of) functions for which IIfllp = ( / I/(x)I
PdX) lip < 00,
have been mentioned before in Chapter 6 (also, see the Bibliographical Notes). In particular, space L I = L I (R) consists of all absolutely integrable functions on the real axis.
7.2 Timefrequency localization and the uncertainty principle Consider a (perhaps complexvalued) signal I(t) such that (1)
The quantity 1I (t) 12 can be thought of as the signal's "mass" density and describes its distribution in time. If the signal I (t) is square integrable but (1) is not satisfied then one can always normalize it by considering I(t)/(j I/(t)1 2dt)1/2. In this context, the quantity
can be interpreted as the location in time of the signal's "center of gravity," or its mean location. For the purposes of this section, and without loss of generality, we
7.2. Timefrequency localization and the uncertainty principle
191
will assume that its mean location is at 0 or, in other words, that J tlf(t)1 2dt = O. In this case, the quantity
(2) measures the average square deviation from the mean time location, or the degree of localization of the signal around its mean in the time domain. On the other hand, the Fourier transform
f(w) = 1 21l'
f
.
f(t)e/(lJt dt
displays no direct information about the signal's time localization, but has explicit information about its frequency localization. The square of its modulus Ij(w)1 2 is the frequency domain counterpart of the time density If(t)12. Note that, by Parseval's formula (3.2.11),
so that 21l' 1j (w) 12 can be viewed as the signal's normalized density in the frequency domain. Assume (again, without loss of generality) that the mean frequency
Then the quantity (3)
measures the mean square deviation from the mean frequency location, or the degree of localization of the signal in the frequency domain. The uncertainty principle asserts that there exists a lower bound on the simultaneous localization of the signal in time and frequency domains. More precisely, it states that
(4)
whenever the variances 0"2[11 and 0"2[j] are well defined. Note the universal constant 1/4. To see why the uncertainty principle holds true, consider the integral I(x) =
f
Ixtf(t)
+ !'(t)1 2dt
0
(5)
192
Chapter 7. Uncertainty principle and wavelet transforms
where x is a real parameter. Then, since
Ixtf(t)
+ f'(t)1 2 =
(xtf
+ f')(xt/* + (f')*)
we get that
The first integral in (6) is just equal to
U
2 [f] (by definition (2». The second integral is
since tlfl2 decays to zero at ±oo in view of the assumption u 2[f] < the third integral is equal to
00.
Finally,
because of Parseval's formula (3.2.11) and the fact that the Fourier transform of f' is i(Oj«(O). As a result, the integral
(7) This is a quadratic polynomial in variable x and, in view of (5), it is nonnegative for all values of x. As such, it has a nonpositive discriminant
which immediately yields the uncertainty principle (4).
Remark 1. The Heisenberg uncertainty principle in quantum mechanics. The (3D version of the) above uncertainty principle concerning timefrequency localization has a celebrated interpretation in quantum mechanics, where the principle asserts that the position and the momentum of a particle cannot be simultaneously measured with arbitrary accuracy. Indeed, in quantum mechanics the particle is represented by a complex wave function f(x), where If(x)1 2 is the probability density of its position in space. The observables are represented by operators A on wave functions; the mean value of the observable is
f
(Af)(x)/*(x)dx.
7.3. Windowed Fourier transform
193
The position observable is represented by a multiplication by variable (vector) x and the momentum observable is represented by the operation of differentiation a/ax. However, via the Fourier transform, the latter also becomes an operation of multiplication but by an independent variable (vector) w in the frequency domain. Thus the uncertainty principle (4) gives the universal lower bound for the product of variances of the probability distributions of the position and of the momentum. In the threedimensional space, and in the physical units, the lower bound 1/4 in (4) .has to be replaced by a different mathematical constant multiplied by a universal physical constant called the Planck constant. The employed above probabilistic concepts of means and variances will be further studied in Chapter 13.
Remark 2. One can check that the equality in the uncertainty principle (4) obtains only for the Gaussian function / (t) = 11: 1/4 exp( _t 2 /2). Thus the optimal simultaneous time and frequency localization is attained for a Gaussianshaped signal.
7.3 Windowed Fourier transform 7.3.1. Forward windowed Fourier transform. The uncertainty principle discussed in Section 7.2 is a basic law of mathematics and it is impossible to fool nature by measuring the frequency of the incoming signal with an arbitrary precision in a finite time interval. Moreover, for most of the signals we have to deal with in practical problems, such as speech, musical sounds, radar signals, the situation is often much worse than the basic uncertainty inequality permits and u (f)u (/)
» 1.
(1)
Nevertheless, it is often possible to process these signals in such a way that, without violation of the uncertainty principle, one can obtain information about the signal's "current" frequency and its time evolution. These various practical signal processing methods are adapted to different kinds of signals and pursue different goals. In this section we will take a look at one of these methods called the win dowed Fourier transform which is closest perhaps to the spirit of the usual Fourier transform. In what follows, to better grasp the mechanisms behind the windowed Fourier transform, it will be instructive to test them on a sample signal that we will call the simplest tune. Mathematically, it is described by the real part of the complex function
/(t) = exp(i 4>(t»,
(2)
194
Chapter 7. Uncertainty principle and wavelet transforms
where
(t) = wot
Q .
+ V
(3)
SID(vt)
is the signal phase. The simplest tune is plotted on Fig. 7.3.1.
Ref(t)
t
v FIGURE 7.3.1 Plot of the simplest tune in case of wo = 10v and fJ = Q/v = 5. It is customary to say in the theoretical physics context that the simplest tune has the instantaneous frequency (admittedly, an oxymoron)
d(t) dt
Winsr(t) =   =
W()
+ Q cos(vt),
(4)
which oscillates with period T = 2Jr / v between its high value Wo + Q and low value W()  Q. By contrast with a theoretician, an experimenter has to deal not with mathematical formulas but with real signals and his job is to come up with a signal processing method that will discover the existence of frequency oscillations in the simplest tune. The mathematical tool that is helpful in this situation is called the windowed Fourier transform which is just the usual Fourier transform
f(w, r)
= 2Jr1
f
.
f(t)g(t  r)e uut dt
(5)
of the timewindowed signal f(t)g(t  r), where g(t) is the windowing function that usually is chosen to have value equal to 1 in a vicinity of the origin t = 0 (say,
7.3. Windowed Fourier transform
195
inside an interval of length),,), and that either vanishes or has values very close to 0 outside this neighborhood. This windowing function property will assure effective timelocalization. Usually, one defines the windowing function g(t) via a windowing shapefunction go(x) of a dimensionless variable x and the formula g(t) = go(t/),,),
(6)
where)" is a scaling parameter. Some typical examples of normalized (1\ gO 1\ windowing shape functions are (see Fig. 7.3.2):
= 1)
(a) Finite memory window
+ 1) 
X(x);
(7a)
go(x) = 2x(x)exp(2x);
(7b)
go(x) = X(x
(b) Relaxation window
(c) Gaussian window (7c)
Shift'l' centers the window at different locations on the timeaxis t. If /(t) is a timedependent signal and processing is performed in the realtime then 'l' is just the current time of the experiment and the timewindow g(t) has to satisfy the causality principle, i.e., g(t) == 0 for t > O. So, in this case, the finite memory and relaxation windows are appropriate but the Gaussian window is not. If the whole signal is recorded before processing, or the variable t has other interpretation (e.g., space, or angle variable) then the experimenter has more freedom in selecting the windowing shape function, and very often the Gaussian window is a good candidate. 7.3.2. Frequency localization. The timewindow g(t) was designed to separate well the timelocalized pieces of duration)" of the incoming signal /(t). Luckily, it turns out that the Fourier image of the timewindow g(t) can help in frequency localization. To see how this happens let us express the original signal/ (t) through its Fourier transform: /(t) =
f
j(w')e ialt dw',
(8)
and substitute it into the righthand side of (5). Note that, in the case of the simplest tune (2), j(w) exists only in the distributional sense. The change of the integration
Chapter 7. Uncertainty principle and wavelet transfonns
196
a
'Ix 1 o 1 2
2
b
2
1
o
1
2
c
1
2
1
o
1
2
FIGURE 7.3.2 Examples of windowing shape functions. (a) Finite memory window; (b) Relaxation window; (c) Gaussian window.
197
7.3. Windowed Fourier transform
order gives that j(w, t")
= eiaJ'r
f
j(w')'g(w  w')e iah: dw'.
(9)
Remarkably, except for the nonessential factor in front of the integral, this expression looks like the symmetric counterpart of (5) in the frequency domain. Now, the role of the signal is played by its Fourier image j(w) and the timewindow has been replaced by the frequencywindow g(w). The uncertainty principle (7.2.4) tells us that if the effective duration of the timewindow is J.. then one can expect the effective width of the frequencywindow to be of order at least 1/J... In terms of the dimensionless window shapes 8o(x) and go(y) where, similarly to (6), g(w) = J..go(J..w),
(10)
both 8o(X) and go(y) have to have a similar effective widths'" 1. However, the actual situation is a bit more complicated than the above juggling of the uncertainty principle may indicate. When the engineers talk about effectively localized frequencywindow, they think about the compact support of the frequencywindowing shape function go(y) or at least about its rapid decay outside a finite frequency band. However, we know from the properties of the Fourier transform that it is impossible for both the function and its Fourier transform to have compact supports. Furthermore, the frequency windowing shape function will decay rapidly for Iyl > 1 only if the time windowing shape function is smooth. This fact eliminates time windowing shape functions (7a) and (7b), which have good timelocalization properties, as good candidates for good frequency localization by their Fourier transforms. Abrupt truncations in them introduce discontinuities of the first kind which slow the decay of their Fourier transforms. For example, the modulus of Fourier image of the relaxation window (7b)
Igo(y)l=
1r
kz 4+ y
(11)
decays to zero slowly as Igo(y)1 '" 1/(lrlyl), (y + 00). So, to achieve better frequency localization one has to take smoother windowing shape functions like, for example, functions described in Section 4.3.
Example 1. Compact timewindow, powerlaw decay of the frequency windOw. Take the windowing shape function 8o(X) =
+ 2) 
x (X)] sin2 (lr;).
(12)
Chapter 7. Uncertainty principle and wavelet transforms
198
corresponding to function (4.3.19), normalized appropriately and shifted to satisfy the causality principle. Its frequency counterpart ·4 _ go(y) = e'Y
siny
(13)
3 y(rr2 _ y2)
decays as l/lyl3, faster than (11), which produces tolerable frequency localization while preserving perfect time localization. The power law of the frequency windowing shape (13) decay was caused by hidden discontinuities (in the second derivative) of the time windowing shape (12). • Example 2. Gaussian time and frequency windows; Gabor transform. Since the Fourier image of a Gaussian time windowing shape gives a Gaussian frequency windowing shape, in this case we have excellent localization in both time and frequency domains. Indeed, if gO (x) is given by (7c) then go(y)
=
1
v2rr
=
rr 1/ 4
v2rr
(14)
exp( _y2 12).
The extra factor 1/.,fii is the result of our asymmetric, but physical definition (3.1.1) of the Fourier transform. In mathematics, for esthetic reasons, one often prefers a symmetric definition of the Fourier transform and its inverse:
1 f(y) =.,fii In this case, go(y)
f
.
f(x)e Ixy dx,
1 f(x) =.,fii
f .
f(x)e 'xy dy.
(15)
== go(Y). The windowed Fourier transform (16)
based on the Gaussian window is called the Gabor transform in honor of the physicist who introduced it for studying quantummechanical problems. • 7.3.3. Energy density in the timefrequency domain. In applied problems the quantity of interest is usually not the complex function (w, r) itself but its squared modulus 2 (17) E(w, r) = 2rrlf(w, r)1 .
1
It follows from (5) that E(w, r)
=
f f dt
dfJ e iw8 f*(t)g*(t  r)f(t
+ fJ)g(t + fJ 
r).
(18)
199
7.3. Windowed Fourier transfonn
For the sake of symmetry between the time and frequency domains we permit both functions I(t) and g(t) to be complexvalued. Integrating the above equality over all w, we get (19) Observe that the three unwieldy integrals on the righthand side were reduced to an elegant single integral by noticing first that
21T
f
e iw8 dw
= l)(O),
(20)
and then using the probing property of the Dirac delta to get rid of another integral. After integration of (19) over all r we have (21)
where the norms on the righthand side are the Hilbertian L2norms introduced in (7.1.18). Since we assumed at the beginning of this section that the time windowing shape function go is normalized, i.e., IIgoll = 1, the squared L2norm of the time windowing function itself (22) i.e., it is equal to the duration of the timewindow. Remembering that the squared norm
11/112 =
f I/(t)1
2 dt
represents the energy Ef of the original signal I(t), formula (21) implies the following energetic relation: Ef =
f f dw
dr E(w, r).
Thus the function E (w, r) /}.. has a physical interpretation as the joint frequencytime density of the signal's energy.
7.3.4. Mean frequency and standard deviation of the windowed Fourier transform. Recall that the windowed Fourier transform was introduced earlier in this section to track (with a precision determined by duration).. of the time window) the time revolution of the "instantaneous frequency" Winst(r) of signal I(t). The latter was sufficiently clearly defined for the simplest tune signal, but for the general
200
Chapter 7. Uncertainty principle and wavelet transforms
situation we need a more rigorous definition. A good, and analytically convenient definition is the mean value w('l') = e('l')
where e(r)
= 1/
f
f
wE(w, r)dw,
1/(t)1 2 Ig(t  r)1 2 dt
(23)
(24)
is the normalizing constant, although some physicists would perhaps prefer to use, as the definition of the instantaneous frequency at time r, the value Wmax = Wmax ( r) which maximizes the joint energy density E(w, r). However, before advising the reader to go ahead and apply the above definition in research problems, let us step back and see what happens in a typical situation where /(t) and g(t) are real functions. Then, in accordance with the Fourier transform properties, !(w, r) = !*(w, r), so the energy density E (w, r) is an even function of w for each r and, necessarily, w (r) == O. In this context, the above notion of the "instantaneous frequency" is useless. The situation is different and more promising for analytic signals /(t) = A(t) exp(i (t»,
(25)
where A (t) and (t) are, respectively, the signal's realvalued amplitude and phase. In physical and engineering applications, the real signal / (t) is often a narrowband process which can be written in the form /(t)
= A(t) cos(W()t + rp(t»,
(26)
where A(t) and rp(t) are slowly varying in comparison to cos(W()t). The "simplest tune" signal introduced earlier in this section is narrowband if W()>>
n,
and
W()
»
v.
(27)
In such cases we will consider "approximately analytic" signals, replacing in (25) the exact amplitude and phase by the amplitude A(t) and phase (t) = W()t
of a narrowband process.
+ rp(t)
(28)
201
7.3. Windowed Fourier transfonn
Let us calculate the current frequency w( r) (23) of an analytic signal (25). Multiplying (18) by w, integrating, and keeping in mind that the differentiation of (20) with respect to () gives
21T we obtain w(r) = ic(r)
f
f
weicuB dw = i8'«()),
f*(t)g*(t  r).!!...[!(t)g(t  r)]dt. dt
Substitute in this formula the expression (25) for the analytic signal and take into account that the window g(t) is a realvalued function to arrive at the formula w(r)
= Winst(r) =
where
f
Winst(t)P(t; r)dt,
d!l>(t) Winst(t) =  dt
(29)
(30)
is the "instantaneous frequency" of the analytic signal, and (31)
is the normalized power of the signal taking into account the window's weight Ig(t  r)1 2 . If I/(t)1 2 is constant, as in the case of the simplest tune signal (24), then the power 1 2 (t  r) = go 1 2  . P(t; r) = g
(t r)
)..
)..)..
For).. 0, the power function P(t; r) weakly converges to 8(t  r) and Winst(r) converges to the instantaneous frequency Wi nst ( r ). Unfortunately, the above conclusion does not imply that the windowed Fourier 0, the precise measurement of the signal's intransform permits, in the limit ).. stantaneous frequency. Actually, the accuracy of such measurement is determined by the frequency deviation a (r) = ,J D( r ), where D(r) = c(r)
f
(w  w(r))2 E(w, r) dw.
(32)
Simple algebra shows that D(r)
= w2(r) 
(w(r))2,
(33)
202
Chapter 7. Uncertainty principle and wavelet transforms
where the second frequency moment (34)
Calculations similar to those that brought us to the expression (29) for the current frequency give
or, after substitution of the analytic signal (25), (35) Here, in agreement with notation from formula (29) (36)
Do(r) = c(r)
f
[:t If(t)g(t
r)lr
dt.
(37)
Finally, substituting (35) into (33) we obtain (38) where (39) Both terms on the righthand side of (38) have an obvious physical interpretation. The first term (39) takes into account the error in determining the instantaneous frequency caused by averaging over the timewindow of duration A, and it converges to 0 as A + O. On the other hand, the second term, (37), blows up to 00 as A + 0 and has a more fundamental nature related to the uncertainty principle considered in Section 7.2.
Example 3. Current frequency for the simplest tune seen through a Gaussian window. Let us consider in some detail the behavior of the current frequency Winst(r) (29), and the competition between two components of the current frequency deviation D(r) (38), in the case of the simplest tune signal (4) and the Gaussian windowing shape (7c).
203
7.3. Windowed Fourier transform
Taking into account that If(t)1 2 = 1 for the simplest tune, elementary calculations give () Winst S
= WO +
I"'t
a 2 /4 COSS,
(40a) (40b)
and D
0=
1 2A 2 '
(40c)
with the dimensionless parameters
a = VA, Fig. 7.3.3 shows, in cases when a
a
o 
0

S
= vr.
(41)
1 and a = 2, the (dimensionless) time
b
o 
' S
FIGURE 7.3.3 Time dependence of (a) the current frequency, and (b) the frequency deviation, for the simplest tune signal for two values, 1 and 2, of the dimensionless parameter a. dependence of the current frequency Winst (r) measured by the windowed Fourier transform (9), and the instantaneous frequency deviation ainst(r) = (Dinst(r»1/2 due to timewindow averaging. Recall that for a 0, the current frequency Winst(r) coincides with the instantaneous frequency Winst(r). The amplitude of Winst(r) is smaller in the case a = 2 than in the case a = 1 which is a consequence of the smoothing action
=
Chapter 7. Uncertainty principle and wavelet transforms
204
of timeaveraging. Fig. 7.3.3(b) shows that the instantaneous frequency deviation (1inst(r) has peaks at the times when the instantaneous frequency changes quickly, and valleys when it changes slowly. To evaluate the efficiency of instantaneous frequency measurement via the windowed Fourier transform we will consider the ratio of the full frequency dispersion (38) at s = 0 to its limit value DOC = 0 2 /2 for A + 00:
p
2(0)
= D(O) = _1_ + (1 ea2/2)2 Doo
4a 2fJ2
(42)
'
where the new dimensionless parameter
fJ =O/v.
(43)
Note that s = 0 is in a sense the "best" case since at that time the instantaneous frequency changes slowly and D(r) has a minimum. The graphs of dependence of p(O) on parameter a are shown in Fig. 7.3.4 for several values of parameter fJ.
p(O)
2
1
a*
1
2
FIGURE 7.3.4
The graphs of dependence of p(O) on parameter a for several values of parameter fJ. As we explained before, the blowup of the graphs to infinity as a + 0 is due to the uncertainty principle effects which guarantee that the window of a shorter duration will lead to greater indeterminacy of the measured frequency. As a increases, the uncertainty principle effects become negligibly small, but the instantaneous frequency measurement error due to the timeaveraging, increases.
205
7.3. Windowed Fourier transform
The value a* for which P attains its minimal value Pmin determines the optimal duration A* = 2a* /v of the timewindow. It is clear from Fig. 7.3.4 that only if {3 » 1 does Pmin « 1 and, as a result, the windowed Fourier transform algorithm is capable of accurately tracking the instantaneous frequency. The latter values correspond to slow and/or large changes of the instantaneous frequency. • To conclude our windowed Fourier analysis of the simplest tune, a word of caution is necessary. Conclusions based on simple integral characteristics of the joint frequencytime density E(w, 1') (17), such as the mean frequency (23) and deviation (32), can be sometimes misleading. They are quite coarse and lose a lot of information contained in the joint density. So, it is useful to indicate their area of applicability which is luckily possible for the simplest tune signal in view of the relatively simple structure of its windowed Fourier transform jew, 1'). First, note the formula
i asinb = where
=.!..
In(a)
L 00
In(a)e inb ,
(44)
n=oo
(1f cos (a sin x _ nx)dx,
10
1r
(45)
are Bessel functions of integer order n which will be encountered often throughout the remainder of this book. Substituting a = Q/v and b = vt in (44) and multiplying it by e iwot we will obtain for the simplest tune signal (2) the formula
L 00
f(t) =
In(Q/v)ei(wo+nv)t.
(46)
n=oo
Hence, f(t) has the following distributional Fourier image:
L 00
jew) =
I n (Q/v)13(w 
wo 
nv).
(47)
n=oo
Applying it to the definition (9) of the windowed Fourier transform we get, in view of (10), that jew, 1')
=
L 00
v n=oo
I n ({3)go(a(y  n»ei(yn)s,
(48)
where y = (w  wo) / v. It follows from (48) that our integral characteristicsbased analysis of the joint density E(w, 1') is certainly not applicable for a » 1 when
206
Chapter 7. Uncertainty principle and wavelet transfonns
the righthand side of (48) collapses to a sum of separate peaks and E(w, r) (see (17» becomes a polymodal function of variable w.
E
E
a
b
o E
c
FIGURE 7.3.5 Plotsofsimplesttune'sjointdensity E(w, r)incaseof,B = 10 and (a) a = a* = 0.27, (b) a = 1, (c) a = 2. The unimodaJity disappears and multimodaJity sets in at a 1. 7.3.5. Inverse windowed Fourier transform. As for any other integral transform, the question whether there is enough information contained in the windowed Fourier transform j(w, r) to recover from it the original function f(t) is of paramount importance. To answer this question let us keep in mind that the windowed Fourier transform of f(t) is nothing but the ordinary Fourier transform ofthe windowed function f(t)g(t  r). Hence, applying the usual formula (3.2.4) for the inverse Fourier transform we immediately get f(t)g(t  r) =
f
j(w, r)e iwt dw.
(49)
This identity permits the recovery of values of function f (t) only inside the window g(t  r), or in practical terms, only where g(t  r) is not too small.
207
7.3. Windowed Fourier transfonn
To remove this limitation let us select a function h*(t), multiply both sides of (49) by its shift by 'l' and integrate them over all 'l'. The result is Af(t) = / dr / dwj(w, 'l')h*(t  r)i wt ,
(SO)
A = / h*(9)g(9)d9.
(SI)
where
If the auxiliary function h(9) is chosen so that the constant A =f:. 0, then the desired inverse windowed Fourier transform formula takes the form f(t)
1/ / 
=A
.
dwf(w, 'l')h*(t  'l')e'wt •
d'l'
(S2)
Clearly, the above inverse formula is not unique as it depends on the choice of function h(9). For example, if we take h(t 
= a(t 
'l')
'l' 
'l'max),
where 'l'max is the time when g(t) has its maximum value, i.e., g(t) :5 gmax = g('l'max),
then A = gmax and the inverse formula takes the form f(t) = 1 / gmax
. t dw. f(w, t  'l'max)e'W
(S3)
The above nonuniqueness of the inverse formulas is caused by the obvious re dundancy ofthe windowed Fourier transform, which maps function f(t) of a single variable into function j(w, 'l') of two variables. Sometimes, however, such an overdetermination is useful. For instance, if one has to reconstruct the entire signal f(t) from an incomplete information about its transform, the overdetermination present in the windowed Fourier transform can be helpful. Among all the possible inverse formulas (S2) one can try to find the optimal one in the sense that it would maximize the value of coefficient A (important in numerical computations) among all the auxiliary functions h(9), such that
IIhll
=
IIgll·
(S4)
Chapter 7. Uncertainty principle and wavelet transforms
208
We assume that the window g(t) is given and is of finite energy, i.e., both g, h E L2. In terms of Section 7.1, A is just the inner product (h, g) and, by the Schwartz inequality Thus, obviously, the greatest possible value of A = IIgll 2 is attained if h(t) = g(t), and the optimal (in the above sense) inverse formula is
f(t) =
1// 
IIgll 2
dr:
.
dwf(w, r:)g*(t  r:)e ,wt ,
or, in view of (6) and (22),
. . 1 (tr:)/ ,wt f(t)= / dr:)".gO )..dwf(w,r:)e
(55)
The asterisk was dropped because the windowing shape function go(x) was realvalued. 7.3.6. From windows to wavelets. Although windowing was adopted in this section as our favored method, it is only fair to take now the last parting look at the windowed Fourier transform 1 / f(w, r:) = 21l"
f(t)go
(t 
r:) e 1wt . dt, )..
(56)
to assess impartially its merits and shortcomings. First of all, note that the righthand side of (56) contains three free parameters r:, w, and ).., but only two of themthe current time r: and the frequency ware variables. The remaining parameter, the window duration ).., is usually assumed to be constant. This line of thinking is tied to the intuition that the windowed Fourier transform is just a parameterized version of the regular Fourier transform introduced for the purpose of tracking instantaneous frequencies of the signal. Such motivation is, however, also the source of limitations of the windowed Fourier transform and keeping).. constant, restricts our ability to simultaneously analyze the timefrequency properties of the signal. To see more precisely what we mean, let us consider the signal f(t) = fo(tja) of a given shape fo(x) but with variable width governed by the parameter a. Suppose that both fo(x) and go(x) are well localized in the vicinity of x o. Then, for a « ).., approximately
=
j(w, r:)
j(w)go(r:j)..).
209
7.3. Windowed Fourier transform
This means that, for a « J.., the windowed Fourier transform really measures the signal's ordinary Fourier image and is not very good at doing its localizationintimeandfrequency job. The window simply becomes too broad to be of any value. In the opposite case, when the window becomes very narrow (J.. «a), j(w, r)
f(r)g(w)eWT:,
the windowed Fourier transform does an excellent job at timefrequency localization but fails to provide any information about its spectral properties. So, the windowed Fourier transform has a limited applicability field and seems to be most suitable for timefrequency analysis of narrowband signals with phase (t) subject to strong but slow nonlinear timeevolution. The example here is the simplest tune (24) for Q » v (wo » Q). However, most of the signals of interest to scientists and engineers, from cardiograms and seismograms to stock market quotations and turbulent velocity fields, do not resemble simple tunes (Fig. 7.3.1) and have a much richer structure which often includes appearance of the wide range of scales. Tools more flexible than windowed Fourier transform are necessary for their satisfactory analysis, and it is clear from the very beginning that the scale parameter J.. has to be treated as one of the primary variables. This leads to a suggestion of the new signal processing algorithm described by f(J.., r) = A(J..) A
f
f(t)1/I* (tr) J.. dt
which is called the continuous wavelet transform and which takes into account both the location and the scale properties via parameters rand J.., respectively. The shape function 1/1 (x) of the wavelet transform kernel is usually called the mother wavelet. In contrast with the windowed Fourier transform where the window's shape gO plays a minor role, the choice of the mother wavelet is of utmost importance and we will devote a lot of attention to it in the following sections. The reader has probably noticed already that the notion of frequency has been lost in the process. This is not accidental and abandoning the frequency paradigm (closely tied to selecting the mother wavelet containing trigonometric functions e iwt ) in favor of the scale paradigm (which permits full flexibility in selecting the mother wavelet) turns out to be not a weakness but the main strength of the wavelet transforms. The wavelet transforms can be tuned to the peculiarities of a signal we have to work with. If the signal is narrowband then we can take an oscillating wavelet 1/I(x) = e iQx go(x) giving rise to a wavelet transform similar to the windowed Fourier transform but with different emphasis. The wavelet transform acts like a microscope, narrowing the visible area of the signal with the growth of the wavelets frequency w = Q/J...
Chapter 7. Uncertainty principle and wavelet transforms
210
7.4 Continuous wavelet transforms 7.4.1. Definition and properties of continuous wavelet transform. In this section we take a general look at the continuous wavelet transform both theoretically and as it relates to physical and engineering problems. Mathematical questions concerning particular wavelet systems will be dealt with in the last three sections of this chapter. The continuous wavelet image of signal I(t) is defined by
10., r) = A()') A
f
l(t)'I/I* (tr) ). dt,
(1)
where '1/1 (x) is a certain function called the mother wavelet and ). and r are called, respectively, the scale variable and the location variable. Function A()') will be specified later. Note that, to distinguish it from the Fourier transform
I(w)
=
1 21r
f
.
l(t)e 1wt dt,
(2)
the continuous wavelet transform will be denoted by applying a "hat" / to the original function I. Observe that the mother wavelet 'I/I(x) plays the role of complex exponentials exp(iz) in the Fourier transform (then, also A = 1/21r). Varying in (2) the frequency w by compressing or dilating the complex exponential function, we obtain, after integration, a new function j(w) which represents the complex amplitude of the corresponding harmonic component of the original signal:
I(t) =
f
j(w)e iwt dt.
(3)
Consequently, the Fourier image j(w) measures the contribution of different harmonies to the in general nonharmonic signal I(t). A similar compression and dilation of the mother wavelet is accomplished for the continuous wavelet transform by the scaling parameter).. Its exact analog for the Fourier transform is the period T = 21r/lwl. In a sense, one can interpret the value of the continuous wavelet transform /()., r) as a measure of the contribution of the rescaled by ). mother wavelet 'I/I«t  r)/).) to the signal I(t) . The coefficient A()') can be selected arbitrarily as to magnify or reduce sensitivity of the transform to different scales. However, very often it is simply selected as A()') = 1/../i, (4)
7.4. Continuous wavelet transfonn so that
f(J..., r) A
1 =..fi.
211
f
f(t)1/I* (tr) J... dt.
(5)
This choice guarantees that the arbitrary rescaling of the mother wavelet preserves the mother wavelet's L2norm. Indeed,
I
1/1*
C r)
r
=
f 11/12C r) I = f \1/1 dt
2 (x) \ dx
= 1\1/1 (x) 1\2 .
(6) One could say that with this choice of A(J...) all the scales carry equal weight. As we already mentioned in the previous section, for all its great features discussed at length in Chapter 3, the Fourier transform has from the point of view of a physicist one essential shortcoming: its "mother wavelet" exp(iz) has unbounded support. As a result, based on information contained in the Fourier image jew) it is difficult to assess where signal f(t) (or its special features) is located on the t axis and where it is equal to O. In particular, this type of information is totally lost in the "spectral density" \j(w)\2 of the distribution of harmonic components over the frequency w axis. That drawback will be removed in the continuous wavelet transform by selecting a localized mother wavelet 1/I(z) which decays rapidly to zero as Z + ±oo. Consequently, in the continuous wavelet transform, in addition to the scale parameter J..., there appears another primary parameterthe location shift r. Varying it we can track the time t evolution of the "events."
Example 1. Wavelet transform expressed via the Fourier transform. To complete the general picture one should note that the continuous wavelet transform can be, obviously, expressed in terms of the Fourier images of the original function and the mother wavelet:
(7) In particular, substituting the distributional Fourier image of function 1/I(z) = exp(iz), we immediately get
if,(wJ...) = 8(wJ...  1) = 8(w  1/J...)/J..., and setting A(J...) = 1/27r we obtain
1(J..., r)
= j(l/J...) exp(ir /J...).
(8)
•
Example 2. Morlet wavelets. The often encountered in practical application mother wavelet 1/I(z) = eiQzcp(z), (9)
Chapter 7. Uncertainty principle and wavelet transforms
212
with the Gaussian windowing function lP(Z) = exp( _Z2 /2),
(10)
is traditionally called the complexvalued Morlet wavelet (the plot of its real part, for Q = 10, is shown in Fig. 7.4.1). As a result, the Fourier image of l/I(z) is also Gaussian: = _1_ exp (_ (w  Q)2). (11) ../2rr 2 Recall that the Gaussian shape of the windowing function is the minimizer in the uncertainty principle (7.2.4) and, consequently, it optimizes the joint resolution in time t and frequency w. Indeed, the continuous wavelet image
= A(J..)
J
/(t) exp
(
r)2)
Q (t i T(t  r) 2J..2
dt
(12)
contains information about the original (not too fast increasing) function /(t) in
FIGURE 7.4.1 The plot of the real part of the Morlet mother wavelet for Q=10. the window of effective length '" u[lP] = J../../2. Expressing, as in (7), j(J.., r) through the Fourier images of the analyzed functions we get
7.4. Continuous wavelet transform
213
This means that /0., 'l') depends on the values of the Fourier image j(w) in the frequency band of width u [t] = 1/}"./i centered at the frequency Q= Q/}...
(14)
In other words, /(}.., 'l') supplies information about the spectral properties of the original function with resolution 1/}.../i. The arbitrary parameter Q entering in the definition (9) of the Morlet wavelet could be called the efficiency factor of the Morlet wavelet since the quantity Q/27r is of the order of the number of Morlet wavelet's periods contained in its window. • Just as the first automobiles of the last century took inspiration from and mimicked the horsedrawn carriages, and only later developed their own identity, the wavelets underwent a similar evolution which started with their identity as "im proved" versions of the Fourier transform and only gradually developed into being recognized for their own outstanding capabilities. These capabilities, still far from being fully tapped, are related to the fact that the mathematical theory of wavelets, as we will see later on, imposes very few restrictions on the choice of the mother wavelet's shape. We will illustrate them on concrete applications in the rest of this section. One of the powerful applications of the continuous wavelet transform is the study of open and hidden singularities in the incoming signal f(t). Usually, the singularities are caused by physical (biological, economic, etc.) laws, whose validity the experimenter is trying to confirm, or come from the existence of the sharp boundaries between the regions where the process f(t) evolves smoothly. The mother wavelets that are useful in this context are quite unlike the Morlet wavelet (910). Example 3. Mexican hat wavelet. Differentiation can bring to the surface function's hidden singularities. For this reason one often selects mother wavelets so that the corresponding continuous wavelet transform converges, for }.. 0, to a desired derivative of the function being analyzed. One of such examples is the Mexican hat (15)
which is just the second derivative of the Gaussian function (10). Its Fourier image is (16)
Chapter 7. Uncertainty principle and wavelet transfonns
214
FIGURE 7.4.2
The Mexican hat mother wavelet. Substituting (15) into (1), and integrating by parts twice, we get
I().., 7:) = A()"»)" 2 A
It is customary to select
I (t qJ
7:) )..
d 22 I(t) dt. dt
A()") = Ij)..3./2ii,
so that j().., 7:) converges, for).. + 0, to exactly
(17)
(18)
1"(7:).
•
Another property of the continuous wavelet transform essential to understand its mechanism is based on the Schwartz inequality (7.1.19)
II
I(t)g*(t) dtl2
I
I/(t)1 2 dt
I
Ig(t)1 2 dt,
(19)
which applied to the function
g(t) =
1 (t7:) ../I1/! )..
(20)
yields the inequality (21) The inequality provides an upper bound on possible values of the modulus of the continuous wavelet transform (5) of 1 (t). Let us assume, without loss of generality,
215
7.4. Continuous wavelet transform
6
4 2
0
2
4
6
FIGURE 7.4.3 The greyscale plot of the wavelet image of f(t) = exp( It I) in the case of the Mexican hat mother wavelet. The horizontal axis represents the Tvariable and the verticalAvariable. The greyscale level changes from black to white as the values of the wavelet image increase. The black oval spot in the lower middle portion ofthe plot is a consequence of the singularity of the original function's second derivative at t = O. that both the signal and the mother wavelet are normalized so that II!II = 111/111 = 1. It is clear that the maximum values are achieved, and the inequality (21) becomes an equality, if the original function f(t) is equal, for certain A = AO and T = TO, to the wavelet f(t) =
(t 
1 TO) . 1/1 
...;;.:0
AO
(22)
Informally, we can say that the continuous wavelet transform is best tuned to, or resonates with signals that have shapes similar to that of the mother wavelet. Note that the more complexstructured the mother wavelet (20) and the resonating signal (22) are, the more pronounced the above resonance property of the corresponding continuous wavelet transform is. To make things a bit more formallet us define the signal as complexstructured if its time (7.2.2) and frequency (7.2.3) localizations satisfy the "strong" uncertainty principle: (23)
216
Chapter 7. Uncertainty principle and wavelet transforms
Example 4. Complexstructured signal. Let us consider signal f(t) whose Fourier image is the familiar Gaussian function
_
Jrl/4
({J}
f(w) =   exp   ( 1 .j2Jrft 2ft2
+ iy)
),
(24)
where y is a real number and the constant ft (with the dimension of frequency) has the meaning of effective width of the Fourier image. The coefficient in front of the exponential function has been selected so that the normalization condition
is satisfied. With the help of the integral formula (3.2.3), we compute the inverse Fourier transform f(t) = Jrl/4.j1
+ iy exp
(t + 2(1
2ft2) iy) .
(25)
In tum, using the integral formula
we find the frequency (7.2.3) and time (7.2.2) localizations of the complexvalued signal (25): 2 
a [f]
2
ft = 2'
212 2ft2 (l y ).
a [f]
=
+
(26)
Substituting these expressions into (25) we get the following condition for the signal f(t) (25) to be complexstructured:
y»1.
(27)
•
Remark 1. To better see reasons why signal (25) turned out to be complexstructured let us write the complexvalued Fourier image i(w) of an arbitrary signal f(t) in the exponential form i(w) = A(w) exp(i (w»,
(28)
where A(w) = li(w)1 is the nonnegative amplitude and (w)the real phase of the complex Fourier image i(w). The amplitude and phase of the Fourier image
217
7.4. Continuous wavelet transfonn (24) of signal (25) are A(w)
1l'1/4 ) = exp (w2  , 2
,J21l' /.L
yw2 cI>(w) = 2/.L2.
2/.L
(29)
The complex structure of signal (25) was conditioned on the fast nonlinear variation of the phase of the Fourier image (24) as a function of w. Indeed, according to (3), the signal can be written in the form
00), the value of Employing the stationary phase method, asymptotically (y signal f(t) at a given instant t is determined by the integral contribution in the small neighborhood of the stationary point, in our case n = 2t /.L 2 / y. Substituting here, instead of n, the effective width /.L of the Fourier image we shall find the effective duration of the complexstructured signal: (T /.L
»
1).
Remark 2. The approximate estimate of the signal (25) duration obtained above via the stationary phase method may seem unnecessary at first sight since we already know the exact form of the signal and the exact formula for its time localization: (30) Nevertheless, the above argument has a heuristic value, emphasizing the principal role of the phase in complexstructure signal formation. It also shows a universal method of calculation of its form and duration. Example 5. Complexstructured mother wavelet. As another example of mother with the complexstructured signal wavelet let us take function y,(z) f(z) (25). The continuous wavelet image f().., 1') of function f(t) to which the mother wavelet is perfectly tuned is K()", 1')
1 =,JI
f
f(t)/*
(t1') )..
dt.
(31)
Recall that the form (5) of the continuous wavelet transform selected here guarantees that, for any).., the normalization condition
218
Chapter 7. Uncertainty principle and wavelet transforms
is satisfied. Notice that we also introduced special notation K (A, 1') for the special continuous wavelet image of the mother wavelet itself. Function K (A, 1') is sometimes called the wideband ambiguity function of the mother wavelet and it plays an important role in wavelet theory. In terms of the Fourier images
K(A, 1') =
2lr,JI!
j(W)j*(wA)eiwT: dw,
(32)
so that substituting (24) we obtain
K(A, 1') =
1IT! ;v ; exp (1 Zp wJL22 + iw1' )
dw,
(33)
where (34)
Finally, evaluation of the integral (33) gives
The above function has a maximum at following dependence on 1 :
l'
= 0, and its modulus square has the
It is natural to interpret function I (A) as a sort of resonance curve which characterizes the response efficiency of the continuous wavelet transform as a function of the scale parameter 1. Fig. 7.4.4 shows graphs of function I (A) for signals of different complexity, as measured by parameter y. It is clear from the illustrations that the resonance is best emphasized for large values of y, that is for signals of large complexity. The maximal value of I (1) is achieved for A = 1. It is related to the fact that for A = 1 function (31) becomes the autocorrelation function
K(1') =
!
f(t)f*(t  1')dt
of the original signal. The autocorrelation function has some remarkable properties. In particular, it transforms any signal, however complex, into a simple signal whose Fourier image, (36)
7.4. Continuous wavelet transform
219
I 1
0.8 0.6
2
1
3
4
5 A
FIGURE 7.4.4 Graphs of function I (A) for different complexitystructure of signals as measured by y. is real and nonnegative, with the phase (w) == O. In electrical engineering one often says that all the harmonics of the autocorrelation function K (t) have identical phases. The autocorrelation function achieves its maximum at r = 0 and decays relatively rapidly as Irl increases. In particular, it is easy to see that its localization properties are determined by u[K] =
1
v:; ,
ILv2
so that, in view of (30), it is clear that for y » 1 the autocorrelation function K(r) = K(A = 1, r) is much better localized on the raxis than the original signal f(t) (25) on the taxis. •
7.4.2. Inversion of the continuous wavelet transform. As for any other integral transform the basic question is: Does the continuous wavelet image j (A, r) contain sufficient information permitting recovery of the original function f(t)? In more practical terms: Does there exist an inversion formula for the continuous wavelet transform? To answer these questions let us multiply equality (1) by r)/A) and integrate over all r. The result is the auxiliary integral
y,«8 
I(A,8) =
f
f(A, A
r)y, (8r) A dr.
(37)
Chapter 7. Uncertainty principle and wavelet transforms
220
Equivalently, I(A,B)=A(A)! dtf(t)!
(38)
It is easy to see that the inner integral can be expressed via the autocorrelation function (35) K(z) =
!
(39)
1/I(s)1/I*(s  z)ds
of the mother wavelet as follows:
! (
t1') A
1/1*
1/1 (B1') A
d1' = AK (Bt) A .
(40)
As a result, I(A,B)=AA(A)!
dt.
(41)
To solve this integral equation for f(t) let us multiply (41) by a function B(A), to be selected later, and integrate over all A:
10
00
I (A, B)B(A) dA =
where g(s) =
10
00
!
f(t)g(B  t) dt,
K(S/A)C(A) dA,
(42)
(43)
and C(A) = AA(A)B(A).
(44)
Clearly, the righthand side of (42) would be reduced to f (B), thus solving equation (41) for function f(t) if g(s)
= 1000 K(S/A)C(A)dA = a(s).
(45)
Let us find C(A) for which the distributional equation (45) is satisfied. Remembering that the Fourier image of the autocorrelation function K(z) is 2:7r1..fr(ev)1 2 , we get the equation (46)
7.4. Continuous wavelet transfonn
221
equivalent to equation (45). To eliminate the dependence of the above integral on lJ) we shall select C(J..) so that
1/ DJ...
J..C(J..) =
(47)
In this case, (46) becomes (48) where
(49) is the normalizing constant that can be calculated from (48) by introducing the new variable of integration" = lJ)J.. to get (50)
Putting together (44), (47) and (49) we get that
1 B(J..) = DJ..3 A(J..) ,
(51)
so that, from (51) and (42),
2. (JO I (J.., 9) dJ.. D 10 J..3A(J..)
_ / 9  ().
Substituting expression (37) for I(J.., t) we finally obtain the inverse continuous wavelet transform
1 /(t) = D
10roo J..3dJ.. A(J..)
f
(t  1') d1' /(J.., 1')t/F J.. . A
(52)
In particular, if the continuous wavelet transform is defined by (45), then the inversion formula takes the form
1 /(t) = D
dJ.. 10roo J..2.ji
f
d7: /(J.., 7:)t/F A
(t7:) J.. .
However, the above inversion formulas require several caveats.
(53)
Chapter 7. Uncertainty principle and wavelet transfonns
222
Remark 3. First of all we have to admit that in the process of making our calculations transparent we cheated quite a bit. The observant reader would have noticed that the passage from (48) to (50) is justified only if is an even function. For that reason formulas (5253) are valid only for twosided mother wavelets, as mother wavelets with even square modulus Fourier image are called. To this class belong all the purely realvalued mother wavelets such as the Mexican hat (1516). On the other hand, the complexvalued Modet wavelet (910) is not of this type. For that reason mathematicians often work with onesided mother wavelets whose Fourier image is
== 0,
w :::
o.
(54)
For such mother wavelets, instead of (46) we have the equality g(w)
10 = 27r D 0
00
dA A

lo/(WA)1 2 
1 = x(w). 27r
(55)
To explain its consequences let us express the righthand side of (42) in terms of the Fourier images j (w) and g(w ):
f
I (A, (})B(A) dA = 27r
f
j(w)g(w)e iw8 dw.
Substituting here (37), (51) and (55), we arrive at the relation that replaces equality (52) for onesided mother wavelets:
f
.()
1 f(w)x (w)e'W dw = D
f
10roo A3dA A(A)
A
dT f(A, T)o/
T)
(0 A
.
As we have shown in Section 6.6, the Fourier integral on the lefthand side is, up to coefficient 1/2, equal to the analytic signal
10
2 (OOdA F(t) = D A3 A(A)
f
dT f(A, T)o/ A
(OT) A
(56)
corresponding to the original signal f(t). Remembering that the real part of the analytic signal coincides with f(t), we arrive at the inversion formula for the continuous wavelet transform for onesided mother wavelets:
10
2 (OOdA f(t) = D Re A3 A(A)
f
dT f(A, T)o/ A
((}T) A .
(57)
7.4. Continuous wavelet transform
223
Example 6. Poisson wavelets. As an example of onesided mother wavelets consider l/Im(Z) = (1  iZ)ml, (58) m >0, which are called Poisson wavelets. Their Fourier images (59)
can be calculated by means of the residues method to be

l/Im(w)
= r(m1+ 1) w m ew x(w).
(60)
Poisson wavelets can be used to identify open and hidden singularities of signal f (t) and, for m = 2, like the Mexican hat, in the search for edges between different regimes of the original function f (t). Indeed, for m = 2, the Poisson wavelet
1 dZ
zz
l/Iz(z) = 2d
Its real part Re l/Iz(z) =
1
1 .• IZ
1 dZ
(61)
1
2 dzz 1 + zZ
has a shape similar to that of the Mexican hat and possesses, for A differentiating properties.
0, the same •
Remark 4. The above derivation indicates that the necessary condition for the existence of the inversion formulas is finiteness of coefficient D (49), or equivalently, the inequality
f
<
00.
(62)
Mother wavelets satisfying condition (62) are called admissible wavelets. The complexvalued Modet wavelet (910) is not admissible since its Fourier image does not vanish for w = 0 and, consequently, the integral (62) diverges. Nevertheless, in practice this is not a serious obstacle. First of all, the inversion formula is not always needed, and, secondly, for sufficiently large values of the "goodness" parameter Q, the Fourier image of the Modet wavelet takes a very small value at w = 0, and it is not difficult to adjust it a little bit to make it admissible.
Remark 5. It follows from condition (62) that the Fourier images of admissible wavelets satisfy condition = 0) = 0, which is equivalent to the condition
f
l/I(z)dz = O.
(63)
224
Chapter 7. Uncertainty principle and wavelet transforms
This in tum implies that any admissible wavelet has to have a oscillatory (signchanging) naturethis provides a partial explanation of the term wavelet. Remark 6. The fact that the continuous wavelet image j (). , r) of the function f(t) of a single variable depends on two variables indicates that the continuous
wavelet transform contains redundant information and is overdetermined. One of the consequences of this fact is that the inversion formula is not unique. This is easily seen by multiplying (1) not by the mother wavelet, as was done earlier, but by another function q;
)..
,
and integrating the resulting equation over all r. For the sake of simplicity of the argument let us assume that q;(z) has a onesided Fourier image, that is ip(w) == 0 for w :s O. As a result we get an analog of expression (41), the only difference being that instead of the autocorrelation function (39) one enters the crosscorrelation function (64) K (z) = q;(s )1{1* (s  z) ds
f
whose Fourierimage is 21fip(w)1fr*(w). Replacing 11fr(w) 12 by ip(w) 1fr* (w) in all the preceding formulas, we arrive at an infinite variety of continuous wavelet inversion formulas:
21
f(t) = Re D
00
d).. 3
o )..
A()")
f
(or)
dr f().., r)q;   , A
)..
(65)
which are all well defined as long as (66) The, complex in general, coefficient D entering in formula (65) is equal to
Notice that the condition (66) can be fulfilled not only for admissible mother wavelets but also if ip(w) converges sufficiently fast to zero as w + 0+. This means that the inversion formula (65) remains valid also in cases when the formulas (5253) do not make sense. For example, formula (65) permits recovery of the original function f(t) from the continuous Morlet wavelet image j().., r) (910) if one takes q;(z) to be the Poisson wavelet (58). The question of how to make wavelet transforms more economical and less redundant, while preserving their good scale and time localization properties is
225
7.5. Haar wavelets and multiresolution analysis
a subtle mathematical problem. For the usual Fourier transform and the Fourier series the lack of overdetermination, and the uniqueness of the inverse Fourier transform (or Fourier coefficient sequences) is guaranteed by the Hilbert space L 2 orthogonality of complex exponentials (or trigonometric functions) on the interval [0, 21l'], that is, by the condition that if
m#n.
The mathematically difficult task of constructing orthogonal wavelet systems will be discussed at some length in the next three sections.
7.5 Haar wavelets and multiresolution analysis In this section we will take a look at a special (one can say digital) series representation for realvalued signals in terms of the socalled Haar wavelets. This idealized system provides a good easy introduction to the concepts of wavelet transforms and multiresolution analysis. Each term of the expansion will provide information about both the time and the frequency localization of the signal. The Haar wavelets will be obtained from a single prototypea mother waveletby translations in time and frequency, although the explicit shift in frequency will be
1
",(1)
'II (t) 1,3
0.5 1
3
t
0.5 1

FIGURE 7.5.1 The Haar mother wavelet, and a wavelet of order (1,3). replaced by a more natural in this case dilation (rescaling, stretching) in time. This
Chapter 7. Uncertainty principle and wavelet transfonns
226
will guarantee that all the wavelets have the same shape. To eliminate redundancy and overdetermination, we will make the wavelet system orthogonal. The Haar mother wavelet is defined as follows
1,
l/I(t) = { 1,
0,
°
for t < 1/2; for 1/2 x < 1; otherwise.
(1)
The Haar wavelet
of order (m, n), m, n = ... , 1,0,1, ... , is obtained by rescaling (dilating or compressing) the time in the mother wavelet l/I(t) by a factor of 2m and then translating the resulting wavelet by an integer n multiplicity of 2 m • The dilation makes the wavelet l/Im,n(t) fit in the interval of length 2 m , and the translation places its support finally in the interval [2 m n, 2 m (n + 1)] (see Fig. 7.5.1). We will call parameter mthe level of resolution of the wavelet, and parameter nthe location parameter of the wavelet. Then the number 2 m can be seen as its resolution, and 2mnas its location. The coefficient 2m/2 in the definition (2) was selected to make all the Haar wavelets normalized in L2(R), that is, to compensate for the dilation operation to guarantee that
(3) It turns out that:
The system of Haar wavelets
is orthogonal, that is (l/Ij,k, l/Im,n)
=
f
m, n = ...  2, 1, 0, 1,2, ...
(4)
= 0,
(5)
l/Ij,k(t)l/Im,n (t) dt
if (j, k)
and complete in L 2(R). The latter means that any function L 2 convergent representation
f
=
00
f
oF (m, n), E
L2(R) has an
00
L L
m=oon=oo
wm,nl/lm,n
(6)
7.5. Haar wavelets and multiresolution analysis
227
where, in view of the orthonormality, the expansion coefficients Wm,n
= wm,n[f] = 2 and coefficients ak (by (7.5.11), N = 2 was necessary and sufficient for Haar wavelets). Then the mother wavelet can be selected to be N
t(t) = L(1)kaNkq;(2t  k), k=O
and the corresponding wavelet system can be built with its help via formula (7.5.2). Such an approach was suggested by Ingrid DAUBECHIES in 1988, and the resulting wavelets are called Daubechies wavelets. Conceptually, the above construction is a clearcut generalization of the construction of Haar wavelets from the scaling function I provided in the previous section. However, for N > 2, the selection of coefficients ak becomes highly nontrivial. Also, as a rule, the smoother one wants the wavelets one wants, the larger the N one has to take. Below, we provide a sketch of the relatively simple construction of continuous Daubechies wavelets which is due to David POLLEN (1992). Their scaling function q;(t) satisfies the scaling relation q;(t) = aq;(2t)
+ (1 
a)q;(2t  1)
where
a=
+ (1 
a)q;(2t  2)
1 +.J3 4
+ aq;(2t 
3),
(1)
(2)
and where, for real numbers of the form at + f3.J3 with ( dyadic) rational at, f3, the overline indicates the "conjugation" operation
The support of the resulting q;(t) is contained in the interval [0, 3] and, additionally,
L 00
q;(k) = 1.
(3)
k=oo
Assume that there exists a scaling function q;(t) supported by [0,3] and satisfying (1) and (3) for integer values of the argument t. The scaling relation (1) written
7.6. Continuous Daubechies' wavelets
233
for t = 0, 1,2,3 becomes a matrix equation
CP(O») ( cp(l) cp(2) cp(3)
=
(a 0 0 0) (CP(O») 0 a0 0 a (1  a) 1  a a _0_ 0 (1  a) 1  a
cp(l) cp(2) cp(3)
which, in view of condition (3), has exactly one solution: cp(O) = 0,
cp(l)
=
1+y'3 2
'
1 y'3
=
cp(2)
2
'
cp(3) = O.
Starting with these prescribed values and using the scaling relation (1) one can produce values of the scaling function cp(t) for any dyadic rational t. For example, cp(lj2) = 2 +4y'3,
cp(3j 2) = 0,
cp(5j2) = 2  y'3, 4
and so on. The values of cp(t) for dyadic t are clearly of the form a + fJy'3 with dyadic a and fJ. One can also prove (see the Bibliographical Notes) that they also satisfy two extended partition of unity (see also (3» formulas
L 00
cp(t  k) = 1
k=oo
and 00
(3 2y'3 +k) cp(tk)=t.
Since the support of cp(t) is contained in [0,3] the above properties also give the interval translation properties for dyadic t E [0,1]: 2cp(t) + cp(t + 1) = t +
1+y'3 2 '
2cp(t + 2) + cp(t + 1) = t + cp(t)  cp(t + 2) = t
+
3y'3 2 '
1 +y'3 2
.
Chapter 7. Uncertainty principle and wavelet transforms
234
Combining them with the scaling relation (1) gives the scaling relations for dyadic [0,1]: O+t) = aq;(t); q; ( 2
t E
1+t) q; ( 2 = aq;(t)
+ at + 2+v'3 4 ;
2+t) _ q; ( 2 = aq;(1 +t) +at 3+t) _ q; ( 2 = aq;(1
q;
t)
4+ ( 2
+ t) 
v'3 + 4;
at
1 + 4";
_+ 3  2v'3
= aq;(2 + t)  at
4
(4) ;
5 +t) = aq;(2 + t). q; ( 2
Compared with the original scaling relation (1), they have a clear advantage: the values of q;(t) at the next resolution level depend only on one value at the previous resolution level (instead of four in (1». The above formulas form a basis for the following recursive construction of the continuous version of the scaling function on the whole interval [0,3]. Start with function 80(t) which is equal to q;(t) at integers 0,1,2,3, and which linearly interpolates q; inbetween these integers. Clearly, 80(t) is continuous. In the next step, form 81 (t) at the second resolution level by applying the (righthand sides of) scaling relations (4) to 80. More precisely, for t E [0, 1], define O+t) = a80(t); 81 ( 2
1+t) 81 ( 2 = a80(t) 81 (
2+t)
2
3+t)
81 ( 2
81
t)
+ at + 2+v'3 4 ;
_ v'3 = a80(1 +t) +at + 4; _ = a80(1
+ t) 
4+ ( 2= a80(2 + t) 
at
1 + 4";
_+ 3  2v'3
at
4
;
7.6. Continuous Daubechies' wavelets
gl (
5;
t)
235
= ago(2 + t).
=
Outside [0, 3] set gl (t) O. Function gl (t) is continuous and coincides with q>(t) at dyadic points with resolution 2 1 (inbetween, it again provides a linear interpolation). Continuing this procedure we obtain a sequence gn of continuous, piecewise linear functions (zero outside [0,3]) which agree with q>(t) at dyadic points of the form k2 n •
1.25 1 0.75 0.5
......
; 0.25 :
0.5
"
1
3
0.25 FIGURE 7.6.1
Values of the Daubechies' scaling function computed at dyadic points t = n . 26 ,0 ::: t ::: 3, via the scaling relation (1). Notice that functions Ign (t) I ::: 3 for all n = 1, 2, ... , and since 0 ::: (see (2», we get that
la I ::: a
< 1
max Igk(t)  gk+j(t) I ::: a k max Igo(t)  gj(t)1 ::: 6a k . t
t
Hence the sequence of functions gk(t) satisfies uniformly the Cauchy condition, and the limit q>(t)
= lim gn(t). n+oo
is a continuous function. This is the scaling function we were searching for. Remark 1. Note that the scaling function q>(t) is not differentiable because _ I· a i q>(1) _ Ii (2a)j (1)Iim q>(2i)  . q>(0) _ I·1m q>(2i) .  lID . m q> 
i+oo
2J
since 2a > 1 and q>(0)
i+oo
# O.
2J
j+oo
2J
j+OO
00,
236
Chapter 7. Uncertainty principle and wavelet transforms
With some additional work one can now establish that
f
q>(t) dt
= 1,
and that the integer translations of q>(t) form an orthonormal system, that is
1.5
l
.il
:v'
1
0.5
__ i
___0.....
...
..
/ 1.5'.
0.5
.
2
:#
3
.,.. ..
'" ."
1
:
.
2.5
.. .....
FIGURE 7.6.2 Values of the Daubechies' mother wavelet computed at dyadic points t = n . 26 , 0 ::::: t ::::: 3, via the formula (5). Following the general scheme explained in detail for the Haar wavelets in Section 7.5, we can now define the mother wavelet ""(t) via equality ""(t)
= iiq>(2t) + (1 
a)q>(2x  1)  (1  a)q>(2x  2)
f
and check that
+ aq>(2x 
""(t) dt = O.
The integer shifts of the mother wavelet are orthonormal, that is
f
""(t)""(t  k) dt
= {O, 1,
=F 0; If k = O.
3), (5)
7.7. Wavelets and distributions
237
Moreover, the scaling function ({J and the mother wavelet t/F are orthogonal as well, that is
f
({J(t)t/F(t  k)dt = 0.
Thus, again, by an argument similar to that used for the Haar wavelets, the set of Daubechies wavelets t/Fm,n(t) = 2m/ 2 t/F(2m  n),
m, n = ... 1,0,1 ... ,
forms an orthonormal complete basis in L 2 (R), and so does the set of functions ({In(t) = ({J(t  n), t/Fm.n(t)
= 2m/ 2 t/F(2m 
n),
n
= ... , 1,0,1, ... ,
m = ... 0,1,2 ... ,
n
= ... , 1,0,1, ....
7.7 Wavelets and distributions The scaling relation (t) = 21/2 Lak(2t  k)
(1)
k
can also have distributional, rather than just function solutions (following our convention we denote distributions by capital letters). In the most trivial case when (1) has only one nonzero term, sayao = ../2, such a solution may be guessed immediately: (t)
= 8(t),
since 8(t) = 28(21). However, in the case of the scaling relation (1) with finitely many (but at least two) nonzero coefficients ak, the scaling distribution can no longer be a linear combination of Dirac deltas, that is, of the form n
(t) = L
Ck8(t  tk)'
(2)
k=O
This can be seen as follows. Suppose that the nonzero coefficients are ao and aI, and perhaps some others. Since the support of must be contained in the interval
238
Chapter 7. Uncertainty principle and wavelet transforms
[0, 1], we may take 0 S to < t1 < ... < tn S 1 in (2). The scaling relation (1) then forces equation n 1 2 n [ Ck ( t (;CkO(t  tk) = 2/ {; aoZo
tk ) "2
1) + ....]
Ck o ( t  2tk + +a1 Z
Comparing coefficients for the same Dirac deltas on either side of the above equation, we get in the case to = 0 that Co = 2
1/2
Co
aoZ '
and since ao =1= 21/2 we get that Co = 0 . If to > 0 then
so that again Co = O. By the same argument we get that C1 = C2 = ... = Cn = 0 because tk/2 < tk < (tk + 1)/2. So, there are no solutions of the form (2). The situation is better if infinitely many nonzero coefficients are permitted in the scaling relation (1), and we will indicate some avenues that can be pursued in such a case. One option is to seek a solution E S' of (1) which is a distribution with compact support. Then its Fourier transform (w) is an analytic function in the entire complex plane C, and the scaling relation (1) translates into the following relation for (w): (w) =
r
1/ 2
L ak exp( iwk/2) (w/2). k
We, however, will not follow this route, and instead will construct the scaling distribution (and the corresponding wavelets) by demanding that its integer translates form an orthonormal basis at the zero resolution level. A mathematical aside: orthogonality of distributions. 2To make a rigorous discussion of the orthogonality of distributions possible we have to select an inner product of (at least some) distributions. One such possibility is the inner product (T, S)
=
2JT
f
T(w)S*(w) dw, w2 +1
2This material may be skipped by the first time reader.
(3)
239
7.7. Wavelets and distributions
which is well defined for all the distribution in the Sobolev space
H. :=
{f e S' : 11/11
2
=
f
dw < 00 } ,
(4)
which is a subspace of the space of tempered distributions S'. Oearly, the Dirac delta B(t) e H. since its Fourier transform is identically equal to 1/21r. Now, our job is to construct a multiresolution analysis of the Sobolev space H., and the first step is to find a scaling distribution 4>, the integer translates thereof would be a orthonormal basis at the zero resolution level, that is for
Vo
= {f e H. : supp feZ}.
Distributions in Vo are of the form T(t) =
L
CkB(t  k),
(5)
k
Any distribution in H. can be then approximated by partial sums of dilations of (5). The first difficulty one encounters is that the usual dilations do not preserve the norm 11.11 in H.. So to give ourselves some leeway, we will allow use of other inner products in H., and introduce the family of inner products
(T S) , a
= 2. 21r
f
T(w)S*(w) d
w2 + a 2
w,
(6)
parameterized by parameter a > o. They all generate norms 1I.lIa equivalent to the original norm 11.11 = 11·111, and the orthogonality notions for all of them are equivalent. Their role will become clear later on. If 4> is of the form (5), then the orthonormality condition gives that
Thus Bok are Fourier coefficients (in L 2 (R» of a function appearing as a fraction in the last integral, which therefore must be equal (in L 2 (R» to its Fourier series
240
Chapter 7. Uncertainty principle and wavelet transfonns
Since ci> has period 21r, the desired scaling function in Vo C
:i..
_ [
... (w, a} 
L n
1i has the Fourier transform
1/2
1
2
(w+21rn) +a
2
_ [2a(COSha  COSW}]1/2 • sinha
]

(7)
That the sum of the above infinite series is as indicated can be seen as follows (see also formula (8.7.7). By expanding exp(ixt) on [0,1] in the Fourier series, and evaluating it at t 0, we get that
=
L
1
00
n=oo
x
so that a(w,a} :=
L n
i
= 4a Ln
(w
1 cos(x /2}
+ 21rn = 2 sin(x/2}
1
(w + 4"'n}2
l 
+ ai}/2 + 21rn
+a2
l
L n
(w  ai}/2
+ 2",n
i [COS(W + ia}/4 cos(w  ia}/4] sin(w + ia}/4  sin(w  ia}/4
= Sa
sin(ia/2)
1
sinha/2
= 4", cosa/2 
cosw/2·
The Fourier transform ci> given by (7) cannot be simply inverted. However, it turns out that the construction of the corresponding mother wavelet distribution 111 leads to an invertible Fourier transform. Indeed, since
1I1(t}
=L
dk ..ti (that is, orthogonal to all integer translations of the Dirac delta) we obtain the following conditions on the Fourier transform Iji:
7.7. Wavelets and distributions
241
In view of the scaling relation (8), If, has period rr, and both series can be simplified by separating odd and even terms, which leads to equations
iii (w)a(w, a) + liI(w + rr)a(w + 21r, a) = O. Their easily identifiable solution is ,i,( ) _ iw/2 ( ... w,a e
a(w + 2rr) ) 1/2 a(w, a) (a(w + 2rr, a) + a(w, a»
. 2 =e1w/ (cosha/2 
cosw/2)
( 4a
.sinha
)1/2 •
(taking just the positive root does not work here as it did in (7), since it does not satisfy the second of the above pair of equations).
For ex = 2 m (and this is just what we need for the construction of a wavelet basis in 11) the easily evaluated inverse transform of "'(ro, ex) gives the following formula for the mother wavelet w(t, 2 m ) = C [2coSh(Tm1)8(t 1/2)  8(t)  8(t 1)] for a certain constant C. Now, we are finally ready for the construction of the orthonormal wavelet basis in the Sobolev distribution space 11, using the mother wavelet 'II as the starting point. Here, the advantage of using the family of norms 1I.lIa becomes obvious. As we observed earlier, the classical rescaling
does not preserve the norm 11.11, but
Therefore, the correct rescaling in 11 which produces the first level resolution wavelet is and the complete multiresolution wavelet basis for 11 is provided by the integer translation of the scaling distribution n(t) = (t  n),
n = ... , 1,0,1, ... ,
242
Chapter 7. Uncertainty principle and wavelet transforms
and the wavelets Wm,n(t) = 2 3m / 2 W(2 mt  n, 2 m ),
n = ... , 1,0,1, ... ,
m = 0, 1, ... ,
or, alternatively, by the wavelet system
complemented by linear combinations of the Dirac delta
Remark 1. One can also expand distributions in weakly convergent series of smooth function wavelets. In this case, the typical result is that a tempered distribution which is the rth derivative of a function (measure) of polynomial growth has a wavelet expansion with the coefficients an = O(lnl k), for some integer k. More precisely, assume that the scaling function q; is in Sr, that is it has r continuous and rapidly decreasing derivatives. In other words, 1q;(k)(t)1 ::::: C p ,k(1
+ Itl)P,
k = 0, 1, ... , r,
pEN, t E R.
The spaces Sr contain the space S of rapidly decreasing test functions. Then the mother wavelet 1/1 is also in Sr, and for any tempered distribution T E S; (of order r) the expansion with respect to q;(t  n) and 1/I(t  n) exists and the coefficients an = O(lnlk) for some integer k. As a result, for T E S;, we have the usual wavelet expansion
where
and where the convergence is in
S; or, alternatively,
T(t) = I>nq;(t  n) + n
L
bm,n1/lm,n(t).
7.B. Exercises
243
7.8 Exercises Function spaces:
=J
1. Prove the inequality P If(t)g(t)1 dt I/fl/I/gl/ which was used, among others, in the proof of the uncertainty principle in Section 7.2. Compare with the Schwartz inequality (7.1.9). 2. Prove the triangle inequality
L2.
1/
f + g 1/
1/
f 1/ + 1/ g 1/ for the functional Hilbert space
Windowed Fourier transform: 3. Let j(w, r) be the windowed Fourier transform of the signal f(t). Denote by f'(t). Express j' in terms of
j'(w, r) the windowed Fourier transform of the derivative j. 4. Signal x(t) is a solution of the differential equation dx(t)
;It
+ hx(t) =
f(t),
where f(t) is a signal with known windowed Fourier image j(w, r) and h is a (real or complex) constant. Express the windowed Fourier image of x(t) in terms of j(w, r). 5. Let j(w, r) be the windowed Fourier image of signal f(t) (i.e. f(t) j(w, r». Find the windowed Fourier images of signals (a) f(t)e iWOI , and (b) f(t + (}). 6. Find the windowed Fourier transform
j
(w, r) of function
f
(t) = e VI.
7. Utilizing results of the previous exercises, find the windowed Fourier image of the signal f(t) = e VI coswot. 8. Assume that the values of the windowed Fourier image j (w, r) of signal f (t) are known outside the interval r E [0, T] only. Is it possible, on the basis of this incomplete information, to recover values of signal f (t) for all values of t?
Wavelets: 9. Provide the expression for coefficient D in (7.4.49) if one employs the more popular in mathematical literature definition (7.3.15) of the Fourier transform. 10. The distributional relation 1
1
00
K
(S)  d)" =8(s)
Do)..
)..2
(1)
played" the principal role in derivation of formulas for the inverse continuous wavelet transform. In practice it is impossible to carry out the integration all the way to ).. = o. Show that the regularized integral
1100
=D
K
(S) i d)" )..2
Chapter 7. Uncertainty principle and wavelet transforms
244
converges weakly to the Dirac delta as E + O. 11. Find an explicit formula for function cI>(x) = «6  x 2 )/(8.,fii» exp( _x 2 /4) (see answer to the Exercise 10) in the case of the Mexican hat mother wavelet (7.4.15). 12. Obtain, in the case of twosided mother wavelets, a formula connecting
1/12, analogous to the Parseval formula for the ordinary Fourier transform.
1/12 and
10(>..,
13. Denote by r) the continuous wavelet image of signal I(at), a > 0, compressed (a > 1) or dilated (a < 1) in comparison with the original signal I(t). Find out how (A, r) is related to the continuous wavelet image of signal I(t) itself in the case of wavelet transform definition (7.4.5).
10
14. Find I(A, r) (7.4.1) for the selfsimilar signal I(t) = Itla.
15. Let I(t) = (5  9t
+ 4t 2 )/(5 12t + 8t 2 )3.
Find numerically and graphically the Haar wavelet expansion of I(t) with resolution level coarser than 0 and finer than 6. Graph the resolution level n contents of I(t) for n = 0, I, ... ,6. Use your computer in order to estimate numerically the maximum error of your approximation. 16. Use your computer and the defining scaling relations to produce numerical values of the Daubechies scaling function and mother wavelet at the dyadic points up to resolution level 6.
ChapterS Summation of Divergent Series and Integrals
The theory of distributions has a flavor similar to the theory of summation of divergent series and integrals and, as we have seen in Chapter 6, is closely related to the theory of singular integrals. A generic alternating series 1  1 + 1  1 + ... =
L ±1
is a good example here. It seems to make no sense to assign a specific value to this infinite sum. Nevertheless, mathematicians have produced certain reasonable rules of summation that assign to it value 1/2. Such an assignment is in complete agreement with the intuition of physicists who encounter similar series. In this chapter, we will see how one can sum this, or even more strange, divergent series and integrals. To gain a better insight into the essence of this problem, let us begin with elementary examples and recall basic notions and theorems of the ordinary theory of convergent infinite series.
S.l
Zeno's "paradox" and convergence of infinite series
8.1.1. Geometric series. Recall the celebrated 4th century B.S. Achillesandtortoise "paradox" due to the Greek philosopher Zeno. At time t = 0 Achilles is at point x = 0 and the tortoise at x = 1. Achilles begins to chase the tortoise with velocity +1, and the tortoise begins to run away with velocity 0 < v < 1. To catch the tortoise Achilles has to first reach point x = 1. This will happen at time ro = 1. By that time the tortoise will have moved by distance vro = v. To cover that distance Achilles needs time rl = v, during which the tortoise will have moved further to the right by distance vri = v 2 (see Fig 8.1.1). The pattern continues ad infinitum, seemingly indefinitely delaying the moment © Springer Nature Switzerland AG 2018 A. I. Saichev and W. Woyczynski, Distributions in the Physical and Engineering Sciences, Volume 1, Applied and Numerical Harmonic Analysis, https://doi.org/10.1007/9783319979588_8
246
Chapter 8. Summation of divergent series and integrals x
11(1 v)
FIGURE 8.1.1 A graph representing paths of Achilles chasing the tortoise. They will meet at the point of intersection of the corresponding straight lines. when Achilles catches the tortoise. Of course, time t needed by Achilles to catch the tortoise is the sum of an infinite series t = iO
where
im
+ i1 + i2 + ... ,
= v m , m = 0, 1, 2, ... , so that t =
1+v
+ v 2 + v 3 + ... =
L 00
vm ,
(1)
m=O
is the sum of a geometric series. Leaving aside philosophical significance of Zeno's "paradox" concerning the nature of space and time, let us underline the basic mathematical difficulty en countered in formula (1): we are trying to perform infinitely many mathematical operations, each of which takes a certain amount of time. At a naive level, this seems to be an impossible task if we want to do it exactly. An insight into how this difficulty can be overcome is gained by a look at the problem from the physical viewpoint. To find time t when Achilles catches the tortoise, it suffices to solve the system of equations of motion for both of them:
x = vt x
+ 1,
= t.
(2a) (2b)
8.1. Zeno's "paradox" and convergence of infinite series
The solution is
1 t=,
247
(3)
1v
which, together with formula (1), gives
Lv (Xl
m
m=O
1 =.
(4)
1v
Equation (4) attaches to the infinite series on its lefthand side, a number from the analytic expression on its righthand side. Note that the latter has a well defined mathematical and physical meaning, even if the lefthand side does not form a convergent series. For example, for v = 1, equation (4) states that 1  1 + 1  1 + ... =
1
2'
(5)
an equality anticipated in the preamble to this chapter. The question is: How can we provide a mathematical justification of this type of identities? The answer is given by the theory of summability of divergent series which formalizes procedures that are consistent with the physical principle of infinitesimal relaxation. On the other hand, we'll be able to elucidate rigorous mathematical constructs of summability theory by looking at their physical roots. Observe that series (1) which appears on the lefthand side of equality (4) arises in an attempt to solve equations of motion (23) by the method of consecutive approximations in parameter v. 8.1.2. Criteria for convergence. The above example illustrates the basic difference between finite sums (6)
and the infinite series
(7) The value of a finite sum can always be explicitly computed (at least in principle), whereas it is not always possible to assign a numerical value to an infinite series. However, for some series, the sequence of partial sums
(8) converges, as n
00,
to a finite limit S = lim Sn, n+(Xl
(9)
248
Chapter 8. Summation of divergent series and integrals
and then it is natural to think about S as the sum of the infinite series, the value of which is being approximated by computable finite partial sums. Such series are called convergent. One learns in calculus that it is possible to find out whether the series converges or not without actually computing the limit of its partial sums. One such approach is based on the socalled Cauchy criterion which states:
Series S converges if and only if for any given number E > 0, one can find a integer N = N(E) such that, for any n N and arbitrary positive integer k = 1,2, ... , (10)
If the series converges, then its remainder
L 00
Rn = S  Sn =
am
m=n+l is a well defined number, and its absolute value measures the error of approximation of the infinite series S by its finite partial sum Sn. In particular, if n N (E) then IRnl < E.
Example 1. Geometric series. Let us return to the geometric series
L qm = 1 + q + q2 + .... 00
S=
m=O
Its partial sums Sn obviously satisfy identity
from which we immediately get that 1 qn+1 Sn = "
1q
Hence, the partial sums sequence Sn converges if and only if converge, its limit
Iq I <
1. If it does
1 1 q
S= lim Sn =   . n+oo
The remainder Rn is also easily computable:
1 qn+1 Rn =    Sn =   .
1q
1q
•
S.l. Zeno's ''paradox'' and convergence ofinfinite series
249
A comparison of terms of an arbitrary series with corresponding terms of the geometric series gives a handy sufficient condition of convergence of general series.
If, for all sufficiently large m, we have that lam 111m < q < 1, then the series L::=o am converges. If a series does not converge then we call it a divergent series. One obvious example is the geometric series with Iq I 1.
Example 2. Harmonic series. Consider the harmonic series 00 1 1 1 H= L   = 1 + 2 ++ .... m=om +1 3
As we have observed in Section 4.1, its partial sums Hn satisfy an asymptotic relation (n + 00) (11) Hn '" In(n + 1) + y, which immediately proves the divergence of the harmonic series. Note, however, that if we alternate signs of the harmonic series to obtain
•
00 (_l)m 1 1 L = L   = l   +   ... , m=O m + 1 2 3
then the latter series converges. This is due to the general phenomenon which makes all the alternating series of the form 00
D = L(l)mcm , m=O
converge. The corresponding formal statement is known in calculus as the Leibniz Theorem. The alternating harmonic series is an example of a convergent series Lm am for which the series of absolute values Lm lam I diverges. Series Lm am which converges together with the series Lm lam I of absolute values is called absolutely convergent. Obviously, any geometric series with Iq I < 1 converges absolutely. 8.1.3. Conditional and absolute convergence. Summation of convergent infinite series is an associative operation. This means that if we group the terms of a convergent series into finite blocks, then the sum of the whole series is equal to the sum of the series consisting of these blocks. More precisely, the associativity means that if we select an increasing sequence of integers
o < ko < kl
< k2 < k3 < ... < km < ...
250
Chapter 8. Summation of divergent series and integrals
and denote the finite sums
sums in finite blocks generated by the above sequence, then the series 00
S'= m=O
also converges and S = S' . This immediately follows from the fact that the ...• • ... } forms a subsequence of the original sequence of partial sums {So' sequence of partial sums {So. S1. ...• Sn •.. .}. In contrast, summation of convergent series is not necessarily a commutative operation. A series formed by a permutation (perhaps infinite) of terms of a convergent series may converge to a different limit, or may even diverge. If a series converges to the same limit after any permutation of its terms, then we call it unconditionally convergent. It turns out that the necessary and sufficient condition for the unconditional convergence of a series is its absolute convergence. If a series converges, but not absolutely, then in view of the above statement, one says that it converges conditionally. For a conditionally convergent series, the Riemann Theorem asserts that any number is the sum of a certain permutation of the original series. This striking theorem will not be proved here (see Bibliography for the proof), but we will illustrate it by a concrete example. Example 3. A conditionally convergent series. Select positive integers p and q, and consider a permutation of the alternating harmonic series in which q negative terms follow p positive terms:
,
1
11
1
1
1
L =1++ ... +... ++ ... +4p 1  .... 3 2p  1 2 2q 2p + 1 Consider partial sums L containing the first n (p + q) terms ofthe rearranged series. Since finite summation is commutative, npl
= ] ; 2m Now,
1
+1 
nql
1
] ; 2(m + 1)
S.l. Zeno's "paradox" and convergence of infinite series
and
1
nql
1 nql
2(m+1) =
251
1
z]; m+1'
which implies that the sums in the above expression for can be expressed through partial sums of the harmonic series. In view of the asymptotic relation (11), as n 00, 2npl
1
L "" In(2np) +y, m +1
m=O
npl
1
L "" In(np) + y, m+1
m=O
so that,
,
1
Ln(p+q) "" In(2np)  Zln(np ) 
Hence, finally, = In
1
Zln(nq ),
(n
(2m·
00).
(12)
The reader can complete the example by proving that not only subsequence but also the entire series L' converges to the same limit. It is clear that, if we choose different integers p' and q', thus selecting a different rearrangement of the alternating harmonic series, then we'll get different limits as long as p / q i p' / q'. In particular, for p = q = 1, we get the well known formula
1 L = 1 
2
1 3
+   ... =
In2
but for p = 2, q = 1, we obtain that
, 11111 L = 1+    +  +   32574
+ ... =
3  In 2. 2
•
The above example gives a taste of the proof of the Riemann Theorem. Although the latter seems to contradict common sense, it is sometimes usedalthough without crediting Riemann by dishonest financiers taking new loans to repay growing old debts. The same argument explains illusory prosperity of declining nations which issue new paper money to cover inflationary budgets.
252
Chapter 8. Summation of divergent series and integrals
8.1.4. Functional series and uniform convergence. In what follows we will often work with infinite series
L am (x), 00
S(x) =
(13)
m=O
where the terms are not numbers but functions. For such series the notion of pointwise convergence, that is the notion of convergence for every point x separately, is not sufficient. For that reason one introduces a stronger notion of uniform convergence which will be discussed below. Assume that series (13) converges for each x E [a, b]. Then, by the definition of pointwise convergence, for any x E [a, b] and for any given number e > 0, one can find an integer N = N(e, x) such that, for any n N(e, x) and for an arbitrary positive integer k = 1,2, ... , ISn+k(X)  Sn(x)1 < e. If we can make the selection of N = N(e, x) independent of x E [a, b] (one says "select N uniformly over [a, b]"), then the series is said to converge uniformly over [a, b]. Notice that the condition of uniform convergence is equivalent to the condition that n(e) = max N(e, x) < 00. xe[a.b]
The importance of uniform convergence becomes obvious when we try to investigate properties of functions am (x) that are inherited by the sum S (x) of the whole series. It is easy to see that, in general, the continuity of am (x) does not imply the continuity of S(x) if the series (13) converges pointwise. Indeed, if Sn (x) = x n , x E [0, 1], then although the corresponding terms am(x) = Sm(x)  Sml(X) = xm  x m 1 are continuous, the limit S(x) which equals 0 for 0 ::: x < 1, and 1 for x = 1, is a discontinuous function. However, the uniform convergence of the series precludes the above situation, and we have the following theorem: If functions am (x), m = 0,1,2, ... , are continuous on [a, b] and the series Lm am (x) converges uniformly on [a, b1 then S(x) is continuous on that interval.
Similarly, pointwise convergence is not sufficient to permit differentiation and integration of the functional series term by term or, in other words, interchanging the order of infinite summation with operations of differentiation and integration. However, in the presence of uniform convergence we have the following two useful theorems: If functions am (x), m = 0,1,2, ... , are integrable on [a, b] and the series Lm am (x) converges uniformly on [a, b1 then
8.2. Summation of divergent series
253
Assume that functions am (x), m = 0,1,2, ... , are defined on [a, b] and are continuously differentiable on (a, b). Then, if the series Lm am (xo) converges for an xo E [a, b1 and the series of derivatives Lm a:" (x) converges uniformly on (a, b), then
for any x E (a, b).
Checking the uniform convergence of a functional series is not always easy but there exist several criteria that can be helpful. One of the most useful is the following Weierstrass criterion.
°
Let Cm > be a sequence of numbers such that the series Lm Cm converges. If functions am (x), m = 0,1,2, ... , are defined on [a, b1 and satisfy, for all m = 0, 1,2, ... , inequalities
then the series Lm am (x) converges uniformly on [a, b]. Example 4. An application of the Weierstrass criterion. Let Iq I < 1. The above Weierstrass criterion, combined with properties of the geometric series, immediately yields the uniform convergence of the functional series 00
Lqmcos(mx) m=O
on the whole real line.
8.2 Summation of divergent series In this section we will study the main topic of this chapter: methods of summa bility of divergent series (1)
We have encountered divergent series before. One example was the series L ±1, which can be obtained by specifying am = (_1)m in expression (1). We also discussed physical situations, where divergent series arise. The solution of equa tions of motion (8.1.2), with v = 1, by the method of successive approximation,
Chapter 8. Summation of divergent series and integrals
254
provides such an example. In mathematics, the study of divergent series is also motivated by other reasons as for example to deal with the fact that the product of two convergent series can be a divergent series. Quite reasonably, one would like to assign to this product a value equal to the product of sums of the two factor series. The basic idea, underlying some of the most useful methods of summation of divergent series, is contained in the relatively innocentlookingAbel's Theorem: If series
am converges to sum S then, for any 0 < q < 1, the power series 00
(2)
S(q) = Lamqm m=O
also converges and, additionally, lim S(q) = S.
(3)
However, sometimes it happens that the series L am diverges, but the power series (2) converges for any 0 < q < 1, and the sums S(q) have a finite limit as q + 1. In such a situation, S can be defined by formula (3), and is called the generalized sum of series (1) in the PoissonAbel sense. Often, the above approach is just called the Abel method of summation, although Abel himself never worked on the summability of divergent series. Example 1. Summing L ±1 by the Abel method. For the series auxiliary power series (2) has the sum S(q)
00
1
m=O
1 +q
= L(I)mqm = ,
Iql
L ±1, the
< 1.
As q + 1  0, the above sum S(q) has a finite limit 1/2. Thus, the Abel method leads to the familiar formula (8.1.5)
1  1 + 1  1 + ... = 1/2.
(4)
•
Another approach, known as the Cesaro method of generalized summation, depends on the formation of arithmetic averages Ao = So,
SI Al = So+ 2 , ... ,
So + SI + ... + Sn An = .....:.......:..., ... n+l
(5)
8.3. Tiring Achilles and the principle of infinitesimal relaxation
255
of partial sums of the original series (1), and on the investigation of their convergence. If the sequence Sn itself converges, then the sequence of Cesaro averages An always converges to the same limit. But sometimes An converges even though Sn diverges. In this case, we shall say that the original series (1) is Cesaro summable and the limit of the sequence An is called the generalized sum of (1) in the Cesaro
sense. Example 2. Summing L ±I by the Cesaro method. In the generic example of the series
L
±I, the partial sums are S2n A
_ k+I +I,
2k 2k
= 1,
S2n+1
A2k+1
= O. Consequently,
1 =2
so that Ak 1/2. Thus, the Cesaro method leads to the same generalized sum as • the Abel method of summation. There exists a large number of other summability methods. To make them useful, mathematicians usually demand that they satisfy regularity and linearity conditions. By definition, a summability method is called regular if it assigns to an already convergent series a generalized sum equal to its usual sum S. The method is called linear, if a generalized sum for the series L(pam +qbm) is equal to pS + qT, provided series Lam and Lbm have generalized sums S and T, respectively. Note, that the regularity of the Abel method follows from the Abel Theorem. Its linearity is obvious. It is also not difficult to prove regularity and linearity of the Cesaro method. Naturally, not all series have finite generalized sums. For example, series 1 + 1 + 1 + . . . has infinite Abel and Cesaro generalized sums, that is, they are not Abel or Cesaro summable. Summable series with alternating signs are sometimes called semiconvergent series.
8.3
Tiring Achilles and the principle of infinitesimal relaxation
As we already noticed, properties of finite and infinite sums can be drastically different. Just recall the Riemann Theorem which implies that the addition in conditionally convergent infinite series is not commutative whereas in finite sums it is. For semiconvergent series addition is not associative either; such series are sensitive not only to rearrangements of their terms, but also to their grouping. For example, for the series ±I, we have
L
(1  1) + (1  1) + (1  1) + ...
= 0 + 0 + o... = 0,
(Ia)
Chapter 8. Summation of divergent series and integrals
256
but 1  (1  1)  (1  1)  ...
=1 
0  0  ...
= 1.
(lb)
A comparison of this result with the result of Abel (or Cesaro) summation of the same series (see (B.1.5» may give you second thoughts about the summability methods introduced in the preceding section. One can argue that (la) and (lb) taken together, although ambiguous, make more sense than the answer 1/2 obtained before, since one would like the sum of integers to be an integer as well. Similar misgivings can help motivate the search for additional summability methods based on physical arguments. To illustrate what we have in mind we shall return to the Achilles and tortoise example. However, this time we will closely watch Achilles' consecutive steps along the xaxis, rather than just the overall time he needs to catch the tortoise. Since Achilles' speed is equal to 1, his total displacement is given by series (B.1.1) with t replaced by x:
=L 00
x
vm
= 1 + v + v 2 + ...
(2)
m=O
Let us imagine the chase as a series of physical stages, each consisting of Achilles reaching the previous position of the tortoise and resting for an instant (as Zeno prescribed originally), before embarking on the next stage. In contrast to formula (B. 1. 1), series (2) has a physical meaning also for negative velocities v of the tortoise. If v < 0, the tortoise always runs towards Achilles. For example, if v = 1, Achilles reaches point x = 1 at the end of the first stage having met the tortoise on his way and, at the same time, the tortoise halts at x = o. Then, Achilles turns around, and at the end of the second stage finds himself at point x = 0, while the tortoise has returned to point x = 1, and so on. The graph of Achilles' motion is presented in Fig. B.3.1a. Now let us strip Achilles of his semigod status and assume, realistically, that he is getting tired chasing the tortoise back and forth. As his, now human, strength is sapped, after each tum he slows down by a factor of q. The graph of Achilles' motion is presented on Fig. 8.3.1b, and his full displacement is given by the sum of an absolutely convergent series
x
1 = L..,.,(1)m q m = . m=O
1 +q
Note that in the case q + 1  0 of infinitesimally slow relaxation of Achilles' strength, the above sum is equal to 1/2; the infinitesimally tiring Achilles will end his chase in the middle of the interval originally separating him from the tortoise. The calculation of Achilles' total displacement, provided above in the case of infinitesimal relaxation of his speed, can be applied to an arbitrary series (B.2.1)
8.3. Tiring Achilles and the principle of infinitesimal relaxation
1
257
x
.5
v
5
V
V V V
10
x
15
20
V 25
t
b
0.5
x

1
0.5
c

25
FIGURE 8.3.1
The graphs of motion of tiring Achilles.
in which case it coincides with the Abel generalized summation procedure (8.2.2) and (8.2.3). The Cesaro method can also be "supported" by a similar physical argument. Suppose that the never tiring Achilles covers the distance am in each stage, but his tiring alter ego decreases his speed linearly (as opposed to the exponential decay seen in the above example related to the Abel method). In this case
Vm
m n+1
=1
Chapter 8. Summation of divergent series and integrals
258
and he comes to a full stop only after n steps. In this case, his total displacement
xn =
m=O
( m) = + + ... + Sn
am 1  n +1
So
S1
(3)
n +1
is equal to the arithmetic mean (8.2.5) of partial sums of series (8.2.1). For am = (_l)m and n 00 (the case of infinitesimal relaxation of Achilles' speed), we see that Xn 1/2. In other words, the Cesarean Achilles will eventually collapse from exhaustion along the Abelian Achilles.
8.4 Achilles chasing the tortoise in presence of head winds The principle of infinitesimal relaxation provides a physical interpretation for the equality 11 + 11 ... = 1/2 (8.1.5), and also for the general framework of Abel and Cesaro summation methods. However, a more detailed analysis of that principle reveals that the ambiguity in assigning values to the series 11 +11 ... and to other semiconvergent series, has not been entirely removed. In the following example, a regular linear method of generalized summation assigns to the series 11 + 11 ... values 1/2 (as in (8.2.4», 0, or 1 (as in (8.3.1», depending on the selection of infinitesimal relaxation rates. Example 1. Summation with variable relaxation rates. Achilles chases the tortoise in asymmetric conditions that slow him down at different rates depending on whether he is running up or down the xaxis (think about the wind blowing in the negative direction along the xaxis). As a result, Achilles' speed up the xaxis decreases by a factor of p at each stage, and his speed down decreases (increases) by a factor q, p =f. q. During the mth stage his speed V2k
= (qp)k,
V2k+1
= p(qpl,
k
= 0, 1,2, ....
(1)
Thus, the original divergent series 1  1 + 1  1 ... is replaced by a series 00
x = L(l)m vm ,
(2)
m=O
°
which, in view of the Cauchy criterion, converges absolutely for < pq < 1. Since the terms of an absolutely convergent series can be rearranged without changing their sum, 00
00
x = L(qpl p L(qpl, k=O
k=O
8.4. Achilles chasing tortoise in presence of head winds
259
where we grouped together the terms of series (2) with even and odd indices, respectively. Formula (8.1.4) for the sum of a geometric series yields 1p
x=. 1qp
Now, following the principle of infinitesimal relaxation, we let qp + 1  O. That can be accomplished in different ways, and it turns out that the final answer is sensitive to the choice of rates at which q and p converge to 1. If, for instance, p = qT, then
and in the limit q + 1, applying L'Hospital's rule we get that
r x=.
r+1
(3)
For r = 1, that is when p = q and the conditions of running to the right and to the left are the same, we get that x = 1/2the result obtained by the Abel and Cesaro methods. However, if we let r ::/= 1 vary then the limiting x can take other values as well. To see this, it is useful to distinguish between the cases of twosided and onesided relaxations. Ifr > Oandq < 1,thenp = qT < 1andtheAchilles'speeddecreases independently of whether he moves to the right or to the left. In this case we say that a twosided relaxation takes place, and the Achilles' total displacement x, given by formula (3), can take any value from the interval (0,1). For r + 0, we have x + 0 which corresponds (see (8.3.1a» to the rearrangement (1 1) + (1  1) + ... = 0 of the series 1  1 + 1  1 + .... For r + 00, we have x + 1 which corresponds (see (8.3.1b» to the rearrangement 1  (1  1)  (1  1)  ... = 1. Note that the condition qp = q1+T < 1, equivalent to the absolute convergence of series (1), can also be fulfilled in the case when p > 1 and p < 11q. This is the model of onesided relaxation corresponding to the values r < 0 and Achilles running in presence of head winds blowing in the direction of the negative x axis. Then, the limiting x < O. The result has a transparent physical meaning. Running to the left Achilles goes further than running to the right. The accumulating drift of infinitesimal displacements leads to the full displacement which can be anywhere in the interval (00,0). In this way the limit value x = 00 can be interpreted as achievable as a result of a perturbation of the rearrangement 1(11)  (11) +... of the series 1  1 + 1  1 ... , in which each expression in parentheses is replaced by a very small number e, so that e + e + ... = 00. •
260
8.5
Chapter 8. Summation of divergent series and integrals
Separation of scales condition
Examples provided in Section 8.4 demonstrated the existence of regular, linear methods of summation satisfying the principle of infinitesimal relaxation, which give sums different than those provided by the Abel and Cesaro methods. In other words, the principle of infinitesimal relaxation by itself does not produce a unique sum for a semiconvergent series. In this context it seems desirable to seek additional natural conditions which would eliminate the above ambiguities. In this section we introduce the separation ofscales condition which will restrict unwanted arbitrariness in the generalized methods of summation discussed so far. Roughly speaking, the condition requires that the relaxation of semiconvergent series terms should proceed at a much slower rate than the rate of internal oscillations around zero of the series itself. Without explicitly formulating the scales separation condition we will proceed to discuss it on a revealing example. Consider the (fJsummation method which associates with the series L am an auxiliary series
=L
00
Srp(8)
a m (fJ(8m),
(1)
m=O
where (fJ(x), x > 0, is a certain function. If series (1) converges for 8 > 0 and its sums have a finite limit (2)
then Srp is taken as a generalized sum of series L am. To make the (fJsummation method work, function (fJ(x) will be assumed to be continuous and satisfy the following conditions: • "Sufficiently" many derivatives (fJ'(x), (fJ"(x), ••. , are continuous and integrable over the halfline (0, (0); • cp(O) = 1; • For "sufficiently" large n = 1,2, ... , lim xn(fJ(x)
x.+oo
= O.
(3)
The conditions were deliberately stated in a somewhat vague form, but their specific role will become clear later on. The method is obviously linear. Its regularity is assured by the continuity of (fJ(x) at 0 which yields limx.+o (fJ(x) = 1. Function (fJ(x), which completely determines the generalized cpsummation method, describes the relaxation law. Condition (3) demands that cp(x) rapidly decreases to 0 at infinity, thus guaranteeing the relaxation principle. The separation of scales condition is fulfilled: inasmuch 8 + 0, the multiplier cp(8m) varies slower and slower as m increases.
8.5. Separation of scales condition
261
Example 1. Abel and Cesaro methods as rpsummation. Formula (1) generates the Abel and the Cesaro methods as special cases. Indeed, taking rp(x) = exp( x) we arrive at the auxiliary series 00
Srp(l) = L
(4)
am (e&)m ,
m=O
which gives formula (8.2.2) of the Abel method after substituting q = exp( l). On the other hand, if we select rp to be the triangular function rp(x) = {1  x,
0,
for 0 < x < 1; for x 1,
(5)
and replace l) + 0 by a sequence l)n = 1/(n + 1), n + 00, then we arrive at the Cesaro summation method. Either function fulfills the conditions imposed on function rp above and, as a result, Abel and Cesaro methods satisfy the condition of separation of scales. • As promised at the beginning of this section, the condition of separation of scales guarantees uniqueness of the sum of series L ±1. More precisely, the generalized sum (3) of series L ±1 does not depend on the selection of the relaxation function rp(x), as long as it satisfies the above listed conditions. To see this, substitute am = (_1)m in (1), and group the terms in pairs to get 00
Srp(l) = L[rp(2kl)  rp(2kl)
+ l)].
(6)
k=O
In view of the Mean Value Theorem, for a continuously differentiable I(b)  I(a) = I'(c)(b  a)
I, (7)
for a certain c E (a, b). Thus, sum (6) can be rewritten in the form
The above series is an approximate sum, with the partition points Xk = 2kl) +Ck, for the integral of function rp' (x) over interval (0, 00). Hence, the assumed integrability of the derivative rp' (x) implies that
Chapter 8. Summation of divergent series and integrals
262
The factor 1/2 is due to the length of the partition intervals being 2c5. Finally,
Srp
11 = &.0 lim Srp(8) = 2 0
00
q/(x)dx
1 1 = , 1 = rp(x) 2 0 2 00
which proves our assertion that under the scales separation condition all generalized rpsummation methods (13) give the same sum 1  1 + 1 1 ... = 1/2.
Remark 1. Divergent improper integrals. In a calculus course, the usual approach is to approximate integrals over finite intervals by finite sums. Extension of the rpsummation method of approximation to integrals over infinite intervals requires certain precautions, and it would be a good exercise for the reader to provide a rigorous proof here, and find assumptions that have to be imposed on functions rp(x) and rp'(x) to make the approximation work in this case. Remark 2. Rapidly divergent semiconvergent series. Note that our proof of the validity of the generalized summation method (13) for the series L ±1 depended only on the integrability of the first derivative rp' (x) over the interval (0,00). This condition was obviously satisfied by the relaxation function rp(x) = exp(x) of the Abel method, as well as the triangular function (5) of the Cesaro method. However, the stronger the growth of oscillations (as m + 00) of the terms of a semiconvergent series, the stronger smoothness conditions will have to be imposed on the relaxation function rp(x). For that reason the Cesaro method does not work well for series diverging stronger than the series L ±1. Example 2. A rpsummation for a rapidly diverging series. Consider series 00
L(l)mm = 1 +2  3 +4 ... m=1
The corresponding auxiliary series in (1) is 00
Srp(8) = L(1)mmrp(8m) = m=1
where g(x)
1
00
L(1)mg (8m), m=O
= xrp(x). As before, we will pair up the terms of this series to get Srp(8) =
1
00
L[g(2k8)  g(2k8
+ 8)].
(8)
m=O
But in this case, the factor 1/8 makes the Mean Value Theorem useless in determining the behavior of Srp(8) as 8 + O. What is needed is a more precise information
263
8.5. Separation of scales condition
contained in the Taylor formula which states that if function f (x) has a continuous derivative of order n + 1 then, for some C E (a, b), f(x) =
Ln
m=O
f(m)() a (x  a)m m!
+
f(n+1)( ) (n
c (x  a)(n+l).
+ I)!
(9)
For n = 1 we get, in particular, that g(2k8)  g(2k8
for a Ck
E
+ 8) =
g'(2k8)8  g"(ck)8 2/2,
(2M, 2k8 + 8) which in combination with (8) gives (10)
As a simple consequence of Taylor's formula (9)
l
a
a +A
f(x)dx = f(a)l:!.
applied to function f(x)
1
=
g'(x), a
1
+ 2· f'(c)l:!.2,
=
2k8 and l:!.
(2k+2)8
g'(x)dx = g'(2k8)28
2k8
C E
(a, a
+ ll.),
=
28, we get that
+ g"(ek)28 2,
(11)
for some ek E (2k8, (2k + 2)8). This equality can now be used to eliminate g' (2k8) from (10) and arrive at
Observe that, for any 8 and for a suitable selection of ek's and Ck'S, the above formula is exact. The integral
10
00
g'(x)dx
= 1:=0 xq>(x)
in view of assumptions on function q>. Integrability of q>"(x) assures then the convergence, as 8 + 0, of the remaining two sums to the corresponding integrals,
264
Chapter 8. Summation of divergent series and integrals
so that finally
.
hm Srp(8)
8+0
11
=
4
00
0
1, 1
= g (0) = .
g " (x)dx
4
4
Thus the principle of infinitesimal relaxation combined with the separation of scales condition permitted us to arrive at a striking result
1 2+3 4+
... =
1/4.
(12)
•
In some physical situations the separation of scales condition is not satisfied. That was the case in the analysis of Achilles chasing tortoise in the presence of head winds. Another, and more convincing example of violation of the scales separation condition will be encountered later on when we study the summability of semiconvergent integrals. However, in physical applications, such cases are rare and the condition is widely applied.
8.6
Series of complex exponentials
Many problems of the theory of divergent (and also convergent) series can be solved by a study of the complex exponential series 00
L e imz ,
(1)
m=O
where z = x
+ iy.
For y > 0 the series converges absolutely and imz m=O
1
= 1 eiz'
In addition, for y + 0+ , the expression on the right gives the Abel summation formula of a divergent series
Loo m=O
1e imx    1 e ix '
(2)
8.6. Series of complex exponentials
265
which, for x = 7r, reduces to the familiar equality 1  1 + 1  1... = 1/2. Separating the real and the imaginary parts in (2), we arrive at the well known formulas of generalized summation: if x# ±27rn, n = 0,1,2, ... ,
L 00
cosmx = 1/2
(3a)
m=O
and
. m=O
x
1
(3b)
= cot. 2 2
If x = ±27rn, n = 0, 1, 2, ... , then the series (2) becomes a realvalued series 1 + 1 + 1 + ... which diverges to +00, and formulas (3) can be complemented by formulas
L 00
cosmx = +00,
m=O
00
LSinmx =0,
x
= ±27rn, n = 0, 1,2, ....
m=O
A more detailed investigation of the character of singularities of series (2) in the neighborhood of these points will be pursued at the beginning of the next section. For y > 0, the absolute convergence of the series obtained by termbyterm differentiation of series (1) justifies formulas
for n = 1,2, .... If, for y 0+, the righthand sides of these equalities converge to finite limits, then we will take them as generalized sums of the corresponding divergent series on the lefthand side. Thus, putting x = 7r, we obtain that
where
2: + 2:1 tanh (Y) 2 .
1 1 f(y) = 1 + eY =
The Taylor expansion of the hyperbolic tangent function is of the form 1
(Y)
2: tanh 2
=
6
22k  1
(2k)! B2kY
2kl
,
lyl <
7r,
Chapter 8. Summation of divergent series and integrals
266
where coefficients Bn entering in the above formula are called Bernoulli numbers and can be determined from the Taylor expansion
x x   = eX  1 2
x x + ooth= 2
2
Bk
n
L..J  x . k=l k!
In particular,
B6
1
= 42'
Bs
1
=  30'
5 691 BlO = 66' B12 =  2730' B14
7
= 6""
Hence, 00
L(1)mm 2k
= 0,
(4)
m=O
For k = 1, (4) gives equality 1  2 + 3  4 + ... we get that
= 1/4
(8.5.12), and for k
=2
m 1 3 1 L.,,(1) m =18+2764+"'=8'
m=l
Differentiation of series (1) led to the generalized summation formulas (4). The integration of the similar series (5)
also gives rise to useful relations. Multiplying both sides of (5) by a function! (y) and integrating them termbyterm with respect to y over (0, 00), we get
L." m=l where F(t) =
F(m)ei(ml)z =
J: !(y)etYdy.
foo
100
!(y)d! , eY  e'X
8.6. Series of complex exponentials
267
In particular, fot' f(y) = ys1, we get that F(t) = r(s)t S, where r(s) is the gamma function (4.4.3). Hence, we obtain the equality
1
L mS e 00
1
00 1 y s1dY = r(s) 0 eY _ eix .
i(m1)x
m=1
(6)
For x = 0, the above expression gives the well known formula for the Riemann Zeta function s(s) =
and, for x =
1f,
,?; 00
1 1 =S m r(s)
1
00
yS1dy
0
eY
1

,?;
s _
00
(_1)m1 _ _ 1_
 r(s)
mS
yS1dy
(00
10
eY + l '
If we substitute s = 1 in (6), then we get
m=1
1 imx eix m 
1
00
0
where we introduced a new variable integral, we finally obtain that 1 imx
e
s> 1,
we get
'1( )
m=1 m
,
z
1
ix
dy eY  eix 
= 1
s > O.
11
dz 1eiz Z
eixeY • After evaluation of the
1
= In(1  e ) = In 2 2(1  cos x)
i + (1f 2
x),
0< x < 21f.
To conclude this section, we derive another useful formula by integrating (5) with respect to y over (p, 00), p > 0, and then substitute x = O. This gives that
L
00
m=1
mp
_e
m
= In(1 + coth(p /2» 
In 2,
p > O.
For small p this series can be interpreted as "quasiharmonic" in which the contribution of large terms is damped by the exponential multiplier. In particular, for p 0 we get the asymptotic formula
L _e00
m=1
mp
m
'" In(1/p),
Chapter 8. Summation of divergent series and integrals
268
which indicates that the main contribution to the sum is made by the first N terms.
1/P
8.7 Periodic Dirac deltas In this section we will consider the infinite series of Dirac deltas and the functional series
= 1+2 L
00
U(x, y)
emy cos(mx)
= 2Re
L eimz 1, 00
z
= x + iy,
(1)
m=O
m=l
which will play a key role in our analysis. For y > 0, the latter converges, and following Section 8.6, one can find its sum to be _ sinhy U(x,y ) coshy  cosx
(2)
Note that for y 0+ the above function gives the Abel generalized sum of the divergent trigonometric series
U(x)
= U(x, y = 0+) =
L cos(mx). 00
1+ 2
(3)
m=l
The same equality may be rewritten in the complex form to get
L m=oo 00
U(x) = U(x, 0+) =
=
e imx .
=
(4)
Setting x 1r in (3), we get a numerical series U(1r) 1  2 + 2  2 + ... , and we already know (see (8.2.4)) that its Abel sum is U(1r) = 1  1 = O. The same answer U (x) == 0 is obtained for any other x if we put formally y = 0 in (2). However, a closer inspection of the righthand side of equality (4) indicates that for x = 2n1r, n = 0, ±1, ±2, ... , the limit U(2n1r) = +1 + 1 + 1 + ... = 00. The situation is elucidated by Fig. 8.7.1, where the graph of function U(x, y) is shown for small y > O. It hugs the xaxis everywhere except for small neighborhoods of points x = 2n1r, where it has sharp peaks.
269
8.7. Periodic Dirac deltas
U(x,y) 10 8 6
O 0, by its very structure, 4> (x, y) is an even, continuous and periodic function of x with period 2rr, which can be expanded into a Fourier series
+L
00
4>(x, y) = co
em cos(mx).
(6)
m=l
The zeroth coefficient
1
Co = 2 rr
l1r 4> (x , y)dx = 21 l1r L ( 1r rr 1r n=oo 00
X 
2ydx
2
rrn)
2
+ y2 .
Mter changes of variables, and gluing the integrals together, we obtain a formula containing a single integral that can be explicitly evaluated to get 1 2rr
co=
f
2ydx
x
22=1.
+y
Chapter 8. Summation of divergent series and integrals
270
Similarly, em
= 1 l1f 7r
1f
ct>(x,y)cosmxdx
= 1 7r
f
2ycos(mx)dx 2
X
+y
2
= 2e my .
Hence, series (6) coincides with series (1), so that, for any y > 0, we have equality ct>(x, y) == U(x, y) or, equivalently,
 = 1+2 L sinh y
00
coshy  COSX
m=l
emy cos(mx)
=
L n=oo (x 00
2y 27rn)2
+ y2
.
(7)
Now, if we let y + 0+, the last series weakly converges to a periodic Dirac delta distribution
L
2y
00
lim
y+O+ n=oo (x  27rn)
2
+y
2
= 27r
L 00
n=oo
8(x  27rn),
and, in view of (4) and (7), we obtain a distributional Abel summation equality
L 00
L 00
eimx = 27r
(8)
8(x  27rn),
n=oo
m=oo
where the lefthand side is the Fourier series representation of the periodic Dirac delta with period 27r. An obvious extension of this formula for the case of an arbitrary period 27r / il has the form (9)
which should be compared with the already familiar formula (3.3.3) for the Fourier image of the Dirac delta. Equality (7), used here as a tool in deriving distributional formulas (8) and (9), can also be used as a summation tool for more ordinary series often appearing in physical applications. For example, setting x = 0 and x = 1l', we obtain that
n=l
4y (27rn)2
+ y2
y
= coth 2 
2
00
4y
" 7r2(2n Y' f:r.
y
1)2 + y2  tanh2·
B.B. Poisson summation formula
8.8
271
Poisson summation formula
A limited supply of elementary and special functions leads to a situation in which analytic solutions of many physical and engineering problems can only be written with the help of series of elementary or special functions. For example, the wellknown method of separation of variables in partial differential equations leads to solutions representable in the form of functional series, and the situation is similar for solutions obtained by the method of successive approximations. Often it turns out that a series obtained in this way converges poorly or does not converge at all. In such cases it is desirable to find a transformation accelerating convergence of that series, or an outright analytic expression for its sum. Sometimes, this goal can be achieved through the Poisson summation formula (1)
which immediately follows from (8.7.9) by multiplying both sides by f(x) and integrating them over the whole xaxis. Observe the main feature offormula (1). The slower the function f(x) on the righthand side varies, the faster the Fourier image j(w) on the lefthand side decays to zero. This means that the more terms one needs on the righthand side for good approximation of the infinite series, the fewer terms of the series are necessary on the lefthand side for accurate computation of its sum. So the Poisson summation formula is capable of transforming poorly convergent series into rapidly convergent ones, and many of its applications rest on the above phenomenon. Relying on properties (3.1.3) and (3.2.5) of the Fourier transform and on formula (1) we can rewrite the Poisson formula in the form
2:n:!1 m=oo L f(mt..)exp(im!1s) = n=_oo!1 L f s + 2:n:n) . 00
00
_ (
(2)
The lefthand side of (2) represents the discrete Fourier transform of function f (t), and the righthand side expresses it in terms of the ordinary Fourier image j(w). Hence formula (2) is useful for interpreting results of the computer implementation of the discrete Fourier transform.
Example 1. Poisson summation formula for series of rational functions. Consider the series
Chapter 8. Summation of divergent series and integrals
272
of rational functions which is often encountered in physical applications. A use of the Poisson summation formula gives that
S(ex)
=
L 00
j(m, ex)
m=oo
=
L 00
(4)
f(2rrn, ex),
n=oo
where
and
f(t, ex) =
f
j(w,ex)cos(wt)dw
in view of the evenness of j. The above integral can be evaluated by the method of residues or by checking the tables of integrals. The result is
f(t, ex)
=
(y cos(yt) ,B sin(YltD),
(5)
where
,B = ex8,
y=exJl8 2 ,
A substitution of (5) into (4) gives
S(ex) = 2rr[y "2 y
+Y L 00
e 21T,Bn cos(2rryn)  ,B
n=l
L e00
21T ,Bn
sin(2rryn) ] ,
n=l
which is recognizable as a familiar series of complex exponentials. The computation of its sums reduces to
=L 00
Q(.l.., v)
n=l
exp[(.l..
+ iv)n] =
1
exp(.l..
. + IV) 
1
for the sum of a geometric series. Since, clearly,
y [Y"2 + yRe Q(2rr,B, 2rry) + ,BIm Q(2rr,B, 2rry) ],
S(ex) = 2rr
an explicit calculation of the real and the imaginary parts of function Q(.l.., v) gives
S(ex)
sinh(2rr,B) = rry y cosh(2rr,B) 
,B sin(2rry) cos(2rry)
.
(6)
8.9. Summation of divergent geometric series
273
This, and related formulas, can also be found by other methods, but the Poisson formula provides, as a rule, the fastest path to the goal. The typical graph of S(a) (6) is pictured in Fig. 8.8.1.
S 14 12 10 8 6
4 2 10
12
a
FIGURE 8.8.1
Graph of series (3) in case of 8 = 0.05 evaluated with the help of formula (6). The larger a, the more terms in (3) are needed to maintain the validity of the result.
8.9 Summation of divergent geometric series The Poisson summation formula was used in the Example 1 of the previous Section 8.8 to explicitly sum a nontrivial but absolutely convergent series; later on we will see its applications to accelerate the convergence of already convergent series. In this section, however, we will utilize it to solve a more exotic summability problem for the everywhere divergent series
S(z) = 1 + 2
L cos(mz), 00
(1)
m=l
where z is a complex variable, and where the Abel method and thus even more so the Cesaro method, fail.
274
Chapter 8. Summation of divergent series and integrals
Let us form an auxiliary perturbed series
+2 L e 00
S(z, e) = 1
_
2
L 00
sm cos(mz) =
exp(em 2 + imz),
(2)
m=oo
m=l
which, obviously, absolutely converges for any e > O. If the limit S(z) =
lim S(z, e)
s+O+
exists for some z, then we will take it as a generalized sum of the divergent series (1). To find the set of z's for which the above limit exists we will transform series (2) by means of the Poisson summation formula
L 00
L 00
/(m) = 2:rr
(3)
!(2:rrn) ,
n=oo
m=oo
with /(t) = exp( et 2 + itz). The lefthand side of (3) coincides with series (2), and the righthand side contains function
/(w) = 1
2:rr
f
exp(et 2 + it(z  w»dt =
Hence, S(z, e) =
exp
It means that for any
4e
[(Z  42:rrn)2] . L exp e 2",:rr e n=oo 1
00
Notice that
I [
[(ZW)2] .
1 2", :rr e
(z  2:rrn)2]
4e
I
=exp
[y2  (x4
2:rrn)2]
e
(4)
.
z = x + iy from the set (blackened out in Fig.
G = {z E C : Iyl < Ix  2:rrnl,
8.9.1)
n = 0, ±1, ±2, ... ,
all the terms of series (4) converge to zero as e + 0+. Correspondingly, it is easy to prove that for z E G we have S(z) = lims+o+ S(z, e) = O. Thus, we demonstrated that
+ 2 L cos(mz) = 00
S(z) = 1
m=l
0,
ZE
G.
8.9. Summation of divergent geometric series
275
y
x
FIGURE 8.9.1 Region G in the complex plane.
In particular, for x =
1f,
we get that (Xl
1 + 2 L(I)m cosh(my) = 0,
Iyl
<
7r.
m=l
Substituting q
= exp y, the above equality can be rewritten in the form (Xl
L(I)m[qm
+ (l/q)m]
= 1,
(5)
m=O
Forq> 1,
converges absolutely. In view of the linearity and regularity of our summation method, we obtain the following striking extension of the familiar formula for the sum of a geometric series: m m
q
m=O
1 = , 1 +q
276
Chapter 8. Summation of divergent series and integrals
8.10 Shannon's sampling theorem In this section, we will obtain another useful result with the help of the Poisson summation formula. Shannon's sampling theorem, which is of importance in information theory, will follow as a particular case. Consider a function g(t) representing the uniformly convergent on the entire realaxis infinite series
L 00
g(t) =
l(m/:1)l{F(t  md).
(1)
m=oo
Here I(t) and l{F(t) are known functions with Fourier images i(w) and ;;'(w), respectively. Let us apply the Fourier transform (3.1.1) to both sides of (1). In view of the uniform convergence of the series, the integration and infinite summation operations can be interchanged so that, taking into account the formula (3.1.3a)
21l'
f
l{F(t  md)e iwt dt = ;;'(w)eiwmll.,
we have
L 00
g(w) = ;;'(w)
l(m/:1)e iwm ll..
m=oo
(2)
Finally, transforming the sum in (2) by means of the Poisson formula (8.8.2), we get
21l' 
g(w) = t;l{F(w)
(
1
w
+ 21l'n) .
(3)
Now, let i(w) be a function with compact support, identically equal to zero for
(4) Then the supports of the summands in (3) have empty intersections and, for a given frequency w, only one term of the infinite series is different from zero. In particular, for Iwl ::: 1l' / d, equality (3) is equivalent to the equality
21l' 

g(w) = t;l{F(w)/(w).
(5)
B.10. Shannon's sampling theorem
277
Assume additionally that if,(w) is of the form 
fl.
l/Io(w) = 21r n(fl.w/rr),
(6)
where the rectangular function n(v) = X(v
+ 1) 
X(v 1).
(7)
Then, equality (5) is valid for any w and assumes the form
(8)
g(w) = j(w).
The equality of the Fourier images implies the equality of the functions themselves. Hence, substituting in (1) the inverse Fourier image
sin1:
rrt
l/Io(t) =  ,
1:=
1:
(9)
of the rectangular function (6), we arrive at the equality /(t) =
f
m=(X)
/(mfl.) sin «rr t/fl.) rrm), (rrt / fl.)  rrm
(10)
which expresses the contents of Shannon's sampling theorem:
Assume that the Fourier image of function /(t) vanishes outside the interval Iwl fl.. Then, for any t, the values of /(t) are completely determined via formula (10) by the values of this function at discrete time instants tm = mfl..
(11)
The theorem has numerous applications in physics and information theory which will not be discussed here. However, we will take a closer look at some of its modifications and extensions which will permit us to grasp the meaning of this result at a deeper level. . First of all, note that if if,(w) is taken to be a onesided rectangular function (as opposed to the symmetric function (6» 
l/I(w)
= rrfl. [X(w) 
X(w rr/ fl.)].
(12)
278
Chapter 8. Summation of divergent series and integrals
Then the righthand side of (8) becomes 8(W) = 2j(w)X(w).
As we observed in Section 6.6, the righthand side is the Fourier image of an analytic signal F(t) corresponding to function f(t). The original function corresponding to the Fourier image (12) is
1/I(t) =
1].
1l'lt
Substituting this expression into (1), we find an explicit formula for the complex function F(t) =  i1l'm) 1,
f
m=oo
representing the analytic signal corresponding to function f(t). Let us rewrite formula (3), replacing
21l'
8(W) = 81/1(w) L
f(w
m=oo
+
(13)
and observe that if (14)
and j(w), as·before, identically vanishes for Iwl then the supports of summands in the series (13) are separated by gaps of length  II As a result, to pass from (13) to (3), it suffices to select 1fr(w) from a broad class of functions for Iwl for Iwl for 1l' I
1fr(w) = { 0, arbitrary,
1l' I <
Iwl
 II M; < 1l'(21  II M.
(15)
The schematic plot of one of possible Fourier images 1fr (w) for which the equality (8) is valid is shown in Fig. 8.10.1. Taking one of the Fourier images (15) and calculating the corresponding original function 1/I(t) we arrive at a formula more general than (to):
L 00
f(t) =
m=oo

(16)
8.10. Shannon's sampling theorem
279
FIGURE S.10.1 The plot of one of possible Fourier images lfr (w) for which the equality (S) is valid. The triangles symbolize here the summands in the series (13). In particular, killing the Fourier image (15) for widely used form of Shannon's formula
f(t) = e
f:
m=oo
Iwl
f(m8) sin «rr t/ ll.)  errm) , (rrt/ll.)errm
rr / ll. we arrive at the most
e = 8/ ll. < 1.
(17)
It would seem that formulas (16) and (17) have, in comparison with the simplest formula (to), the shortcoming of requiring the knowledge of values f (t) at densely distributed time instants tm = 8m. However, the indicated drawback is partly amended by the great flexibility in the selection of function 1/I(t). The freedom to impose values of the Fourier image lfr (w) in the intervals rr < ll.
Iwl
< rr 
(2 1)   8 ll.
(18)
can be utilized to improve the speed of convergence of the series appearing on the righthand side of the generalized Shannon's theorem. Indeed, the experience gathered by the reader while studying the Fourier transform's asymptotics in Chapter 4 suggests that to achieve that goal one should choose lfr(w) decaying to zero inside the intervals (18) as smoothly as possible. An infinitely differentiable (in the classical sense) lfr(w) would be ideal. In this case, the corresponding original function 1/1 (t) would decay for It I + 00 faster than any power function 1/ tn. As a result, it may turn out that in computing values of f(t) with a given accuracy fewer
Chapter 8. Summation of divergent series and integrals
280
terms are needed in the series (16) than in (10). An example of a damped Fourier image 1fr(w) and the corresponding rapidly decaying function 1/I(t) is provided in Exercises at the end of this chapter. In electrical engineering applications one often deals with narrowband / (t), that is, with functions whose Fourier images are concentrated in a narrow neighborhood of specific frequencies ±wo and vanish outside the intervals
Iw  wol
<
n,
Iw + wol
<
n.
Using the simplest Shannon formula (10) we should impose values of function /(t) in intervals of length (19)
1
because the Fourier image (w ) of a narrowband function / (t) is identically equal to zero only for frequencies satisfying condition
analogous to (4). At the same time it is intuitively clear that for a narrowband signal for which n « wo, one should have a more adequate version of formula (16), with sufficiently large intervals B between the readings: B 21l' / n » 1l' / wo. Let us derive it. To begin with, note that the compactly supported Fourier image of a narrowband signal can be represented in the form of two components
I(w)
= I+(w) + I(w),
(20)
concentrated in the neighborhoods of the central frequencies +WO and wo, rein the sum spectively. Choose B so that the central frequencies of components (13) would be located halfway inbetween adjacent central frequencies of the comTo accomplish this it is necessary that for some positive integer 1 we ponents have the equality 21l' 1l' wo +  I = WO  .
1+
1.
B
Consequently,
B
B = 1l'(1 + 1/2)/WO.
(21)
Select the value of I in such a way that, for any w, only one of the summands in (13) is different from zero. To achieve this it is sufficient to demand that B satisfies the inequality 1l' /2B > Hence by (21) it follows that
n.
(22)
B.11. Divergent integrals
281
It remains to select as {fr(w) a function that would tum equality (13) into g(w)
= i+(w) + i(w).
By analogy with (15), it is clear that here it is sufficient to let {fr(w) =
1
for Iw±wol w; for Iw±wol rr/8  Q; forQ< Iw±wol Q, that is selecting

/) [(wwo) n + n (w+WO)]
1/I'0(w) = 21r
where to be
Q
Q
,
(24)
n (v) is given by (7), we can calculate the corresponding original function 1/I'0(t) = 28 sin(Qt) cos(wot).
rrt
(25)
Substituting it into (16), we arrive at the simplest variant of the narrowband Shannon's theorem: 28
f(t) = 
LJ f(m8)
sin(Qt  m/)Q)
t  m8
rr m=oo
cos(wot  mw(8).
(26)
Note that, if f(t) is a real narrowband process, that is if Q « WO then, without violating inequality (22), one can select I » 1, and make the distance 8 (21) between readings tm much larger than tl. (19), as in the case of the standard Shannon's theorem.
S.ll Divergent integrals Summability problem for divergent integrals is close in spirit to that for divergent infinite series. Let us demonstrate this using the Fourier transform X(w) =  1
1
2rr 0
00
• eUJJtdt
282
Chapter 8. Summation of divergent series and integrals
of the Heaviside function X (t) as an example. It is a typical divergerit integral. Evaluating it by the infinitesimal relaxation method gives
1
00
o
1
00
. elwtdt = lim
a+O+ 0
. elwtatdt = .1 IW +0
(1)
and illustrates an application of the Abel summation method to divergent integrals. The generalized sum (1) of the above divergent integral is quite stable with respect to a wide class of regularizing functions. Indeed, let us replace the above exponential regularizing function exp( at) by an arbitrary, absolutely integrable and continuously differentiable function I(at), and consider the integral
(2) The change of variable x = at transforms (2) into the integral
11
a 0
00
. e IPX I(x)dx,
where we also introduced a new parameter p = wJa. If a + 0 then p + 00. Hence, to evaluate the above integral, we can employ results of Chapter 4 on the asymptotic behavior of Fourier transforms of discontinuous functions. In particular, by analogy with (4.3.3), we have that
11

a
0
00
1
. e IPX I(x)dx '" ./(0)
lap
1
= :1(0), IW
a + 0,
W
:;C o.
(3)
The regularity assumption for our summation method requires that 1 (0) = 1. So, for W :;C 0, the summation result is the same as in (1). The above summation method obviously satisfies the separation of scales condition (as discussed in the context of series summability methods). However, one can easily produce examples of methods which violate it and thus can potentially give nonunique answers.
Example 1. Divergent integral with nonunique generalized sum. Consider the divergent integral
10
00
sin x dx.
(4)
Summation of this integral following the prescription given in (2) gives value 1, which can be also found by formally writing the integral as a series E an with terms an =
l
lT (n+l)
lTn
sinxdx = 2(1)n,
B.11. Divergent integrals so that
10
00
283
sin x dx = 2  2 + 2  2 + ... = 1.
However, if one considers a method of summation based on the integral
[(a) =
10
00
eax sin x [1 + 2aq sinx]dx,
(5)
one obtains a result different from 1. Indeed, the integral (5) can be easily evaluated to give
1
/(a) =  2  
a +1
4q a +4
+ 2·
As a
0, the integrand in (5) converges uniformly on bounded sets of x's to sin xthe integrand in (4). However, for any q ¥= 0, lim /(a) = 1 +q
a+O
¥= 1,
and for different values of parameter q, one gets different summation results. This is related to the fact that the relaxation function e ax [1 + 2aq sin x] is not of the form !(ax) and violates the principle of separation of scales.
8.12 Exercises 1. Find the sum of the series
L emy sin(mx), 00
S(x, Y) =
Y>
o.
m=l
2. Find the sum of the divergent series 00
S = Lmsin(mx). m=l
3. Using the Poisson summation formula find the functional action of the distribution 00
E = Lsin(mx). m=l
284
Chapter 8. Summation of divergent series and integrals
4. The analysis of waves propagating in resonators and waveguides leads to series of the following type:
L 00
S(w, M =
f(mMexp(img(w»,
m=OO
where f (x) is an absolutely integrable function such that f (0) = 1. Transform the above series by means of the Poisson summation formula and find its weak limit for /). + o. 5. Using formula (8.7.8), transform the series P(x) =
L eXP(i(m. a (x) ). m
where the summation is extended to all vectors m with integer components (m1, m2, m3), and y = a(x) is an infinitely differentiable vector field which provides a onetoone mapping of the xspace into the yspace such that the Jacobian J(x) = laa/axl of the transformation is everywhere positive and continuous. The inner product appearing under the summation is (a(x) . m) = a1 (x)m1 + a2(x)m2 + a3(x)m3. 6. The 3D Poisson summation formula turns out to be useful in solid state physics and, in particular, in the study of crystal properties. Using the result from Exercise 5, derive the righthand side of the formula if the lefthand side is 00
F(x) =
L
ml,m2,m3=oo
where a1, a2, a3 are three not coplanar vectors, x = (Xl, x2, X3) in a certain Cartesian coordinate system, and f(x) is an absolutely integrable function with 3D Fourier image
j(k).
7. Calculate the discrete Fourier transform of the rectangular function n(t) = X(t
+ 1) 
(1)
X(t  1).
8. Find the discrete Fourier transform of the function f(t) = n(t)cos4 (1l't/2)
and, then, compare it with its usual Fourier image (4.3.22). 9. Find the discrete Fourier image of the periodic function where M is a positive integer. 10. Find the sum of the functional series S (/).)
sin(m/).)
= m=1
with the help of the Poisson summation formula.
m
f (t), with period T
= M /).,
285
B.11. Divergent integrals 11. The functions IN(t) =
t
m=l
=
are partial sums of the Fourier series of the periodic function 1 (t) (1f  t) 12+1f It I (21f) J (see solution to Exercise 10). Compare functions IN(t) and I(t), and investigate the asymptotic behavior (for N + 00) of IN(t) in the vicinity ofthe discontinuities t = 21fn of I(t). 12. The Gibbs phenomenon discovered in Exercise 11 is undesirable in many physical and engineering applications. Find the method of summation of the first N terms of the series sin(mt) IN(t)
=
N)m
m=l
which, at the continuity points of I(t) from Exercise 11, would guarantee convergence of IN(t) to I(t), and would avoid the Gibbs phenomenon at the discontinuity points of I(t).
13. Derive a formula analogous to the formula (8.9.26) for the analytic signal F(t) corresponding to the narrowband signal I(t), whose Fourier image is identically zero for Iw±wol::: n. 14. Construct a function 1/I(t) entering into the Shannon series (8.9.16), and possessing a Fourier image of the type (8.9.15). Use common sense. 15. Find the maximal distance" between readings in the narrowband Shannon's formula (8.9.25) for n = wo/lO. 16. Suppose that I(t) is a narrowband function with the Fourier image vanishing for Iw ± wol ::: n = wo/lO. Find function 1/I(t) entering into the generalized Shannon formula (8.9.16), which decays sufficiently rapidly as It I + 00. Utilize results of the Exercises 14 and 15.
Appendix A Answers and Solutions
A.l
Chapter 1. Definitions and operations
1. (a) (./1i/2)8'(x); (b) «rr/2)8'(x); (c) The "zero distribution" 0; (d) X(x). 2.
f
f(x)dx = 0,
f
F(x)dx = 1, where F(x) =
f(y)dy.
3. (a) x'(ax) = 8(x) and is independent of a. (b)
L 00
x'(e Ax sin ax) =
(1)"8 (x _
n=oo
:n)
and is independent of A.
4. f'(x) = 8(x  1)  8(x + 1), and
lim f'(x/s)/s2
8+0+
= U'(x).
5. y(x) = A8(x  1) + B8'(x  1) + C8(x). where A, B, C are arbitrary constants. 6. Taking into account the multiplier probing property, we have
=
8. 8a 8(g(x)  a)IVg(x) I. The action of this distribution as a functional on an arbitrary test function t/J e V(R3 ) is 8a [t/Jl =
f
8(g(x) a)IVg(x)It/J(x)d3 x =
© Springer Nature Switzerland AG 2018 A. I. Saichev and W. Woyczynski, Distributions in the Physical and Engineering Sciences, Volume 1, Applied and Numerical Harmonic Analysis, https://doi.org/10.1007/9783319979588
1
t/Jdu.
288
Answers and solutions
9.
I[Vg1 (x)
8l =
10.
J IVy,(x) 18(y, (x) 
X
2
IT 8(gk (x) 
V g2(X)] I
ak)·
k=1
e)d2 x.
11. P[ cP] is equal to the flux of the vector field surface y,(x) = e: P[t/>] =
1cP·
cP inside the region bounded by a level nda,
where n is the interior unit normal vector to the level surface a. The interior points x are those for which y, (x) > e.
A.2 Chapter 2. Basic applications Ordinary differential equations
1. ji + yy
+ (J)2y =
(b
+ ya)8(t) + a8(t).
2. y(O) = 1, Y(O) = y.
3. y(t) = (1/2) sin Itl. 4. y(t) = X (t)(a
+ (b 
a) cost).
5. ji + y = 8(t). 6. y = Aeht
+ BX(t)eht •
Wave equation 7. The Green function G(x, t) is a solution of the following initial value problem:
82 G
a2 G
at 2 = e2 ax 2 '
G(x, t = 0) = 0,
8
at G(x, t
= 0) = 8(x).
Using the D' Alembert formula we obtain that G(x, t)
= ;eX(t)( X(x + et) 
X(x  et)).
Hence the solution of a nonhomogeneous wave equation is of the form
1 1t
u(x, t) = 2c
dT
dy f(y, T).
00
8. u(x, t) = g(x) the spatial variable.
* 8G(x, t)/8t + h(x) * G(x, t), where * denotes the convolution in
289
Chapter 2. Basic applications Continuity equation 9.
f
p(x, t) = C(x, t) =
Po(y)8(x  yeBt)dy = e gt Po(xe gt ),
f
Co(y)8(y  xegt)dy = Co(xe gt ).
For g > 0, as time increases, the density at each point x decreases to 0, and the concentration asymptotically becomes homogeneous and converges everywhere to a constant. IT the total mass of the passive traces is equal to m, then for g < 0 and t .... (Xl the density weakly converges to the distribution m8 (x) and the concentration to the zero distribution.
J
10. p(x, t) = Po(y)8(x  y  gyt)dy = (1/1 + gt)Po«x/l + gt». For g > 0 the density converges everywhere to zero but much slower than in the previous problem. However, for g < 0 the density becomes singular in finite time t* l/g.
=
=
J
=
11. p(x, t) PoV 8(x  v.  g.2/2)d. Po(1 + ag/v2)1/2. One can see a similar phenomenon watching the stream of water flowing out of a faucet. The further it is from the faucet, the thinner it gets. 12. The equation in question has the integral p(x)v(x) = const = PoV. Eliminating the time t from the equations of motion x
we get v(x)
= vt + gt 2/2,
v(t) = v
+ gt,
= (v 2 + 2gx)1/2. Thus, p(x) = Pov/(v2 + 2gx)I/2.
13. In this case the x axis is simply reversed upwards and we are not solving the "rain" problem but a "fountain" problem where the droplets move upwards. In this case the density of droplets is p(x) = 2Po ( 1  ag/v 2)
1/2
.
=
if x < v 2/2g, and it is 0 if x > v 2/2g. The infinite singularity of the density as x .... v 2 /2g  0 describes the effect of concentration of droplets at the point where rising droplets reach their highest elevation. By the way, explain why an "extra" coefficient 2 appeared in the numerator of the above formula. 14. The continuity equation
+ 0\11 !!!... _ 0\11 !!!... = 0 at OX2 OXI OX1 OX2 has the integral of motion p = g(\II (Xl, X2», where g(\II) is an arbitrary function. IT we Op
additionally impose a physical requirement that g 0 then we obtain a class of densities of the passive tracer that will be independent of time. 15. Since in this example the dimension of the Dirac delta is [8]
1
= [t]/[x 2], we get that
[a] = [m]/[t]. Also, m = a dl/lvl, where the integration is performed on the contour (or contours) given by equation \II(x) = b.
290
Answers and solutions
16. The Liouville equation for the particle density is
a!  + (v· V x )! + (V v ·g(x, v)f) at
=
o.
17. In this case the Liouville equation is
a! 8i + (v· V x )! = h(V v · vf). Its Green function is G(x, v,Y,
U,
t) = 8(x  Y  (ul h)(l  e T»8(v  ue T).
"C = ht introduced above is dimensionless time. Using the Green function we can write the solution in the form
!(x, v, t)
=
f
!o(y, u)8(x  Y  (ul h)(l  e T»8(v  ueT)d3y d 3u
= e3T !o(x  (vi h)(eT  1), veT).
18. p(x, t)
= J !(x, v, t)d3v = my3(t)w(xy(t», where y(t) = heTl(eT 
19. lim,....oo !(x, v, t)
1).
= mh 3w(xh)8(v).
20. Density p(x, t) satisfies the continuity equation and velocity v(x, t) satisfies the (nonlinear) Riemann equation
av 8i + (v· V)v = O. 21. Replacing the substantial derivative by its expression in Eulerian coordinates we obtain the Riemann equation with external forcing acting on particles in the hydrodynamic flow
av at
22. v(x, t) =
+ (v . V)v = g(x, v, t).
Jv!(x, v, t)d3vI p(x, t).
23. The analogue of formula (2.4.11) in the ID case is
_ ( ) =!!.
A
Po x
usg
(x  f3(il:lS») I
.
1=00
For l:l + 0 and smooth functions g(z) and f3(y), the sum converges to the integral
 (x ) = !!.I Po
f
g
(x  f3(y») I
d y.
Chapter 2. Basic applications
291
Now, observe that for I + 0, the function g (x / I) / I weakly converges to the Dirac delta 8(x). In this fashion the sought limit function is Po(x) = p /
8(x  f3(y»dy = pa'(x),
where y = a(x) is the function inverse to the function x = f3(y). In our case, for x > anda(x) = (x 1 + ../x 2 +s2)/2, _
Po(x) =
IXI)
p (
'2
°
1 + ../x2 + s2
°
°
.
Remark. The condition ll. = 0(/) for I + permitted us to correctly choose the order of the limit passages: first let ll. + for a positive I, and then I + 0.
Pragmatic approach
24. JL = (lal 25.
f =
(a
+ If3D/I2af3l, if one assumes that Dirac delta is even.
+ 13)/2.
26. y(o_) = y(O+), .y(0_) = eYy(O+). 27. y(O_) = eYy(O+), y'(O_) = y'(O+). 28. (a) R(y) = 8(y) q,(x) dx + x (y)q,(y); (b) R(y) = x(lyl 1)q,«lyll); (c) R(y) = 8(y) q,(x) dx + q,(lyl + 1). 29. The generalized particle density is p(x, t) = m(t)8(x)
where m (t)
+ p(x, t),
= Po.;w; is the mass of the cluster of particles sticking at the origin, and 1 p(x, t) = Po2 (1
+
Ixl
../x2 +wt
)
=
=
is the continuous component of the density. It has a minimum p(O, t) Po/2 at x 0, and increases to its original value Po as Ix I + 00 where the particles do not move. The mass conservation law is here reduced to the requirement that the mass of the cluster is equal to the deficit of the mass of unglued particles: m(t)
= /(PoP(X,t)]dX.
A direct substitution leads to the equality
to
10
dx
../x2
+ l(x +.JX2+1)
= 1,
one of the standard definite integral that can be found in the tables of integrals or via Mathematica. Here it was derived via a "physical" argument.
292
A.3
Answers and solutions
Chapter 3. Fourier transform
1. We have f+(w) = 1 2rr
1
00
0
' ) tdt = 1 .1 . e (Y+IW 2rrlw+y
2. In view of the properties of the Fourier transform the answer can be obtained from the answer to the previous exercise: 1 1 f(w) = 2Ref+(w) =  2 2' rr w + y 3. j(w) = (1/2y)e Y1wl •
4. j(w)
= 8(w) 
(y/2)e Y1wl •
5. The answer is
2:
00 1 f(w) = 1 e (.IWY ) n = _1 . 2rr n=O 2rr 1  exp(iw  y)
6. We have
fo(w)

1
sinw ' cosw
= 1m f(w) = 4rr cosh y 
sinh y fe(w) = Re f(w) = 1 [ 4rr cosh y  cos w 7. r(y) = coth(y /2)
00
for y
+ 1] .
O.
8. The Fourier image of the derivative is iwj(w) = (1/2)(e iW  e iw ). Hence, by the shift theorem, f'(t) = rr[8(tl)8(t+ 1)]. Integrating t, and assuming that f( 00) = 0 we get that f(t) = rrTI(t), where TI(t) = 1 for It I :::: 1, and = 0 for It I > 1.
9. The second derivative f"(t) = 2(8(t+ 1)+8(t 1)  TI(t)). Using the answer to the previous problem we get that its Fourier transform is (2/rr)(cosw  sinw/w). Therefore j(w) = (l/w2)(2/rr)(sinw/w  cosw). 10. The third derivative flll(t) = 8(t
+ 2) 
28(t
+ 1) + 28(t 1) 
The corresponding Fourier image F'(w) =
 2sin(w»).
Multiplying this expression by iw3 we finally obtain 2sinw(1  cosw) f(w) = 3' rrw
8(t  2).
293
Chapter 3. Fourier transform 11. The second derivative
L 00
f/,(t) =
l:!.2f(nl:!.)8(t  nl:!.),
n=oo
where we used the standard notation l:!.2f(nl:!.) = f«n
+ 1)l:J.) 
for the secondorder difference of function

2f(nl:!.)
+ f«n 
1)l:!.)
f. The corresponding Fourier image
1
.
t:.
fi(w) =  23r l:!.w2 L.J l:!.2f(nl:!.)e"J} n. n=oo
Regrouping terms of the above series, we arrive at
which is more convenient for calculations.
12. First, observe that the sought Fourier image is related to the Fourier image of function h(t) via the formula j(w) = 23rlh(w)1 2 • Calculate first the Fourier image of function h(t). Its second derivative is h" = 8(t)  8(t  (})  (}8'(t  ()).
Thus the Fourier image of function h(t) is 
(}2
[
'0
h(w) = 23r02 e' (1
where 0 =
or, finally
w(}
+ in) 
]
1 ,
is a new dimensionless argument. In this way,
_
(}4
f(w) = 23r04 [0 2 + 2(1  cos 0  0 sin 0)] .
13. f(4)(t) = 8(t)  8(1tl (})  (}8'(ltl ())  (}28"(t). 14. Apply the Fourier transform to both sides of the recurrence relation to obtain 
f(w; p) =
p(p 1) 2
p  w
2
f(w; P  2).
In the case of even p = 2n, the previous identity implies that


f(w; 2n) = 2n!f(w; 0)
nn m=l
1 2 w2. 4m 
294
Answers and solutions
Substituting f(w; 0)
sin(Jrw/2) = Jr1 17r/2 cos(wt) dt = , 0 JrW
we obtain sin(Jrw/2) f(w; 2n) = 2n! JrW
+ 1,
Similarly, for odd p = 2n

f(w; 2n
nn m=l
2
 1) nn + 1) = (2n + 1)!f(w;
1 f(w;I)=Jr
2·
1
m=l (2m
A substitution
1
4m w
+ 1)
2
 w
2.
17r/2 cost cos(wt)dt = cos(Jrw/2)
2'
Jr(1  w )
0
gives, finally,
+ 1) = (2n + I)! cos(Jrw/2) nn
f(w; 2n
A.4
Jr
m=O (2m
1
+ 1)
2
 w
2·
Chapter 4. Asymptotics of Fourier transforms
1. f(x) '" x 2 (x
+
2. f(x) '" x/3 (x 3. f(x) '" x/3 (x
0).
+
0), f(x) '" l/x (x + 00).
+ 0).
4. It is clear that the root x(a) + 0 as a + 00. Hence we can use the asymptotic formula cotx = l/x  x/3 + O(x 3 ) and replace the original equation by a simpler quadratic equation 2 3 3 (1) x  ax +  =0
2
2
which gives the following asymptotic behavior of the root: (a + 00).
(2)
Remark. A related question of finding positive roots of the transcendental equation
x = tan x also arises in several mathematical physics problems.
(3)
Chapter 4. Asymptotics of Fourier transforms
295
x,tan(x)
X FIGURE A.4.1 Graphs of the functions marked on the xaxis.
x
and tan x. The first two roots, Xl, x2, of equation (3) are
The roots xn (see Fig. A4.1) can be expressed via the solution x (a) of the equation (1) as follows: where
=
an
rr(2n
2
+ 1)
Hence, in view of (2), the corresponding approximate values of the roots of (3) are Xn
an
4
+4
J
a2 n
(4)
3
The larger n the more accurate the approximation is. However, even for n = 1, the approximate value Xl 4.493397 given by (4) differs from the true solution by less than 1.2 . 105 • So in many practical situations expressions (4) are as good as exact analytic solutions of the equation (3). 5. Taking logarithms of both sides to replace the product by a sum, we obtain that In f(N) = :L:=lln(l+anl N). The boundednessofsequence {an} implies that ani N + o as N + 00, uniformly in n. This permits use of asymptotics In(1 + x) ..... x, for each term separately, and adding them up. As a result, we get In f (N) ..... (II N) an. Hence, f(N) ..... exp[(I/N)
:L:=l
:L:=l an].
6. The answer is
D(1 +
= exp
(1'
cp(r)dr) (1 R),
296
Answers and solutions
o<
R < !l
l'
q12('r) dr:.
We tum your special attention to this formula as it forms a basis of numerous scientific (Physical, chemical and biological) laws.
7. Integral Sew) decays, as w .... 00, more rapidly than any power of w. Integral C(w) has a power asymptotics C(w) '" (a/w2) (w .... 00). 8. jew) = (6a/7rW4) + O(w 10 ). 9. We have
n!A. 7r few) '" 7rw"+1 sm(wr:  '2 n )
(Iwl .... 00).
The trigonometric factor in the above formula describes, as physicists say, the interference contributions from hidden discontinuities of the function at two points t = ±r:.
10. J(w) '" (3/w2)cosw. 11. jew) '" sinw/(i7rw2), (w .... 00).
12. The answer is
r(p + 1) exp (.7r ) f(w) '" 27rlwIJ:l+l ''2(1 P)'Slgnw,
(w .... 00).
13. One obtains sin Iwl  cos Iwl few) '" 2.Jil'Wflwl '
A.5
(Iwl .... 00).
Chapter 5. Stationary phase and related methods
1. Apply the general asymptotic formula (5.2.3). In our case, take x = kp .... 00, pet) = cosh(t), f(t) = 1/27r. There exists only one stationary point r: = 0 where pl/(O) = 1. Consequently, there remains only one term in the formula (5.2.3), and that term has to be multiplied by 1/2 because the stationary point coincides with the lower limit of integration. Hence, G(p) '" 
1 . exp[Ikp), "j8mkp .
(kp .... 00),
an expression equivalent with (9.4.3).
2. Initially, it is easier to find asymptotics of the complex function Dw(z) described in Chapter 4. To make use of it notice that Dw(z) (5.6.1) is the Fourier transform of function f(t/» = 2eizsinif>[X(t/»  X(t/>  7r»),
Chapter 5. Stationary phase and related methods which has two jumps: Lf(O)l relation (4.3.3), we get
297
= 2 and Lf(Jr)l = 2. 1
.
Dw(z)   . (1  e'W1r), JrCUI
So, in view of the asymptotic
(cu + 00).
Now, separate the real and imaginary parts of this expression to obtain the sought asymptotics of Anger and Weber functions:
1.W _ sin(cuJr) , CUJr
EW _ 1 cosCUJr ,
(CU
CUJr
+ 00).
3. Again it is simpler to find initially the asymptotics of Dw(z). Separation of it real and imaginary parts is straightforward. The asymptotics is obtained by the stationary phase method. There is only one stationary point rp* = Jr /2. Therefore, with help of formula (5.3.1), we get
4. In our case, calculation of integral (5.6.1) in the Fresnel approximation is reduced to an application of approximate equality
where, as in Exercise 3, the stationary point rp* = Jr /2, and to the extension of the domain of integration to the entire line. Thus, in the Fresnel approximation, the exact integral (5.6.1) is replaced by the approximate expression
The last integral can be evaluated in closed form using the standard formula (3.2.3) which gives Dw(z)
Notice that as
z +
V{2 ;zi exp [( i z  2
CUJr)
2
iCU ] + 2z '
(z
»
1).
00 the above expression tends to the expression obtained in Exercise
3. 5. Let us write integral (5.6.1) in the form
(2) where p(rp) = rp  psinrp,
(3)
which is more convenient for asymptotic analysis. In our case (0 < p < 1) function p(rp) is strictly monotone so that there exists the unique and strictly monotone and smooth
298
Answers and solutions
inverse function cP = q(s). Consequently, it is possible to write the integral (2) in the familiar Fourier transform form:
(4) where !(s)
= q'(s)[X(s)  Xes n)].
(5)
This function has two jumps
L!(On
= 2/ P'(O) = 2/(1 
Thus, in view of (4.3.3), we get Dw(pw) _ nWI
p),
L!(n)l
(_1__ 1 p
= 2/ P'(n) = 2/(1 + p).
e iW1r )
,
1+P
(W _
00).
(6)
Notice that this relation is a natural generalization of the answer to Exercise 2 (in the case z < W (p < 1)) which follows from (6) by taking p _ O. 6. Observe that in the case p = 1, the inverse function cp = q(s) is no longer smooth over the entire interval cp E [0, n]. Indeed, in the vicinity of cp = 0 the original function P(cp) and its derivative P'(cp) have the following asymptotic behavior:
=
Consequently, cp(s)  (6s)1/3 and the function !(s) (5) has at s 0 a singularity of the order 16 )1/3 (s _ 0+), !(s) ( 9s 2 Asa  1,
=
where
A = (16/9)1/3,
a
= 1/3.
Thus, the asymptotics of Dw(w) is described by (4.6.2), i.e., r(a) Dw(w)  A 2n(iw)a'
(w 
00).
Inserting the numerical values of constants A and a and utilizing the symmetrization formula r(a)r(1  a) = n/ sin(na) one gets Dw(w) 
2 ) 1/3 .j3 _ i ( 9w "J273
,
(w _ 00).
(7)
=
7. In this case function P(cp) (3) has a simple stationary point cp* arccos(1/ p). For convenience, express p via an auxiliary variable 9: p 1/ cos 9, 0 < 9 < n/2. Then p* = 9 and elementary calculations yield
=
P(cp·)
=9 
tan 9,
P"(cp*)
= tan 9.
Chapter 5. Stationary phase and related methods
299
A substitution of these quantities into (5.3.1) gives D(Apw) '"
1 . exp [I (w(tan9  9) ",/4)]. pcos9 ../2rrwtan9
Noticing that wtan9 =
J Z2 
9 = arctan
w2 ,
= 1,
(w
4
(0).
(8)
Jz2/w2 1,
it is easy to see that (8) is a natural generalization of the answer to Exercise 3 in the case z > w (p > 1).
8. Transform (5.6.3) applying the following "physical" argument: if t represents the time and wthe frequency, then it is natural to analyze the integral (5.6.3) in the dimensionless variable of integration u wt and rewrite the former as
=
1
I(w) =
(9)
'\I 2rrw
The deliberately separated factor 1/../21rw has the dimensionality of the original integral and (10)
is a dimensionless function of a dimensionless variable y = v=Ju+y 2 Ieadsto A(y)
where (y) =
i1 rr
°o·2
Y
e IV dv
.fWT. The substitution
= e'Y• 2 (y),
=
1
M"! '\121
(11)
C(v 2/rry) + i S(v 2/rry)
is the complex Fresnel integral expressed through the Fresnel sine and cosine integrals discussed in Sections 5.34. The discussions of Section 5.4 indicate that ( ) '" { 1/..tiI, y 1/(iy..[iii)eiy2 ,
(y
4
0);
(y
4
(0).
Consequently, it follows from (911) that I (w) '" { 1/(2../rriw), 1/(2rriwv'L),
(wr
4
(WT 4
0); (0).
(12)
Having arrived this far, we already have noticed the remarkable fact that the above formula contains the main asymptotics of our integral: for T = 0 (5.6.5) and for T > 0 (5.6.4). Although, for any T > 0, the formula (16) implies the asymptotic power law with a = 1 (w 4 (0), for w « T, we already witness the appearance of the intermediate asymptotic power law with a = 1/2 (5.6.5). For T 4 0, the region of frequencies where the power law with a = 1/2 obtains, expands towards large w's, and for T = 0 the power law with a = 1/2 is valid everywhere.
300
Answers and solutions
9. Introduce a new function (13)
which is well defined for any real x > 1. For x = n, n = 1, 2, ... , it coincides with the factorial n!. Let us rewrite the integral (13) in the form suitable for asymptotic analysis by introducing a new variable of integration T such that t = XT. Then x! = x(x+l)
1
00
exp[xP(T)]dT,
where P(T) = T In T.
This function has a unique minimum at of (5.5.3),
= 1, where P(l) = P"(l) = 1. Hence, in view
T
(x + 00),
which, in particular, gives the Stirling formula. Note that the asymptotic Stirling formula gives a good approximation of the factorial for any finite n. For 1! = 1 we get a decent approximation 0.9221, and even for (1/2)! = ..fii/2 0.8862, the Stirling formula gives 0.7602. In a certain sense one can claim that the Stirling's formula is most precise for x = 1. The relative error decreases with the growth of n but the absolute error increases!
10. Introduce in the integral (5.6.9) a new variable of integration y
=x
 v(x, t)t.
(14a)
Taking into account equation (5.6.6) satisfied by the solution v, the old variable of integration is expressed in terms of the new variable via the equality
x = P(y, t) := y
+ vo(y)t,
(14b)
and the Fourier integral takes the form V(K, t) =
2:n:
f
vo(y)eiKP(y,t)dP(y, t).
Integrating by parts we arrive at a more convenient for asymptotic analysis expression ii(K, t) = _1._
2:n:1K
f
v' (y)eiKP(y,t)dy. 0
(15)
Here, and below, the prime denotes differentiation with respect to y. As long as 0 :::: t < l/u, the function P(y, t) is a strictly monotone smooth function, with P'(y, t) > 0 for any y, and the Fourier image (15) decays to zero, as K + 00, more rapidly than any power function. However, at the time t = l/u, there appears on the yaxis a point y = z where P'(z, l/u) = O. It means that the behavior of P(y, l/u) irI the vicinity of this point
Chapter 5. Stationary phase and related methods
301
might have the power asymptotics of the integral (15) as K 00. To find this asymptotics expand P(y, l/u) in the Taylor series in powers of (y  z): P(y, l/u) = P(z, l/u)
1 " (z, l/u)(y + '2P
1 11/ (z, l/u)(y  z) 3 + ... z) 2 + 6"P
The term of order 1 has disappeared since, in view of the assumption of this exercise, = u and P'(z, l/u) = O. Moreover, since z is the minimum of function vo(z), we also have P"(z. l/u) = O. Suppose that
Then the above Taylor expansion gives the asymptotic relation P(y, l/u)  P(z. l/u) '" P . (y 
d.
(y + z).
Since (y  Z)3 is an odd function of (y  z) we have to use the asymptotic formula (5.1.9) and obtain that
_ l/u) '"  ur(4/3) . P (z, 1 /u)] Re V(K. . exp [IK 2:7r1
(p.IK 4)1/3 ,
(K
00).
(16)
In the concrete case (5.6.10), where
z=O,
u=l,
P(z,l) =0,
P=l,
we get (K + 00).
Remark. Physicists and engineers usually do not work with the complex Fourier image of the signal but with the real and nonnegative spectral density of the signal's energy
J
In this fashion, its integral over the entire K axis gives the "total energy" v 2 (x , t) dx of the signal. The corresponding asymptotics of the energy density at the time t l/u is then (K + 00).
=
302
A.6
Answers and solutions
Chapter 6. Singular integrals and fractal calculus
Principal value 1. Transforming the original partial differential equation by Fourier transfonn in the space and time coordinates x and t produces the algebraic equation 2 2 iwu(k  w ) = 21r few).
Its solution u(k, w)
i = few) 41r
[1
1]
   . k  w k +w
First let us calculate the inverse Fourier image in x u(x, w) =
f
u(k, w)eikx dk
of each of the two tenns in the brackets and write u(x, w) = u_(x, w)  u+(x, w),
where i u±(x, w) = 41r few)
f
eikx dk k ±w .
To obtain u_(x, w) observe that in view of the causality principle the frequency wappearing inside the integral should be replaced by w  iO. The resulting integral is then evaluated with the help of (6.2.4) which gives
f f .+ = f = f = eikx dk k  w 10
PV
eikx dk k w
 
. i1re'ox.
To calculate the above principal value integral note that
PV
f
eikxdk kw
. e'ox PV
e
isx s
ds
. e,oxi
sin(sx)    ds s
. = i1re'QX sign (x).
Substituting the righthand side of this equation into the preceding expression we obtain
f .+ = eikxdk
k w
So u_(x,w)
Similar calculations give
10
. '21rie' OX x(x).
1 . = "if(w)e'OXx(x).
Chapter 6. Singular integrals and fractal calculus
303
Therefore u(x, w) =
j(w) exp( iwlxl).
Finally, taking the inverse Fourier transform in w we obtain
1
2. f(f Ixl).
u(x, f) =
As expected, the obtained solution satisfies the radiation condition which, in this case, means that the waves generated by a point source at the origin should run away from the origin. Hilbert transform 2. 1/I(f)
= (1/:7r) In I(,r 
t}/fl.
3. The problem can be solved by passing to the Fourier images of the corresponding functions. We know that iP(w) = exp( lwl1')/21'. According to (6.5.2), 1/I(f) = 21m
1
00
Hence, 1/I(f) = 1'(f2 4. 1/I(f) = ie mt sign (n),
 sin Inlf. If fP(f)
n # o.
(1)
iP(w)i wt dw.
f
+ 1'2) .
In particular, if fP(f) = cos nf, then 1/I(f) =
= sin Inlf, then 1/1 (f) = cos nf.
5. If one remembers that the Fourier image iP(w) = 21)X(W + v)  X(w  v)),
then the simplest way to find 1/I(f) is to utilize the formula (1) which holds true for any real function fP(f). Thus, 1
l
v
•
1/I(f)=SInwfdw= v 0
cos vf  1 vf
.
Analytic signals 7. Recall that the Fourier image of the original function is
= 411' (elwnIT
+ e1w+nIT).
The corresponding Fourier image of the analytic signal (6.6.1) is equal to
= Utilizing the inverse Fourier transform t(f)
=
21
00
dw
Answers and solutions
304
yields
8. Two cases should be considered separately:
vt e,'n = sin __ ..t, vt
0> v,
and
(e ivt cosOt 1), 0< v. Ivt Notice that in the first case the analytic signal is obtained from the original function just by replacing cos Ot by emt . This is a particular case of a more general result that is the subject of Exercise 8. In the case 0 = v both expressions coincide. =
9. Proof' It is evident that 1 '0 = 2"[/(t)  ig(t)]e' t
1 + 2"[/(t) + ig(t)]e''0t.
Now it follows from the assumptions that the Fourier image of the second summand is equal to zero for w > 0, hence the Fourier image of the first summand coincides with onehalf of the Fourier image of the analytic signal. As a consequence, the analytic signal is twice the first summand. Remark 1. Recall that the imaginary part of the analytic signal (6.6.2) is equal to minus This implies the following corollary to the above the Hilbert transform (6.5.1) of result: If qJ = / cos Ot + g sin Ot then its Hilbert transform is "'" = g cos Ot  / sin Ot. Remark 2. Signals with finitesupport Fourier image seldom appear in electrical engineering applications. However, for narrowband signals the replacement of the actual analytic signal by the expression (6.9.3) gives a rather good approximation. For example, the signal from Exercise 6 has an unbounded Fourier image but is narrowband if 0. » 1. In this case, it is easy to see that the approximate expression
is very close to
10. Function g(t) has the Fourier image g(w) = iB'(w) supported by the single point
= O. Hence, according to Exercise 8. = itemt . Generally, if in formula (6.11.1) /(t) = Pn(t), g(t) = Qm(t), are polynomials of degrees nand m, respectively, then
w
because the Fourier transforms of these polynomials have a onepoint support w = O. 11. The corresponding analytic signal is
= A(t)ei«l>(t) = Alem}t  iA2e m2t .
Chapter 6. Singular integrals and fractal calculus Consequently, A(t) =
J
Ai +
305
+ 2AIA2 sin(02 
01)t.
Assuming, for the sake of concretness, that Al > A2 (p = A2/ Al < 1), we get A.(
'¥
n ( pCOS(02  01)t ) t) = .'It  a r c t a n . , 1 + p sm(02  01)t
and
12. For convenience, let us replace the original function by its complex twin = X (t)e iflr and recall that the sought function 11(t) is equal to minus Hilbert transform (6.6.2) of the original function. So the corresponding complex twin 11c(t)
1 = pv 7r
f
ds s t
1 = eiflr PV 7r
1
i9
00
Or
e dO. 0
Here we have introduced the new integration variable 0 = O(s  t). After splitting the integral into real and imaginary parts one gets 11c(t) = e iflr
where Ci(x) = 
1
00
x
coso   dO,
0
 i
Si(Ot»)] ,
. l
x> 0,
SI (x)
=
0
x
sinO
0 dO,
are, respectively, the integral cosine and sine functions introduced before. Clearly, the imaginary part of 11c(t) coincides with the sought function 11(t) so that finally 11(t) =
.!.7r Ci (lOtI) sin Ot 
(!2 + .!.
Si (Ot») cos Ot.
7r
Notice that, in a sense, the corresponding analytic signal violates the causality principle. Indeed, at t < 0, when the original signal is identically equal to zero, the analytic signal along with 11(t), is aheady nonzero. 13. The Fourier image (4.4.10) of sponding analytic signal is equal to
(4.4.9) has aheady been found. So, the corre
(2) The last integral obviously reduces to integral (4.4.8) which we will rewrite in another, more suitable for our purposes, form (0 < {3 < 1):
1
00
o
wfJleiwr dw = f({3)1tIfJ x { I
t > 0 ,1ft < O.
Answers and solutions
306
FIGURE A.6.1 Graphs of (a)
and (b) 71 as functions of the dimensionless variable x = Ot.
Taking P = 1  ot and substituting the result into (2) we obtain r()_ Itl a  1 .. t  . x SID7rot
{Sin7rot+iCOS7rot,
.
I,
ift>O; if t < 0.
14. It follows from the Hilbert transform property (c) of Exercise 5 that
= W) + where the asterisk means the complex conjugate and previous exercise. Thus
is the analytic signal from the
15. Similar to the sought function is a singular distribution. In order to find it, let us multiply the formal equality 71(t)
= 7r1 PV
f
8'(s) ds S 
t
Chapter 6. Singular integrals and fractal calculus
307
by a test function cp(l) and integrate it over all II's. Mter simple distributiontheoretical transformations we get that
f
cp(t)'1(t) dt
= rr1 'PV
f
cp'(s) ds. s t
The last functional has been described in detail in Section 6.7 and the results of that section imply that '1(t) =
'PV
(t; ).
16. Applying the Fourier transform to the last equality we get an algebraic equation CU
>
o.
(3)
It has no ordinary continuous nonzero solutions ;j, (cu), but there are distributional singular solutions like solution (1.5.5) of the equation (1.5.6). Indeed, if CUn is a root of the "dispersion equation" i = e;wT. then 8(cu  cun ) is a solution of equation (3). Clearly, CUn
=
(i + 2rrn).
so that we obtain cp(t) = exp
(iif)
n = 0.1.2•...•
exp
t).
(4)
where an are arbitrary constants. In other words, we have found solutions
'"2 Tt) t(t).
.rr cp(t) = exp (
(5)
where t (t) is an arbitrary periodical analytic signal with period T. Furthermore, since Ii is a real operator, both Re cp(t) and 1m cp(t) are also solutions to this exercise. The simplest example here is obtained for t(t) == 1. Then cp(t) = em,. Recp = cosOt. Imcp = sin Ot. 0 = rr /2T. A little more complicated example is obtained setting an = an. la I < 1. Then, forinstance, ""( cos(Ot)  a cos(30t) Re'l' t) = 2 • 1 + a  2a cos(40t)
Fractal calculus 17. As it often happens, it is simpler to derive a more general formula involving the integral
(6) This is a linear functional of g(t) satisfying the causality principle and as such it has an integral representation x(t) =
h(t.s)g(s)ds
308
Answers and solutions
the kernel thereof can be found by substituting g(t)
= l5(t 
s) in (6):
In particular
18. We shall utilize the identity g(b) =
1
g(u) BBu X(u  b) duo
Substituting it in the original integral and changing the order of integration we get
1=
1 [I . i '!.
g(u) BBu
a(xt. ... , xn)X (u  b(Xl, ... , xn) )dXl ••. dXn] duo
Note that the inner ntuple integral is an integral of function a over (not necessarily connected) domain Cu = An Bu created by intersection of domains A and Bu, where the latter is the set of points satisfying the inequality
Denote that integral by 1/I(u) = I·'!·
[u
a(xl, ..• ,Xn)dxl ... dxn.
So, the desired formula (also called the Catalan formula) has the form
1= Ig(U).!!.1/I(U)dU. du
If b(Xl, ... , xn) is a bounded continuous function with minimum m and maximum M then the Catalan formula simplifies to
I =
A.7
i
M
g(u) d1/l(u).
Chapter 7. Uncertainty principle and wavelet transform
Function spaces:
1. Consider an auxiliary function l(x) = j O.
m=l
Selecting an even function I(t), with values at t = ml!. coinciding with those of the function 1 in the above series, the equality can be rewritten in the form
1 1 S(l!.) = 2/(0) + 2
L 00
I(m l!.) ,
m=oo
which is more suitable for an application of the Poisson summation formula. The latter gives
Since the Fourier image i(w) of an even function I(t) is also even, the above equality can be rewritten in the form
S(l!.)
=
rr1 1(0)  2/(0)
L
2rroo_ I(Zrrnj l!.).
+ t;
n=l
319
Chapter 8. Summation of divergent series and integrals
If j(w) has a compact support, then the series on the right has only finitely many terms different from zero. In our case it is convenient to choose f(t)
sint = 1:1, t

f(w)
1:1 = 2[X(w + 1) 
X(w 1)].
Substituting these function in the preceding equality we get
1 $(1:1) = 2(11' x)
00
+11' LX(1 21rn/I:1), n=l
so that finally
f
m=l
sin(ml:1) = !(1r  1:1) + 11' m 2
21r
'
where LxJ is the largest integer x. Note that the expression on the left hand side represent the Fourier series of the periodic function appearing on the righhand side. 11. Counting on the curiosity of the reader we take a look at a more general problem of the Fourier integral (11) where, in the case of a 21rperiodic function f(t), the Fourier image appearing inside the integral
L 00
j(w) =
(12)
 m).
m=OO
If we replace (11) by the equality fN(t) =
f
(13)
j(w)h(w, N)e ifUt dw
where h(w) is the rectangular function, for example, h(w, N)
= n(w/N),
n(t) = X(t
+ 1) 
X(t 1),
(13a)
then the series (12) retains only the needed number 2N + 1 of terms. The corresponding original function is (13b) h(t, N) = 2sin(Nt)/t. Using formula (3.2.6), we can rewrite the equality (13), with the help of the convolution fN(t)
= 21r1 f(t) * h(t, N).
(14)
By (3.3.7), the function h(t, N)/21r weakly converges, as N _ 00, to the Dirac delta. So, if f(t) is sufficiently smooth, then fN(t) converges pointwise (and even uniformly) to f(t). However if, as is the case in this exercise, the original function f(t) is only piecewise continuous, then it is necessary to study in more detail the asymptotic behavior (for N » 1) of the convolution integral (14) in the vicinity of discontinuities of the first kind.
Answers and solutions
320
Assume that the piecewise continuous function f(t) has a jump at t = lfl f(T + 0)  f(T  0). Hence, asymptotically, as t + T,
=
 0)
f(t) '"
+ f(T + 0)) + lfl sign (t 
T
of size
T).
Substituting this expression into (13), we find that the behavior of fN(t) in the vicinity of the discontinuity point is described by the following asymptotics: fN(t) '"
 0)
tT
+ f(T + 0)] + lfl Si (x),
X
="N'
(15)
where the integral sine function Si (z) =
 dx. Iooz sinx x
(16)
Observe certain features of the asymptotic formula (15). First of all, Si (0) = 0, which means that at the discontinuity point t = T of function f(t), the truncated Fourier integral (13) (and in our case, the series in Exercise 11) converges to the arithmetic average of the onesided limits of the function. Secondly, if we remove the already analized first summand on the righthand side of (15), place the discontinuity at the origin t 0 (T = 0), and rewrite (15) in the form
=
fN(t) '" lflg(tIN),
g(x) =
..!.7r Si(x)
then a new, important in physics and engineering phenomenon can be observed. The odd function g(x) entering in the above formula, normalized by the size lfl of the jump, describes fN(t) behavior in the vicinity of the discontinuity. For x + ±oo, we have g(x) + ±1/2. However, since the integrand in (16) changes sign, the approach of g(x) to the limit is not monotone. In particular, as is clear from the the graph of function g(x) shown on Fig. A.8.2, g(x) assumes the maximal value g(7r) 0.59 for x 7r (t 7rIN). This means that, for arbitrarily large N, at the distance 7r INto the left and to the right of the jump point of function f(t), the graph of function fN(t) has sharp up and down excursions. This anomalous behavior of function fN(t) in the neighborhood of jump points of function f(t) is called the Gibbs phenomenon.
=
=
12. As before, we shall rely on formula (14). Observe that the Gibbs phenomenon was caused by the fact that function h(t, N) entering (14) changed signs. So the phenomenon can be avoided if we select a nonnegative function h(t, N) such that its Fourier image has a compact support. The latter is needed so that only finitely terms of the series from Exercise 11 are different from zero. So let us consider function h(t, N) = A 4Sin;2(Nt)
(17)
as a candidate. Coefficient A will be selected later from a normalization condition, which, as is clear from (15), takes the form
f
h(t, N)dt = 1.
321
Chapter 8. Summation of divergent series and integrals
g(x) 0.6 0.5 0.4 0.3 0.2 0.1 0
5
10
15
20
25
30
X
FIGURE A.S.2 Graph ofthe function (1/7r) Si (x). First of all, let us check that the Fourier image of h has a compact support. In view of (3.2.10), the Fourier image of the square ofthe original function is
h(w, N) = AfI(w/ N)
* fI(w/ N).
The convolution of two functions with compact support also has compact support. In our case h(w, N) 2AN { 1 Iwl/2N, for Iwl 2N; (19) 0, for Iwl > 2N. With this information it is easy to find the correct norming A, which has to be such 1. Thus A 1/2N. It is not difficult to show that the above function that h(O, N) h(t, N)/27r weakly converges to the Dirac delta and thus guarantees the convergence of fN(t) to f(t) at the continnuity points of the latter. On the other hand, the positivity of function (18) removes the Gibbs phenomenon. Finally, it should be mentioned that function (19) substituted in (12) preserves in the sum (13) 4N  3 terms: 2N  1 on the left and on the right of the terms m = O. To keep only N terms of both sides of m = 0 one has to replace N by (N + 1) /2 in the preceding formulas. As a result we obtain that
=
=
=
h(
w,
N) _ { 1  Iwl/(N  0,
+ 1),
for for
Iwl Iwl
2N;
> 2N;
h t N = _4_sin2«N + l)t/2) (, ) N+1 t2 The corresponding sum in the formulation of Exercise 12 assumes then the familiar shape fN(t) =
E N
m=l
(
m ) sin(mt) 1   N +1 m
(20)
322
Answers and solutions
of the Cesaro sum for the series from Exercise 11. Fig. A.8.3 shows the graphs of functions fN(t) sin(mt)/m obtained by the simple summation of the first N 25 terms of the Fourier series, and by the Cesaro summation method. The plots clearly shows how the Cesaro method helps to avoid the Gibbs phenomenon.
= r::=1
=
f
(t)
1t
FIGURE A.8.3
t
r::=1
TbegrapbsoffunctionsfN(t) = sin(mt)/mobtainedbythesimplesummatioD ortbe first N = 25 terms of the Fourier series, and by the Cesaro summation method. 13. It sufficies to select

8
1/Io(w) = TI«w  WO)/O)· 1C
The inverse Fourier transform 1/Io(t)
. 215 = e"·'Ot  sin(nt), 1C
substituted into (8.9.16), gives F(t)
= 2c5 eiwot 1C
m=oo
f(mc5) sin(nt  mc5n) eimwo6 tmc5
(21)
Recall that one of the features of the analytic signal that is attractive for the engineers is that it uniquely determines its amplitude A(t) and phase q,(t) so that f(t)
= A(t) cos(wot + q,(t».
To see this it is sufficient to write the complex signal in the form F(t)
= A(t) exp(iwot + q,(t»,
Chapter 8. Summation of divergent series and integrals
323
separate the real and imaginary parts F(t) = le(t) cos(wot)
+ ils(t) sin(wot),
and compare this equality with with the previous one to obtain the slowly evolving in time envelope of the narrowband signal. 14. Initially, let us construct an appropriate class of sufficiently smooth Fourier images
{fr(w) such that the corresponding origianal function 1/I(t) decays to zero (as It I 4 (0) faster than 1/10 (t) (8.9.6), and such thatthe series (8.9.16) converges faster than the standard
Shannon's series (8.9.17). Our experience with the Fourier transform suggests that it is useful to write {fr (w) in the form of the convolution 1/I(w) = 1/Io(w)
l::!.
J.t * ;P(J.tw/rc) rc
of the rectangular function l::!., where {fro(w) is given by the equality (8.9.6), and of the "enveloping" it compactsupport function 1.t;P(l.tw/rc)/rc. The constant I.t is determined from the equality rc/I.t =  1/l::!.). If ;p(v), q:{r) being the original function, vanishes identically outside the interval Iv I ::: 1 and satisfies norming condition qJ(O) = 1, then {fr (w) is a function of type (8.10.5). With such choice of {fr (w), the function of interest sin '[' 1/I(t) = qJ(l::!.'['/I.t), '[' = rct/l::!.. l::!. '[' As ;p(v) one can take, for example, the sufficiently smooth function (0) = { (4/3) cos4 (rcO/2), qJ 0,
for 101 < 1; for 101 > 1;
the graph thereof is shown of Fig. 4.3.3b. Then, in view of (4.3.22), qJ('[')
sin '[' = 4rc 4 '['('['2 _ rc2)('['2 . _ 4rc2)
15. First, let us find the largest value of I which determines the length of the interval between readings (8.9.21). In our case, it follows from (8.9.22) that I < 9/2. Therefore, we choose I = 4. Thus, by (8.9.22), = 9rc /2WO. 16. Taking advantage of the freedom to choose arbitrarily values of the Fourier image Iw ± wol <  0, the widths thereof are p = wo/45 in
{fr(w) in the intervals 0 <
our case, let us smooth out the Fourier image {fro(w) (8.9.24) with the help of convolution with the function {fr(w) = {fro(w)
* ;P(l.tw/rc),
I.t
= 9Orc/wo.
Calculate the inverse Fourier image with the help of (8.9.25) to obtain
= !!:.qJ(rct / I.t) sin(Ot) cos(wot). rc rct = 9rc/2wo, I.t = 9Orc/WO.
1/I(t)
Recall, that 0
= WO/10,
AppendixB Bibliographical Notes
The history of distribution theory and its applications in physics and engineering goes back to [1] O. HEAVISIDE, On operators in mathematical physics, Proc. Royal Soc. London, 52(1893),504529, and 54 (1894),105143, and [2] P. DIRAC, The physical interpretation of the quantum dynamics, Proc. Royal Soc. A, London, 113(19267), 621641. A major step towards the rigorous theory and its application to weak solutions of partial differential equations was made in the 1930s by [3] J. LERAY, Sur Ie mouvement d'un liquide visquex emplissant l' espace, Acta Mathematica 63 (1934), 193248. [4]
R. COURANT R., D. HILBERT, Methoden der Mathematischen Physik, Springer, Berlin 1937.
[5] S. SOBOLEV, Sur une theoreme de 1'anayse fonctionelle, Matemat iceski Sbornik 4 (1938),471496. The theory obtained its final elegant mathematical form (including the locally convex linear topological spaces formalism) in a classic treatise of [6] L. SCHWARTZ, Theorie des distributions, vol. 1(1950), volII (1951), Publications de 1'Institut de Mathematique de L'Universite de Strasbourg, which reads well even today. In its modem mathematical depth and richness the distribution theory and its application to Fourier analysis and differential equations can be studied from many sources starting with massive multivolume works by [7] I.M. GELFAND et al. Generalized functions, 6 volumes, Moscow, Nauka 19591966. © Springer Nature Switzerland AG 2018 A. I. Saichev and W. Woyczynski, Distributions in the Physical and Engineering Sciences, Volume 1, Applied and Numerical Harmonic Analysis, https://doi.org/10.1007/9783319979588
Bibliographical notes
326
and [8] L. HORMANDER, The Analysis of Linear Partial Differential Op erators, 4 volumes, Springer, BerlinHeildelbergNew YorkTokyo 19831985, to smaller, one volume monographs from research oriented [9] E.M. STEIN, G. WEISS, Introduction to Fourier Analysis on Eu clidean Spaces, Princeton University Press 1971, [10] L.R. VOLEVICH, S.G. GINDIKIN, Generalized Functions and Con volution Equations, Moscow, Nauka 1994, to, textbook style [11] R. STRICHARTZ, A Guide to Distribution Theory and Fourier Trans form, CRC Press, Boca Raton 1994, [12] V.S. VLADIMIROV, Equations of Mathematical Physics, Moscow, Nauka 1981, An elementary, but rigorous, construction of distributions based on the notion of equivalent sequences was developed by
[13] J. MIKUSINSKI, R. SIKORSKI, The Elementary Theory ofDistribu tions, I (1957), II (1961), PWN, Warsaw. [14] P. ANTOSIK, J. MIKUSINSKI, R. SIKORSKI, Generalized Func tions, the Sequential Approach, Elsevier Scientific, Amsterdam 1973. The applications of distribution theory have appeared in uncountable physical and engineering books and papers. As far as more recent, applied oriented textbooks are concerned, which have some affinity to our book, we would like to quote [15] F. CONSTANTINESCu,Distributions and Their Applications in Physics, Pergamonn Press, Oxford 1980. [16] T. SCHUCKER, Distributions, Fourier Transforms and Some ofTheir Applications to Physics, World Scientific, Singapore 1991. which however, have a different spirit and do not cover some of the modem areas covered by our book. The classics on Fourier integrals are [17] S. BOCHNER, Vorlesungen aber Fouriersche Integrale, Akademische Verlag, Leipzig 1932, [18] E. C. TITCHMARSCH, Introduction to the Theory ofFourier Integrals, Clarendon Press, Oxford 1937, with numerous modem books on the subject, including the above mentioned monograph [9] and elegant expositions by
Bibliographical notes
327
[19] H. BREMERMAN, Complex Variables andFourier Transforms, AddisonWesley, Reading, Mass. 1965, [20] H. DYM, H.P. McKEAN, Fourier Series and Integrals, Academic Press, New York 1972, [21] T. W. KORNER, Fourier Analysis, Cambridge University Press 1988. The asymptotic problems (including the method of stationary phase) discussed in this book are mostly classical. The well known reference is e.g. [22] N.G. DE BRUIJN,Asymptotic Methods in Analysis, NorthHolland, Amsterdam 1958, with the newer reference being [23] M.B. FEDORYUK, Asymptotics, Integrals, Series, Nauka, Moscow 1987. The special functions have a rich literature including [24] H. BATEMAN, A, ERDELYI, Higher Transcendental Functions, 2 volumes, McGrawHill, New York 1963, [25] F.W.J.OLvER,AsymptoticsandSpeciaIFunctions,NewYork1974, and their connections with harmonic analysis on groups are explained, e.g., in [26] N. VA. VILENKIN, Special Functions and Group Representations, Nauka, Moscow 1965.
As always, it is handy to keep around [27] JAHNKEEMDE, Tables ofHigher Functions, TeubnerVerlag, Leipzig 1960, [28] LS. RYZHIK, LM. GRADSHTEYN, Tables ofIntegrals, Sums, Series and Products, FM, Moscow 1963, [29] M. ABRAMOWITZ, LA. STEGUN,HandbookofMathematicaIFunc tions, National Bureau of Standards, 1964 although the role of such compendia has been recently diminished with introduction of computer symbolic manipulation software such as Mathe11Ultica and Maple. A good source on the mathematical theory of singular integrals is [30] E.M. STEIN, Singular Integrals and Differentiability Property of Functions, Princeton University Press 1970, with vast literature spread through mathematical journals. The well known sources on fractal calculus are [31] A.H. ZEMANIAN, GeneralizedIntegral Transformations, Interscience, New York 1968,
328
Bibliographical notes [32] K.B. OLDHAM, J. SPANIER, The Fractional Calculus. Theory and applications of Differentiation and Integration to Arbitrary Order, Academic Press, San Diego 1974, [33] A.C. McBRIDE, Fractional Calculus and Integral transforms of Generalized Functions, Pitman, London 1979,
and more recent advances in the area can be gleaned from the collection of research papers [34] A.C. McBRIDE, G.F. ROACH, Editors, Fractional Calculus, Research Notes in Mathematics, Pitman, Boston 1985. The wavelets have obtained recently several excellent expositions, and the reader can benefit from consulting [35] Y. MEYER, Wavelets and Operators, Cambridge University Press 1992, [36] I. DAUBECHIES, Ten Lectures on Wavelets, SIAM, Philadelphia 1992. [37] G. KAISER,A Friendly Guide to Wavelets, BirkhiiuserBoston 1994. This is a very active research area and some new results have appeared in the following volumes of articles [38] I. DAUBECHIES, Editor, Different Perspectives on Wavelets, American Mathematical Society, Providence, R.I. 1993, [39] C.K. CHUI, Editor, Wavelets: A Tutorial in Theory and Applications, Academic Press, New York 1992. The latter contains articles by D. Pollen on construction of Daubechies wavelets and G.G. Walter on wavelets and distributions which were used in Chapter 7. An engineering perspective can be found in [40] A. COHEN, R.D. RYAN, Wavelets and Digital Signal Processing, Chapman and Hall, New York 1995. The classic text on divergent series is [41] G.H. HARDY,DivergentSeries, Clarendon Press, Oxford 1949, but the problem has broader implications and connections with asymptotic expansions and functional analytic questions concerning infinite matrix operators, see e.g. [42] R.B. DINGLE, Asymptotic Expansions, Academic Press, New York 1973, [43] I.J. MADDOX ,InjiniteMatricesofOperators, SpringerVerlag, Berlin 1973. The modem viewpoint is presented in [44] B. SHAWYER, B. WATSON, Borel's Methods ofSummability: The ory andApplications, Clarendon Press, Oxford 1994.
Bibliographical notes
329
Finally, an exhaustive discussion of the Shannon's sampling theorem and related interpolation problems can be found in
[45] R.J.
MARKS II, Introduction to Shannon's Sampling and Interpola tion Theory, Springerverlag, Berlin 1991.
Index
Abel method of summation, 254 Abramovitz, M., 177,327 Achilles, 245 tiring, 255 Anger function, 146 Antosik, P., 326 approximation Born, 61 diffusion,91 Fresnel, 14, 141 geometric optics, 142 weak,1O asymptotic equivalence, 93 asymptotic expansion, 103 principal term of, 104 asymptotics general scheme, 137ff of fractal relaxation, 178 autocorrrelation function, 218 Bateman, H., 327 beads on a string, 37 beetle on a string, 96 Bernoulli numbers, 266 Bessel function, 102, 204 beta function, 169 Bochner, S., 326 Bremerman, H., 327 camel passing through needle's eye, 94 Catalan formula, 308 Cauchy criterion, 248 Cauchy density, 12 Cauchy formula, 157, 166
Cauchy integral principal value of, 152 causality principle, 39, 153 Cesaro method of summation, 254, 261 characteristics, 51ff, 86 Chui, C.K., 328 Cohen, A., 328 Constantinescu, F., 326 continuity equation, 70 in Lagrangian coordinates, 53 of continuous medium, 47 of single particle, 44 convergence absolute, 250 conditional, 250 uniform, 8, 101 weak, 10 convolution, 27, 80 coordinates, Eulerian, 50 Lagrangian, 50 Courant, R., 325 d' Alembert solution, 44 Daubechies, I., 328 Daubechies wavelet, 231ff scaling relation, 232 Dawson integral, 151 de Bruijn, N.G., 327 decibels, 107 deltafunction, 3 derivative fractal 17Off substantial, 52
© Springer Nature Switzerland AG 2018 A. I. Saichev and W. Woyczynski, Distributions in the Physical and Engineering Sciences, Volume 1, Applied and Numerical Harmonic Analysis, https://doi.org/10.1007/9783319979588
332
differential equations fractal, 175 ordinary, 39 with singular coefficients, 60 with timedependent coefficients, 41 diffusion approximation, 88 dilation, 225 Dingle, R.B., 328 dipole, singular, 28, 38 Dirac, P., 6, 325 Diracdelta, 6, 9 analytic representation of, 160 on lines, 28ff on R n , 28ff on surfaces, 28ff periodic, 268 selfsimilarity of, 26 with nonmonotonic argument, 61 Dirichlet integral, 83 discontinuities of the first kind, 10 1 of the second kind, 112, 131 dispersion relation, 162 distribution convolution of, 27 derivative of, 15 direct product of, 65 elementary construction of, 326 integral of, 20 multiplication of, 18 nonlinear transformation of, 63 of composite argument, 24, 61 on finite interval, 58 orthogonality of, 238 parameterdependent, 64 principal value, 150ft' regular, 9 singular, 9 supersingular, 65 support of, 10 tempered, 82 vectorvalued, 35 dual space, 6 Duhamel integral, 41 Dym, H., 327 electrodynamics, xii
Index
energy density, 198 total,80 equation characteristic, 51, 86 divergence form of, 46 Liouville, 290 Riemann, 147 wave,70 equivalence relation (asymptotic), 95 Er&lyi, A., 325 Euler constant, 97, 131 integral, 115 Eulerian coordinates, 50 Fedoryuk, M.B., 327 filtering, 11 0 Fourier image, 75 Fourier transform, 75 asymptotics, 93ff discrete, 271 generalized, 81 inverse,78 inverse windowed, 206 of convolution, 80 of derivative, 78 of fractal derivative, 174 of Heaviside function, 125, 163 of power function, 112, 123 of smooth function, 78 redundancy of windowed, 207 windowed, 193ff, 309 fractal (fractional) differential equation, 176 differentiation, 170ft', 307 causality, 174 nonlocal character, 174 scale invariance of, 174 integration, 166ff, 307 iteration of, 169 Laplacian, 174 relaxation, 175 frequency, 77 angular, 75 current, 202 instantaneous, 194 localization, 191
333
Index
mean, 199 relaxation, 175 standard deviation of, 199 Fresnel approximation, 141ff function( s), Anger, 146 autocorrelation, 219 beta, 169 Bessel, 100,205 Cauchy,12 error, 178 Gamma, 112, 115 Gaussian, 11 Green's, 39, 49 Heaviside, 21 infinitely differentiable, 7 jumpof,22 locally integrable, 9 Lorentz, 12 MittagLeffler, 177 of the order not greater than ..., 94 of the order smaller than ... , 94 of the same order, 94 power,80 pulse,I09 rapidly decreasing, 82 Riemann's zet, 267 scaling, 228 stream,57 transfer, 39 unit step, 21 Weber, 146 wideband ambiguity, 218 windowing, 195 with compact support, 7 functional analysis, 189 continuous, 7, 8 linear, 5, 8 Gabor transform, 198 Gelfand, I.M., 325 geometric optics approximation, 142 Gibbs phenomenon, 285,320 Gindikin, 5.G., 326 Gradshteyn, I.M., 327 Green's function, 39,49
Haar wavelet, 226ff location, 227 orthogonality, 227 resolution level, 227 scaling relation, 231 selfsirnilarity, 231 Hardy, G.H., 328 harmonic oscillator, 41 with damping, 41 Heaviside, 0., 325 Heisenberg uncertainty principle, 192 Hermitian property, 186 Hilbert space, 183 completeness of, 189 transform, 160ff, 303 Hilbert, D., 325 Holder inequality, 170 Horrnander, L., 326 incompressible medium, 55 initialvalue problem for ID wave equation, 43 integral collision, 84 cosine, 130, 305 Dawson, 151 Dirichlet, 81, 85 divergent, 262ff, 309 Duhamel,41 Euler,115 fractal, 166 Fresnel sine and cosine, 142 Lebesgue, 187 principal value of, 125 scattering, 84 sine, 305, 3201 singular, 149ff isomorphism, 54 JahnkeEmde, 327 Jordan Lemma, 116 jumps of the first kind, 101 Keiser, G., 328 kernel,6 Korner, T.W., 327
334
Lagrangian coordinates, 50 Laplace operator, 56 Laplace's method, 145 Lebesgue integral, 188 spaces LP, 190 Leibnitz formula, 17ff lemma RiemannLebesgue, 98, 132 Leray, J., 325 linear combination, 17 linear functional, 6, 9 Liouville equation, 289 Lipschitz condition, 150 loglog scales, 105 Maddox, LJ., 328 majorant, 103 Marks II, R.J., 329 mass conservation law, 46 McBride, A.C., 328 McKean, H.P., 327 medium incompressible, 55 Meyer, Y., 328 Mikusinski, J., 326 multiresolution analysis, 225 norm, 184 observables momentum, 193 position, 193 Oldham, K.B., 328 Olver, F.W.J., 327 operator Laplace, 56 linear, 77 Parceval equality, 80, 192 particles sticky,66 partition of unity, 233 passive tracer, concentration, 54 density, 54 phase space, 46, 84
Index
transition, 147 Poisson summation formula, 271ff, 313ff Pollen, D., 232, 328 potential scalar, 56 vector, 56 powertype behavior, 107 principal value of Cauchy integral, 152 of integral, 126, 149 principle causality, 40, 153 of asymptotic attenuation, 157 of infinitesimal relaxation, 156, 255ff superposition, 77 uncertainty, 183ff, 190 probing property, 4 multiplier, 19 recursive method, 175 redundancy, 207,223 regime multistream, 67 singlestream, 66 relaxation asymptotics, 178 fractal in frequency, 179 fractal in time, 179 Kohlrausch, 179 of order 112, 177 remainder term, 103 renormalization, 127 Riemann equation, 147 Riemann Theorem, 250 RiemannLebesgue Lemma, 98, 132 Riemann's zeta function, 267 Riesz Theorem, 167 Roach, G.F., 328 Ryan, R.D., 328 Ryzhik, I.S., 327 Schucker, T., 326 Schwartz inequality, 169, 185, 189, 207, 215 Schwartz, L., 8, 325 selfsimilarity of Direc delta, 26 of function, 5
335
Index
of Haar wavelets, 231 separation of scales, 260ff series Abel's summability of, 254, 261 Cesaro's summability of, 254, 261 convergence of, 245ff absolute, 249 Cauchy criterion, 248 conditional, 249 pointwise, 252 unconditional, 250 uniform, 252 Weierstrass criterion, 253 divergent, 245ff, 312 geometric, 273ff summability of, 253 functional, 252 tpsummability of, 260ff geometric, 245, 248 generalized sum of, 253 harmonic, 96, 249 of complex exponentials, 264ff Shannon's sampling theorem, 276ff, 320 Shawyer, B., 328 signals analytic, 162, 303 and systems, xii cornplexstructured, 215 input,39 narrowband, 162 output,39 smoothing and filtering of, 110 timewindowed, 195 Sikorski, R., 326 simplest tune, 193 smoothing, 109 Sobolev, S., 8, 325 SobolevSchwartz space, 8 solution d' Alembert, 44 fundamental, 39 space,
S,82 S',82 Sobolev, 238 SobolevSchwartz, 8 Spanier, J., 328 stationary phase method, 140 accuracy of, 143 stationary point, 140 steepest descent method, 145 Stegun, I.A., 327 Stein, E.M., 325, 327 Stirling formula, 147 stream function, 57 Strichartz, R., 326 substantial derivative, 52
V,8 V',8
Taylor formula, 264 test function, 5 theorem Abel's, 254 Cauchy's, 115, 157 Fubini's, 65 fundamental of vector calculus, 56 Gauss divergence, 47 Leibniz, 249 mean value, 261 Riemann , 250 Shannon sampling, 276ff Weierstrass, 99 timefrequency localization, 190 Titchmarsch, B.C., 326 topology dual,32 weak, 32 tortoise, 245 transfer function, 39 transform Fourier,77 Gabor, 198 Hilbert, 160,303 transport equation, 84ff phenomena, xii triangle inequality, 185
Hilbert, 183 ff inner product, 184 linear topological, 31 phase,46
uncertainty principle, 183ff Heisenberg, 192 strong, 215
336
Vilenkin, N. Ya., 327 Vladimirov, V.S., 326 Volevich, L.R., 326 Walter, G.G., 326 Watson, B., 328 wave monochromatic, 153 propagation, xi ID equation, 43 wavelet admissible, 223 Daubechies wavelet, 231ff scaling relation, 232 Haar, 226ft" location, 226 orthogonality, 226 resolution level, 226 scaling relation, 231 selfsimilarity, 231 image,210 Mexican hat, 213 Morlet,211 efficiency factor of, 212 mother, 210ft" complexstructured, 217 onesided, 221 twosided, 221 Poisson, 222 transform, 210ft", 310 inversion, 219 wavenumber, 154 Weber function, 146 Weiss, G., 327 window finite memory, 195 relaxation, 195 Gaussian, 195 Zemanian, A.H., 327 Zeno's paradox, 245
Index
When life gives you a hundred reasons to cry, show life that you have a thousand reasons to smile
© Copyright 2015  2023 AZPDF.TIPS  All rights reserved.