Construction of Hilbert Transform Pairs of Wavelet Bases and Gabor-like Transforms

Construction of Hilb ert T ransform P airs of W a v elet Bases and Gab or-lik e T ransforms Kunal Nara yan Chaudh ury and Mic hael Unser ∗ Octob er 24, 2018 Abstract W e prop ose a no v el metho d for constructing Hilbert transform (HT) pairs of w av elet bases based on a fundamen tal approximation-theoretic c haracterization of scaling functions—the B-spline factorization theorem. In particular, starting from w ell-lo calized scaling functions, we construct HT pairs of biorthogonal wa velet bases of L 2 ( R ) by relating the corre- sp onding wa velet ﬁlters via a discrete form of the contin uous HT ﬁlter. As a concrete application of this metho dology , we iden tify HT pairs of spline w av elets of a sp eciﬁc ﬂa vor, which are then combined to realize a family of complex w av elets that resem ble the optimally-lo calized Gabor function for suﬃcien tly large orders. Analytic wa velets, derived from the complexiﬁcation of HT wa velet pairs, exhibit a one-sided spectrum. Based on the tensor-product of such analytic wa velets, and, in eﬀect, by appropriately combining four separa- ble biorthogonal wa v elet bases of L 2 ( R 2 ), we then discuss a methodology for constructing 2D directional-selectiv e complex wa velets. In particu- lar, analogous to the HT corresp ondence b et ween the comp onen ts of the 1D counterpart, we relate the real and imaginary comp onen ts of these complex w av elets using a multi-dimensional extension of the HT—the di- rectional HT. Next, we construct a family of complex spline wa v elets that resem ble the directional Gab or functions proposed b y Daugman. Finally , w e present an eﬃcien t FFT-based ﬁlterbank algorithm for implementing the asso ciated complex wa velet transform. 1 INTR ODUCTION The dual-tree complex w av elet transform (DT- C WT) is a recent enhancemen t of the con ven tional discrete wa v elet transform (DWT) that has gained increasing p opularit y as a signal pro cessing to ol. The transform was originally introduced ∗ Corresponding Author: Kunal Naray an Chaudhury . The authors are with the Biomed- ical Imaging Group, Ecole Polytechnique F ´ ed´ erale de Lausanne (EPFL), Station-17, CH- 1015 Lausanne VD, Switzerland. F ax: +41 21 693 37 01, e-mail: { kunal.c haudhury , michael.unser } @epﬂ.c h. This work w as supp orted b y the Swiss National Science F oundation under gran t 200020-109415. 1 b y Kingsbury [17, 19] to circumv en t the shift-v ariance of the decimated DWT, and in volv ed t wo DWT channels in parallel with the corresp onding wa velets forming a quadrature pair. In particular, Kingsbury realized the quadrature relation b y in terp olating the lowpass ﬁlters of one D WT “mid-wa y” betw een the lo wpass ﬁlters of the other DWT. Moreo ver, based on appropriate combinations of separable w a velets, he extended the dual-tree construction to tw o-dimensions, where the corresp onding transform, b esides impro ving on the shift-inv ariance of the 2D D WT, exhibits b etter direction selectivity as w ell. There is no w go od evidence that the transform tends to p erform better than its real coun terpart in a v ariet y of applications such as suc h as deconv olution [10], denoising [29], and texture analysis [16]. The crucial observ ation that the dual-tree wa v elets inv olv ed in Kingsbury’s construction form an approximate HT pair was made by Selesnic k [28, 26]. He also demonstrated that a particular phase relation b etw een the lowpass (reﬁne- men t) ﬁlters of the tw o c hannels resulted in the desired HT corresp ondence. This link consequently transp osed the problem of designing diﬀeren t ﬂa v ors of dual-tree wa velets to that of identifying new HT pairs of wa v elets. Indeed, fol- lo wing this remark able connection, several new paradigms and extensions hav e b een prop osed: design of HT pairs of biorthogonal wa velet bases [40], alterna- tiv e frameworks for complex non-redundant transforms [12], and the M-band extension [8], to name a few. 1.1 Motiv ation The deploymen t of complex signal represen tations for the determination of in- stan taneous amplitude and frequency is classical [13, 38]. Gab or and Ville [13, 39] prop osed to unambiguously deﬁne them using the concept of the an- alytic signal —a unique complex-v alued signal representation speciﬁed using the HT. Sp eciﬁcally , the analytic signal s a ( x ) = s ( x ) + j H s ( x ) corresp ond- ing to a real-v alued signal s ( x ) ( H denotes the HT op erator), w as used to stipulate the instantaneous amplitude and phase via the p olar representation s a ( x ) = A ( x )e j φ ( x ) . In particular, this represen tation allows one to retriev e the time-v arying amplitude and frequency of an AM-FM signal of the form s ( x ) = A ( x ) cos  2 π R x 0 ν ( τ ) dτ + ξ 0  via the estimates A est ( x ) = | s a ( x ) | and ν est ( x ) = (2 π ) − 1 dφ ( x ) /dx , assuming A ( x ) to b e slo wly-v arying compared to ν ( x ). The analytic signal has b ecome an imp ortan t complex-v alued representa- tion in signal pro cessing, especially in applications such as phase and frequency mo dulation, speech recognition and pro cessing of seismic data. These concepts ha ve also b een transp osed to the m ulti-dimensional setting: the lo cal frequency has b een used as a measure of lo cal signal scale; structures such as lines and edges hav e b een distinguished using the lo cal phase; and the lo cal amplitude and phase ha ve b een used for edge detection and for texture and ﬁngerprin t analysis [25]. The adv antage of viewing the dual-tree wa velets as a HT pair is that we can mak e a direct connection with the formalism of analytic signals. Indeed, if w e transp ose the abov e concept to the w av elet domain and consider the input signal 2 to b e lo cally of the AM-FM form, we obtain a resp onse where the lo cal energy of the signal is encoded in the magnitude of the wa v elet coeﬃcients, while the relativ e displacement is captured by the phase. In fact, this turns out to b e the fundamen tal reason for the superiority of the DT- C WT ov er con ven tional real-v alued transforms whose resp onse is necessarily oscillating. 1.2 Our Con tribution In this con tribution, w e in v oke the B-spline factorization theorem [37]—a fun- damen tal sp ectral factorization result—along with certain fractional B-spline calculus [36], to construct HT pairs of biorthogonal wa velets from well-localized scaling functions. In particular, we do so b y relating the corresp onding w av elets ﬁlters via a discrete version of the con tinuous HT ﬁlter. Next, we iden tify a family of analytic spline w av elets, of increasing v anishing momen ts and regularity , that asymptotically con verge to Gab or-like functions [13]. As far as the implementation is concerned, unlik e Kingsbury’s sc heme that uses diﬀerent ﬁlters for diﬀerent stages (often with ﬁlter-swapping b et ween the dual-trees), our implemen tation uses the same set of ﬁlters at all stages of the ﬁlterbank decomp osition. Notably , we use an appropriate pair of pro jection ﬁlters for coherent signal analysis which, in turn, allows us to identify a discrete coun terpart of the analytic wa v elet—the so-called analytic wavelet ﬁlter that exhibits a one-sided sp ectrum. The construction is then extended to t wo-dimensions through appropriate tensor-pro ducts of the one-dimensional analytic w a velets. In particular, we con- struct a family of directional complex wa v elets that resemble the directional Gab or functions prop osed by Daugman [9] for suﬃcien tly large orders. More- o ver, we also relate the real and imaginary comp onen ts of the complex wa v elets using the directional HT—a m ultidimensional extension of the HT—that pro- vides further insight into the directional-selectivity of the dual-tree wa velets. 1.3 Organization of the Paper W e b egin b y recalling certain fundamental deﬁnitions and prop erties p ertain- ing to the HT and the fractional B-splines in § 2. W e characterize the action of the HT op erator on B-splines in § 3, which, along with the B-spline factor- ization theorem, is used to prop ose a formalism for constructing HT pairs of biorthogonal w av elet bases in § 4. The implemen tation asp ects are discussed in § 5. As a concrete application, w e construct the Gab or-lik e wa v elets in § 6. In § 7, directional complex wa velets are constructed by appropriately combining the wa velets corresp onding to certain separable multiresolution analyses; the highligh t of this section is the construction of 2D Gab or-like spline wa velets. The implementation asp ects of the corresp onding 2D Gabor-like transform are pro vided in § 8, b efore concluding with § 9. 3 2 PRELIMINARIES W e b egin b y introducing sp eciﬁc op erators and functions that play a ma jor role in the sequel follow ed by a discussion of their relev ant prop erties. In what follo ws, we use ˆ f ( ω ) = R R d f ( x )e − j x T ω d x to denote the F ourier transform of a function f ( x ) on R d ( d > 1), with x T ω b eing the usual inner-pro duct on R d . W e also frequen tly use the notations f ( · − s ) and f ( λ · ), corresponding to some s in R d and λ > 0, to denote the function obtained by translating (resp. dilating) f ( x ) by s (resp. λ ). W e denote the Kroneck er-delta sequence b y δ [ n ]: its v alue is 1 at n = 0, and is zero at all other integers. 2.1 Hilb ert T ransform and W a v elets The Hilb ert transform, that generalizes the notion of the quadrature transfor- mation cos( ω 0 x ) 7→ sin( ω 0 x ) b ey ond pure sinusoids [5], forms the cornerstone of this pap er. F rom a signal-pro cessing p erspective, the HT can b e interpreted as a ﬁltering op eration in which the amplitude of the frequency components is left unchanged, while their phase is altered by ± π / 2 depending on the sign of the frequency . Mathematically , the HT of a suﬃciently w ell-b eha ved function is deﬁned using a singular integral transform [2, 30]. Ho wev er, in the context of ﬁnite- energy signals, it admits a particularly straightforw ard form ulation based on the F ourier transform on L 2 ( R ). In particular, the Hilb ert transform on L 2 ( R ) is characterized by the equiv alence H f ( x ) F ← → − j sign( ω ) b f ( ω ) (1) where the m ultiplier sign( ω ) is deﬁned as ω / | ω | for non-zero ω , and as zero at ω = 0. Based on the ab o ve deﬁnition 1 , and the prop erties of the F ourier transform on L 2 ( R ), the following prop erties of the HT can b e readily derived: • Line arity and T r anslation-invarianc e : It is a linear and translation-inv ariant op erator; that is, it acts as a conv olution op erator. • Dilation-Invarianc e : It commutes with dilations: H{ f ( λ · ) } ( x ) = ( H f )( λx )for all λ > 0. • Anti-Symmetry : It an ti-commutes with the ﬂip op eration f T ( x ) = f ( − x ), so that ( H f T )( x ) = − ( H f ) T ( x ); thus the HT of a symmetric function is necessarily anti-symmetric. • Unitary (Isometric) Natur e : It acts as a unitary op erator on L 2 ( R ), so that hH f , H g i = h f , g i for all f and g , where h· , ·i denotes the usual inner- pro duct on L 2 ( R ). Equiv alently , this means that the inv erse HT op erator is given by its adjoint: H − 1 = H ∗ . 1 The deﬁnition can also b e extended to temp ered distributions such as the Dirac delta and the sin usoid [2, § 2.5]. 4 It is w ell-known that HT of a wa v elet is also a w av elet. The implication of the simultaneous in v ariance to dilations and translations is that the HT of a dilated-translated wa velet is a wa v elet, dilated and translated by the same amoun t: H{ ψ ( λx − s ) } = ( H ψ )( λx − s ). Moreov er, an immediate consequence of the unitary prop ert y is that the HT op erator maps a basis into a basis: if { ψ n } form a (Riesz) wa v elet basis of L 2 ( R ), then so do es {H ψ n } . It even preserves biorthogonal wa velet bases of L 2 ( R ): if { ψ n } and { ˜ ψ m } form a biorthogonal w av elet basis of L 2 ( R ), satisfying the duality criteria h ψ n , ˜ ψ m i = δ [ m − n ], then using the same unitary prop ert y , we hav e hH ψ n , H ˜ ψ m i = h ψ n , ˜ ψ m i = δ [ m − n ] (2) so that {H ψ n } and {H ˜ ψ m } form a biorthogonal w av elet basis of L 2 ( R ) as well. It is exactly the abov e inv ariance prop erties that make the construction of HT pair of wa v elet bases of L 2 ( R ) feasible. Unfortunately , the HT exhibits certain inheren t pathologies in the context of multiresolution analyses and wa velets. The impulse resp onse of the HT, H δ ( x ) = 1 /π x (in the sense of distributions), clearly indicates the non-lo cal nature of the op erator. This has t wo serious implications: ( i ) the HT of a compactly-supp orted scaling function/w a velet is no longer of ﬁnite supp ort; ( ii ) the HT-transformed function has a O (1 / | x | )-deca y in general, and hence is not in tegrable; and ( iii ) the (anti-symmetric) HT suppresses the dc-comp onen t of symmetric scaling functions that is essential for fulﬁlling the partition-of-unity criterion. Therefore, the HT of a scaling function is not a v alid scaling function, and cannot b e used to sp ecify a m ultiresolution analysis in the sense of Mallat and Meyer [22, 23]. Next, we recall the notion of an analytic signal that generalizes the phasor transformation transformation cos( ω 0 x ) 7→ exp( j ω 0 x ) to ﬁnite-energy signals using the HT as the quadrature transformation. In general, the analytic signal f a ( x ) asso ciated with a real-v alued signal f ( x ) is deﬁned as the complex-v alued signal f a ( x ) = f ( x ) + j H f ( x ) . (3) In particular, f a ( x ) = exp( j ω 0 x ) when f ( x ) = cos( ω 0 x ). Imp ortan tly , note that the F ourier transform of the analytic signal ev aluates to b f a ( ω ) = (1 + sign( ω )) b f ( ω ), so that b f a ( ω ) v anishes for all negativ e frequencies. It is exactly this one-sided sp ectrum that makes the analytic signal particularly interesting in signal pro- cessing [13]; we exploit this property for constructing directional w av elets in § 7. 2.2 F ractional B-spline Multiresolution The family of fractional B-splines [4]—fractional extensions of the polynomial B-splines—will pla y a key role in the sequel. In particular, we recall that the fractional B-spline β α τ ( x ), corresp onding to a degree α ∈ R + 0 and a shift τ ∈ R , 5 is sp eciﬁed by its F ourier transform β α τ ( x ) F ← →  1 − e − j ω j ω  α +1 2 + τ  1 − e j ω − j ω  α +1 2 − τ . (4) The parameters α and τ control the width and the av erage group dela y of the scaling function resp ectiv ely . In particular, when τ = ( α + 1) / 2, the fractional B-spline β α τ ( x ) corresp onds to the causal B-spline β α + ( x ) deﬁned in [36]. The fractional B-splines, in general, do not ha ve a compact supp ort (except for in teger degrees); how ever, their O (1 / | x | α +2 ) decay ensures their inclusion in L 1 ( R ) ∩ L 2 ( R ). Another relev ant prop ert y that will b e in vok ed frequently is that the shift τ inﬂuences only the phase of the F ourier transform; that is, | b β α τ ( ω ) | is indep endent of τ . The fundamental role pla yed by fractional B-splines in this pap er is, ho wev er, based on the fact they satisfy certain admissibility criteria [36, 4] needed to generate a v alid multiresolution of L 2 ( R ): ( C1 ) The appro ximation space V 0 = span ` 2 { β α τ ( · − k ) } k ∈ Z admits a stable Riesz basis. ( C2 ) There exists an in tegrable sequence h α τ [ k ] (reﬁnement ﬁlter) such that the t wo-scale relation 1 2 β α τ  x 2  = X k ∈ Z h α τ [ k ] β α τ ( x − k ) (5) holds. In particular, the transfer function of the reﬁnemen t ﬁlter is sp eciﬁed by H α τ (e j ω ) = 1 2 α +1  1 + e j ω  α +1 2 − τ  1 + e − j ω  α +1 2 + τ . (6) ( C3 ) P artition of unit y: The in teger-translates of β α τ ( x ) can repro duce the unit y function. W e brieﬂy discuss the signiﬁcance of these admissibility conditions. The criterion ( C1 ) ensures a stable and unique represen tation of functions in V 0 using co eﬃcien ts from ` 2 ( Z ); equiv alently , this also signiﬁes that the transfer function of the auto correlation (Gram) ﬁlter, A α (e j ω ) = P k ∈ Z | b β α τ ( ω + 2 π k ) | 2 , is uniformly b ounded from ab o ve, and aw ay from zero [21]. On the other hand, ( C2 ) implies the inclusion of β α τ ( x/ 2) in V 0 , which, in turn, allows one to deﬁne a hierarchical embedding of approximation spaces {V j } j ∈ Z that is key to the m ultiresolution structure of the asso ciated wa v elet transform. Finally , the tech- nical condition ( C3 ) ensures that the multiresolution {V j } is dense in L 2 ( R ): arbitrarily close appro ximations of functions in L 2 ( R ) can b e achiev ed using elemen ts from {V j } . 3 HILBER T TRANSF ORM AND B-SPLINES It turns out that the action of the HT on B-splines can b e eﬀectively character- ized in terms of certain fractional ﬁnite-diﬀerence (FD) operators. In particular, 6 corresp onding to an order α ∈ R + 0 and shift τ ∈ R , we consider the op erator ∆ α τ deﬁned on L 2 ( R ) by ∆ α τ f ( x ) F ← → D α τ (e j ω ) b f ( ω ) , (7) where D α τ (e j ω ) =  1 − e − j ω  α 2 + τ  1 − e j ω  α 2 − τ . One recov ers the conv en tional n -th order FD op erator by setting α = n and τ = n/ 2. Since the op erator has a p eriodic frequency resp onse, one can asso ciate with it a digital ﬁlter d α τ [ k ] through the corresp ondence ∆ α τ f ( x ) = P d α τ [ k ] f ( x − k ) . The FD op erator that is esp ecially relev ant for our purp ose is the zeroth-order op erator ∆ 0 − 1 / 2 (henceforth, w e simply denote it by ∆). The corresp onding frequency resp onse D (e j ω ) = D 0 − 1 / 2 (e j ω ) reduces 2 to D (e j ω ) = − j sign( ω )e − j ω 2 for ω ∈ ( − π , π ) , (8) signifying that D (e j ω ) is in L 2 (( − π , π )); the corresp onding ﬁlter co eﬃcien ts d [ k ] = d 0 − 1 / 2 [ k ] in ` 2 ( Z ) are then sp eciﬁed 3 b y d [ k ] = 1 2 π Z π − π D (e j ω )e j kω d ω = 1 2 π Z π − π − j sign( ω )e j ( k +1 / 2) ω d ω = 1 π ( k + 1 / 2) , k ∈ Z . (9) Th us, similar to the HT op erator, ∆ is also unitary , and the corresp onding ﬁlter d [ k ] can b e in terpreted as a discrete form of the contin uous HT ﬁlter 1 /π x . In particular, we can relate the action of the HT on the B-splines solely in terms of this ﬁlter. Indeed, it can easily b e seen that the F ourier transform of the B-spline can b e factorized as b β α τ ( ω ) = ( j ω ) 1 / 2 ( − j ω ) − 1 / 2 D (e j ω ) b β α τ +1 / 2 ( ω ) , (10) whic h, along with the identit y ( j ω ) 1 2 ( − j ω ) − 1 2 = j sign( ω ), results in the equiv- alence H β α τ ( x ) F ← → − j sign( ω ) b β α τ ( ω ) = − j sign( ω ) · j sign( ω ) D (e j ω ) b β α τ +1 / 2 ( ω ) = D (e j ω ) b β α τ +1 / 2 ( ω ) F ← → ∆ β α τ +1 / 2 ( x ) , (11) that establishes the desired result: 2 W e sp ecify the fractional p ow er of a complex num b er z by z γ = | z | γ e j γ arg( z ) corre- sponding to the principal argument | arg( z ) | < π . On this principal branch, the identit y ( z 1 z 2 ) γ = z γ 1 z γ 2 holds only if arg( z 1 ) + arg ( z 1 ) ∈ ( − π , π ) [31, Chapter 3]. 3 The in verse F ourier transform ov er the principal perio d ( − π, π ) is inv oked. 7 Prop osition 3.1 The HT of a fr actional B-spline c an b e expr esse d as H β α τ ( x ) = X k ∈ Z 1 π ( k + 1 / 2) β α τ +1 / 2 ( x − k ) . (12) In particular, the digital ﬁlter d [ k ] acts as a unitary conv olution op erator on L 2 ( R ) when applied to functions, and as a discrete ﬁlter on ` 2 ( Z ) when applied to sequences. The theoretical diﬃculty with the HT stems from the fact that its frequency resp onse has a singularity at ω = 0 , which results in a p oor decay of the transformed output. The remark able feature of (12) is that w e ha ve b een able to express the slowly decaying HT as a linear combination of the b etter-behav ed B-splines. Sp eciﬁcally , the sequence d [ k ] decays only as O (1 / | k | ), whereas β α τ +1 / 2 ( x ) decays as O (1 / | x | α +2 ). Th us, by expressing the HT using shifted B-splines as in (12), we hav e, in eﬀect, mov ed the singularity on to the digital ﬁlter. In the sequel, we shall apply this ﬁlter to the w av elets where its eﬀect is m uch more inno cuous since b ψ ( ω ) = 0 around the origin. Half-Delay Filters : As remark ed earlier, the shift parameter τ only aﬀects the phase of the F ourier transform of the fractional B-spline and the corresp ond- ing reﬁnement ﬁlter [4]. In particular, based on the factorization H α τ + 1 2 (e j ω ) = 1 2 α +1  1 + e j ω  α +1 2 − ( τ + 1 2 )  1 + e − j ω  α +1 2 +( τ + 1 2 ) = (1 + e j ω ) − 1 / 2 (1 + e − j ω ) 1 / 2 H α τ (e j ω ) = e − j ω 2 H α τ (e j ω ) , for ω ∈ ( − π , π ) , w e arrive at the follo wing result: Prop osition 3.2 The spline r eﬁnement ﬁlters h α τ [ k ] and h α τ +1 / 2 [ k ] ar e “half- sample” shifte d versions of one another in the sense that H α τ +1 / 2 (e j ω ) = e − j ω 2 H α τ (e j ω ) (13) for al l ω in ( − π , π ) . Indeed, if we consider the bandlimited function h α τ ( x ) = P h α τ [ k ]sinc( x − k ) that satisﬁes the constraint h α τ ( x ) | x = k = h α τ [ k ], then we ha ve, as a consequence of (13), the relation h α τ +1 / 2 [ k ] = h α τ ( k − 1 / 2): each ﬁlter provides the bandlimited in terp olation of the other mid-wa y b et ween its samples. Finally , we mak e a note of the fact that the abov e reﬁnement ﬁlters can also b e related through a c onjugate-mirr or e d version of the FD ﬁlter: H α τ +1 / 2 (e j ω ) = D ( − e − j ω ) H α τ (e j ω ) . (14) 8 4 HT P AIR OF W A VELET BASES Before stating the main results, we recall the approximation-theoretic notion of appr oximation or der , and a fundamental sp ectral factorization result in volving B-splines. Appr oximation Or der : Scaling functions pla y a fundamental role in wa v elet theory . The technical criteria for a v alid scaling function w as discussed earlier in the context of B-splines (cf. § 2.2). Next we recall the fundamental notion of order for a scaling function that characterizes its approximation p o wer [37]. A scaling function ϕ ( x ) is said to hav e an appro ximation order γ if and only if there exists a p ositiv e constant C suc h that for all elements of the Sob olev space W γ 2 ( R ), of order γ , we hav e the estimate || f − P a f || 6 C a γ || ∂ γ f || . (15) Here P a denotes the pro jection operator from W γ 2 ( R ) on to the appro ximation subspace span ` 2 { ϕ ( · /a − k ) } k ∈ Z , and ∂ γ denotes the (distributional) deriv ative of order γ . In other w ords, the approximation order provides a c haracterization of the rate of decay of the approximation error for suﬃciently regular functions as a function of the scale. It turns out that, akin to their p olynomial counterparts, the order of frac- tional B-splines is entirely controlled b y their degree [36, 4]; in particular, we ha ve γ = α + 1. Equiv alen tly , this signiﬁes that any p olynomial of degree 6 d α e can b e repro duced by the set { β α τ ( · − k ) } , which is crucial for capturing the lo wpass information in images is concerned. Char acterizazion of Sc aling F unctions : A fundamental result in wa v elet the- ory is that it is alwa ys p ossible to express a v alid scaling function as a conv olu- tion b et ween an fractional B-spline and a distribution [37]. The original result in [37] inv olves causal B-splines; how ev er, the result can readily b e extended to the more general fractional B-splines since the shift parameter τ do es not inﬂuence the order of the scaling function. Indeed, note the theorem in [37] asserts that H (e j ω ) is the reﬁnement ﬁlter of a v alid scaling function (cf. § 2.2) of order α + 1 if and only if it can b e factorized as H (e j ω ) = β α + spline part z }| {  1 + e − j ω 2  α +1 distributional part z }| { Q (e j ω ) , (16) where Q (e j ω ) is stable: | Q (e j ω ) | < C < + ∞ for all ω . Rewriting (16) in terms of a ( α, τ ) B-spline reﬁnemen t ﬁlter, w e then hav e the follo wing equiv alent represen tation: H (e j ω ) = β α τ spline part z }| {  1 + e − j ω 2  α +1 2 + τ  1 + e j ω 2  α +1 2 − τ distributional part z }| { P (e j ω ) (17) 9 with P (e j ω ) = e − j ω ( α +1 2 − τ ) Q (e j ω )for ω ∈ ( − π , π ). Note that P (e j ω ) is stable, with | P (e j ω ) | < C < + ∞ for all ω . That is, H (e j ω ) is the reﬁnement ﬁlter of a v alid scaling function of order α + 1 if and only if it admits a stable factorization as in (17). W e then arrives at the following extension: Theorem 4.1 ( B-spline F actorization ) A valid sc aling function ϕ ( x ) is of or- der α + 1 if and only if its F ourier tr ansform c an b e factorize d as b ϕ ( ω ) = b β α τ ( ω ) b ϕ 0 ( ω ) (18) for some τ ∈ R , wher e b ϕ 0 ( ω ) is a function of ω that is b ounde d on every c omp act interval, and e quals unity at the origin. In the signal domain, this corresp onds to a w ell-deﬁned conv olution ϕ ( x ) = ( β α τ ∗ ϕ 0 )( x ) b etw een a B-spline and the temp ered distribution ϕ 0 . The crux of the abov e result is that it is the constituent B-spline that is solely resp onsible for the approximation prop ert y , and other desirable features of the scaling function [37]. 4.1 Construction of HT Pairs of W a velets In what follows, w e use the notation f j,k ( x ), corresp onding to a function f ( x ), and in tegers j and k , to denote the (normalized) dilated-translated function 2 j / 2 f (2 j x − k ). The HT of a wa velet is also a wa velet in a well-deﬁned sense. In particular, if ψ ( x ) is a wa v elet whose dilations-translations { ψ j,k } form a Riesz basis of L 2 ( R ), then H ψ ( x ) is also a v alid wa velet with {H ψ j,k } constituting a Riesz basis of L 2 ( R ). As remarked earlier, this follows from the fundamental in v ariance prop erties of the HT. W e now establish a formalism for constructing the HT of a given w av elet ψ ( x ). In particular, if ϕ ( x ) be the asso ciated scaling function, say of order α + 1, and g [ k ] b e the generating wa v elet ﬁlter, then we hav e the relation ψ ( x/ 2) = X g [ k ] ϕ ( x − k ) . (19) F ollowing Theorem 4.1, let us factorize ϕ ( x ) as ϕ ( x ) = ( β α τ ∗ ϕ 0 )( x ) (20) corresp onding to some real τ . Then, consider the scaling function ϕ 0 ( x ), of the same order, speciﬁed b y ϕ 0 ( x ) = ( β α τ +1 / 2 ∗ ϕ 0 )( x ). Let ψ 0 ( x ) b e any arbitrary w av elet, corresp onding to the multiresolution analysis associated with ϕ 0 ( x ), that is sp eciﬁed by ψ 0 ( x/ 2) = P g 0 [ k ] ϕ 0 ( x − k ). W e then hav e the following necessary and suﬃcien t condition for the desired HT correspondence in terms of the discrete HT ﬁlter d [ k ] (see § 11.1 for a pro of ): Theorem 4.2 ( HT Pair of W a velets ) The wavelets ψ ( x ) and ψ 0 ( x ) have the c orr esp ondenc e ψ 0 ( x ) = H ψ ( x ) if and only if g 0 [ k ] = ( d ∗ g )[ k ] . 10 Moreo ver, the construction has the following c haracteristics: • Both ϕ ( x ) and ϕ 0 ( x ) hav e the same Riesz b ounds and the same deca y , • The reﬁnemen t ﬁlters H (e j ω ) and H 0 (e j ω ) corresp onding to ϕ ( x ) and ϕ 0 ( x ) resp ectiv ely , are related as H 0 (e j ω ) = e − j ω 2 H (e j ω ) for all ω in ( − π , π ). The equality of the Riesz b ounds follows from the observ ation that the auto cor- relation ﬁlters of ϕ ( x ) and ϕ 0 ( x ) are identical. Indeed, we hav e a [ k ] = h ϕ , ϕ ( · − k ) i = ( β 2 α +1 0 ∗ ϕ 0 ∗ ϕ T 0 )( k ) , a 0 [ k ] = h ϕ 0 , ϕ 0 ( · − k ) i = ( β 2 α +1 0 ∗ ϕ 0 ∗ ϕ T 0 )( k ) . The assertion regarding the decay is based on the observ ation that b oth β α τ ( x ) and β α τ +1 / 2 ( x ) hav e the same decay . Finally , using (13) and (17), we can relate the transfer functions on ( − π , π ) as follows H 0 (e j ω ) = H α τ +1 / 2 (e j ω ) P (e j ω ) = e − j ω 2 H α τ (e j ω ) P (e j ω ) = e − j ω 2 H (e j ω ) , where P (e j ω ) denotes the transfer function of the ﬁlter asso ciated with the distribution ϕ 0 . R emark : Note that although H ψ ( x ) is unique, the scaling function ϕ 0 ( x ) and the corresp onding ﬁlter g 0 [ k ] generating H ψ ( x ) are by no means unique. F or instance, the particular choice ϕ 0 ( x ) = H ϕ ( x ) and g 0 ≡ g is suﬃcient to ensure that ψ 0 ( x ) = H ψ ( x ). Moreov er, if ϕ 0 ( x ) and g 0 [ k ] generate the w av elet ψ 0 ( x/ 2) = P g 0 [ k ] ϕ 0 ( x − k ) such that ψ 0 ( x ) = H ψ ( x ), then so do ϕ 0 eq ( x ) = P r [ k ] ϕ 0 ( x − k ) and g 0 eq [ k ] = ( g 0 ∗ r inv )[ k ]. Here the ﬁlter r [ k ] is such that 0 < | P k r [ k ]e − j ωk | < + ∞ for all ω so that the con v olutional in verse r inv [ k ] is w ell-deﬁned. The condition g 0 [ k ] = ( d ∗ g )[ k ] is b oth necessary and suﬃcien t only for our preferred choice of the scaling function ϕ 0 ( x ) = ( β α τ +1 / 2 ∗ ϕ 0 )( x ). This particular c hoice of the scaling function against the more direct choice H ϕ ( x ) is justiﬁed on the following grounds: • The function ϕ 0 ( x ) is well-localized with b etter deca y properties than H ϕ ( x ); the latter is not ev en integrable in general (e.g., the Harr scaling function), • The scaling function ϕ 0 ( x ) satisﬁes the partition-of-unity requirement, whereas H ϕ ( x ) is not a v alid scaling function since d H ϕ (0) is not nec- essarily unity . F or example, if ϕ ( x ) is symmetric and H ϕ ( x ) is in tegrable, then we ha ve d H ϕ (0) = R ( H ϕ )( x ) dx = 0 following the fact that H ϕ ( x ) is an ti-symmetric. 11 4.2 HT P airs of Biorthogonal W av elets A biorthogonal wa velet basis of L 2 ( R ), corresp onding to the dual-primal scaling function pair ( ϕ , ˜ ϕ ) of order ( N + 1 , ˜ N + 1), in volv es the nested multiresolution { 0 } ⊂ · · · ⊂ V − 1 ⊂ V 0 ⊂ V 1 ⊂ · · · ⊂ L 2 ( R ) , and its dual { 0 } ⊂ · · · ⊂ ˜ V − 1 ⊂ ˜ V 0 ⊂ ˜ V 1 ⊂ · · · ⊂ L 2 ( R ) , where the approximation subspace V j (resp. ˜ V j ) is generated b y the translations of ϕ j, 0 ( x ) (resp. ˜ ϕ j, 0 ( x )) [21]. Let ( ψ , ˜ ψ ) b e the wa v elets associated with these m ultiresolutions, which, along with their dilated-translated copies, enco de the residual signal—the diﬀerence of the signal appro ximations in successive subspaces. In particular, the wa v elet ψ j, 0 ( x ) (resp. ˜ ψ j, 0 ( x )) and its translates span the complemen tary space W j = V j  V j − 1 (resp. ˜ W j = ˜ V j  ˜ V j − 1 ). The crucial aspect of the construction is that the dilated-translated ensemble ψ j,k ( x ) and ψ j 0 ,k 0 ( x ) form a dual basis of L 2 ( R ), i.e., they satisfy the biorthogonality criteria h ψ j,k , ˜ ψ j 0 ,k 0 i = δ [ j − j 0 , k − k 0 ]. The expansion of a ﬁnite-energy signal f ( x ) in terms of this biorthogonal basis is then given by f ( x ) = X ( j,k ) ∈ Z 2 h f , ˜ ψ j,k i ψ j,k ( x ) . (21) In other w ords, the w a velets { ˜ ψ j,k ( x ) } and { ψ j,k ( x ) } , in terpreted as the analysis and syn thesis wa velets resp ectiv ely , together constitute a biorthognal wa velet basis of L 2 ( R ). In particular, let ˜ ϕ ( x ) and ϕ ( x ) b e the scaling functions, of order ˜ N + 1 and N + 1 resp ectively , asso ciated with a given biorthogonal wa v elet basis, with asso ciated w av elets ˜ ψ ( x/ 2) = X ˜ g [ k ] ˜ ϕ ( x − k ) , ψ ( x/ 2) = X g [ k ] ϕ ( x − k ) . No w, let ˜ ϕ ( x ) = ( β ˜ N ˜ τ ∗ ˜ ϕ 0 )( x ) and ϕ ( x ) = ( β N τ ∗ ϕ 0 )( x ) b e the resp ectiv e fac- torizations of ˜ ϕ ( x ) and ϕ ( x ). Consider the scaling functions ˜ ϕ 0 ( x ) = ( β ˜ N ˜ τ +1 / 2 ∗ ˜ ϕ 0 )( x ) and ϕ 0 ( x ) = ( β N τ +1 / 2 ∗ ϕ 0 )( x ), with asso ciated w av elets sp eciﬁed by ˜ ψ 0 ( x/ 2) = X ˜ g 0 [ k ] ˜ ϕ 0 ( x − k ) , ψ 0 ( x/ 2) = X g 0 [ k ] ϕ 0 ( x − k ) . Then the following result comes as a direct consequence of Theorem (4.2). Corollary 4.3 ( HT Pair of Biorthogonal W av elets ) The fol lowing ar e e quiva- lent: • The primal and dual wavelets form HT p airs, ˜ ψ 0 ( x ) = H ˜ ψ ( x ) and ψ 0 ( x ) = H ψ ( x ) , and { ˜ ψ 0 j,k ( x ) } and { ψ 0 j 0 ,k 0 ( x ) } to gether c onstitute a biortho gonal wavelet b asis of L 2 ( R ) . 12 • The discr ete HT c orr esp ondenc es ˜ g 0 [ k ] = ( d ∗ ˜ g )[ k ] and g 0 [ k ] = ( d ∗ g )[ k ] hold. The ab o ve construction also exhibits the following prop erties: • The tw o biorthogonal systems ha ve the same order, and the same Riesz b ounds. • If the pair ( ˜ ϕ , ϕ ) satisfy the biorthogonality relation, then so do ( ˜ ϕ 0 , ϕ 0 ). Indeed, using the identit y H β α τ +1 / 2 ( x ) = − ∆ 0 1 / 2 β α τ ( x ), w e can express the inner-pro duct h ˜ ϕ 0 , ϕ 0 ( · − k ) i as h β ˜ N ˜ τ +1 / 2 ∗ ˜ ϕ 0 , ( β N τ +1 / 2 ∗ ϕ 0 )( · − k ) i = hH  β ˜ N ˜ τ +1 / 2 ∗ ˜ ϕ 0  , H  β N τ +1 / 2 ∗ ϕ 0  ( · − k ) i = h− ∆ 0 1 / 2  β ˜ N ˜ τ ∗ ˜ ϕ 0  , − ∆ 0 1 / 2  β N τ ∗ ϕ 0  ( · − k ) i = h ˜ ϕ , ϕ ( · − k ) i , whic h establishes the assertion. • The lowpass ﬁlters on b oth the analysis and syn thesis side are “half- sample” shifted versions of one another, and are related via the mo dulation of the discrete HT ﬁlter: ˜ H ( z − 1 ) = D ( − z − 1 ) ˜ H 0 ( z − 1 ) , H 0 ( z ) = D ( − z − 1 ) H ( z ) . In particular, the ﬁlter is “half-sample” delay ed on the analysis side, whereas on the synthesis side the ﬁlter has a “half-sample” adv ance. • The highpass ﬁlters on b oth the analysis and syn thesis side are related through the FD ﬁlter as ˜ G ( z − 1 ) = D ( z ) ˜ G 0 ( z − 1 ) , G 0 ( z ) = D ( z ) G ( z ) . • If the analysis and synthesis ﬁlters of the original biorthogonal system satisfy the PR conditions G ( z − 1 ) ˜ G ( z ) + H ( z − 1 ) ˜ H ( z ) = 1 , G ( z − 1 ) ˜ G ( − z ) + H ( z − 1 ) ˜ H ( − z ) = 0 , then so do the ﬁlters of the HT pair. Indeed, since D ( z ) D ( z − 1 ) = 1, w e ha ve G 0 ( z − 1 ) ˜ G 0 ( z ) + H 0 ( z − 1 ) ˜ H 0 ( z ) = D ( z − 1 ) D ( z ) G ( z − 1 ) ˜ G 0 ( z ) + D ( − z ) D ( − z − 1 ) H ( z − 1 ) ˜ H ( z ) = G ( z − 1 ) ˜ G ( z ) + H ( z − 1 ) ˜ H ( z ) . 13 Similarly , G 0 ( z − 1 ) ˜ G 0 ( − z ) + H 0 ( z − 1 ) ˜ H 0 ( − z ) = 0 . Note that ab o v e prop erties relate to a common theme: the unitary nature of the op erators H and ∆ inv olv ed in the wa velet and the ﬁlterbank construction, resp ectiv ely . 5 1 D IMPLEMENT A TION Signal Pr e-ﬁltering : In order to implement the DT- C WT, w e need to emplo y t wo parallel wa velet decompositions corresp onding to the w av elets ψ ( x ) and ψ 0 ( x ). Moreov er, to hav e a coherent signal analysis—same input applied to b oth wa velet branches—w e need to pro ject the input signal f ( x ) separately on to V ( ϕ ) and V ( ϕ 0 ) before applying the resp ective DWTs. In particular, giv en a ﬁnite-energy input signal f ( x ), we consider its orthogonal pro jection f 0 ( x ) = P c 0 [ k ] ϕ ( x − k ) on to the space V ( ϕ ). The J -lev el wa v elet decomp osition of the signal f 0 ( x ) is subsequently given by f 0 ( x ) = X k ∈ Z c J [ k ] ϕ J,k ( x ) + X 1 6 j 6 J, k ∈ Z d j [ k ] ψ j,k ( x ) , (22) where the w av elet co eﬃcien ts d j [ k ], and the coarse appro ximation co eﬃcien ts c J [ k ] are recursively deriv ed from the pro jection co eﬃcien ts c 0 [ k ] using Mallat’s ﬁlterbank algorithm [22]. Ho wev er, in practice one has access only to the discrete samples of the input signal f ( x ); let { f [ k ] } k ∈ Z b e such (uniform) signal samples. It turns out that b y assuming the input signal f ( x ) to bandlimited, a particularly simple digital ﬁltering algorithm for computing the pro jection co eﬃcien ts is obtained: c 0 [ k ] = ( f ∗ p )[ k ] , (23) where the frequency resp onse P (e j ω ) of the digital ﬁlter p [ k ] is uniquely sp eciﬁed b y the restriction P (e j ω ) = b ˚ ϕ ( ω ) for ω ∈ ( − π , π ) (deriv ation details in § 11.3). As for the second branc h, the input signal is pro jected onto the corresp ond- ing appro ximation space V ( ϕ 0 ): the same t yp e of pre-ﬁltering is applied with an appropriate modiﬁcation of the frequency response, i.e., c ˚ ϕ 0 ( ω ) is used instead of b ˚ ϕ ( ω ). T o implemen t the ﬁlters for ﬁnite input signals, we use a FFT-based algorithm, similar to the one used in [3] for implemen ting the DWT ﬁlters. A nalysis & R e c onstruction : T o simplify the notation, w e shall henceforth use matrix notation to represen t the linear transformations asso ciated with the discrete DT- C WT. F or instance, corresp onding to an input signal f ∈ R N , the least-square pro jections are speciﬁed by c L 0 = Pf and c 0 L 0 = P 0 f , where P and P 0 are the N × N circulant matrices corresp onding to the tw o pre-ﬁlters. Let ( ˜ h, ˜ g , h, g ) and ( ˜ h 0 , ˜ g 0 , h 0 , g 0 ) b e the set of p erfect-reconstruction ﬁlters asso ciated with the biorthogonal systems ( ˜ ϕ , ˜ ψ , ϕ , ψ ) and ( ˜ ϕ 0 , ˜ ψ 0 , ϕ 0 , ψ 0 ) re- sp ectiv ely . The lo wpass { ( c L i , c 0 L i ) } and the highpass { ( c H i , c 0 H i ) } subbands at 14 successiv e lev els i = 1 , . . . , J are then giv en b y the recursiv e ﬁlterbank decom- p ositions: c L i = F ˜ h c L i − 1 , c H i = F ˜ g c L i − 1 c 0 L i = F ˜ h 0 c 0 L i − 1 , c 0 H i = F ˜ g 0 c 0 L i − 1 , (24) where F ˜ h and F ˜ g (resp. F ˜ h 0 and F ˜ g 0 ) denote the comp osition of the do wn- sampling matrix and the D WT matrix representing the lo wpass and highpass analysis ﬁlters of the ﬁrst (resp. second) channel. The complex wa velet sub- bands w 1 , . . . , w J are then sp eciﬁed by w i = c H i + j c 0 H i . In fact, the analysis can b e summarized b y the single frame op eration T : f 7→  c L J , c 0 L J , w 1 , . . . , w J  . (25) from a lo wer-dimensional space to a higher-dimensional space: dim( T f ) > dim( f ). In sev eral signal processing applications (e.g., denoising) one also needs to p erform an in verse transform, that is, reconstruct the denoised signal from the pro cessed complex wa velet coeﬃcients. Since T is realized through the concate- nation of bases, it is injective: T f = T ´ f only if f = ´ f ; how ever, as result of the redundancy , T exhibits non-unique left-inv erses. In our case, we use a simple left-in verse: T † :  c L J , c 0 L J , w 1 , . . . , w J  7→ f = 1 2 ( P − 1 c L 0 + P 0− 1 c 0 L 0 ) , where ( c L 0 , c 0 L 0 ) are obtained via the recursion c L i = F h c L i +1 + F g Re ( w i +1 ) , c 0 L i = F h 0 c 0 L i +1 + F g 0 Im ( w i +1 ) , (26) for i = J − 1 , . . . , 0 ( Re ( z ) and Im ( z ) denote the real and imaginary comp onen ts of z respectively). Here F h and F g (resp. F h 0 and F g 0 ) represent the comp osition of the DWT matrix corresp onding to the lowpass and highpass synthesis ﬁlters of the ﬁrst (resp. second) channel and the upsampling matrix. In short, the ab o ve in version op eration essentially amounts to inv erting the t wo parallel transforms and av eraging the inv erses. R emark : The role play ed by the tw o pro jection ﬁlters P (e j ω ) and P 0 (e j ω ) is critical as far as the issue of analyticity is concerned. Note that, while the analytic w av elet has an exact one-sided F ourier transform (b y construction), the corresp onding complex wa v elet ﬁlter ˜ G (e j ω ) + j ˜ G 0 (e j ω ) do es not inherit this prop ert y naturally; it is only the combination of the pro jection and wa velet ﬁlters, P a (e j ω ) = P (e j ω ) ˜ G (e j ω ) + j P 0 (e j ω ) ˜ G 0 (e j ω ) , that exhibits this prop ert y: P a (e j ω ) = 0 , for ω ∈ ( − π, 0]. Figure 1 sho ws the one-sided magnitude resp onse of the ﬁlter. 15 − 3 − 2 − 1 0 1 2 3 0 0.5 1 1.5 2 2.5 ANGULAR FREQUENCY (RADIANS) MAGNITUDE OF THE TRANSFER FUNCTION TRANSFER FUNCTION OF THE ANALYTIC FILTER TRANSFER FUNCTION OF THE REFERENCE FILTER Figure 1: T ransfer function of the analytic wa v elet ﬁlter P a (e j ω ). 6 GABOR-LIKE W A VELETS The “quantum law” for information—the principle that the join t time-frequency domain of signals is quan tized, and that the joint time-frequency supp ort of signals alwa ys exceed a certain minimal area—w as enunciated in signal theory b y Dennis Gab or [13]. He also identiﬁed the fact that the family of Gaussian- mo dulated complex exp onen tials (and their translates) provide the b est trade-oﬀ in the sense of Heisenberg’s uncertaint y principle. The canonical Gab or transform analyzes a signal using the set of “optimally- lo calized” Gab or atoms: g m,n ( x ) = 1 √ 2 π T 1 exp  − 1 2 T 2 1 ( x − mT )  e j n Ω( x − mT ) (27) generated via the mo dulations-translations of a Gaussian-mo dulated complex exp onen tial pulse [13, 1]. In particular, this paradigm in volv es the analysis of a ﬁnite-energy signal f ( x ) using the discrete sequence of pro jections c m,n = h f , g m,n i corresp onding to diﬀeren t mo dulations and translations ( m, n ) ∈ Z 2 . Note that the Gab or atoms hav e a ﬁxed size (sp eciﬁed b y the width T 1 of the Gaussian window), and hence the asso ciated transform essentially results in a “ﬁxed-windo w” analysis of the signal. Moreov er, the analysis functions { g m,n ( x ) } form a frame [21] and not a basis of L 2 ( R ); consequen tly , the recon- struction pro cess inv olving the dual frame is often computationally exp ensiv e and/or unstable [1]. 16 6.1 Analytic Gab or-lik e W av elets As a concrete application of the ideas developed in § 4, w e no w construct a family of analytic spline w av elets that asymptotically conv erge to Gab or-lik e functions. In particular, w e consider the family of semi-orthogonal B-spline w av elets that are b etter lo calized in space than their orthonormal counterparts, and that exhibit remark able joint time-frequency lo calization prop erties [34]. In particular, consider the multiresolution in § 2.2, generated by the fractional B-spline β α τ ( x ). The transfer function of the wa velet ﬁlter that generates the so-called B-spline w av elet [35] associated with this multiresolution is speciﬁed b y G α τ (e j ω ) = e j ω A α ( − e j ω ) H α τ ( − e − j ω ) . (28) W e denote the wa velet by ψ α τ ( x ). The dual multiresolution (resp. w av elet) is sp eciﬁed by the unique dual-spline function ˚ β α τ ( x ) (resp. dual wa velet, denoted b y ˜ ψ α τ ( x )). F ollowing Corollary 4.3, it can b e shown (pro of pro vided in § 11.2) that the family of B-spline w av elets { ψ α τ ( x ) } τ ∈ R , of a ﬁxed order α , and their duals are closed with resp ect to the HT: Prop osition 6.1 (HT Pair of B-spline W av elets) The HT of a B-spline (r esp. dual-spline) wavelet is a B-spline (r esp. dual-spline) wavelet of same or der, but with a diﬀer ent shift: H ψ α τ ( x ) = ψ α τ +1 / 2 ( x ) , H ˜ ψ α τ ( x ) = ˜ ψ α τ +1 / 2 ( x ) . (29) The imp ortance of this result is that it allows us to identify the analytic B-spline wa velet of degree α and shift τ : Ψ α τ ( x ) = ψ α τ ( x ) + j ψ α τ +1 / 2 ( x ) . (30) In the sequel (cf. § 7.3), w e shall make particular use of this analytic spline w av elet. W e would, ho wev er, like to highligh t a diﬀerent aspect: the remark- able fact that the w av elet Ψ α τ ( x ) resembles the celebrated Gab or functions for suﬃcien tly large α . Indeed, it was shown in [34] that the B-spline wa v elets asymptotically conv erge to the real part of the Gabor function; by appropri- ately mo difying the pro of in [34], the follo wing asymptotic conv ergence can established: ψ α τ ( x ) ∼ α → + ∞ M exp  − ( x − 1 / 2) 2 2 σ 2  cos  ω 0 x − ω 0 2 − π τ  . (31) where M = 2 M α +1 0 ∆ ω 0 / p 2 π ( α + 1); σ = √ α + 1 / ∆ ω 0 , with M 0 = 0 . 670 , ω 0 = − 5 . 142 and ∆ ω 0 = 2 . 670. W e recall that the asymptotic notation f α ( x ) ∼ g α ( x ) signiﬁes that f α ( x ) /g α ( x ) → 1 as α → + ∞ for all x . Immediately , we hav e the follo wing result: 17 Prop osition 6.2 (Gab or-lik e W av elet) The c omplex B-spline wavelet Ψ α τ ( x ) r esembles the Gab or function for suﬃciently lar ge α : Ψ α τ ( x ) ∼ M exp  − ( x − 1 / 2) 2 2 σ 2  e j ( ω 0 x − ω 0 2 − π τ ) . (32) The ab o ve conv ergence happ ens quite rapidly . F or instance, w e hav e ob- serv ed that the joint time-frequency resolution of the complex cubic B-spline w av elet ( α = 3) is already within 3% of the limit speciﬁed b y the uncertain ty principle. Fig. 2 depicts the complex w av elets generated using HT pair of B- spline wa v elets; the wa velets b ecomes more Gab or-lik e as the degree increases. Also sho wn in the ﬁgure is the magnitude env elop e | Ψ α τ ( x ) | of the complex w av elet which closely resembles the well-localized Gaussian windo w of the Ga- b or function. F rom a practical viewp oin t, this means that one could use the non-redundan t and n umerically stable multiresolution spline transforms to ap- pro ximate the Gab or analysis. R emark : While the B-spline wa v elets tend to be optimally lo calized in space, w e ha ve already observ ed that they are not orthogonal to their translates. The reconstruction therefore requires the use of some complementary dual functions. The ﬂip side is that these dual-spline wa v elets ha ve a comparativ ely po or spatial lo calization, that deteriorates as the degree increases. This is evident in Fig. 3, whic h shows quadrature pairs of such w av elets of diﬀerent degrees. How ev er, w e should emphasize that the dual (syn thesis) wa v elets hav e the same math- ematical rate of deca y as their analysis counterpart, and that the associated reconstruction algorithm is fast and numerically stable. 6.2 Gab or-lik e T ransform The dual-tree Gab or-lik e transform is based on the analytic B-spline w av elet Ψ α 0 ( x ) = ψ α 0 ( x ) + j ψ α 1 / 2 ( x ), where the degree α is suﬃciently large (the c hoice τ = 0 is arbitrary). The analysis and syn thesis ﬁlters for the ﬁrst and second c hannel are as sp eciﬁed b elo w [35]: • First channel: ˜ H ( z ) = 2 − ( α +1) (1 + z ) α +1 2 (1 + z − 1 ) α +1 2 , ˜ G ( z ) = z A α ( − z ) H ( − z − 1 ) , H ( z ) = ˜ H ( z ) A α ( z ) / A α ( z 2 ) , G ( z ) = ˜ G ( z ) / A α ( z 2 ) A α ( − z ) . (33) 18 !"#$%&'() ! *$%+,&(-./&%&0*("1(2&34&&5678 ( ( 9&.%(:.40 ;#.3+,.4<(:.40 ="2>%>* !"#$%&'() ! *$%+,&(-./&%&0*("1(2&34&&5678 (a) α =3 (b) α =6 Figure 2: HT pairs of B-spline wa velets. In either case, Blue (solid line): ψ α 0 ( x ), Red (broken line): ψ α 1 / 2 ( x ), Black (solid line): | ψ α 0 ( x ) + ψ α 1 / 2 ( x ) | 19 (a) α =3 (b) α =6 Figure 3: HT pairs of dual-spline w av elets. In either case, Blue (solid line): ψ α 0 ( x ), Red (broken line): ψ α 1 / 2 ( x ). • Second channel: ˜ H 0 ( z ) = 2 − ( α +1) (1 + z ) α 2 (1 + z − 1 ) α 2 +1 , ˜ G 0 ( z ) = z A α ( − z ) H 0 ( − z − 1 ) , H 0 ( z ) = ˜ H 0 ( z ) A α ( z ) / A α ( z 2 ) , G 0 ( z ) = ˜ G 0 ( z ) / A α ( z 2 ) α A ( − z ) . (34) The DT- C WT corresp onding to this Gab or-like w a velet w ould then result in the analysis of the input signal f ( x ) in terms of the sequence of multiscale pro jections h f , √ 2 m Ψ α 0 (2 m · − k ) i onto the (normalized) dilated-translated tem- plates of the Gab or-lik e wa velet Ψ α 0 ( x ). Note that here the Gab or-lik e wa velet is used for analysis, whereas its dual is used for synthesis. The corresp ond- ing DWTs ( F ˜ h , F ˜ g , F h , F g ) and ( F ˜ h 0 , F ˜ g 0 , F h 0 , F g 0 ) are eﬃcien tly implemen ted using a practical FFT-based algorithm, outlined in [3]. This metho d is exact despite the inﬁnite support of the underlying w av elets, and ac hieves p erfect- reconstruction up to a very high accuracy . The pre-ﬁlters, P and P 0 , are also implemen ted in a similar fashion. An added adv antage of the frequency domain implemen tation is that the execution time is indep endent of the order of the spatial ﬁlters. Moreov er, the ﬁlters in (33) need to be pre-computed once and for all in order to apply the transform to diﬀeren t signals (of a ﬁxed length). 20 7 BIV ARIA TE EXTENSION Next, based on ideas similar to those of Kingsbury [27], w e construct 2D com- plex wa velets, and 2D Gab or-lik e wa velets in particular, using a tensor-pro duct approac h. Moreov er, we also relate the real and imaginary components of the complex wa velets using a multi-dimensional extension of the HT. Sep ar able Biortho gonal Wavelet Basis : Biorthogonal wa velet bases of L 2 ( R ) can b e combined to construct a biorthogonal wa v elet basis of L 2 ( R 2 ). The underlying principle used to construct such a basis using tensor-pro ducts is as follo ws [21]: Theorem 7.1 L et ( ψ p , ˜ ψ p ) b e the primal and dual wavelets of a biortho gonal wavelet b asis of L 2 ( R ) , with c orr esp onding sc aling functions ( ϕ p , ˜ ϕ p ) . Simi- larly, let ( ψ q , ˜ ψ q ) c onstitute another biortho gonal wavelet b asis with c orr esp ond- ing sc aling functions ( ϕ q , ˜ ϕ q ) . Consider the fol lowing sep ar able wavelets and their duals ψ 1 ( x ) = ϕ p ( x ) ψ q ( y ) , ˜ ψ 1 ( x ) = ˜ ϕ p ( x ) ˜ ψ q ( y ) , ψ 2 ( x ) = ψ p ( x ) ϕ q ( y ) , ˜ ψ 2 ( x ) = ˜ ψ p ( x ) ˜ ϕ q ( y ) , ψ 3 ( x ) = ψ p ( x ) ψ q ( y ) . ˜ ψ 3 ( x ) = ˜ ψ p ( x ) ˜ ψ q ( y ) . (35) Then the dilation-tr anslations of ( ψ 1 ( x ) , ψ 2 ( x ) , ψ 3 ( x )) and ( ˜ ψ 1 ( x ) , ˜ ψ 2 ( x ) , ˜ ψ 3 ( x )) to gether c onstitute a biortho gonal wavelet b asis of L 2 ( R 2 ) . The functions ψ 1 ( x ) , ψ 2 ( x ) and ψ 3 ( x ) are p opularly referred to as the ‘low-high’ (LH), ‘high-low’ (HL) and ‘high-high’ (HH) w av elets, resp ectively , to emphasize the directions along whic h the lowpass scaling function and the highpass wa v elet op erate (here x = ( x, y ) denote the spatial coordinates). Note that the primal and dual appro ximation spaces for the ab ov e construction are V ( ϕ p ) ⊗ V ( ϕ q ) and V ( ˜ ϕ p ) ⊗ V ( ˜ ϕ q ) respectively , where V ( ϕ p ) ⊗ V ( ϕ q ) denotes the subspace span { ϕ p ( · − m ) ϕ q ( · − n ) } ( m,n ) ∈ Z 2 . 7.1 W av elet Construction A drawbac k of 2D separable wa v elets is their preferen tial resp onse to horizontal and v ertical features. Fig. 4 sho ws the three separable w av elets arising from the separable construction. The pulsation of the LH and HL w a velets are orien ted along the directions along which the constituent 1D wa v elets op erate. How ever, the HH wa velet, with its constituent 1D wa velets op erating along orthogonal directions, do es not exhibit orientation purely along one direction; instead it sho ws a chec kerboard app earance with simultaneous pulsation along the diago- nal directions. This is exactly where the analytic wa velet ψ a = ψ + j H { ψ } , with its one- sided frequency spectrum comes to the rescue: if instead of employing separable 21 w av elets of the form ψ ( x ) ψ ( y ), complex w av elets of the form ψ a ( x ) ψ a ( y ) are used, then the corresponding spectrum b ψ a ( ω x ) b ψ a ( ω y ) will ha ve only one pass- band, and consequently the real w av elets Re ( ψ a ( x ) ψ a ( y )) and Im ( ψ a ( x ) ψ a ( y )) will indeed b e orien ted. The motiv ation then is to use HT pairs of 1D biorthogonal wa v elets to con- struct orien ted 2D wa v elets. In particular, w e do so b y appropriately combin- ing four separable biorthogonal wa v elet bases using Theorem (7.1). T o b egin with, we immediately identify the tw o scaling functions ϕ p ( x ) = ϕ ( x ) and ϕ q ( x ) = ϕ 0 ( x ), associated with the analytic w av elet ψ a ( x ) = ψ ( x ) + j ψ 0 ( x ), where ψ 0 ( x ) = H ψ ( x ). This naturally leads to the p ossibilit y of four separable biorthogonal wa velet bases corresponding to the following p ossible choices of appro ximation spaces: V ( ϕ ) ⊗ V ( ϕ ) , V ( ϕ ) ⊗ V ( ϕ 0 ) , V ( ϕ 0 ) ⊗ V ( ϕ ) and V ( ϕ 0 ) ⊗ V ( ϕ 0 ). In fact, as will b e demonstrated shortly , w e will employ all of these to obtain a balanced construction. First, w e iden tify the separable wa velets corresponding to the four scaling spaces: ψ 1 ( x ) = ϕ ( x ) ψ ( y ) , ψ 4 ( x ) = ϕ ( x ) ψ 0 ( y ) , ψ 2 ( x ) = ψ ( x ) ϕ ( y ) , ψ 5 ( x ) = ψ ( x ) ϕ 0 ( y ) , ψ 3 ( x ) = ψ ( x ) ψ ( y ) , ψ 6 ( x ) = ψ ( x ) ψ 0 ( y ) , ψ 7 ( x ) = ϕ 0 ( x ) ψ ( y ) , ψ 10 ( x ) = ϕ 0 ( x ) ψ 0 ( y ) , ψ 8 ( x ) = ψ 0 ( x ) ϕ ( y ) , ψ 11 ( x ) = ψ 0 ( x ) ϕ 0 ( y ) , ψ 9 ( x ) = ψ 0 ( x ) ψ ( y ) , ψ 12 ( x ) = ψ 0 ( x ) ψ 0 ( y ) . (36) The corresp onding dual w av elets ˜ ψ ` ( x ) are sp eciﬁed iden tically except that the dual wa velets are used instead of the primal ones. Finally , by judiciously using the one-sided sp ectrum of the analytic wa velet ψ a ( x ) = ψ ( x ) + j ψ 0 ( x ), and b y com bining the four separable w av elet bases (36), we arrive at the follo wing w av elet sp eciﬁcations: Ψ 1 ( x ) = ψ a ( x ) ϕ ( y ) = ψ 2 ( x ) + j ψ 8 ( x ) , Ψ 2 ( x ) = ψ a ( x ) ϕ 0 ( y ) = ψ 5 ( x ) + j ψ 11 ( x ) , Ψ 3 ( x ) = ϕ ( x ) ψ a ( y ) = ψ 1 ( x ) + j ψ 4 ( x ) , Ψ 4 ( x ) = ϕ 0 ( x ) ψ a ( y ) = ψ 7 ( x ) + j ψ 10 ( x ) , Ψ 5 ( x ) = 1 √ 2 ψ a ( x ) ψ a ( y ) =  ψ 3 ( x ) − ψ 12 ( x ) √ 2  + j  ψ 6 ( x ) + ψ 9 ( x ) √ 2  , Ψ 6 ( x ) = 1 √ 2 ψ ∗ a ( x ) ψ a ( y ) =  ψ 3 ( x ) + ψ 12 ( x ) √ 2  + j  ψ 6 ( x ) − ψ 9 ( x ) √ 2  . (37) 22 D WT : Check er -Boar d Eff ect Figure 4: W a velets associated with the separable basis. The ﬁgure sho ws the LH, HL and HH wa v elets in the space domain. The dual complex w av elets, ˜ Ψ k ( x ), are speciﬁed in an identical fashion using the dual wa v elets ˜ ψ ` ( x ). Imp ortan tly , the ab o ve construction is c omplete in the sense that it inv olves all the 4 × 3 = 12 separable w av elets of the four parallel m ultiresolutions. The factor 1 / √ 2 ensures normalization: the real and imaginary comp onen ts of the six complex wa v elets hav e the same norm. 7.2 Directional Selectivit y and Shift-In v ariance A real wa v elet has a bandpass sp ectrum that is symmetric w.r.t. to the origin. As a result, for the w av elet to b e oriented, it is necessary that its sp ectrum b e bandpass only along one preferen tial direction. W e claim that the real and imaginary components of the ab o ve complex wa v elets are oriented along the primal directions θ 1 = θ 2 = 0, θ 3 = θ 4 = π/ 2 , θ 5 = π/ 4, and θ 6 = 3 π/ 4, resp ectiv ely . Indeed, it is easily seen that the supp ort of b Ψ 1 ( ω ) and b Ψ 2 ( ω ) is restricted to the half-plane { ( ω x , ω y ) : ω x > 0 } , since their F ourier transform can b e written as b Ψ k ( ω ) = (1 + sign( ω x )) \ Re (Ψ k )( ω ). As it is necessary for the real functions Re (Ψ k ) and Im (Ψ k ) to hav e symmetric passbands, the claim ab out their orientation along the horizontal direction then follo ws immediately . The orientation of the comp onen ts of the wa velets Ψ 3 ( x ) and Ψ 4 ( x ) along the v ertical direction follows from a similar argumen t. As far as the wa velet Ψ 5 ( x ) is concerned, note that b Ψ 5 ( ω ) = (1 + sign( ω x )) (1 + sign( ω y )) b ψ ( ω x ) b ψ ( ω y ). As a consequence, the supp ort of b Ψ 5 ( ω ) is restricted to the quadrant { ( ω x , ω y ) : ω x > 0 , ω y > 0 } . The symmetry requirements on the spectrums of Re (Ψ 5 )( x ) and Im (Ψ 5 )( x ) then establish their orientation along π / 4. A similar argument establishes the orientation of the real comp onen ts Re (Ψ 6 )( x ) and Im (Ψ 6 )( x ) along 3 π / 4. The ab ov e-men tioned directional prop erties allude to some kind of analytic c haracterization of the complex wa velets. Indeed, akin to the 1D counterpart, it turns out that the comp onen ts of the ab o v e complex wa velets can also b e re- lated via a multi-dimensional extension of the HT that provides further insights in to the orientations of the wa velets. In particular, w e consider the following directional version of the HT [14]: H θ f ( x ) F ← → − j sign( ω T u θ ) ˆ f ( ω ) , (38) sp eciﬁed by the unit vector u θ = (cos θ , sin θ ) pointing in the direction 0 6 23 Figure 5: 2D Gab or-lik e W av elets. Left: Real comp onen t of the six complex w av elets, Right: Magnitude en velope of the six complex wa velets. The diago- nally placed wa v elets are identical, they are used twice to balance the represen- tation. θ < π . That is, the directional HT is p erformed with resp ect to the half-spaces { ω : ω T u θ > 0 } and { ω : ω T u θ < 0 } sp eciﬁed by the v ector u θ , and it maps the directional cosine cos( ω T θ x ) in to the directional sine sin( ω T θ x ). Based on the w av elet deﬁnitions (37), the following corresp ondences (for a pro of see § 11.4) can then b e deriv ed: Prop osition 7.2 The r e al and imaginary c omp onents of the c omplex wavelets Ψ k ( x ) form dir e ctional HT p airs. In p articular: Im (Ψ k ( x )) = H θ k Re (Ψ k ( x )) , 1 6 k 6 6 . (39) A signiﬁcant problem with the decimated DWT is that the critical down- sampling makes it shift-v ariant. The redundancy of the dual-tree transform has b een successfully exploited for partially mitigating this shift-v ariance problem [17, 19]. Our design further mitigates this shift-v ariance problem by using a ﬁner sub-sampling scheme in the 0 ◦ and 90 ◦ directions. Observe that Ψ 1 ( x ) ≈ Ψ 2 ( x, y − 1 / 2) owing to the fact that β α τ ( x ) ≈ β α τ +1 / 2 ( x + 1 / 2). These wa velets pro vide a ﬁner sampling in the y -direction. Similarly , the v ertical wa velets give us a ﬁner sampling in the x -direction. 7.3 Gab or-lik e W av elets Daugman generalized the Gab or function to the following 2D form G ( x ) = 1 2 π σ x σ y e −  ( x − x 0 ) 2 2 σ 2 x + ( y − y 0 ) 2 2 σ 2 y  e j ( ξ 0 x + ν 0 y ) , (40) 24 in volving the mo dulation of an elliptic Gaussian using a directional plane-w av e, to model the receptiv e ﬁelds of the orien tation-selective simple cells in the visual cortex [9]. W e are particularly interested in the dual-tree wa v elets, denoted by G 1 ( x ; α, τ ) , . . . , G 6 ( x ; α, τ ), deriv ed from the quadrature B-spline wa velets ψ p ( x ) = ψ α τ ( x ) and ψ q ( x ) = ψ α τ +1 / 2 ( x ). These complex wa v elets inherit the asymptotic prop erties of the constituen t spline functions. Indeed, by appropriately mo difying the pro of in [34], it can b e shown that β α τ ( x ) ∼ 1 √ 2 π σ exp  − ( x − τ ) 2 2 σ 2  (41) for suﬃcien tly large α , where σ = √ α + 1 / 2 √ 3. This, combined with (32), then results in the following asymptotic characterization: Prop osition 7.3 ( 2D Gabor-like W av elets ) The c omplex wavelets G k ( x ; α, τ ) r esemble the 2 D Gab or functions for suﬃciently lar ge α : G 1 ( x ; α, τ ) ∼ M 1 e − „ ( x − 1 / 2) 2 σ 2 1 + ( y − τ ) 2 σ 2 2 « e j ( ω 0 x − ω 0 2 − π τ ) , G 2 ( x ; α, τ ) ∼ M 1 e − „ ( x − 1 / 2) 2 σ 2 1 + ( y − τ − 1 / 2) 2 σ 2 2 « e j ( ω 0 x − ω 0 2 − π τ ) , G 3 ( x ; α, τ ) ∼ M 1 e − „ ( x − τ ) 2 σ 2 2 + ( y − 1 / 2) 2 σ 2 1 « e j ( ω 0 y − ω 0 2 − π τ ) , G 4 ( x ; α, τ ) ∼ M 1 e − „ ( x − τ − 1 / 2) 2 σ 2 2 + ( y − 1 / 2) 2 σ 2 1 « e j ( ω 0 y − ω 0 2 − π τ ) , G 5 ( x ; α, τ ) ∼ M 2 e − „ ( x − 1 / 2) 2 σ 2 1 + ( y − 1 / 2) 2 σ 2 1 « e j ( ω 0 ( x + y ) − ω 0 − 2 πτ ) , G 6 ( x ; α, τ ) ∼ M 2 e − „ ( x − 1 / 2) 2 σ 2 1 + ( y − 1 / 2) 2 σ 2 1 « e j ω 0 ( y − x ) , (42) wher e M 1 = 2 √ 3 M α +1 0 ∆ ω 0 /π ( α + 1) ; M 2 = 2 M 2( α +1) 0 ∆ ω 2 0 /π ( α + 1) ; σ 1 = √ α + 1 / ∆ ω 0 ; and σ 2 = p ( α + 1) / 6 . W e call the w a velets “Gab or-lik e” since they form approximates of 2D Gab or functions similar to the ones prop osed b y Daugman (40). The dual-tree trans- form (cf. § 8) corresp onding to a sp eciﬁc family of such Gab or-lik e wa v elets (ﬁxed α and τ ) results in a multiresolution, directional analysis of the input im- age f ( x ) in terms of the sequence of pro jections h f , 2 i G k (2 i x − m ; α, τ ) i . Fig. 5 sho ws the 2D Gab or wa velets corresp onding to α = 6 and τ = 0. The ensem ble sho ws the mo dulus | G k ( x ; 6 , 0) | and the real comp onen t Re ( G k ( x ; 6 , 0)) of the six complex wa v elets; the former shows the pulsations of the directional plane w av es, whereas the latter shows the elliptical Gaussian env elop es. 25 7.4 Discussion Before moving on to the implemen tation, we digress brieﬂy to discuss certain k ey asp ects of our construction: • Dir e ctionality : The six complex w av elets in Kingsbury’s DT- C WT sc heme are oriented along the directions: ± 15 ◦ , ± 45 ◦ , and ± 75 ◦ [19]. Though we use similar separable building blo c ks in our approac h, our w av elets are orien ted along the four principal directions: 0 , π / 4 , π / 2 and 3 π / 4. The added redundancy along the horizontal and v ertical directions yields b etter shift-in v ariance along these directions. Alternativ ely , we could also ha ve applied Kinsbury’s construction to obtain Gab or-lik e wa v elets orientated along ± 15 ◦ , ± 45 ◦ , and ± 75 ◦ . • L o c alization Vs. F r ame Bounds : In this pap er, we placed emphasis on time-frequency lo calization, and were able to construct new wa velets that con verge to Gab or-lik e functions. These basis functions should pro ve use- ful for image analysis tasks such as extraction of AM-FM information and texture analysis. How ev er, the price to pa y for this improv ed lo cal- ization is that the asso ciated transform—in contrast with the transforms constructed by Kingsbury et al. [27, 18]—is no longer tight, and con- sequen tly requires a diﬀerent set of reconstruction ﬁlters. Nevertheless, the tightness of the frame-bounds—a desirable property for image pro- cessing applications such as denoising and compression—can, in principle, also b e achiev ed within our prop osed framework by replacing the B-spline w av elets with the orthonormal ones (Battle-L ´ emarie wa velets). • Analytic Pr op erties : Our metho d of construction tak es a primary wa velet transform and obtains an exact HT pair us ing a simple unitary mapping. The consequence is that all fundamental appro ximation-theoretic prop er- ties of contin uous-domain wa v elets, such as v anishing momen ts and regu- larit y , are automatically preserv ed, and that the asso ciated ﬁlters inherit an exact one-sided resp onse. W e also obtain an explicit space-domain expression for the Gab or-lik e wa v elets. • Multidimensional HT pr op erties : The directional HT corresp ondences (39) for our complex w av elets follo ws as a direct consequence of the tensor- pro duct construction. W e w ould ho w ever like to note that there exist other m ultidimensional extensions of the HT as well: the “single-orthant” exten- sion of Hahn [15] in volving the b oundary distribution of analytic functions; the “h yp ercomplex” extension due to B ¨ ulo w et al. [6]; the “monogenic” signal due to F elsb erg et al. [11]; and the spiral-phase quadrature trans- form of Larkin et al. [20]. The last tw o in the list are closely related to the Riesz transform of classical harmonic analysis [32]. Design of directional w av elets based on these and other alternative extensions are a promising topic of research [7, 24]. 26 Complex W av elet Coefﬁcients Input Signal f P er mutations 2D D WT Projection Filters Orthonormal T ransforms P 1 P 2 P 3 P 4 F 4 F 3 F 2 F 1 w Π Λ R , Λ I Figure 6: Blo c k Diagram of the 2D Complex W av elet T ransform. 8 2 D IMPLEMENT A TION Pr e-ﬁltering : The input signal has to b e pro jected onto eac h of the four sep- arable spaces of the form V ( ϕ p ) ⊗ V ( ϕ q ) before initiating the m ultiresolution decomp ositions. As in the 1D setting, the orthogonal pro jection is achiev ed in a separable fashion using an appropriate pre-ﬁlter along eac h dimension. In par- ticular, if { f [ k ] } k ∈ Z 2 b e the uniform samples of a bandlimited input signal f ( x ), then the pro jection co eﬃcients are given by c LL 0 [ k ] = ( f ∗ p )[ k ], where the sep- arable pre-ﬁlter p [ k ] is speciﬁed by P p [ k 1 , k 2 ]e − j ( k 1 ω x + k 2 ω y ) = b ˚ ϕ r ( ω x ) b ˚ ϕ s ( ω y ) for ω = ( ω x , ω y ) in ( − π , π ) 2 . In general, there would b e four such pro jections c LL 0 ( n ) = f ∗ p n , 1 6 n 6 4, corresp onding to the 2D pre-ﬁlters p 1 , . . . , p 4 asso ciated with the four appro x- imation spaces. Note that the ﬁlters can b e implemen ted eﬃcien tly through successiv e 1D ﬁltering along either dimension. A nalysis : W e consider the implementation asp ects for a ﬁnite input signal f ∈ R M × N . The transform, corresponding to the complex wa velets (37), in volv es four separable DWTs with diﬀerent ﬁlters applied along the x and y directions (cf. T able 8 for the list of ﬁlters), and result in four subbands at each decom- p osition level. Sp eciﬁcally , let c LL i ( n ) , c LH i ( n ) , c H L i ( n ) and c H H i ( n ) , 1 6 n 6 4 , denote the low-lo w, lo w-high, high-low and high-high subbands, resp ectiv ely , of the four D WT decomp ositions at resolution i = 1 , . . . , J . The low-lo w subbands c LL 0 ( n ) are identiﬁed as the four set of pre-ﬁltered signals c LL 0 ( n ) = P n f , with P n b eing the (blo c k) circulant matrices asso ciated with the 2D pre-ﬁlters. The 27 coarser subbands at levels i = 1 , . . . , J are then given by c LL i ( n ) = F n ( ˜ h x , ˜ h y ) c LL i − 1 ( n ) c LH i ( n ) = F n ( ˜ h x , ˜ g y ) c LL i − 1 ( n ) c H L i ( n ) = F n ( ˜ g x , ˜ h y ) c LL i − 1 ( n ) c H H i ( n ) = F n ( ˜ g x , ˜ g y ) c LL i − 1 ( n ) , (43) where F n ( q x , q y ) denotes the comp osition of the n th D WT matrix (employing analysis ﬁlters q x and q y in the x -direction and y -direction), and the do wnsam- pling matrix. The complex subbands w i = ( w 1 i , . . . , w 6 i ) , 1 6 i 6 J, are sp eciﬁed b y w i = Λ R ζ i + j Λ I ξ i , where ζ i =  c H L i (1) , c H L i (2) , c LH i (1) , c LH i (3) , c H H i (1) , c H H i (4)  , ξ i =  c H L i (3) , c H L i (4) , c LH i (2) , c LH i (4) , c H H i (2) , c H H i (3)  , (44) are obtained through a particular p erm utation of the 12 highpass subbands; and the blo c k matrices Λ R and Λ I are sp eciﬁed as Λ R = 1 √ 2         √ 2 I 0 0 0 0 0 0 √ 2 I 0 0 0 0 0 0 √ 2 I 0 0 0 0 0 0 √ 2 I 0 0 0 0 0 0 I − I 0 0 0 0 I I         , Λ I = 1 √ 2         √ 2 I 0 0 0 0 0 0 √ 2 I 0 0 0 0 0 0 √ 2 I 0 0 0 0 0 0 √ 2 I 0 0 0 0 0 0 I I 0 0 0 0 I − I         . (45) In short, the transform can b e formally summarized via the frame op eration T : f 7→ ( c J (1) , . . . , c J (4) , w 1 , . . . , w J ) (46) in volving the sequence of transformations: pro jections P 1 , . . . , P 4 ; discrete wa velet transforms F 1 , . . . , F 4 ; p erm utation Π, and orthonormal transformations Λ R and Λ I . Figure 6 provides a sc hematic of these sequence of transformations. R e c onstruction : Note that the p erm utation Π : { c LH i ( n ) , c H L i ( n ) , c H H i ( n ) } 1 6 n 6 4 7→ ( ζ i , ξ i ) (47) in volv ed in (44) is in vertible, and that the matrices Λ R and Λ I are orthonormal, with corresponding in v erses giv en b y Λ T R and Λ T I resp ectiv ely . Starting with the 28 Six Directional Subbands Input Images Figure 7: Directional decomp osition (one-level) of a syn thetic image (Octagon) and a natural image (Cameraman) using the Gab or-like transform. Ordering of the subbands in either case: First column : | w 1 | , | w 2 | ; Second column : | w 3 | , | w 4 | ; Third column : | w 5 | and | w 6 | . 29 complex w av elet subbands w 1 , . . . , w J , the highpass subbands c b i (1) , . . . , c b i (4) corresp onding to the bands b = H L, LH , and H H , are then computed from the v ectors ζ i = Λ T R Re ( w i ) and ξ i = Λ T I Im ( w i ), via the p erm utation Π − 1 at lev els i = 1 , . . . , J . These, along with the lowpass subbands c LL J (1) , . . . , c LL J (4), are then used to reconstruct the pro jected signals c LL 0 (1) , . . . , c LL 0 (4) using the recursion c LL i ( n ) = F n ( h x , h y ) c LL i +1 ( n ) + F n ( h x , g y ) c LH i +1 ( n ) , + F n ( g x , h y ) c H L i +1 ( n ) + F n ( g x , g y ) c H H i +1 ( n ) (48) for i = J − 1 , . . . , 0. Here, F n ( m x , m y ) represen ts the comp osition of the up- sampling matrix and the synthesis matrix corresponding to the n th D WT, with ﬁlters m x and m y in the x -direction and y -direction, respectively , as sp eciﬁed in T able 8. The input signal samples are ﬁnally reco vered as f = 1 / 4 P 4 n =1 P − 1 n c LL 0 ( n ). 2 D Gab or-like T r ansform : The Gab or-lik e transform is based on the ana- lytic B-spline w av elets sp eciﬁed in § 7.3, where the complex subbands w k i [ m ] represen t the directional decomp ositions of the input image along the four pri- mal directions using the optimally-lo calized Gab or-like wa v elets G k ( x ; α, τ ) at diﬀeren t resolutions. The ﬁlterbank analysis (43) and synthesis (48) op erations are implemented in a separable fashion using the 1D spline D WT ﬁlters sp eciﬁed in (33) and (34). Fig. 7 shows the magnitude resp onse of the six complex wa v elet subands obtained by applying our Gab or-like transform to a syn thetic and a natural im- age. In particular, the wa velet subbands corresp onding to the syn thetic image, with directional edges along 0 , π / 4 , π / 2 and 3 π/ 4, highligh t the directional- selectivit y of the transform. The simulation w as carried out in MA TLAB 7 . 5 on a Macin tosh 2 . 66 GHz Intel dual-core system. The a verage execution time for one-lev el wa velet analysis and reconstruction (including pre- and p ost-ﬁltering) of a 512 × 512 image is 1 . 2 seconds, and the reconstruction error is of the order of 10 − 16 . SUBMITTED T O IEEE TRANSA CTIONS ON SIGN AL PR OCESSING 24 T ABLE I 2 DC O M P L E X W A V E L E T T R A N S F O R M F I L T E R S Index ( n ) 2 D Scaling Spaces Analysis Filters Synthesis Filters x -direction y -direction x -direction y -direction 1 V ( ϕ ) ⊗ V ( ϕ ) e h, e g e h, e g h, g h, g 2 V ( ϕ ) ⊗ V ( ϕ ! ) e h, e g e h ! , e g ! h, g h ! ,g ! 3 V ( ϕ ! ) ⊗ V ( ϕ ) e h ! , e g ! e h, e g h ! ,g ! h, g 4 V ( ϕ ! ) ⊗ V ( ϕ ! ) e h ! , e g ! e h ! , e g ! h ! ,g ! h ! ,g ! V I I I . I M P L E M E N T A T I O N The implementation of the 2 D comple x w a v elet transform in v olv es four separable D WTs in parallel. The comple x w a v elet subbands at each stage of the decomposition are then obtained via simple linear combinations of the subbands of the four D WTs in accordance with the deﬁnitions of the six comple x w a v elets (34). A. Pr e-ﬁltering Similar to the 1 D case, the input signal needs to be projected onto each of the four separable space of the form V ( ϕ r ) ⊗ V ( ϕ s ) , before we be gin the multiresolution analysis. The orthogonal projection is achie v ed e xactly as in the 1 D case, e xcept that dif ferent projection ﬁlters are used along the tw o dimensions. Speciﬁcally , let { f [ k ] } k ∈ Z 2 be the input signal samples. Assuming a bandlimited model for the input signal, the projection coef ﬁcients c LL 0 [ k ] are then gi v en by c LL 0 = f ∗ p , where the pre-ﬁlter { p [ k ] } k ∈ Z 2 is speciﬁed by the DTFT ! p [ k , l ] e − j ( k ω x + l ω y ) = ! ( m,n ) ∈ Z 2 ν ( ω x +2 π m, ω y +2 π n ) , with ν ( ω x , ω y ) := " ˚ ϕ r ( ω x ) " ˚ ϕ s ( ω y ) , for ( ω x , ω y ) ∈ ( − π , π ] 2 and 0 else. As in the 1 D case, ˚ ϕ r and ˚ ϕ s are the (unique) duals of ϕ r and ϕ s , respecti v ely . In particular , corresponding to the each of the four scaling spaces, there w ould be four such projections c LL 0 ( n )= f ∗ p n , 1 ! n ! 4 , with p n being the corresponding 2 D pre-ﬁlters. Note that the abo v e 2 D ﬁltering can be ef ﬁciently implemented by successi v e 1 D ﬁltering operations along the tw o dimensions. B. Analysis Similar to the 1 D case, we consider the ﬁnite input signal f [ k 1 ,k 2 ] , 0 ! k 1 ! M − 1 , 0 ! k 2 ! N − 1 and the its corresponding v ectorial representation f := ( f [0 , 0] , . . . , f [ M − 1 ,N − 1]) . The projected signals are then gi v en by c LL 0 ( n )= P n f , 1 ! n ! 4 , with P n being the (block) circulant matrices corresponding to the 2 D pre-ﬁlters. Figure 8: Analysis and synthesis ﬁlters corresp onding to the four multiresolu- tions. 30 9 CONCLUDING REMARKS The primary ob jective of this con tribution was to com bine the attractive fea- tures of Gab or analyses and multiresolution wa velet transforms into a single theoretical framework, and to pro vide a fast algorithm for the same. Sp eciﬁ- cally , w e prop osed a formalism for constructing exact HT pairs of biorthogonal w av elets based on (i) the B-spline factorization theorem, and (ii) a natural dis- cretization of the contin uous HT ﬁlter identiﬁed via the action of the HT on fractional B-splines. Based on this methodology , analytic wa v elets resembling the Gab or function w ere then designed using HT pair of B-spline wa v elets. W e then extended our scheme to 2D: starting from HT pair of 1D biorthogo- nal wa v elet basis, w e constructed directional complex wa velets b y appropriately com bining four separable biorthogonal wa velet bases. In particular, we related the real and imaginary comp onents of the complex wa v elets using a directional extension of the HT. The particular family of wa velets constructed using B- splines was shown to resem ble the directional Gab or function family prop osed b y Daugman. Finally , we demonstrated how the discrete Gab or-like transforms could b e implemented using fast FFT-based ﬁlterbank algorithms. 10 A CKNO WLEDGEMENTS The authors would like to thank Dr. T. Blu for sharing his researc h ﬁndings (the asymptotic form (31) in particular) on fractional splines, and Dr. P . Th´ evenaz and Dr. C. S. Seelamantula for pro ofreading the man uscript. 11 APPENDIX 11.1 Pro of of Theorem 4.2: W e b egin with the following sequence of equiv alences H ψ ( x/ 2) = X k ∈ Z g [ k ] H ϕ ( x − k ) = X k ∈ Z g [ k ]( H β α τ ∗ ϕ 0 )( x − k ) = X k ∈ Z g [ k ] X n ∈ Z d [ n ]( β α τ +1 / 2 ∗ ϕ 0 )( · − n ) ! ( x − k ) = X m ∈ Z ( g ∗ d )[ m ] ϕ 0 ( x − m ) , based on (12), and the linearity , asso ciativit y , and commutativit y of the under- lying conv olution op erators. The suﬃciency part of the theorem then follo ws immediately: if g 0 [ k ] = ( g ∗ d )[ k ], then ψ 0 ( x ) = H ψ ( x ). 31 Con versely , let ψ 0 ( x ) = H ψ ( x ), so that P g 0 [ k ] ϕ 0 ( x − k ) = P ( g ∗ d )[ k ] ϕ 0 ( x − k ). Now, since { ϕ 0 ( · − n ) } forms a Riesz basis of the subspace V ( ϕ 0 ) = span ` 2 { ϕ 0 ( · − k ) } k ∈ Z , every element in V ( ϕ 0 ) necessarily has a unique rep- resen tation. Hence, g 0 [ k ] = ( g ∗ d )[ k ]. 11.2 Pro of of Prop osition 6.1: The primal scaling functions can b e trivially factorized: ϕ ( x ) = ( β α τ ∗ ϕ 0 )( x ) and ϕ 0 ( x ) = ( β α τ +1 / 2 ∗ ϕ 0 )( x ), where ϕ 0 is the Dirac delta distribution. Similarly , the dual scaling functions can b e factorized as ˜ ϕ ( x ) = ( β α τ ∗ ˜ ϕ 0 )( x ) and ˜ ϕ 0 ( x ) = ( β α τ +1 / 2 ∗ ˜ ϕ 0 )( x ), where ˜ ϕ 0 = P q α [ k ] δ ( · − k ) with P q α [ k ]e − j ωk = 1 / A α (e j ω ). Note that in the latter case we hav e particularly used the fact that A α (e j ω ), and hence q α [ k ], are indep enden t of τ . The prop osition then follows from Corollary (4.3) since the wa velet ﬁlters satisfy the suﬃciency conditions: ˜ g 0 [ k ] = ( d ∗ ˜ g )[ k ] and g 0 [ k ] = ( d ∗ g )[ k ], resp ec- tiv ely . Indeed, from (14) and (28), we hav e G 0 (e j ω ) = e j ω A α ( − e j ω ) H α τ +1 / 2 ( − e − j ω ) = e j ω A α ( − e j ω ) D (e j ω ) H α τ ( − e − j ω ) = D (e j ω ) G (e j ω ) . The other condition ˜ G 0 (e j ω ) = D (e j ω ) ˜ G (e j ω ) can b e similarly derived. 11.3 Deriv ation of Equation (23) : It is well-kno wn that the least-square appro ximation op erator P V ( ϕ ) : L 2 ( R ) → V ( ϕ ), deﬁned by P V ( ϕ ) f = arg min f 0 ∈ V ( ϕ ) k f − f 0 k (49) giv es the orthogonal pro jection of f ( x ) on to V ( ϕ ). The solution to the ab ov e problem is explicitly given by P V ( ϕ ) f ( x ) = P c 0 [ k ] ϕ ( x − k ) where the co eﬃ- cien ts are sp eciﬁed by c 0 [ k ] = h f , ˚ ϕ ( · − k ) i . Here ˚ ϕ ( x ) denotes the dual of ϕ ( x ) that satisﬁes the biorthogonalit y criterion h ϕ , ˚ ϕ ( · − n ) i = δ [ n ]. Moreov er, under the constrain t that ˚ ϕ ( x ) ∈ V ( ϕ ), w e recov er a unique dual that is sp eciﬁed b y the F ourier transform b ˚ ϕ ( ω ) = b ϕ ( ω ) / P | b ϕ ( ω + 2 π k ) | 2 [33]. Next, using the Poisson summation form ula, w e derive the expression C 0 (e j ω ) = P n ∈ Z \ ( f ∗ ˚ ϕ T )( ω + 2 π n ) for the (discrete) F ourier transform of c 0 [ k ]. The ban- dlimited mo del f ( x ) = P f [ k ] sinc( x − k ) ﬁnally results in the simpliﬁcation C 0 (e j ω ) = X n ∈ Z \ ( f ∗ ˚ ϕ T )( ω + 2 π n ) = F (e j ω ) X n ∈ Z rect  ω + 2 π n 2 π  b ˚ ϕ ( ω + 2 π n ) = F (e j ω ) P (e j ω ) , (50) 32 where P (e j ω ) equals b ˚ ϕ ( ω ) on ( − π , π ), and F (e j ω ) is the F ourier transform of f [ k ]. 11.4 Pro of of Prop osition 7.2: W e establish the corresp ondence for the w av elets Ψ 1 ( x ) and Ψ 5 ( x ) (the rest can b e derived similarly). The corresp ondence for the former is direct: H 0 Re (Ψ 1 ( x )) = H 0 { ψ ( x ) } ϕ ( y ) = ψ 0 ( x ) ϕ ( y ) = Im (Ψ 1 ( x )). Next, note that the F ourier transforms of Re (Ψ 5 ( x )) and Im (Ψ 5 ( x )) can b e written as \ Re (Ψ 5 )( ω ) = 1 √ 2  1 + sign( ω x )sign( ω y )  ˆ ψ ( ω x ) ˆ ψ ( ω y ) , and \ Im (Ψ 5 )( ω ) = − j √ 2 (sign  ω x ) + sign( ω y )  ˆ ψ ( ω x ) ˆ ψ ( ω y ) . (51) The corresp ondence Im (Ψ 5 ( x )) = H π / 4 Re (Ψ 5 ( x )) then follo ws from the iden- tit y (sign( ω x ) + sign( ω y )) = sign( ω x + ω y )  1 + sign( ω x )sign( ω y )  . References [1] M. J. Bastiaans, Gab or’s exp ansion of a signal into Gaussian elementary signals , Pro c. IEEE 68 (1980), 538–539. [2] J. J. Benedetto, Harmonic Analysis and Applic ations , CRC Press, 1996. [3] T. Blu and M. Unser, The fr actional spline wavelet tr ansform: Deﬁnition and implementation , Pro c. IEEE International Conference on Acoustics, Sp eec h, and Signal Pro cessing (ICASSP’00) (2000), 512–515. [4] , A c omplete family of sc aling functions: The ( α, τ ) -fr actional splines , Pro c. IEEE In ternational Conference on Acoustics, Sp eech, and Signal Pro cessing (ICASSP’03) VI (2003), 421–424. [5] R. Bracew ell, The Fourier Tr ansform and Its Applic ations , McGra w-Hill, 1986. [6] T. B ¨ ulo w and G. Sommer, Hyp er c omplex signals - A novel extension of the analytic signal to the multidimensional c ase , IEEE T rans. Signal Pro cess. 49(11) (2001), 2844–2852. [7] W. L. Chan, H. Choi, and R.G. Baraniuk, Coher ent multisc ale image pr o- c essing using dual-tr e e quaternion wavelets , IEEE T rans. Image Pro cess. 17 (2008), 1069–1082. [8] C. Chaux, L. Duv al, and J.C. Pesquet, Image analysis using a dual-tr e e M-b and wavelet tr ansform , IEEE T rans. Image Process. 15 (2006), no. 8, 2397–2412. 33 [9] J. G. Daugman, Two-dimensional sp e ctr al analysis of c ortic al r e c eptive ﬁeld pr oﬁle , Vision Researc h 20 (1980), 847–856. [10] P . F. C. de Riv az and N. G. Kingsbury , Bayesian image de c onvolution and denoising using c omplex wavelets , Pro c. IEEE International Conference on Image Pro cessing (ICIP’01) 2 (2001), 273–276. [11] M. F elsb erg and G. Sommer, The mono genic signals , IEEE T rans. Signal Pro cess. 49(12) (2001), 3136–3144. [12] F. C. A. F ernandes, R. L. C. V an Spaendonck, and C. S. Burrus, A new fr amework for c omplex wavelet tr ansforms , IEEE T rans. Signal Pro cess. 51 (2003), no. 7, 1825–1837, 43. [13] D. Gabor, The ory of c ommunic ation , J. Inst. Elect. Eng. 93 (1946), 429– 457. [14] G. H. Granlund and H. Knutsson, Signal Pr o c essing for Computer Vision , c h. 4, Dordrech t, The Netherlands: Kluw er, 1995. [15] S.L. Hahn, Multidimensional c omplex signals with single-orthant sp e ctr a , Pro c. IEEE 80(8) (1992), 1287–1300. [16] S. Hatip oglu, S.K. Mitra, and N.G. Kingsbury , Image textur e description using c omplex wavelet tr ansform , Pro c. IEEE International Conference on Image Pro cessing (ICIP’00) 2 (2000), 530–533. [17] N. G. Kingsbury , Shift invariant pr op erties of the dual-tr e e c omplex wavelet tr ansform , Pro c. IEEE International Conference on Acoustics, Sp eec h, and Signal Pro cessing (ICASSP’99) (1999). [18] , A dual-tr e e c omplex wavelet tr ansform with impr ove d ortho gonality and symmetry pr op erties , Pro c. IEEE In ternational Conference on Image Pro cessing (ICIP’06) 2 (2000), 375–378. [19] , Complex wavelets for shift invariant analysis and ﬁltering of sig- nals , Journal of Applied and Computational Harmonic Analysis 10 (2001), no. 3, 234–253. [20] K. G. Larkin, D. J. Bone, and M. A. Oldﬁeld, Natur al demo dulation of two-dimensional fringe p atterns. I. Gener al b ackgr ound of the spir al phase quadr atur e tr ansforms , J. Opt. So c. Am. A 18(8) (2001), 1862–1870. [21] S. Mallat, A Wavelet Tour of Signal Pr o c essing , San Diego, CA: Academic Press, 1998. [22] S. G. Mallat, A the ory for multir esolution signal de c omp osition: The wavelet r epr esentation , IEEE T rans. Pattern Anal. Mach. In tell. 11 (1989), no. 7, 674–693. [23] Y. Meyer, Ondelettes et op ´ er ateurs i : Ondelettes , Hermann, P aris, 1990. 34 [24] S.C. Olhede and G. Metik as, The hyp er analytic wavelet tr ansform , Imp erial College Statistics Section T echnical Rep ort TR-06-02 (2008), 1–49. [25] M. S. Pattic his and A. C. Bovik, Analyzing image structur e by multidimen- sional fr e quency mo dulation , IEEE T rans. Pattern Anal. Mach. In tell. 29 (2007), no. 5, 753–766. [26] I. W. Selesnick, The design of appr oximate hilb ert tr ansform p airs of wavelet b ases , IEEE T rans. Signal Pro cess. 50 (2002), no. 2, 1143–1152. [27] I. W. Selesnick, R. G. Baraniuk, and N. C. Kingsbury , The dual-tr e e c om- plex wavelet tr ansform , IEEE Signal Pro cess. Mag. 22 (2005), no. 6, 123– 151. [28] I.W. Selesnick, Hilb ert tr ansform p airs of wavelet b ases , IEEE Signal Pro- cess. Lett. 8 (2001), no. 6, 170–173. [29] L. Sendur and I. W. Selesnic k, Bivariate shrinkage with lo c al varianc e es- timation , IEEE Signal Pro cess. Lett. 9 (2002), no. 12, 438–441. [30] E. M. Stein, Singular Inte gr als and Diﬀer entiability Pr op erty of Functions , Princeton Universit y Press, 1970. [31] E. M. Stein and R. Shak archi, Complex analysis , Princeton Universit y Press, 2003. [32] E. M. Stein and G. W eiss, F ourier Analysis on Euclide an Sp ac es , Princeton Univ ersity Press, 1971. [33] M. Unser, Sampling—50 ye ars after Shannon , Pro c. IEEE 88 (2000), no. 4, 569–587. [34] M. Unser, A. Aldroubi, and M. Eden, On the asymptotic c onver genc e of B-spline wavelets to Gab or functions , IEEE T rans. Inf. Theory 38 (1992), no. 2, 864–872. [35] M. Unser and T. Blu, Construction of fr actional spline wavelet b ases , Pro c. SPIE Conference on Mathematical Imaging: Wav elet Applications in Sig- nal and Image Pro cessing VI I 3813 (1999), 422–431. [36] M. Unser and T. Blu, F r actional splines and wavelets , SIAM Review 42 (2000), no. 1, 43–67. [37] , Wavelet the ory demystiﬁe d , IEEE T rans. Signal Pro cess. 51 (2003), no. 2, 470–483. [38] B. V an der Pol, The fundamental principles of fr e quency mo dulation , Jour- nal IEE 93 (1946), 153–158. [39] J. Ville, The orie et applic ation de la notion de signal analytique , Cables and T ransmissions 93 (1948), no. I I I, 153– 158. 35 [40] R. Y u and H. Ozk aramanli, Hilb ert tr ansform p airs of biortho gonal wavelet b ases , IEEE T rans. Signal Pro cess. 54 (2006), 2119–2125. 36

Construction of Hilbert Transform Pairs of Wavelet Bases and Gabor-like Transforms

Original Paper

Comments & Academic Discussion

Leave a Comment

Original Paper

Related Papers

Comments & Academic Discussion

Leave a Comment