Document text
The Monogenic Synchrosqueezed Wavelet Transform: A tool for
<n! the Decomposition/Demodulation of AM-FM images
o
(N
>
O
o
(N
<
>
00
M. Clausel, T. Oberlin, V. Perrier*
November 22, 2012
Abstract
The synchrosqueezing method aims at decomposing ID functions as superpositions
of a small number of "Intrinsic Modes", supposed to be well separated both in time
and frequency. Based on the unidimensional wavelet transform and its reconstruction
properties, the synchrosqueezing transform provides a powerful representation of multi-
component signals in the time-frequency plane, together with a reconstruction of each
mode.
In this paper, a bidimensional version of the synchrosqueezing transform is defined,
by considering a well-adapted extension of the concept of analytic signal to images: the
monogenic signal. The natural bidimensional counterpart of the notion of Intrinsic Mode
is then the concept of "Intrinsic Monogenic Mode" that we define. Thereafter, we inves-
tigate the properties of its associated Monogenic Wavelet Decomposition. This leads to a
. natural bivariate extension of the Synchrosqueezed Wavelet Transform, for decomposing
and processing multicomponent images. Numerical tests validate the effectiveness of the
method for different examples.
| Keywords: monogenic signal, wavelet transform, directional time-frequency image analysis,
synchrosqueezing .
MSC 2010: 65T60, 92C55, 94A08.
• i-H -
X.
1 Introduction
In the past few years, there has been an increasing interest in representing images using
spatially varying sinusoidal waves. As the understanding of the theory advanced, amplitude-
and frequency-modulation (AM-FM) decompositions have been applied in a large range of
problems: for example motion estimation using a flow-optical method based on an assumption
of phase-invariance [H [3] , reconstruction of breast cancer images [U [5] , texture analysis [6j
E], or ultrasound image segmentation methods [8] based on the concept of monogenic signal
and quadrature filters [SJ. A survey of applications of AM-FM decompositions in medical
imaging can also be found in [lOj .
In each case, the main challenge is to decompose any input image s(x) into a sum of
bidimensional AM-FM harmonics of the form
L L
(1) s(x 1 ,x 2 ) = y^^t) = } y Aj(x 1 ,x 2 ) cos(y»f (g?i, x 2 )) ,
i=i i=i
"University of Grenoble and CNRS, Laboratoire Jean Kuntzmann UMR 5224, Saint Martin d'Heres, France.
1
where An > denotes a slowly-varying amplitude function, tpi denotes the phase, and t =
1, • • • ,L indexes the different AM-FM harmonics. To each phase function, one can associate
an instantaneous frequency vector field defined as uj£ = V<pg. Finding the components S£ from
the bidimensional signal s is called the decomposition problem.
When the decomposition is given, another subsequent problem is to determine for each
component its amplitude, phase and frequency functions involved in ([1]). The amplitudes
can be related to the energy contained at each point of the image, whereas significant texture
variations are captured in the frequency content. For example, in the case of a single compo-
nent, the instantaneous frequency vectors are orthogonal to the isovalue intensity lines of an
image, while their magnitudes provide a measure of local frequency content. The problem of
estimating the amplitudes, phases and instantaneous frequencies of a signal of the form (pQ)
(bidimensional or not) is called the demodulation problem.
In the one dimensional case, these two problems have been widely investigated both from
the theoretical and practical point of view. To address the demodulation problem in the one
dimensional setting, a first and necessary step is to define, in a proper way, the concept of
amplitude, phase and instantaneous frequency of a given signal. To this end, one can use the
well-known Hilbert transform, defined for any / 6 L 2 (R) as
For a.e. t G M, Hf(t) = lim ( - [ ds
e^O y-K J\t- S \ >£ t - S
The analytic signal associated to / is then the complex- valued function F = {1 + \H)f .
Thereafter, one defines amplitude A and phase (p of / as, respectively, the modulus and
argument of the analytic signal F associated to / (the uniqueness of phase being ensured
under smoothness assumptions on / and under some initial conditions of the form <p(0) = ipo).
One then obtains
For a.e. t 6 R, F(t) = A{t)e lLp{t) .
The instantaneous frequency of / is thus the derivative of the phase ip' [llj . The practical
estimation of amplitude, phase and instantaneous frequency was the subject of numerous
studies (see [12] for a survey). The "naive" method, which consists in using the Hilbert trans-
form to estimate A and tp, is known to be numerically unstable [13]. Alternative methods use
the wavelet transform associated to reallocation techniques: in the one dimensional case, the
time-frequency localization of wavelets allows to compute robust estimations of instantaneous
frequency lines in the time-frequency representation: several approaches have been developed
in the 90th, one should cite the method of "wavelet ridges" (see the two pioneer works [14J
and [H>]), the squeezing method introduced in |16J or the reassignment method [17j .
On the other side, the decomposition problem consists in developing decompositions of
any signal into a sum of "well-behaved" AM-FM components, and was widely studied. A
recent attempt in this context is Empirical Mode Decomposition (EMD), initially introduced
by Huang et al. in [18] and later popularized by Flandrin and his co-authors (see e.g. [19]
or [20]).
Basically, EMD is the output of an iterative algorithm which provides an adaptive decom-
position of any signal into several AM-FM components. It is now used in a wide range of
applications including meteorology, structural stability analysis, and medical studies [T2J. In
spite of its simplicity and its efficiency, this method is hard to analyze mathematically, since
it is defined in an empirical way. To tackle this problem, in [21], the authors propose and
analyze an alternative method called "synchrosqueezing" derived from reassignment methods:
2
they introduce a class of functions which can be viewed as the superposition of a reasonable
number of modes, i.e. which read Yle Ae(t)e ilpe ^ , where the amplitude Ag(t) of each mode
is slowly varying with respect to its instantaneous frequency wg(i) = (p'At). For each func-
tion belonging to this class, the synchrosqueezing provides an estimation of this (unknown)
decomposition. Synchrosqueezing presents many similarities with Empirical Mode Decompo-
sition as shown in [21] (see also [22J for a comparison of these two methods or [23] for more
about potential applications). In addition, since synchrosqueezing is based on the wavelet
transform, it can be used to address simultaneously the decomposition and demodulation
problem as done in |21j .
In this paper, we will focus on the two dimensional case and extend the synchrosqueezing
approach to images. Our first step consists in defining a convenient extension of the concept of
analytic signal, in order to properly define the notion of amplitude and phase of an image. Two
main generalizations of the analytic signal to the bidimensional setting have been considered in
the literature: the hypercomplex and the monogenic signal defined respectively in [21] and [9].
Here, we focus on the monogenic setting, since its characteristics lead to easier interpretation.
The monogenic signal associated to an image is defined using the Riesz transform. If we are
given a real valued function / 6 L 2 (R 2 ), one defines its Riesz transform IZf as the following
vector-valued function
where for any i = 1, 2 and a.e. i£l 2 one has
'£<■/(■'■) lim . I - I rrl/fe) dy
\x~ y\>e 1*^ y\ j
Then, the monogenic signal associated to / is
Mf =
f
Kf
The characteristics of the monogenic signal are usually defined using the quaternionic formal-
ism (see Appendix [A] for some background about quaternionic calculus and Appendix [B] for
the main properties of the Riesz transform): the monogenic signal Aif associated to / can
be read
Mf = f + n 1 fi + n 2 fj,
where (1, i, j, k) denotes the canonical basis of the algebra H of the quaternions. In [9j and [25J,
it is proved that the monogenic signal can be written as
Mf = A e^™ 8 = A (cos (p + ng sin ip) with ng = cos 9 i + sin 9 j ,
where A is a positive real valued function and <p and 9 are two real valued functions defined
in a unique way (under suitable assumptions). The function A is called the amplitude of /, (p
its phase and ng its local orientation. The function uj = Vip will be called the instantaneous
frequency of /. To illustrate these definitions, we consider the case where f(x) = Aq cos(k ■ x)
with k = (ki, ^2) an d k\ > 0. Set 9q = Arctan(A;2/A:i). One then has
\ a fsin(k 1 x 1 + k 2 x 2 ) cos 6> \ . k , .
\ v sm(A;iXi + ft 2 x 2 ) smfc'oy \k\
3
and
— A e (k-x)(cos0 i+sin6> j)
Then for all x G M 2
one can set
A(x) = Aq, <p(x) = k ■ x, 9{x) = 9q = Arctan(/c2/^i) •
The aim of the present paper is to define the synchrosqueezing of bidimensional images,
in the monogenic setting, in order to address both the decomposition and the demodula-
tion problems. To this end, we will consider monogenic signals which are superposition of
monogenic waves of the form
where A and 9 are slowly varying functions with respect to ip.
We will consider successively the two different problems of deconvolution and decompo-
sition. As in [21j, our method is based on the continuous wavelet transform. Therefore, in
Section [2J we recall a recent extension of the continuous bidimensional wavelet transform
to the monogenic setting. We then investigate in Section [3J the problem of estimating the
amplitude and instantaneous frequency of a monogenic spectral line of the form ([2]). Subse-
quently, we define in Section H] the bidimensional synchrosqueezing and its application to the
decomposition of multicomponent images. Finally, Section [5] gives some insights on the prac-
tical implementation of bidimensional synchrosqueezing, and provides numerical experiments.
Proofs are postponed to Section [6] to make reading easier.
2 2D-wavelet transform and monogenic wavelet transform
In the one dimensional setting, the synchrosqueezing method is based on the one-dimensional
wavelet analysis: indeed, the wavelet transform is used both to decompose an analytic signal
into several components, and also to estimate the amplitude and the instantaneous frequency
of each component. The main advantage of using the wavelet analysis is to avoid the use of
the discrete Hilbert transform, since there is a link between the wavelet coefficients of the
analytic signal associated to a real function /, and its analytic wavelet coefficients (i.e. the
wavelet coefficients of / among an analytic wavelet).
In the bidimensional case, it is then quite natural to wonder if this approach can be
extended using bidimensional wavelet analysis, replacing the analytic signal by the monogenic
signal associated to our data. Since the monogenic signal can be viewed as a 3D vector field,
it can be componentwise analyzed by the usual wavelet transform: Section 12.11 recalls first
some basics about the bidimensional continuous wavelet transform. Thereafter in Section \2.2\
we recall the main characteristics of the monogenic wavelet analysis, which is an extension of
the usual wavelet analysis introduced in [26] and |27| . We would like to point out that, for
any bidimensional real image /, the wavelet coefficients of the monogenic signal Mf can be
related to the monogenic wavelet coefficients of / (see Proposition 12 . 1 1 for a precise statement).
Then, using monogenic wavelet analysis, one can recover the characteristics of the monogenic
signal associated to our data, without estimating it.
In what follows, for any (a, 6, a) 6 x R x (0, 2ir), we denote by D a the dilation operator,
T;, the translation operator and R a the rotation operator defined on L 2 (R 2 ) by:
(2)
D a f(x)
a- 1 fix/a), T b f(x) = f(x - b), R a f(x) = f{r~ x x)
4
where r a is the usual 2x2 rotation matrix of angle a:
/cos a —sin a
(3) r a =
ysm a cos a
2.1 The usual bidimensional wavelet transform
We briefly recall classical definitions about the bidimensional continuous wavelet transform.
A (real or complex) function ip G L 2 (M. 2 ) is called an admissible wavelet if it satisfies 0:
(4) C^ = (2nf r ifigfu < +oc
As usual (see [28]), the wavelet family {^ a ,^&}( a ,a,6)eR* x(o,2tt)x]R2 is defined by ip a ^ b =
T b R a D a i\). One then defines the wavelet coefficients of any / G L 2 (M 2 ) as follows:
c f (a,a,b) = / f(x) tp a ,a,b(x) dx .
Jk 2
If the wavelet is assumed to be isotropic, the wavelet coefficients do not depend on a. In this
case, we will denote c/(a, b) = c/(a, 0, b) and Tp at b = ip a ,o,b-
Any square integrable function / can be recovered from its wavelet coefficients using the
so-called reconstruction formula:
(5) f(x) = ^- f f f c f (a,a,b) ip a>a ,b(x) ^da db ,
U V ^6eR 2 iaG(0,27r) JaG(0,+oo) a
where this equality stands in L 2 (M 2 ). When the wavelet is assumed to be isotropic, we get a
simpler expression:
f( x ) = TT [ [ c f (a,b) ip a ,b(x) ^dfe .
W JbeK 2 JaG(0,+oo) a
Other reconstruction formulas are available (see e.g. [28]), among which the pointwise recon-
struction formula, obtained by summing over scales only. In the isotropic case, it reads:
2?r f +co . . da . ~ f
(6) f(x) = 7 r- c f (a,x) -5- with = / d£ .
Jo a Jr 2 I? I
As in the one dimensional case [21], this formula will play a special role in the bidimensional
synchrosqueezing.
For any real function / G L 2 (M 2 ), one can also define the wavelet transform of its monogenic
signal F = Mf = f + K x f i + K 2 f j , as:
■ ( c <
c F = c f + c nif 1 + cn 2f j = I c nif
W2//
lr rhe Fourier Transform of if> being defined by: V>(£) = ^ J R2 tp{x)e '^' x dx
5
2.2 The monogenic wavelet transform
We now present the continuous monogenic wavelet transform as denned in We consider
(*\
a real admissible wavelet ip and define ip( M > = A4i[) = I TZiip I the associated monogenic
wavelet. Observe that %fj^ M ^ is also an admissible wavelet (taking values in R 3 ), since each of
its components satisfies relation (j4]). We define below the monogenic wavelet coefficients of a
real function.
Definition 2.1 Let f £ L 2 (R 2 ). The monogenic wavelet coefficients Cj^^a, 6) of f are
defined by:
(7) cf 4) (a,a,b)= [ f(x)^ a M J b (x)dx,
where for any (a, a, b) <G (0, +00) x (0, 2vr) x M 2 , = T b R a D a (Mip).
Remark 2.1 For any (a, a, 6), c^\a,a,b) is a Clifford vector, and thus can be written as
Cj M ^(a, a, b) = Cf(a, a, b) + c'p (a, a,b) i + c'p (a, a, b) j where we denote for i = 1, 2
cy(a,a,b) = / /(x) (T b R a D a Kiip)(x) dx ,
('see Appendix\M for more details).
The next proposition states that there exists an explicit relationship between the mono-
genic coefficients of a real image / and the usual wavelet coefficients of its monogenic signal
Mf.
Proposition 2.1 Let f £ L 2 (M 2 ) and denote by F = Aif the monogenic signal associated
to f. Then for any (a, a, b) G (0, +00) x (0, 2ir) x M 2
(8) cf(ci, a, b)
Proof. We will prove that
C/ M) (°' a ' 6 ) = (J Ci? ( a ' a ' 6 ) '
which is equivalent to (|8j). By definition of the monogenic wavelet transform one has:
c ( f M) (a,a,b) = [ f{x) T b R a D a (MiP)(x) dx .
The translation-, scale-invariance and steerability properties of the Riesz transform recalled
in Propositions IB.31 and IB~4l of Appendix [Bl imply that:
T h R a D a {n^) = r- l n{^ b ) .
6
Hence
(9)
f(x) T b R a D a (K*p)(x) dx = r~ l / f(x) n{i> a ^ b ){x) dx .
Since by Proposition IB. 51 the Riesz transform is a componentwise antisymmetric operator on
L 2 (R 2 ), one deduces that:
(Kf)(x) 1p a ,a,b{x) dx .
Equations ([9]) and ([10]) then directly imply relation
□
We now consider the case an admissible isotropic real-valued wavelet ip. Examples include
the Mexican hat (Laplacian of Gaussian) or the Morse wavelets (see |29j and Section [5]) . In
such a case, the wavelets tp a ,a,b as well as the wavelet coefficients cf(cl, a, b) of the monogenic
signal do not depend on the orientation a. From now, we then denote tp a ,a,b and cp(a,a,b)
respectively by ^> 0) & and cp(a,b). Thus equation ([8]) can be simplified in the following way:
c F (a,b)
-r,
a, a, b)
c f (a,b)
-TV
/ c f (a,b)
Cj\a, q, b)
\cf\a,a,b)
- cos a ci (a, a, b) + sin a (a, a, 6)
y— sin a Cy. (a, a, 6) — cos a ci 2 (a, a, 6) /
In practical situations, it will be easier to consider the case in which a = 0. We then get that
(12)
c F (a,b)
c ( , M) (a,0,6).
To sum up, in the isotropic case, the monogenic wavelet coefficients of the 2D real image / can
be related in a very simple way to the wavelet coefficients of the monogenic signal using f)12|) .
It will be useful in the sequel since when analyzing images, one deals with a real-valued signal
/ and aims at recovering the wavelet coefficients cp(a,6) of its monogenic signal F without
computing F.
3 Wavelet analysis of bidimensional spectral lines
In this section, we tackle the problem of estimating the amplitude and instantaneous fre-
quency of a bidimensional mode of the form A(x)e !fi ^ n ^ x \ In the one dimensional setting, the
"wavelet-ridge" method has proved to be efficient to address this problem |14tll5j. In |26tl29|
it has been extended to the bidimensional setting using monogenic wavelet analysis.
Here we follow an alternative approach consisting in extending synchrosqueezing to the
bidimensional context. We first define precisely the bidimensional modes that we are consid-
ering, and that we will call "Intrinsic Monogenic Mode Function" (IMMF). Then we deduce
the notion of bidimensional frequency vector and state an approximation result using the
monogenic wavelet transform.
Let us first define our concept of bidimensional modes:
7
Definition 3.1 Let e > 0. An Intrinsic Monogenic Mode Function (IMMF) with accuracy e
and smoothness a > is a function F G B^ T l (IR 2 , M) of the form
(13) F(x) = A(x)e^ x)nB ^ with n e(x) = cos(0(x))i + sin(#(x))j ,
where A>0andA,6£ C 1 (M 2 ) n L°° (R 2 ) , <p G C 2 (IR 2 ). The function A is called the amplitude
of F, whereas (p and n$ are called respectively the scalar phase and the local orientation of F.
The function ip is assumed to have bounded derivatives of order 1 and 2:
(14) Vi G {1,2}, < inf \d Xi <p(x)\ < sup \d Xi <p(x)\ < oo .
< oo .
(15) M = max sup d 2 . x . ip(x
Further we assume that A, 9, X7 2 tp are slowly varying functions with respect to V(p, namely
that for i = 1 , 2
(16) \fx G R 2 ,max(\d Xt A(x)\,\d Xt 9(x)\) < e\d Xi <p(x)\ .
and
(17) Vx G M 2 , . max n d*M x ) < ■
ii,«2=li2
)2
Remark 3.1 Observe that, in the sequel, we shall consider non-square integrable functions.
In this case, the two reconstruction formulas |2P and may not hold (see Appendix B
of }30^). We then need an additive assumption to ensure that the integral f + °° cf(o, 6)^
exists to give some sense to the reconstruction formula Moreover, we will estimate the
behavior of
/ \cF{a,b)\ —
for any ao > 0.
The additional hypothesis F G B^ l (IR 2 , H) for a > is equivalent to the following condition
on the wavelet coefficients (see for e.g. ]31^):
f°° da
(18) IIFIId-o- = / sup I cf (a, b)\—z — < +oo .
°°> l Jo bm 2 a a
In such case, we easily deduce that for any ao > 0:
f°° da f°° 1 da
(19) sup/ \c F (a,b)\^ = sup / \c F (a, b)\ — ^— < a a \\F\\^-a
bm 2 Ja a bm 2 Ja a* a 2 ° a oo,i
Moreover, by hypothesis, F G C 1 (M 2 ,EI) and VF is bounded. Using the expression of the
wavelet coefficients then implies that:
\c F (a, b)\ < a 2 ||VF|| LO c [|a^|| £ i , V6 G IR 2
which ensures that for any oq >
f a ° da
(20) sup / \c F (a, b)\— < oo .
beR 2 Jo a
and finally leads to the existence of the integral f + °° cp(a, b)^.
8
Remark 3.2 In ]29^ authors considered spectral lines of the form f{x) = A(x) cos((p(x)) with
A,Vip slowly varying with respect to (p. The authors then approximate the monogenic signal
F associated to f by an IMMF of the form 1113]) with constant amplitude, linear phase and
constant local orientation.
We now focus on the estimation of the amplitude and instantaneous frequency of an IMMF.
To this end, we need some assumptions on the wavelet that we list below:
Assumptions (W)
(i) the wavelet ip is isotropic and real-valued.
(ii) V e W 1 ' 1 (R Z ).
(hi) The hrst three moments of tp and are hnite, namely
sup I a < oo, sup l' a < oo ,
aG{l,2,3} ae{l,2,3}
where for any a >
(21) I a = [ |x| Q |V>(x)|dx , I' a = f \x\ a \Vip(x)\dx .
JR2 JR2
We first deal with the simple case of an IMMF of constant amplitude and phase:
Lemma 3.1 Let F be an IMMF of the form F(x) = ^ e (fc ' x)ns with k = (ki,k 2 ) G R 2 ,
Aq G R+,# G R, and let ip be an isotropic real wavelet belonging to iy 1 ' 1 (R 2 ).
The wavelet transform of F is given by:
(22) c F (a,b) = A a^(ak)(cos(k-b) + sm(k-b)(cos9 i + sin#j)) = a^(ak) (A Q e {k - b)ne
and one has for i = 1,2:
(23) d h c F (a, b) = hn e (af (a*;)) (A e^ n ^ .
The vectors k and ng can then be recovered using the two following relations:
kin e = d bl c F (a,b) x (c^a,^) -1 ,
k2ng = d b2 CF(a,b) x (cj?(a,6)) _1 .
Proof. By definition, since the wavelet is isotropic real, one has:
c F (a, b) = a -1 / F(x)ip ( ] dx = a f F(au + b)ip(u)du
A simple calculation leads to:
(J R2 cos(k ■ (au + b))ij}(u)du \ / cos(k.b)
cos 9 f R2 sin(fc • (au + b))ij}(u)du = AQa^(ak) sin(k.b) cos(9) = aip(ak)F(b) .
sin9 L a sin(k ■ (au + b))ip(u)du J \sm(k.b) sin(0) /
9
where we have used ip(-ak) = ip{ak) since the wavelet is isotropic.
Now, since tp G W 1 ' 1 ^ 2 ) and F G L°°(IR 2 ), we get:
d h c F {a,b) = a' 1 [ F(x)d bi L (—)] dx = -a~ 2 [ F(x)(d Xi ip) (—) dx .
Jr-2 I V a )\ J R2 \ a J
A similar calculation, using d Xi ip(£) = i£i*0(£),V£ = (£1,62)) and tp isotropic, yields:
(— sin(k.b) \
cos(fc.6) cos(6>) = akitp(ak) ngF(b) ,
cos(fc.fr) sin(0) /
where we have used:
ngF(b) = ng(cos(k ■ b) + ng sin(k ■ b)) = — sin(/c ■ b) + rig cos(k ■ b).
This leads to □
Remark that, more generally, if F(x) = J 4 e( fc ' a;+a ) ne , (a £ 1):
(24) c F (a,b) = A a^{ak)e^ b+a ^ = a$(ak) F(b)
(25) 3^0^(0,6) = akiip(ak) ngF(b)
The formulas (|22p and (|23[) can be extended in the general case of an IMMF F of the form f)13[) .
Proposition 3.1 iei -0 a wavelet satisfying assumptions (W) and F an IMMF with accuracy
e and smoothness a > of the form Iil3\) . For all (a, b) G x M 2 , one has
(i)
c F (a, 5) = a$(aVip(b)) (A(b)e^ b)n(b) ) +ea 2 R 1 {a,b) ,
with
(26)
|#iM)l < Ii(V^A(6) + l)|V^(6)|+aJ 2 A(6)(|V^(6)|+V2M)/2
+aI 2 M/2 + a 2 / 3 ^(6)M/6 .
^c F (a,6) = d b Mb)n e(b) (a#(aVp(&))) (yl(6)e^^w) +eaii 2 (a,6) ,
with
\R 2 (a,b)\ < A(b){a\V^(b)\r 2 /2 + a 2 MrjQ)
{ ' + (>/2A(6) + 1) (|V^(6)|I( + aM/ 2 /2) .
Remark 3.3 As a consequence of this proposition, since in the isotropic case c F {a, b) and
Cj M ^ (a, a, b) are related by All]) , the above proposition specifies the results stated in \20$
and ]29$ . Notably, it gives a bound of the approximation error.
10
Proof. See Section EU □
We now define the two following Clifford vectors:
(28) Ai(o,6) = d bl c F (a,b) x (cfM)) -1 ,
A 2 (a,6) = d b2 CF(a,b) x (cF(a,b)) L .
These two vectors will provide an approximation of the vectors c\.(p(b)nQn,\ for an IMMF F
of the form (|13p . More precisely, we can state the following result:
Theorem 3.2 Let F be an IMMF with accuracy e > and smoothness a > of the
form [T3\) . Assume that we are given tp a wavelet satisfying assumptions (W) and consider
(a, b) G W* + x M 2 such that
\cF(a, b)\ >e u for some v G (0, 1/2) .
Then for i = 1, 2:
(29) \Ai(a,b) - d h <p(b)n e{b) \ < e 1 ^ a (\R 2 (a,b)\+ a|fli(o,6)| |Vp(6)|)
where R\, R2 have been defined respectively in [26]) and (27\ ). In addition for any oq > 0,
there exists some Eq > such that, for all (a,b) G (0, ao) x M 2 and any < e < e^:
(30) |Ai(a,6)-9 6 ^(6H (6) | <e v
Proof. See Section E2J □
Since cos(tp(x)) = cos(— <p(x)), there is an ambiguity in the definition of the instantaneous
frequency since ip can be replaced by — (p. To solve this problem, practitioners usually
assume that the first component of the instantaneous frequency is positive (see [10] and
references therein for more details). Under this additional assumption, our method yields an
estimate of the amplitude and of the instantaneous frequency of an IMMF:
Proposition 3.3 The notations and assumptions are those of Theorem \3.2\ In addition, we
assume that:
Vx G M 2 , d xi ip(x) > .
Then for i = 1,2, and (a, b) G (0, ao) x M:
(31) |^^(b) - |A x (a,6)|| < e" .
and
(32) \d b Mb) ~ |A 2 (o,6)| sgn(Re(9 fel c F (a, b) d b2 c F (a,b)))\ < e v .
Proof. See Section E3J □
11
4 Decomposition of a mult i— component image into spectral
lines
In this section we will consider a more general context, where our model is now a superpo-
sition of several IMMFs (introduced in Section [3]), assumed to be well separated, as defined
in the following definition 14.11 We will show in Theorem 14.11 how such functions can be de-
composed into monogenic modes, using synchrosqueezing. Thereafter, for each component,
the associated amplitude and instantaneous frequency will be estimated using the results of
Section [3j
Definition 4.1 A function F defined from R 2 to H is said to be a superposition of well
separated Intrinsic Monogenic Mode Components with smoothness a > up to accuracy
e > and with separation d > 0, if there exists a finite integer L such that
L
(33) F(x) = ^F e (x),
where all the Fi are IMMFs with accuracy e and smoothness o~i > a > of the form i!3\):
Fg(x) = Ag(x) Q Pe ^ n9 t( x ) ; and moreover satisfy for any x, £ > £' and i G {1, 2};
(34) \d Xi <pt(x)\ > \d x .tp£>(x)\ ,
and
(35) \d Xi ipi(x)n edx) - d Xi ip e >(x)n 9i/ ( x) \ > d [\d x .<pt(x)\ + \d Xi ipt> (x)\] .
We denote by A £ a d the class of all functions F satisfying these conditions.
We now define the monogenic synchrosqueezing transform (MS ST) of a function F belonging
to A Eja> d-
Definition 4.2 Assume that we are given ip a wavelet satisfying assumptions (W) such that
supp(V0 C {£ G R 2 ; 1 - A < |f | < 1 + A} ,
with
(36) A < d
2{\ + d)
Consider h G C£°(R,R) such that f R h(x)dx = 1. Let 5 > 0, v G (0,1), and e > 0. For
F G A e>a ,d, the monogenic synchrosqueezing transform of F is defined for any (6, k, n) G
l 2 xl 2 x S 1 fS 1 being the unit sphere of R 2 ) by
(37)
6u f . 1 /fci-Re(n Ai(a,6))\ fk 2 -Re(n A 2 (a,6))\ da
S F ^(b, k, n) = c F (a, b) T2 h { j h { j ,
with
A £)F (b) = {aeR + ; \c F {a,b)\ > e u ] ,
and where A\(a,b), ^(a, b) are defined by equation A28\) .
12
We can now state our main result:
Theorem 4.1 The notations are those of Definition ^. ^ Let F G A £t<Ti d an d v £ (0)1/(2 +
4/cr)). Provided that e > is sufficiently small, the following assertions hold:
• |cp(a, 6)| > e v only if there exists some I such that
(38) (a, b) G Z e = {(a, b), |a|V W (6)| - 1| < A} .
• For eac/i i = {1, • • • , L} and eac/i pair (a, 6) G 2^ /or which |cp(a, 6)| > e", one /ias /or
i = 1,2 and for some C > 0:
|Ai(a,6) - dnpi(b)n et ( b ) \ < Ce v .
• Moreover for any £ G {1, • • • , L}, there exists C > 0, suc/i £/iaf /or any b G M 2
lim
5^0
1
27rC^ J§i JB e {e»,n,b)
S 5 f E {b, k, n)dkdn - A e (b)e w(b)n0 iW
< Ce u .
where we denote Be(e u ,n,b) = {k G R ; maxj |/cjn — 9 bi (pi(b)n0 e ( b \ | < e^}, and where
C^p was introduced in ([6]).
Proof. See Section E3
□
5 Implementation and numerical experiments
This section aims to illustrate the efficiency of the monogenic synchrosqueezing transform,
and its potential applications. We will first describe the computation of the discrete MSST
from a N x N discrete image, then we will show some examples on artificial and real images.
5.1 The discrete monogenic synchrosqueezing transform
The first step consists in the discretization of the variables, therefore one builds discrete
sets A, JC, O for the scales a, the normalized frequencies (&i,&2), and the orientations 9
respectively, the frequency resolution in k being denoted by Note that one usually chooses
a logarithmic scale for a and k, e.g. A = {2^ nv }j = o.-n a -i, so that we have n v coefficients
per octave. Then, we compute the discrete monogenic wavelet transform cf(cl, 6), and the
estimations of the instantaneous frequencies Aj.(a, b) and A 2 (a, b) defined in equation (f2"8j) by
writing for i = 1,2:
d bt c f (a,b) = f(x)d bi tp a , b (x) dx,
the other components d bi CTi 1 f(a,b) and d bi c-]z 2 f{a,b) being computed in the same manner.
The monogenic synchrosqueezing transform consists in a partial reallocation of the mono-
genic wavelet transform according to space, frequency and orientation parameters. Given a
threshold 7, and for all a £ A, (k\, k 2 ) G K? ,6 G O, one simply writes
(39) S F ^{b,ki,k 2 ,n d ) = > .
n„ ^— ' a
(\c F (a,b)\ > 7
\h - Re(Ai(a, b) n e )\ <
\k 2 -Re(A 2 (a,b)n e )\ < ^
13
The term 1( ^ 2 ' comes from the measure ^§ in equation (|37p . as we deal with logarithmic
scales. Note that if we use a logarithmic scale for the frequencies like in [21], the resolution
Afc depends on the value k.
This discrete MSST is relatively easy and cheap to compute, but it is not of great interest
in practice because of the large number of variables which prevents both easy representation
and fast processing (i.e. decomposition, demodulation). To decrease this number of variables,
one can remark that if a superposition of IMMFs satisfies the separation condition of equation
(fM|) . then for any I > £' it also satisfies |V<p^(x)| > |V(/?^(x)|, so that the IMMFs are also
separated in terms of the norm of the frequency vector. This leads us to define the following
discrete isotropic MSST for all a G A and k G fC:
(40) S F>1 (b,k)- l ° s{2) ^ CfM)
a£A S.t.
— 2
V|Ai(a,6)| 2 + |A 2 (a,6)|^
Intuitively, Sf^Q), k p ) contains all the coefficients Ci?(a,6) whose "instantaneous isotropic
frequency" y / \A\(a, b)\ 2 + | A2 (a, 6) | 2 is around the value k. The latter isotropic MSST can
be particularly useful when visualizing the time-frequency distribution, as it only depends
on three scalar parameters (bi,b2,k), that is by far easier to represent. Hence, this version
will be used in the experiments hereafter. Additionally we will use the Morlet wavelet with
parameters \x > and a > 0, defined in the Fourier domain as follows
(41) ^ i(T (0 = exp(-7rV(|e|-^ 2 ).
5.2 Representing a synthetic 3-components signal
This section illustrates the interest in using the MSST to represent multicomponent signals
in 2 dimensions. Let us first define the synthetic 3-components test-signal / = f± + f 2 + f'3
with
( frfay) = e- 10 « a; - a5 ) 2+ ^- ' 5 )- 2 ))sin(10^(x 2 + y 2 + 2(x + 0.2y))
(42) { f 2 (x,y) = 1.2sin(40vr(x + y))
[ / 3 ( x , y ) = cos(2vr(70x + 20x 2 + 50y - 20y 2 - 41xy))
We compute both its monogenic isotropic wavelet transform cp and its isotropic MSST Sf-
We then aim at representing |cf| (resp. |»Sf|), which depends on the three scalar parameters
x, y and a (resp. k). We will here propose two distinct visualizations: either a 3-dimensional
scatter "density" plot, or a 2-dimensional representation for x or y fixed, which looks similar
to a ID continuous wavelet representation. The following Figure Q] shows these different
representations for the wavelet and the MSST of the signal (|42p . We see how the MSST
sharpens the time-scale representation, making it more readable.
14
5.3 Decomposition and Demodulation
Once the MSST has been computed, the key point to separate and reconstruct the modes is
to identify the instantaneous (anisotropic) frequency of each mode \Vipi(b)\ at each point b.
The t th IMMF Fi> is then approximated by F e (b) as follows:
(43) F t (b)= S Fn {b,k).
ke(b)—K<k<k((b)+K
where ke(b) is an appropriate estimate of |V(/?£(6)|. Obtaining this estimate is an involved
problem which will not be discussed here, so we simply use the approach used in [21] based on
a greedy algorithm. Our approximation Ff(b) of the I th IMMF is computed by summing the
coefficients of the synchrosqueezed transform in the vicinity of this instantaneous frequency,
to regularize the solution. In practice the width k of the window is set to To evaluate
the method, we compare the mode F% of the test-signal defined in equation (|42|) to the
corresponding mode Fj> obtained from the discrete MSST after extraction, by computing the
normalized Mean-Squared Error (MSE):
(44) MSE(F e ) = " h ~ f '"
\\m
Each reconstructed IMMF is displayed on Figured where one also provides the corresponding
MSE. It is clear that the modes are reconstructed with very high accuracy, although our
method for ridge extraction is quite simple.
Finally, we aim here at showing that the MSST remains efficient on real images or tex-
tures, even though they can not be considered as multicomponent signals. The following
example shows the image Lenna where we artificially added an oscillating pattern. After a
2D synchrosqueezing transform one is able to extract the oscillating component to roughly
reconstruct the original image. Note that this kind of image was first introduced in |27| to
illustrate the discrete monogenic wavelet transform.
15
(a) (b) (c)
Figure 2: The extracted and reconstructed modes 1, 2 and 3 for the test signal of equation
(|42p . The normalized MSE are respectively 0.03, 0.06 and 0.05. We used the Morlet Wavelet
ifi^o- defined in (|4"Tj) with a = 2 and \i = 1, and removed the border (| at each side for each
dimension) to avoid border effects.
(a) (b) (c)
Figure 3: Illustration of possible applications on real images, (a): Input image, Lenna plus
an oscillating pattern, (b): the extracted mode by the 2D MSST. The parameters are the
same as for Figure [2j The MSE is 0.12. (c): The residual approximates the original Lenna
image.
6 Proofs
6.1 Proof of Proposition 13.11
We prove Proposition 13.11 into several steps.
6.1.1 A preliminary result
We first need a preliminary result:
Lemma 6.1 Assume that F is an IMMF with accuracy e of the form t!3\) . The following
estimates hold for any y,h G M 2 :
(45) \A(y + h)- A(y)\ < e (\h\\V<p{y)\ + \h\ 2 ^
(46) \ cos e{y + h)- C o S e{y)\<e[\h\\V^y)\+\h\ 2 ^j ,
16
(47) |sin% + /0-sin%)| < e (\h\\V<p(y)\ + H'y)
(48) \V<p(y + h)- V V {y)\ < e (\h\\V V {y)\ + \h\ 2 ^j
where M is the constant defined in fTl
Proof. We only prove inequality ()45|) . the others follows in the same way. Observe that:
A(y + h)-A(y)= [ V A(y + th) ■ h dt
Jo
By assumption (|16p . VA is slowly varying with respect to then:
\A(y + h)- A(y)\ < e\h\ I \V<p{y + th)\dt
J
Applying the usual mean value theorem, we deduce that for any t:
\V<p(y + th)\ < \V<p(y)\ + \V<p(y + th)-V<p(y)\
< \V<p(y)\ + t\h\ max ( sup \d 2 ip{y)\ )
< |V^(i/)|+Mt|fe|,
applying Assumption (fl5j) on second order partial derivatives of (p. To sum up, one has
\A(y + h)- A(y)\ < e\h\ J \V<p(y)\dt + eM\h\ 2 J tdt < e (\h\\V<p(y)\ + M^-\ ,
which is the required result. □
6.1.2 Proof of Proposition I3TT1
We now prove Proposition 13.11
Proof of point ([!]) of Proposition 13.11
By definition of the wavelet coefficients of F, one has
c F (a,b)= [ A(x)e^ ne Ma- 1 ip(- — dx
Jm.2 V a J
We now split this sum into three terms
c F {a,b) = a- l A(b)j R2 B^ n ^i>(^)dx
(49) +a~ 1 A{h) J R2 \ e ^ n <>M - ^ x ) n m j ^ [*=±) dx
+CL- 1 J R2 [A(x) - A(b)] e^W0 (2^) dx .
Observe now that
(50) (p(x) = (p(b) + V<p(b) -(x-b)+ [ [V<p(b + t(x - b)) - Vtp(b)] ■ (x - b)dt .
17
We set u = and use equation with k = V(p(b). We then deduce that:
(51) a' 1 I e Mb)+Vrtb)ix-b))n e(b) ^ fEz]L) dx = e <p(b)n e(b) a ^( aV ^(&)) .
Jm. 2 V a )
Combining equations ([75|) . ([50]) and (f5"Tj) implies:
c F (a,6)-a^(aV^(6)) e ^ n °m
= A(b) f e Mb)+V<p(b)-(x-b))n g{b) f e n g(b) fi[Vip(b+t(x-b))-Vip(b)]-(x-b)dt _ A Q -l^
x — b
+A(b)
^tp(x)n g(x) _ <p(x)ng( b)
a~ x ip
x — b
dx
+ / [A(x) - A(6)]e v(:c)n ^)a~V
x — b
dx .
Hence
c F (a,b) - a$(aVip{b)) \A(b) e^ (6 H(&)
< a _1 yl(6) / e [/ 1 [V V (6+t( a; -fe))-V^(6)]-( a; -fe)dt]n e(i , ) _ 1
Jr 2
V'
x — b
a
dx
+a- 1 / \A(x)-A(b)\
x — b
dx .
dx
We now give an upper bound of each of these three terms. We bound the first one using
mean value theorem and the triangular inequality:
a^A(b)
^([V<p(b+t(x-b))-V<p(b)Hx-b))dt]ne (b)
1
x — b
\(y<p(b + t(x - &)) - Vp(&)) • (x - b)\ dt
We now use inequality (jl8"|) with y = b and h = t(x — b). We then obtain:
a -1 A(6)
< ea _1 A(&)
|(V^(6 + t(x - b)) - V(p(b)) -(x-b)\dt
)dt
■4>
\\t\\x - b\ 2 \V<p(b)\ + M\t\ ^ x n 6|S
dx
x — b
a
x — b
a
x — b
a
dx .
dx
dx
x - b\ 2 ,„ „ Jx - b\ 3
6
x — b
dx
where we set u = s — in the last integral.
Let us now give an upper bound of the second term. Since
^(x)n e(x) _ ^(x)n e(b) = sin{lp{x)){ne(x) _
18
one has:
jp{x)n e(x) _ ip(x)ng( b )
< a 1 A(b) / \n e{x) - n e[b) \
x — b
x — b
dx
dx
1 A{b)V2 [\x-b\\V^{b)\ + \x-b\ 2 —\
< ea 2 A(b)V2 [ (\u\\V<p(b)\ + a\u\ 2 — J Mu)\du
JlR2 V 2 /
< eA(b)V2 (a 2 \Vy{b)\h + a 3 ^I 2
x — b
dx
(53)
where we have used inequalities (|46p and (|47p with y = b and h = x — b for the third inequality
and set u = (x — b)/a for the fourth one.
Finally, using the same approach we prove that the last term is bounded by, using (I45p :
(54)
\A(x)-A(b)\
x — b
M
dx < e ( a 2 \V(f(b)\h + a 3 —I 2
Combining inequalities ([52l |53"1 [54^) leads to:
Ka,6) - a4>(aV<p{b))A(b)e^ n ew | < ea 2 A(6) (Y| |Vp(6)| J 2 + yMJ 3 ) + v 72 f|V^(6)|Ji +
+ ea 2 ^|V^(6)|Ji+ayJ 2 ^
which gives (|26|) .
x — b
dx .
Proof of point ([TT]) of Proposition 13.11
Since A G L°°(M 2 ) and ip G W 1 ' 1 ^ 2 ), one has:
a bi c F (a,6) = - / ^(xle^W^Wa- 2 ^^
As below we split the sum into three terms
d bt c F (a,b) = - a - 2 A(b)J R2 e^ n ^d x ^{^)dx
(55) -a~ 2 A{b) f R2 (e^GOM*) - e*w n «( 6 )) <9 Xi V> (^) dx
-a~ 2 f R2 [A(x) - A{b)} ert'fr'Wd^ dx
We again use (f50|) and equation ([24"]) with fc = V(p(b) to write:
(56) - a" 2 y e Mb)+v V (b)ix-b))n e(b) f^-^\ dx = e ^ ne ^ d bi <p{b)n m a^(aV<^(6))
19
Hence we deduce that
d bi c F (^b) - d bi cp(b)n e(b) (a${aV<p(b)j) (A(b)e*® n m
< a- 2 A(b)
+ a - 2 A(b)V2 [
(fo[Vv(b+t(x-b))~\7<p(b)]-(x-b)dt)n e(b) _ 1
x — b
x — b
dx
\ n e(x) - n ${b)\
dx
-2
\A{x)-A(b)\
x — b
dx
A similar approach than in the proof of point (JTJ) leads to the following upper bound:
d hi c F (^b) ~ d b .Mb)n e{b) (a^aV^ft))) (A(b)^ n ^)
< eA{h) (^\ Vlp{b) \^ + ^ M A
+eA(b)V2 (a|Vv9(6)|/;+a 2 ^/^
+e (a\Vtp{b)\l' 1 + a 2 ^-I^
which leads to the estimate (Jn]) of Proposition 13.14 with:
\R 2 (a, 6)| < A(b) C| |V^(6)| l' 2 + ^mA + (flA{b) + l) (\Vip(b)\l[ + a™l' 2
6.2 Proof of Theorem
We can now prove Theorem 13.21 By definition for i = 1, 2:
Ai(a,b) = (d bi c F (a,b)) (cFia^b))' 1 .
Set B = a^(aVip(b)) (A(b)e v ^ n ^. Proposition O implies that:
(57) \c F (a,b) -B\ < ea 2 \Ri(a,b)\ ,
and
(58) \d bi c F (a, b) - d bi <p(b)ng(b)B\ < ea \R 2 (a, b)\ ,
where Ri(a,b) and R 2 (a,b) satisfy inequalities ([26]) and ([2?]) respectively. Then, for i = 1,2:
ki(a,b)-d bi ip{b)n e[b) = [d bi c F (a, b) - d bi ip{b)n eib) B + d bi ip{b)n e(b) (B - c F (a, b))] [cp^b)}^ 1
Using ([57|) and ([58]) implies that:
\ki(a,b) - d bi tp(b)n e{b) \ < e[a\R 2 (a,b)\ + a 2 \Ri(a,b)\\d bi tp{b)\] |cF(a,6)| _1 .
20
By assumption \c F {a, b)\ > e v for all (a, b), which leads to (|29p . The last estimate (|30p follows
from the fact that for any given ao >
sup (o|# 2 (a,6)| +o 2 |i?i(a,6)||V^(6)|) < oo .
(a,b)e(0,oo)xK
Then, since f G (0, 1/2), there exists some £o > depending on ao such that for any < e < £o
a\R 2 (a,b)\ + a 2 \Ri(a,b)\\V<p(b)\ < e 2v ~ x ,
which is equivalent to
e x - v (a|i? 2 (a,6)| +a 2 |i? 1 (a,6)||V V 9(6)|) < e u .
This last inequality and inequality ([29]) clearly implies inequality (f30|) .
6.3 Proof of Proposition 13.31
The proof of Proposition 13.31 relies on Theorem 13.21 and on the following lemma:
Lemma 6.2 Let F an IMMF with accuracy e > 0. Provided that e > is sufficiently small,
the sign of Re(db 1 c F (a, b)db 2 c F (a, b)) is this of db 1 ip(b)db 2 (p(b) .
Proof. Let us first observe that Theorem 13.21 implies that for all (a, b) under consideration:
Ax (a, b) = d b M b ) n e(b) +0{e u )
A 2 (a,b) = d b2 ip{b)n e(b) +0(e v )
Then
(59) Ax (a, 6)A 2 (a, b) = ^(0)^(6) + 0(e"
Hence, for e > sufficiently small, the sign of Ai(a, 6)A 2 (a, b) is this of db 1 (p(b)db 2 <p(b).
To get the required conclusion, we now relate Ai(a, 6)A 2 (a, b) and db l CF(a,b)db 2 CF(a,b).
Let us first remark that the definition of the two Clifford vectors Ai(a, b), A 2 (a, b) from equa-
tion ([28]) implies that :
(60) Re ( Ai(a, 6)A 2 (a, b) ) = Re (db l c F (a, b) (c F (a, b)) l db 2 c F (a,b) (c F (a,b)) 1
We use now the fact that db i c F (a,b) and c F (a,b) are both Clifford vectors. We then apply
equality ([83]) with q = db 2 c F (a,b) and q' = (c F (a, Hence we get that
d b2 c F (a,b) (c F (a,b)) 1 = (cf(o, &)) 1 x d b2 c F (a,b) .
Combining this last equality with (|60p leads to
(61) Re ( Ai(a, 6)A 2 (a, 6) j = Re (c^c^a, 6) db 2 c F (a, b) ) |cf(o, 6)|~ 2
Taking into account (|6T|) , equation [59] reads
Re ( 9^0^(0, 6) db 2 c F (a,b)
21
which implies that Re ydb 1 cp(a, b) db 2 CF{a,b)j and db 1 ip(b)db 2 tp(b) have the same sign for
e > sufficiently small. □
We can now prove Proposition 13.31
Proof of Proposition 13.31
Theorem 13.21 directly implies that for any i = 1,2 and (a, b) under consideration:
||Ai(a,6)|- |^(6)|| <£ v .
Since 9^(^(6) is assumed to be always positive, we get ([3T]) . Since db 2 (p(b) and db 1 ^p(b)db 2 ^p(b)
have the same sign, Lemma 16.21 implies (|32p .
6.4 Proof of Theorem 14.11
The proof of Theorem 14.11 consists in several steps.
Lemma 6.3 Let F be a function belonging to A e ,d- For any (a,b), there can be at most one
£ 6 {1, • • • , L] such that
(62) |a|V<^0)| - 1| < A .
Remark 6.1 Condition Ii6ty) is a necessary condition for having ?Jj(aV(pi(b)) ^ 0. Lemma \6.'3\
then states there is at most one £ such that ip(aV(pe{b)) ^ 0. Then the two following sums,
involved in the estimation o/c^(a,6) and d^CF^a, b):
L
J^Mb) e w(J) "W a$(aVip e (b)) ,
t=\
and
L
J2Mb) e^ n °cW a^{aV^{b))d m {b)n ei{b) ,
e=i
contain at most one term.
Proof. We follow the same line as in |21| . Assume that there exists some (^1,^2) such
that ([62]) is satisfied with £\ ^ £2- One may suppose that l\ < £2- One then has for j = 1,2,
a~ 2 (l-A) 2 <
which directly implies
9<Pli 0)
+
dx 2
< a- 2 (l + A) 2
+
dxo
&Ph (b)
dx\
+
9^2 (b)
dxo
> 2a~ 2 (l - Af
and
&Ph (b)
dx\
+
8x2
9ft! (b)
dx\
Agfa (ft)
8X2
< 4a" 2 A .
22
Further, by assumption (|34|) . for any i = 1,2 one has
dip e2 (b)
dxi
d<Pti 0)
dxi
>
d<pt a (b)
dxi
dtp i2 -i(b)
dxi
> d
d(p h (b)
dxi
+
d<ft 2 -i(b)
dxi
where we applied assumption ([35]) :
1 1 0*0 1 - 1^^2-1(^)11 > \dx^e 2 (x)n ee{x) - d Xi (pi 2 -i(x)ne e ,( x )\ > d [\d Xi <pe 2 (x)\ + \d Xi (pe 2
The usual inequality (u + v) 2 > u 2 + v 2 , valid for any u, v > implies that
d(fe 2 (b)
dxi
dxi
> d
d(p h (b)
dx;
dx.
Summing over i = 1, 2 gets:
4a -2 A > 2da _2 (l - A) 2 > 2da~ 2 (l - 2A)
that is
A > dj [2(1 + <£)]
which leads to a contradiction.
□
Lemma 6.4 Let F be a superposition of IMMF with accuracy e and separation d of the
form For all (a, b) £R*x M 2 one has
c F (a,b) -aJ2 ${aV<pi(b)) (A e {b) e ^ (6)n W
i'=i
< ea 2 ri(a,6)
with
(63)
ri(o, b) = h Eii(V2^(6) + 1) |V^(6)| +af [Eti ^(&) (|V W (6)| + Af*)
(T a freing introduced in 1121]) ). In particular
• If for some ^ G {l,--- , L]
(64) \a\V(pe (b)\ - 1| < A .
f/ien
(65)
(a, 6) - a^(aV^ (6))^ (6) e ^ (b)n % (6) < ea^M)
If for all £ 6 {1, • • • , L}, Condition fl&^[ ) is not satisfied then
(66) |c F (a,6)| < ea 2 ri(a,6) .
23
Proof. By assumption F is of the form
L
e <Pe n i
=1
where for each £, Fj> = A(,e tptni is an IMMF of accuracy e. Proposition 13.11 applied successively
to F\ , ■ ■ ■ ,Fl and the linearity of the continuous wavelet transform imply:
c F (a,b) - a^2i>(aV<p£(b))Ai(b) e
<Pt(b)n g (6 )
< ea 2 ri(a,6) .
The rest of the proof follows from Remark 16.11 which states that under condition (|64p , the
following sum
L
Y^Mb) e^ (6 Kw atj>(aV(pt(b)) ,
reduces to A £o (b) e w ° (fc)n V b) a$(aV(p £o (b)), and if condition (JMI is never satisfied. □
Lemma 6.5 Let F be a superposition of IMMF with accuracy £ and separation d. For all
a, b such that
(67) |a|V^(6)| — 1| < A ,
one has, for any i = 1,2:
(68) d bt c F (a, b) - ad m (b)n ei{b) ^{aV W (b)) (A e (b)e n{b)n ^) < eaT 2 (a, b) ,
with
(69)
T 2 (a, 6) = l'x Eti I V^(6)|( V2^(6) + 1) + a4f ^ti {M t + A e (b)(V2M e + |V^(6)|))
(I' a being introduced in IF?])).
Proof. As in the proof of Lemma 16.4] F is of the form F = Y2i=i A.t^ ptnt > where each
.F) = Aie Vtni is an IMMF of accuracy e. Proposition 13.11 applied successively to F\, ■ ■ ■ ,Fl
and the linearity of the continuous wavelet transform imply:
L
d bi c F (a,b) - a]T dm'(b)n eAb) ii(aVMb)) (if(4)e"'' (^) > (i, ) < eaT 2 {a,b) .
i'=i
Remark 16. II states that under condition (I64p . the following sum
L
^ $M*K,(6)#*V^(&)) (^(6)e w(6)n V<
£'=i
reduces to d m {b)n dl(b) i>{aV W {b)) (^(&)e<«< 6 Kw) . □
24
Lemma 6.6 Let F be a superposition of IMMF with accuracy £ > and separation d. Let
I G {1, • • • , L} . For all a,b such that
(70) \a\V<p e (b)\-l\<A,
one has, for any i = 1,2
(71) \Ha,b) - d bi w{b)n 9e(b) \ < ae 1 -" {r 2 {a,b) + a|V^ (6)|ri(a,6)) .
Proof. The proof follows the same lines than that of Theorem 13. 2^ replacing the two esti-
mates (|57p and (|58p on the wavelet coefficients and their partial derivatives, by (|65p and (|68p
respectively. □
We can now prove Theorem 14.11
Proof of Theorem 14.11
In the proof of Theorem 14.1} the following notations will be needed:
1 - A „ _ 1 + A
^min "
max te{1 ... ;L} sup b€K 2(|V^(6)|)' max min te{lj ... jL} inf 6eK2 (|V^(6)|) '
By assumptions on the phases (pg of the modes of F, a m \ n and a max are both finite and positive.
For any t G {1, • • • ,L}, e > 0, b G R 2 and n G S 1 , one also defines:
Bi(e, n, b) = {k G M 2 ; max - d b .(p e (b)ne e (p)\ < e} .
We now show Theorem 14.11 We set e = e v and fix t G {1, • • • , L}. Let us first prove that
= 2 ^ fA £ , F (b)n{a-,m aXt \A l (a,b)-d bz Mb)no e(b) \<^ a ~ C H a > & ) da •
By definition of Sp", we have
Sp U £ (b, k, n)dkdn
S 1 JB e (e,n,b)
JS 1 JB e (e,n,b) \JA £tF (b) fj[ V * / « /
We first use Fubini Theorem to interchange the order of summations in this integral. Remark
that:
Js 1 JBe{e,n,b) J A E<F {b) fj[ \ ° J
< [ a-\ F {aM I [ r 2 ff h ( *j-MU«,b) n) \
JA EF (b) JS 1 JR 2 ~_i V ° J
^-dfcdn
a z
'A s , F (b)
le,F(6)
/ a 2 |cf(o, b)\da < + oo ,
J A, F tb)
25
where the last inequality is coming from the fact that F is a finite superposition of IMMFs
Fg, whose wavelet coefficients satisfy (|19|) and (|20p (see Remark 13.11 for more details). Then
by Fubini Theorem:
F,e
S 1 JB e (e,n,b)
Sp l ' (b, k, n)dkdn
[ a-*c F {a,b)( [ [ 6-ZT\h( ki - Rei Y a,b)n) )dkdn)da.
2
A : ., i.V) V-/:' 1 J6 ; r. „.h) "
Since the integrand is bounded by:
a~ 2 c F (a,b) H jT (?wfc) ^ 2 n^( ^" Re( t (a ' b) dMn ) ^ 2vr||/ l ||i 1(R) a- 2 |cHa,6)|
which belongs to L l (A £jF (b)), the Dominant Convergence Theorem applies. Hence,
lim / / S^(6, /c, n)d£;dn
8->0J§i JB e (e,n,b)
= [ a"V(a,6) (lim / / 6^ f\ h ( ^^^) dkdn) da
JA e , F (b) \^J § i J Bt (e,n,b) f = i V d J )
2
a- 2 c F (a, b) ( lim / / «T 2 TT J» ( ~ w) ) dfcdn ] da
A £ , F (6) y^Ojgi JB e (£,n,b)nB'(SR,n,b) fj[ V ° J J
where we set B'(5R,n,b) = {k £ M. 2 ; maxj |/cjn — Aj(a, 6)| < with i? the bound of the
support of h. Remark now that, if maxj |Aj(a,6) — db i <pi(b)riQ e ^\ > e, then for 5 sufficiently
small
B e (e,n,b)nB'(SR,n,b) = ,
and the limit above is equal to 0. On the other side, if maxj |Aj(a, b) — d^cpiQijng^l < e, for
5 sufficiently small:
B e (e, n, b) n B'(SR, n, b) = B'(SR, n, b) ,
and
lim / / s - 2 fr h (h Z MHa 1 b) n)\ ^
J§i JB'(5R,n,b) fJi \ 6 /
2
= 2vr,
where the last equality comes from the assumption L /i = 1. One then deduces equality (f72|)
Observe now that for any e > sufficiently small there exists some ao > a m ax such that
(73) 2vr a Q a < e v and sup e x ~ v a (T 2 (a, b) + a|V^(6)|ri(a, b)) < e v .
00,1 (a,fe)e(0,a )xR 2
26
Indeed, by definition of Ti(a,b),T2(a,b) given respectively in equations ([65]) and ([69]) . there
exists Cr > not depending on ao such that
sup a (r 2 (a, b) + a\Vipt{b)\T\(a, b)) < Crag .
(a,6)e(0,a )x]R 2
This inequality implies that to prove (J73J) , it is sufficient to find ao such that
(74) a max < (2vr \\F\\ B -„ e-») 1/a < a < (C^e 2 ^ 1 ) 1 ^ .
Since v < 1/(2 + 4/<r), one deduces that g^ 2 " -1 )/ 4 /e~ u ' a tends to oo when e — > 0. Hence, for
e sufficiently small, there exists some ao > a max satisfying (|74p . From now to the end of the
proof, we fix such ao-
We now split the integral in the right-hand side of (|72p into two terms
lim^o / s i I Be (e,n,b) k > n )dkdn
( 75 ) = 27F X4 e , J ,(6)n{ae(0,ao);max i |A 4 (a,6)-ft iW (6)n tfi( 6 ) |<6} Ci? ( a ' W
+ 2w JA e , F (b)n{a>a ;ma Xi |A i (a,6)-a i , < ^(6)n^ (b) |<g} CF ( a ' ? ')# "
We now deal with each of these two terms successively. Recall that since by assumption F is
a superposition of IMMFs, then inequality (|19p holds for some a > (see Remark 13. ip . We
then get that for any ao satisfying ([73]) :
27r X4 e , F (6)n{o>ao;ma Xi [A i (a,6)-a i) . w (6)n 9 . (i , ) |<g} c H a , & )^
(76) < 2vr/+ oo | CF (a,6)|^<2 7 r||F|| B - :i a -<e.
Let us now consider the integral
I n da
27T / c F (a,b)— ,
J A e>F (b)n{ae(0,ao); max, |Ai(a,6)-S i , i ^(6)n^( i ,) |<e) a
and prove that under assumption ([75])
(77)
A £jF (6)n{a G (0,a ), maxlA^a^)-^^^)^^^! < e} = A £ , F {b)n{a G (0,a ); |o|V^(6)|-l| < A} .
Let us first assume that a G A £t p(b)n{a G (0, ao); |a|V^(6)| — 1| < A}. Since |a| V^(6)| —
1| < A, Lemma 16, 61 and assumption (j73p imply that
max | A; (a, 6) - d bi ^i{b)n dl{b) \ < ae 1 ^ (T 2 (a,b) + a\Vip e (b)\Fi(a, b)) < e ,
that is a G {a G (0,a ), max; |Aj(a,6) - ^^(6)n^(6)| < e}>
Conversely, let us assume that a G A £) i?(?)) H {a G (0, ao); maxj |Aj(a, b) — ^^(6)^(6)1 <
e}. Since a G j4 £j j7(6), |ci?(a,6)| > e. Then by Lemma E3] there exists some t! G {1, ••• ,L}
such that |a|V^/(6)| — 1| < A, otherwise assumption ([73]) and equation ([66]) would lead to a
contradiction. In addition, remark that if £' ^ £, one has
max|Aj(a,6)-<9 6 ^(&)rt^( 6 )| > \d bi ipi(b)n ei{b) -d bi ip e (b)n 8e , {b) \-max |Aj(a, 6) — (6)n^, (6) |
27
As above, Lemma 16.61 and assumption (|73p imply that
max | Aj (a, b) - d bi ip v (b)n 6t ,{b)\ < £ ■
Hence, we get that
max|Ai(a,6) - d bi tp e (b)n ei{b) \ > d [\d b .ip £ (b)\ + \d bt (pe(b)\] - e
> e
provided that e is sufficiently small, namely that:
(78) e< [-[\d bi w{b)\ + \d bi w
i.e. d [|^(6)l + \d bi ip e/ (b)\] > 2e, W,^,Vi
Hence, necessarily, I = £' and then a G A £ ^(b) n {a G (0, ao), |a|V^(6)| — 1| < A}.
In addition, one has
2ir I CF(a,b)a 2 da
' J 4 e , F (6)n{ae(0,a ),|a|V^(6)|-l|<A}
2tt
{ae(0,ao), |a|V<^(6)|
Ci?(a, 6)a 2 da — 2ir I
-1|<A} ./{,
ae(0,a ),|a|V^(b)j-l|<A}\A EiF (6)
cj?(a, 6)a
with |/| < e. Using this last equation and equations (|75|) and (jT6j) . we deduce that
li m J_ f f S S / e {b,k,n)dkdn - yl £ (6)e^ (b)n <W
5^0 (7^ y§i JB e (e,n,b)
<
2tt
+
C^. V{a6(0,o ),|a|V^(6)l-l|<A}
2vr f
y{ae(0,ao),|a|V^(6)|-l]<A}\ J 4 e]F (6)
c F (a, 6)a~ 2 da - A e (b)e Mb)n<> ^
cp(a, b)a~ 2 da + e .
Observe that the domain {a G (0, ao), |a|V(/?£(6)| — 1| < A} is bounded from below by o n
which only depends on F. In addition, if a G" A s f(6), inequality ([66]) is satisfied. Hence
(79)
/■fa
{ae(0,a ), \a\V(pe{b)\-l\<A}\A s ,F(b)
CF(a,b)a 2 da
< e
< e
< e
■CL [ su Pae(o,a )a 2r iM)
su Pag(o,a ) a 2 ri(a,6)
^min
SU Pae(0,a )« 2r i( a ^)
a 2 da
l-u
< Ce,
where C = (a m ; n inf fceK 2 |V^(6)|) _1 . The last inequality comes from assumption ([73]) valid
for £ > sufficiently small. This last inequality and inequality (I65p proved in Lemma 16.41
28
then imply that
<
lim
2vr r
Cij} J{ae(0,a ),\a\V<p e (b)\-l\<A}
2vr
+
J{a£(0,a ), |a|V</^(ft)|-l|<A}
S^(b, k, n)dkdn - A e (b)e w{b)n W
a^(aV^(6))^(6)e^ (6)n <w) a~ 2 da - A t {x)e^ x)ne i^)
[ea 2 Ti(a, b))a~ 2 da
+ (C + l)e
To conclude, observe now that the first term of the right-hand side vanishes. Indeed, the
wavelet is real isotropic and for any b G M 2 , a — > "ifi(aV(pi(b)) is supported in (0, a max ). Since
do > amax, we then have:
do
/ $(aV<p e (b))
J{oe(0,oo),|o|V^(6)|-l|<A} a
il>(aV<p t (b))
da _ J_
a 2n
12 ^
iei s
27T '
where CU has been defined in ([6]). The second term of the last inequality can be bounded
using the same approach than in inequality ()79p . We then get Theorem 14.11
A Quaternionic calculus
In this appendix, we give a short introduction to the algebra of quaternions H. More details
can be found for instance in [32J.
EI is the real 4-D vector space spanned by {l,i,j,k}, i.e. a quaternion is of the form q =
Qo + Qi 1 + °2j + 93k, where the algebra product is defined by i 2 = j 2 = k 2 = — 1 and ij = — ji =
k, jk = -kj = i, ki = -ik = j.
For a given quaternion 5 = 50 + 5ii + ffej + 53k, one defines:
• its real part: Re(g) = qo,
• its vectorial part: Vect(g) = (51,52,93) £
If Re(5) = 0, 5 = 5ii + q2] + 93k is called a pure quaternion.
The conjugate of the quaternion 5 is 5 = 50 — 5ii — 52j — 53k, and its norm (or modulus) is
defined by:
\q\ = Vtfl = \j ql + q\ + ql + <&
For any (5,5') G H 2 , one has:
55' = q' 5 and \qq'\ = \q\ \q'\ .
Additionally, the quaternions are a division algebra, which means that any (non-zero) quater-
nion admits a multiplicative inverse given by:
(80) q- 1 = .
\q\
29
Observe that one has: q 1 = -A™ = (q) . The set of unit quaternion {q £ EI ; |g| = 1} will
be denoted by S 3 .
The exponential function is defined on the quaternion algebra as
+ °° q e
exp :g ^e g = ^-,
which converges since exp(|g|) converges. Using the exponential map, each unit quaternion
of S 3 (\q\ = 1) such that Vect{q) ^ can be written as:
Vect(q)
(81) q = (cos ip+n sin ip) = ^ pn with n = — -, cosp = Re(q) and sin ip = \ Vect(q)\
\Vect{q)\
which is an extension of the complex exponential. Notice that the usual property e^e u = e^ +u
is no more satisfied in general. The polar form of any quaternion q is then given by:
(82) q= \q\ (cos (p + n sine/?) = \q\ e^ 71
where n is a pure unit quaternion: n = ai + 6j + ck with a 2 + b 2 + c 2 = 1, and
If (<p,n) £Rx§ 3 satisfies ®, one says that <p is a scalar argument of q and n a vectorial
orientation of q.
A quaternion of the form q = qo + Qi^ + Q2] with (qo, q±, 52) £R 3 (</4 = 0) is called a Clifford
vector. Observe that if q = ao+aii+a2j and q' = a^+a'^i+a^i with (ao, 01, 02), a 'i> a 2) ^ ^ 3
then
(83) ? = 7xg.
In this case any vectorial orientation of q is of the form:
n = ai + b] with a 2 + b 2 = 1 .
Then there exists 9 £ M such that n = cos i + sin 9 j . This implies that q admits the following
polar form:
(84) q=\q\ (cos ip + sin tp cos 9 i + sin ip sin j) = \q\ e^ (cos ei + sine J) .
If (</?, #) £ M 2 satisfies (|84p . one says that ip is a scalar argument of q and a scalar orientation
of g. Observe that <p and are also Euler angles of the vector \q\\ £t 3 .
w
B Hilbert and Riesz Transform
In this appendix we briefly recall the Hilbert transform and the analytic signal analysis, and
then proceed with a natural extension in two dimension, the Riesz transform, used to define
the 2-D monogenic signal in section [TJ
30
B.l The 1-D Hilbert Transform
We first define the Hilbert transform of a 1-D real valued function.
Definition B.l Let f G L 2 (R, R). The Hilbert transform of f , denoted as Hf is defined in
the Fourier domain as follows: for a.e. £ G R,
ft?(0 = "i sgn(0/(0 •
Definition B.2 The complex analytic signal associated to f G L 2 (R, R) is t/ien
F = f + iHf,
which reads in the Fourier domain as
F=(l+sgn(0)7.
We recall below some properties of the Hilbert transform and the analytic signal.
Proposition B.l The Hilbert transform is an antisymmetric operator on L 2 (R, R):
<nhj 2 >=-<h,nf 2 > .
Proposition B.2 1. f and rif are orthogonal
2. ||F|| 2 = 2||/|| 2
3. The spectrum of the analytic signal is one sided.
B.2 2-D Riesz transform
We begin with the definition of the 2-D Riesz transform of a real valued image.
Definition B.3 Let f G L 2 (R 2 ). The Riesz transform of f , denoted as IZf is the vector
valued function:
where for any i = 1,2, Ttif is defined in Fourier domain as follows: for a.e. £ G R 2 ,
nTf(0 = -i 47(0 •
We now present the key properties of 7Z (|33j.[27j). The first ones concern the invariance with
respect to dilations, translations, and the steerability property (relation with the rotations).
Proposition B.3 The Riesz transform commutes both with the translation, and the dilation
operator: for any f G L 2 (R 2 ) ; a > and b G R 2 , one has:
KD a f = D a lZf and KT b f = T b Kf .
31
Proposition B.4 The Riesz transform is steerable, that is, for any f E L 2 (M?) one has
(85) R e (nf)-r d u(R 9 f) - ^_ shld ni{Ref) + cos en 2 {Ref)) '
where rg is the rotation matrix defined by
Proof. We will prove in the Fourier domain that IZ(Rgf) = rgRg{lZf) , which is equivalent
to equation (185h . Let us first remark that for a.e. (£l 2 ,
(86) (R e f)^) = f( r ^) •
Using the definition of the Riesz transform in Fourier domain, one has, for a.e. £ £ R 2 :
nRd)(0 = r e r e l x f-ij-/( r -^
r e (Kf) (r-^)
r e (Wizf) (0 ,
the last relation coming from (|86p . One then deduces the required result. □
The Riesz transform is also a unitary and componentwise antisymmetric operator on L 2 (M 2 ):
Proposition B.5 For any i £ {1,2}, the i-th component of the Riesz transform IZi is an
antisymmetric operator, namely for all f,g£ L 2 (M. 2 ):
(87) < Kif, g > L 2 (R 2)= - < /, TZig
Since TZ 2 + TZ 2 = —Id, it implies in particular that:
< 7tf,1lg > L 2( R 2 !R 2)=< TZif,TZig >l2( K 2) + < IZ 2 f, 1Z 2 g >l2 (R 2)=< f,g
References
[1] D.J. Fleet and A.D. Jepson. Computation of component image velocity from local phase
information. Int J Comput Vision., 5(1):77-104, 1990.
[2] A. Basarab, H. Liebgott, and P. Delachartre. Analytic estimation of subsample spatial
shift using the phases of multidimensional analytic signals. IEEE Trans. Imag.Proc,
18(2), 2009.
[3] A. Basarab, P. Gueth, H. Liebgott, and P. Delachartre. Phase-based block matching
applied to motion estimation with unconventional beamforming strategies. IEEE Trans-
actions on Ultrasonics, Ferroelectrics, and Frequency Control, 56(5), 2009.
32
[4] M. Elshinawy and M. Chouikha. Using AM-FM image modeling technique in mammo-
grams. Micro-NanoMechatronics and Human Science, 2003 IEEE International Sympo-
sium on., 2, 2003.
[5] J. Elshinawy, M. Zeng, S.C.B. LO, and M. Chouikha. Breast cancer detection in mammo-
gram with AM-FM modeling and gabor filtering, in Proc 7th International Conference
on Signal Processing, 3, 2004.
[6] I. Kokkinos, G. Evangelopoulos, and P. Maragos. Texture analysis and segmentation
using modulation features, generative models, and weighted curve evolution. IEEE Trans
Pattern Anal Mach Inteli, 31(1):142-157, 2009.
[7] P. Tay. AM-FM image analysis using the Hilbert Huang transform. Southwest Symposium
on Image Analysis and Interpretation SSIAI, pages 13-16, 2008.
[8] A. Belaid, D. Boukerroui, Y. Maingourd, and J.F. Lerallut. Phase-based level set seg-
mentation of ultrasound images. IEEE Trans Inf Tech Biomed., 15(1):138-147, 2011.
[9] M. Felsberg and G. Sommer. The monogenic signal. IEEE Transactions on Signal
Processing, 49(12):3136-3144, 2001.
[10] V. Murray, M.S. Pattischis, E.S. Barriga, and P. Soliz. Recent multiscale AM-FM meth-
ods in emerging applications in medical imaging. EURASIP Journal on Advances in
Signal Processing, 23:1-14, 2012.
[11] T. Qian. Analytic signals and harmonic measures. Journal of Mathematic Analysis and
Applications, 314:526-536, 2006.
[12] N.E Huang, Z. Wu, S.R. Long, K..C. Arnold, K.. Blank, and T.W.. Liu. On instantaneous
frequency. Advances in Adaptive Data Analysis, 1:177-229, 2009.
[13] P.V. Rodriguez and M.S. Pattichis. New algorithms for fast and accurate AM-FM
demodulation of digital images. IEEE International Conference on Image Processing,
2005.
[14] R.A. Carmona, W.L. Hwang, and B. Torresani. Characterization of Signals by the Ridges
of Their Wavelet Transforms. IEEE Transactions on Signal Processing, 45(10):2586-
2590, 1997.
[15] R.A. Carmona, W.L. Hwang, and B. Torresani. Multiridge Detection and Time-
Frequency Reconstruction. IEEE Transactions on Signal Processing, 47(2):480-492,
1999.
[16] I. Daubechies and S. Maes. A nonlinear squeezing of the continuous wavelet transform
based on auditory nerve models. A.Aldroubi, M.Unser (Eds.), Wavelets in Medicine and
BiologyCRC Press, 1996.
[17] F. Auger and P. Flandrin. Improving the readability of time-frequency and time-scale
representations by the reassignment method. IEEE Trans. Signal Process., 43(5):1068-
1089, 1995.
33
[18] N.E Huang, Z. Shen, S.R. Long, S.C. Wu, H.H. Shih, Q. Zheng, N.C. Yen, C.C. Tung,
and H.H. Liu. The empirical mode decomposition and the Hilbert spectrum for nonlinear
and non-stationary time series analysis. Proc. R. Soc. A, 454:903-995, 1998.
[19] P. Flandrin, G. Rilling, and P. Goncalves. Empirical mode decomposition as a filter
bank. IEEE Signal Process. Lett, 11(2):112-114, 2004.
[20] G. Rilling and P Flandrin. One or two frequencies? The empirical mode decomposition
answers. IEEE Trans. Signal Process., 56(l):85-95, 2008.
[21] I. Daubechies, J. Lu, and H.-T. Wu. Synchrosqueezed Wavelet Transforms: an Em-
pirical Mode Decomposition-like Tool. Applied and Computational Harmonic Analysis,
30(2):243-261, 2011.
[22] H.T Wu, P. Flandrin, and I. Daubechies. One or two frequencies?The synchrosqueezing
answer. Advances in Adaptive Data Analysis, 3(l-2):29-39, 2011.
[23] S. Meignen, T. Oberlin, and S. Mclaughlin. A new algorithm for multicomponent sig-
nal analysis based on synchrosqueezing: With an application to signal sampling and
denoising. IEEE Transactions on Signal Processing, 60(ll):5787-5798, 2012.
[24] M. Biilow and G. Sommer. A novel approach to the 2D analytic signal. Proc. CAIP99,
F. Solina and A. Leonardis, Eds. Ljubljana, Slovenia, pages 25-32, 1999.
[25] Y. Yang, T. Qian, and F. Sommen. Phase Derivative of Monogenic Signals in Higher
Dimensional Spaces. Complex analysis and operator theory, 2011.
[26] S.C. Olhede and G. Metikas. The Monogenic Wavelet Transform. IEEE Transactions
on Signal Processing, 57(9):3426-3441, 2009.
[27] M. Unser and D. Van De Ville. Multiresolution Monogenic Signal Analysis Using the
Riesz-Laplace Wavelet Transform. IEEE Trans. Imag. Proc, 18(11):2402-2418, 2009.
[28] J.-P. Antoine, R. Murenzi, P. Vandergheynst, and S. Twareque Ali. Two-dimensional
Wavelets and their Relatives. Cambridge University Press, 2004.
[29] S.C. Olhede and G. Metikas. Multiple Multidimensional Morse Wavelets. IEEE Trans-
actions on Signal Processing, 55(3):921-936, 2007.
[30] S. Jaffard, Y. Meyer, and R.D. Ryan. Wavelets: tools for science and technology. SIAM,
2001.
[31] H. Triebel. Interpolation theory, function spaces, differential operators. Amsterdam,
North-Holland, 1978.
[32] A Sudberry. Quaternionic analysis. Math. Proc. Camb. Phil. Soc, 85:199-225, 1979.
[33] E. M. Stein. Singular Integrals Differentiability Properties of Functions. Princeton Uni-
versity Press, New York, second edition, 1970.
34