Complex-valued Wavelets, the Dual Tree, and Hilbert Pairs - Signal ...
Complex-valued Wavelets, the Dual Tree, and Hilbert Pairs: why these lead to shift invariance and directional m-D wavelets? Nick Kingsbury Signal Processing and Communications Laboratory Department of Engineering, University of Cambridge, UK. email: [email protected]
web: www.eng.cam.ac.uk/~ngk
Part 3, UDRC Short Course, 21 March 2013
Nick Kingsbury (Univ. of Cambridge)
Dual-Tree Wavelets and Hilbert pairs
UDRC 2013 (3)
1 / 60
Introduction
Important continuous-time complex-valued wavelets The Gabor function ψ(t) = k e−t
2 /2σ 2
1 2 2 ˆ ψ(ω) = K e− 2 σ (ω−ω0 )
e jω0 t
The Morlet wavelet ψ(t) = k e−t
2 /2σ 2
(e jω0 t − κ0 )
1 2 1 2 2 2 ˆ ψ(ω) = K e− 2 σ (ω−ω0 ) − κ0 e− 2 σ ω
The Cauchy wavelet ψ(t) = k (1 − jβt)−α
ˆ ψ(ω) =
K ω α−1 e−ω/β , ω ≥ 0 0, ω 0 ψg (ω)
(10)
Note that the behaviour of the RHS near ω = 0 is immaterial, because, for wavelets to be admissible bandpass functions, ψˆg (0) = ψˆh (0) = 0. Substituting (8) into (5) and (9) into (7), we get the following expression for this ratio
ψˆh (ω) ψˆg (ω)
=
Q∞
−k ˆ k =2 H0 (2 ω) φh (0) Q∞ −k ˆ π) k =2 G0 (2 ω) φg (0)
e−jmω/2 H0∗ ( ω2 ± π)
e−jmω/2 G0∗ ( ω2 ± "∞ # Y = R0∗ ( ω2 ± π) R0 (2−k ω)
(11)
k =2
where R0 (ω) = H0 (ω)/G0 (ω) and is 2π-periodic. R0 will give the desired relationship between H0 and G0 if (10) and (11) are both satisfied. Nick Kingsbury (Univ. of Cambridge)
Dual-Tree Wavelets and Hilbert pairs
UDRC 2013 (3)
17 / 60
Hilbert-pair delay condition
Solving for R0 (ω) Hence
" R0∗ ( ω2
± π)
∞ Y
# −k
R0 (2
ω) =
k =2
j if ω < 0 −j if ω > 0
(12)
Since the modulus of the RHS of (12) is unity everywhere, and the LHS contains an infinite product of terms in R0 , which will tend to grow or shrink if R0 does not have unit magnitude, we deduce that |R0 (ω)| = 1. Now consider the phase θ(ω) of R0 , by letting R0 (ω) = e jθ(ω) Equating the phases in (12), we require that ( π ∞ X 2 −k ω −θ( 2 ± π) + θ(2 ω) = − π2 k =2 Nick Kingsbury (Univ. of Cambridge)
Dual-Tree Wavelets and Hilbert pairs
(13)
if ω < 0
(14)
if ω > 0 UDRC 2013 (3)
18 / 60
Hilbert-pair delay condition
Deducing the form of θ(ω) Because of the infinite sum in (14), we require θ(ω) → 0 as ω → 0. Hence θ(0) = 0. Since g0 (n) and h0 (n) are purely real and lowpass, R0 (ω) must be conjugate symmetric and so θ(ω) must be a continuous odd function about ω = 0. It can be shown (by Fourier analysis on θ0 (ω)) that any non-linear terms in θ(ω) would prevent (14) from being satisfied, because in (14) the gradient of the first term must cancel out the gradient of the summation terms at all ω 6= 0. Therefore we let θ(ω) = αω for −π < ω < π, where α is a constant. (15) ( Hence α( ω2 + π) if −4π < ω < 0 ω (16) θ( 2 ± π) = α( ω2 − π) if 0 < ω < 4π Note that θ(ω) must be 2π-periodic for |ω| ≥ π, and so it will have discontinuities at ω = ±π, if α is not an integer. These become discontinuities at ω = 0 in θ( ω2 ± π). Nick Kingsbury (Univ. of Cambridge)
Dual-Tree Wavelets and Hilbert pairs
UDRC 2013 (3)
19 / 60
Hilbert-pair delay condition
Typical plots of θ(ω) and terms in equ.(14)
θ(ω) ≈ αω, −π < ω < π and θ(ω) is 2π-periodic −θ( ω2 ± π) θ( ω4 ) θ( ω8 ) −θ( ω2 ± π) +
∞ X
θ(2−k ω)
k =2
−5
−4
−3
−2
−1
0
1
2
3
4
5
ω/π Nick Kingsbury (Univ. of Cambridge)
Dual-Tree Wavelets and Hilbert pairs
UDRC 2013 (3)
20 / 60
Hilbert-pair delay condition
Calculating α Noting that
∞ X
θ(2−k ω) = αω[ 41 +
1 8
+
1 16
+ . . .] =
αω 2
if −4π < ω < 4π,
k =2
and substituting (16) into (14) gives −α( ω2 + π) + and
−
α( ω2
− π) +
αω 2 αω 2
=
π 2
=
− π2
if −4π < ω < 0
(17)
if
(18)
0 < ω < 4π
(17) and (18) are both satisfied if α = − 12 , and therefore H0 (ω) = R0 (ω) = e jθ(ω) = e jαω = e−jω/2 G0 (ω)
for −π < ω < π
(19)
This is the half-sample delay solution that makes ψh (t) the Hilbert transform of ψg (t). Ozkaramanli and Yu (Dec 2005 and June 2006) have shown this solution to be unique and applicable to biorthogonal as well as orthonormal wavelet transforms. Nick Kingsbury (Univ. of Cambridge)
Dual-Tree Wavelets and Hilbert pairs
UDRC 2013 (3)
21 / 60
Hilbert-pair delay condition
Level 4 x0000a G0 Level 3 ↓2 Tree a (3q) x00a G0 Level 2 - ↓2 (3q) x0001a x0a G0 G1 Level 1 - ↓2 - ↓2 (q) (3q) 1 G1 - G0 - ↓ 2 - ↓2 - x 001a (q) (1) G1 - ↓2 - x x0000b 01a x000b H0 (q) 1 - ↓2 - G1 - ↓ 2 - x (q) x00b H0 (0) 1a - ↓2 (q) x0001b x0b H1 - H0 - ↓ 2 - ↓2 (q) (3q) 1 H1 - H0 - ↓ 2 - ↓2 - x (3q) 001b (0) H1 - ↓2 - x (3q) 01b 1 Tree b - H1 - ↓ 2 - x 1b (1) x000a
1-D Input Signal x
Fig. 1: Dual tree of real filters for the Q-shift CWT, giving real and imaginary parts of complex coefficients from tree a and tree b respectively. Figures in brackets indicate the approximate delay for each filter, where q = 1
1 4
sample period. Special level 1
1
filters, G and H , allow for the finite number of levels. Nick Kingsbury (Univ. of Cambridge)
Dual-Tree Wavelets and Hilbert pairs
UDRC 2013 (3)
22 / 60
Hilbert-pair delay condition
Why does the delay condition give shift invariance?
The half-sample delay between the G0 and H0 lowpass filters means that their output samples are uniformly interleaved at all scales, and hence the sample rate is effectively doubled everywhere. The doubled sampling rate is sufficient to virtually eliminate aliasing if filters of 12 or more taps are used. If aliasing is eliminated in the lowpass branch of each 2-band reconstruction block, then it must also be eliminated in the highpass branch (since perfect reconstruction is achieved). Elimination of aliasing means that each subband can be represented by a unique z-transfer function, and hence the filtering is LTI, linear time-invariant (i.e. shift-invariant). At level 1 of a finite dual tree, the delay difference must increase to one sample to compensate for the absence of delay differences at finer levels. Nick Kingsbury (Univ. of Cambridge)
Dual-Tree Wavelets and Hilbert pairs
UDRC 2013 (3)
23 / 60
Q-shift filters
Why do we use Q-shift filters (below level 1)? Half-sample delay difference is obtained with filter delays of 14 and 34 of a sample period (instead of 0 and 21 a sample for our original DT CWT). This is achieved with an asymmetric even-length filter G0 (z) and its time reverse H0 (z) = z −1 G0 (z −1 ). G1 (z) and H1 (z) are the CQFs of these. Due to the asymmetry (like Daubechies filters), these may be designed to give an orthonormal perfect reconstruction wavelet transform in each tree. Tree b filters are the reverse of tree a filters, and reconstruction filters are the reverse of analysis filters, so all filters are from the same orthonormal set, yielding a tight-frame transform. Both trees have the same frequency responses (in magnitude). The combined complex impulse responses are conjugate symmetric about their mid points, even though the separate responses are asymmetric. Hence symmetric extension still works at image edges. At level 1, any DWT filters can be used. Nick Kingsbury (Univ. of Cambridge)
Dual-Tree Wavelets and Hilbert pairs
UDRC 2013 (3)
24 / 60
Q-shift filters
Q-shift DT CWT Basis Functions – Levels 1 to 3 1−D complex wavelets & scaling functions at levels 1 to 3. 1
Level 1: scaling fn.
Tree a
0
Tree b |a + jb|
wavelet −1
Basis functions for adjacent sampling points are shown dotted.
Level 2: scaling fn. −2
wavelet −3
Level 3: scaling fn. −4
wavelet −5
−20
Nick Kingsbury (Univ. of Cambridge)
−15
−10
−5
0 5 sample number
Dual-Tree Wavelets and Hilbert pairs
10
15
20
UDRC 2013 (3)
25 / 60
Q-shift filters
Frequency Responses of 18-tap Q-shift filters Frequency responses of 18−tap Q−shift wavelets (n = 9) 1.6 Wavelets at level:
1.4
4 3
2
1
1.2 Scaling fn. at level 4 1 0.8 0.6 0.4 0.2 0 −8
−6
Nick Kingsbury (Univ. of Cambridge)
−4
−2 0 2 frequency (sample freq / 16) Dual-Tree Wavelets and Hilbert pairs
4
6
8 UDRC 2013 (3)
26 / 60
DT CWT in 2-D
The DT CWT in 2-D – why good directional selectivity? When the DT CWT is applied to 2-D signals (images), it has the following features: It is performed separably, using 2 trees for the rows of the image and 2 trees for the columns – yielding a Quad-Tree structure (4:1 redundancy). The 4 quad-tree components of each coefficient are combined by simple sum and difference operations to yield a pair of complex coefficients. These are part of two separate subbands in adjacent quadrants of the 2-D spectrum. This produces 6 directionally selective subbands at each level of the 2-D DT CWT. The next slide shows basis functions of these subbands at level 4, and compares them with the 3 subbands of a 2-D DWT. The DT CWT is directionally selective because the complex filters can separate positive and negative frequency components in 1-D, and hence separate adjacent quadrants of the 2-D spectrum. Real separable filters cannot do this! Nick Kingsbury (Univ. of Cambridge)
Dual-Tree Wavelets and Hilbert pairs
UDRC 2013 (3)
27 / 60
DT CWT in 2-D
HL2 (z) H0 (z2 ) G0 (z2 )
0.4 0.2
2-D Basis Functions at level 4
0 −0.2 −15 0.1
DT CWT real part
−10
−5
0
5
10
15
0
20
40
60
Level 4 Sc. funcs. 0.05 0
15
45
75
−75
−45
−15(deg)
−60 0.1
−40
−20
Level 4 Wavelets
Magnitude
0.05 0 −0.05
DT CWT imaginary part
−60
Real −40
−20
Imaginary 0
20
40
60
Time (samples)
Real DWT
e jω1 x . e jω2 y = e j(ω1 x+ω2 y ) only one wave! 90
45(?)
0
Basis functions of 2-D Q-shift complex wavelets (top), and of 2-D real wavelet filters (bottom), all illustrated at level 4 of the transforms. The complex wavelets provide 6 directionally selective filters, while real wavelets provide 3 filters, only two of which have a dominant direction. The 1-D bases, from which the 2-D complex bases are derived, are shown to the right. Nick Kingsbury (Univ. of Cambridge)
Dual-Tree Wavelets and Hilbert pairs
UDRC 2013 (3)
28 / 60
DT CWT in 2-D
Frequency Responses of 2-D Q-shift filters at levels 3 and 4 Contours of 2−D DT CWT freq responses at levels 3 & 4 4
3
2
1
Contours shown at −1 dB and −3 dB.
0
−1
−2
−3
−4 −4
Nick Kingsbury (Univ. of Cambridge)
−3
−2
−1 0 1 frequency (sample freq / 32)
Dual-Tree Wavelets and Hilbert pairs
2
3
4
UDRC 2013 (3)
29 / 60
Q-shift filter design
Q-shift Filter Design Requirements - H0 (z) - ↓ 2 - ↑ 2 - G0 (z) ?Y (z) X (z) + 6 - H1 (z) - ↓ 2 - ↑ 2 - G1 (z) Fig 2: 2-band analysis and reconstruction filter bank. H1 (z) = z −1 G0 (−z)
No aliasing:
G1 (z) = zH0 (−z) ;
Perfect reconstruction:
H0 (z)G0 (z) + H0 (−z)G0 (−z) = 2
Orthogonality:
G0 (z) = z −k H0 (z −1 )
Group delay '
1 4
sample period for H0 .
Good smoothness properties when iterated over scale.
Nick Kingsbury (Univ. of Cambridge)
Dual-Tree Wavelets and Hilbert pairs
UDRC 2013 (3)
30 / 60
Q-shift filter design
Filter Design — Delay To get 2n-tap lowpass filters, H0 (z) and G0 (z), with delays:
1 4
and
3 4
sample
Design a 4n-tap symmetric lowpass filter HL2 (z) with half the required bandwidth and a delay of 12 sample; Subsample HL2 (z) by 2:1 to get H0 (z) and G0 (z). 0.6 0.4 0.2 0 −0.2 −15
−10
−5
0
5
10
15
Fig 3: Impulse response of HL2 (z) for n = 6. The H0 and G0 filter taps0.6are shown as circles and crosses respectively. 0.4 Nick Kingsbury (Univ. of Cambridge)
Dual-Tree Wavelets and Hilbert pairs
UDRC 2013 (3)
31 / 60
Q-shift filter design
Filter Design – Perfect Reconstruction (PR) For PR and orthogonality, H0 (z) G0 (z) = H0 (z) H0 (z −1 ) must have no terms in z 2k except the term in z 0 . Therefore H0 (z 2 ) H0 (z −2 ) must have no terms in z 4k except the term in z 0 . But HL2 (z) = H0 (z 2 ) + z −1 H0 (z −2 ) and so HL2 (z) HL2 (z −1 ) = 2 H0 (z 2 ) H0 (z −2 ) + z −1 H02 (z −2 ) + z H02 (z 2 ) | {z } odd powers of z only Therefore HL2 (z) HL2 (z −1 ) must have no terms in z 4k except the term in z 0 . Hence we can include PR as a direct design constraint on HL2 (z) HL2 (z −1 ).
Nick Kingsbury (Univ. of Cambridge)
Dual-Tree Wavelets and Hilbert pairs
UDRC 2013 (3)
32 / 60
Q-shift filter design
Filter Design — Smoothness To obtain smoothness when iterated over many scales: Ensure that the stopband of H0 (z) suppresses energy at frequencies where unwanted passbands appear from subsampled filters operating at coarser scales. Consider the combined frequency response of H0 over just two scales: H0 (z) H0 (z 2 )|z=e jω = H0 (e jω ) H0 (e 2jω ) If the stopband of H0 (e jω ) covers ωs ≤ ω ≤ π, then the unwanted transition band and passband of H0 (e 2jω ) will extend from π − ω2s to π. For H0 (e jω ) to suppress the unwanted bands of H0 (e 2jω ): ωs ≤ π −
Nick Kingsbury (Univ. of Cambridge)
ωs 2
and so ωs ≤
2π 3
Dual-Tree Wavelets and Hilbert pairs
(see fig 4)
UDRC 2013 (3)
33 / 60
Q-shift filter design
Gain of filters for N = 12
i = 1, nzhpi = 1
10 T
0 −10 H0 (z) −20 −30 dB −40
HL2 (z) H0 (z)H0 (z2 )
−50 −60 −70 20
30
40
−80 0
0.2
0.4 0.6 Frequency ( ω / π )
0.8
1
Fig 4: Frequency responses of HL2 (z) (blue), H0 (z) (green), H0 (z) H0 (z 2 ) (magenta), and the gain correction matrix T (red) for n = 6 (12 taps for H0 ). Nick Kingsbury (Univ. of Cambridge)
Dual-Tree Wavelets and Hilbert pairs
UDRC 2013 (3)
34 / 60
Q-shift filter design
Optimization for MSE in the frequency domain We have now reduced the ideal design conditions for the length 4n symmetric lowpass filter HL2 to be: Zero amplitude for all the terms of HL2 (z) HL2 (z −1 ) in z 4k except the term in z 0 , which must be 1 (these are quadratic constraints on coef vector hL2 ); Zero (or near-zero) amplitude of HL2 (e jω ) for the stopband, (these are linear constraints on hL2 ).
π 3
≤ω≤π
If all constraints were linear, the LMS error solution for hL2 could be found using a matrix pseudo-inverse method. Therefore we linearise the problem and iterate. If hL2 at iteration i is hi = hi−1 + ∆hi , then hi ∗ hi = (hi−1 + ∆hi ) ∗ (hi−1 + ∆hi ) = hi−1 ∗ (hi−1 + 2∆hi ) + ∆hi ∗ ∆hi Since ∆hi becomes small as i increases, the final term can be neglected and the convolution (∗) is expressed as a linear function of ∆hi . Nick Kingsbury (Univ. of Cambridge)
Dual-Tree Wavelets and Hilbert pairs
UDRC 2013 (3)
35 / 60
Q-shift filter design
Hence we solve for ∆hi such that: Ci−1 (hi−1 + 2∆hi )
=
[0 . . . 0 1]T
F (hi−1 + ∆hi ) ' [0 . . . 0]T where Ci−1 calculates every 4th term in the convolution with hi−1 , and F evaluates the Fourier transform at M discrete frequencies ω from π3 to π (typically M ' 8n). Note that only one side of the symmetric convolution is needed in the rows of Ci−1 , and the columns of Ci−1 and F can be combined in pairs so that only the first half of the symmetric ∆hi need be solved for. To obtain high accuracy solutions to the PR constraints, we scale the equations in Ci−1 up by βi = 2i to get the following iterative LMS method for ∆hi and then hi : 2βi Ci−1 βi (c − Ci−1 hi−1 ) ∆hi = and hi = hi−1 + ∆hi F −F hi−1 where c = [0 . . . 0 1]T . Nick Kingsbury (Univ. of Cambridge)
Dual-Tree Wavelets and Hilbert pairs
UDRC 2013 (3)
36 / 60
Q-shift filter design
Two final refinements To include transition band effects, we scale rows of F by diagonal matrix Ti , the gain (at iteration i) of H0 (z 2 )/H0 (1) at frequencies corresponding to π3 ≤ ω ≤ π2 in the frequency domain of HL2 (Ti is the red curve in fig 4). To insert predefined zeros in H0 (z) or HL2 (z), we first note that a zero at z = e jπ in H0 will be produced by a pair of zeros at z = e ±jπ/2 in HL2 . We can force zeros in HL2 by forming a convolution matrix Hf such that Hf h0i = hi , where h0i is the coef vector of the filter which represents all the zeros of HL2 that are not predefined, and Hf produces convolution with the predefined zeros. Hence we now solve for ∆h0i and then hi using 2βi Ci−1 βi (c − Ci−1 hi−1 ) 0 Hf ∆hi = with hi = hi−1 + Hf ∆h0i Ti−1 F −Ti−1 F hi−1
Nick Kingsbury (Univ. of Cambridge)
Dual-Tree Wavelets and Hilbert pairs
UDRC 2013 (3)
37 / 60
Q-shift filter design
10 0 −10 −20 n=8 −30 dB −40 −50 −60 n = 16 −70 −80 0
0.2
0.4 0.6 Frequency ( ω / π )
0.8
1
Fig 5: Frequency responses of HL2 (z) for n = 8 (blue), n = 12 (green) and n = 16 (red). Each filter has one predefined zero at ω = π2 and one at ω = π. Nick Kingsbury (Univ. of Cambridge)
Dual-Tree Wavelets and Hilbert pairs
UDRC 2013 (3)
38 / 60
Q-shift filter design
Initialisation To initialise the iterative algorithm when i = 1, we must define h0 and hence C0 and T0 . This is not critical and can be achieved by a simple inverse FFT of an ‘ideal’ lowpass frequency response for HL2 (e jω ) with a root-raised-cosine transition band covering the range π π 6
web: www.eng.cam.ac.uk/~ngk
Part 3, UDRC Short Course, 21 March 2013
Nick Kingsbury (Univ. of Cambridge)
Dual-Tree Wavelets and Hilbert pairs
UDRC 2013 (3)
1 / 60
Introduction
Important continuous-time complex-valued wavelets The Gabor function ψ(t) = k e−t
2 /2σ 2
1 2 2 ˆ ψ(ω) = K e− 2 σ (ω−ω0 )
e jω0 t
The Morlet wavelet ψ(t) = k e−t
2 /2σ 2
(e jω0 t − κ0 )
1 2 1 2 2 2 ˆ ψ(ω) = K e− 2 σ (ω−ω0 ) − κ0 e− 2 σ ω
The Cauchy wavelet ψ(t) = k (1 − jβt)−α
ˆ ψ(ω) =
K ω α−1 e−ω/β , ω ≥ 0 0, ω 0 ψg (ω)
(10)
Note that the behaviour of the RHS near ω = 0 is immaterial, because, for wavelets to be admissible bandpass functions, ψˆg (0) = ψˆh (0) = 0. Substituting (8) into (5) and (9) into (7), we get the following expression for this ratio
ψˆh (ω) ψˆg (ω)
=
Q∞
−k ˆ k =2 H0 (2 ω) φh (0) Q∞ −k ˆ π) k =2 G0 (2 ω) φg (0)
e−jmω/2 H0∗ ( ω2 ± π)
e−jmω/2 G0∗ ( ω2 ± "∞ # Y = R0∗ ( ω2 ± π) R0 (2−k ω)
(11)
k =2
where R0 (ω) = H0 (ω)/G0 (ω) and is 2π-periodic. R0 will give the desired relationship between H0 and G0 if (10) and (11) are both satisfied. Nick Kingsbury (Univ. of Cambridge)
Dual-Tree Wavelets and Hilbert pairs
UDRC 2013 (3)
17 / 60
Hilbert-pair delay condition
Solving for R0 (ω) Hence
" R0∗ ( ω2
± π)
∞ Y
# −k
R0 (2
ω) =
k =2
j if ω < 0 −j if ω > 0
(12)
Since the modulus of the RHS of (12) is unity everywhere, and the LHS contains an infinite product of terms in R0 , which will tend to grow or shrink if R0 does not have unit magnitude, we deduce that |R0 (ω)| = 1. Now consider the phase θ(ω) of R0 , by letting R0 (ω) = e jθ(ω) Equating the phases in (12), we require that ( π ∞ X 2 −k ω −θ( 2 ± π) + θ(2 ω) = − π2 k =2 Nick Kingsbury (Univ. of Cambridge)
Dual-Tree Wavelets and Hilbert pairs
(13)
if ω < 0
(14)
if ω > 0 UDRC 2013 (3)
18 / 60
Hilbert-pair delay condition
Deducing the form of θ(ω) Because of the infinite sum in (14), we require θ(ω) → 0 as ω → 0. Hence θ(0) = 0. Since g0 (n) and h0 (n) are purely real and lowpass, R0 (ω) must be conjugate symmetric and so θ(ω) must be a continuous odd function about ω = 0. It can be shown (by Fourier analysis on θ0 (ω)) that any non-linear terms in θ(ω) would prevent (14) from being satisfied, because in (14) the gradient of the first term must cancel out the gradient of the summation terms at all ω 6= 0. Therefore we let θ(ω) = αω for −π < ω < π, where α is a constant. (15) ( Hence α( ω2 + π) if −4π < ω < 0 ω (16) θ( 2 ± π) = α( ω2 − π) if 0 < ω < 4π Note that θ(ω) must be 2π-periodic for |ω| ≥ π, and so it will have discontinuities at ω = ±π, if α is not an integer. These become discontinuities at ω = 0 in θ( ω2 ± π). Nick Kingsbury (Univ. of Cambridge)
Dual-Tree Wavelets and Hilbert pairs
UDRC 2013 (3)
19 / 60
Hilbert-pair delay condition
Typical plots of θ(ω) and terms in equ.(14)
θ(ω) ≈ αω, −π < ω < π and θ(ω) is 2π-periodic −θ( ω2 ± π) θ( ω4 ) θ( ω8 ) −θ( ω2 ± π) +
∞ X
θ(2−k ω)
k =2
−5
−4
−3
−2
−1
0
1
2
3
4
5
ω/π Nick Kingsbury (Univ. of Cambridge)
Dual-Tree Wavelets and Hilbert pairs
UDRC 2013 (3)
20 / 60
Hilbert-pair delay condition
Calculating α Noting that
∞ X
θ(2−k ω) = αω[ 41 +
1 8
+
1 16
+ . . .] =
αω 2
if −4π < ω < 4π,
k =2
and substituting (16) into (14) gives −α( ω2 + π) + and
−
α( ω2
− π) +
αω 2 αω 2
=
π 2
=
− π2
if −4π < ω < 0
(17)
if
(18)
0 < ω < 4π
(17) and (18) are both satisfied if α = − 12 , and therefore H0 (ω) = R0 (ω) = e jθ(ω) = e jαω = e−jω/2 G0 (ω)
for −π < ω < π
(19)
This is the half-sample delay solution that makes ψh (t) the Hilbert transform of ψg (t). Ozkaramanli and Yu (Dec 2005 and June 2006) have shown this solution to be unique and applicable to biorthogonal as well as orthonormal wavelet transforms. Nick Kingsbury (Univ. of Cambridge)
Dual-Tree Wavelets and Hilbert pairs
UDRC 2013 (3)
21 / 60
Hilbert-pair delay condition
Level 4 x0000a G0 Level 3 ↓2 Tree a (3q) x00a G0 Level 2 - ↓2 (3q) x0001a x0a G0 G1 Level 1 - ↓2 - ↓2 (q) (3q) 1 G1 - G0 - ↓ 2 - ↓2 - x 001a (q) (1) G1 - ↓2 - x x0000b 01a x000b H0 (q) 1 - ↓2 - G1 - ↓ 2 - x (q) x00b H0 (0) 1a - ↓2 (q) x0001b x0b H1 - H0 - ↓ 2 - ↓2 (q) (3q) 1 H1 - H0 - ↓ 2 - ↓2 - x (3q) 001b (0) H1 - ↓2 - x (3q) 01b 1 Tree b - H1 - ↓ 2 - x 1b (1) x000a
1-D Input Signal x
Fig. 1: Dual tree of real filters for the Q-shift CWT, giving real and imaginary parts of complex coefficients from tree a and tree b respectively. Figures in brackets indicate the approximate delay for each filter, where q = 1
1 4
sample period. Special level 1
1
filters, G and H , allow for the finite number of levels. Nick Kingsbury (Univ. of Cambridge)
Dual-Tree Wavelets and Hilbert pairs
UDRC 2013 (3)
22 / 60
Hilbert-pair delay condition
Why does the delay condition give shift invariance?
The half-sample delay between the G0 and H0 lowpass filters means that their output samples are uniformly interleaved at all scales, and hence the sample rate is effectively doubled everywhere. The doubled sampling rate is sufficient to virtually eliminate aliasing if filters of 12 or more taps are used. If aliasing is eliminated in the lowpass branch of each 2-band reconstruction block, then it must also be eliminated in the highpass branch (since perfect reconstruction is achieved). Elimination of aliasing means that each subband can be represented by a unique z-transfer function, and hence the filtering is LTI, linear time-invariant (i.e. shift-invariant). At level 1 of a finite dual tree, the delay difference must increase to one sample to compensate for the absence of delay differences at finer levels. Nick Kingsbury (Univ. of Cambridge)
Dual-Tree Wavelets and Hilbert pairs
UDRC 2013 (3)
23 / 60
Q-shift filters
Why do we use Q-shift filters (below level 1)? Half-sample delay difference is obtained with filter delays of 14 and 34 of a sample period (instead of 0 and 21 a sample for our original DT CWT). This is achieved with an asymmetric even-length filter G0 (z) and its time reverse H0 (z) = z −1 G0 (z −1 ). G1 (z) and H1 (z) are the CQFs of these. Due to the asymmetry (like Daubechies filters), these may be designed to give an orthonormal perfect reconstruction wavelet transform in each tree. Tree b filters are the reverse of tree a filters, and reconstruction filters are the reverse of analysis filters, so all filters are from the same orthonormal set, yielding a tight-frame transform. Both trees have the same frequency responses (in magnitude). The combined complex impulse responses are conjugate symmetric about their mid points, even though the separate responses are asymmetric. Hence symmetric extension still works at image edges. At level 1, any DWT filters can be used. Nick Kingsbury (Univ. of Cambridge)
Dual-Tree Wavelets and Hilbert pairs
UDRC 2013 (3)
24 / 60
Q-shift filters
Q-shift DT CWT Basis Functions – Levels 1 to 3 1−D complex wavelets & scaling functions at levels 1 to 3. 1
Level 1: scaling fn.
Tree a
0
Tree b |a + jb|
wavelet −1
Basis functions for adjacent sampling points are shown dotted.
Level 2: scaling fn. −2
wavelet −3
Level 3: scaling fn. −4
wavelet −5
−20
Nick Kingsbury (Univ. of Cambridge)
−15
−10
−5
0 5 sample number
Dual-Tree Wavelets and Hilbert pairs
10
15
20
UDRC 2013 (3)
25 / 60
Q-shift filters
Frequency Responses of 18-tap Q-shift filters Frequency responses of 18−tap Q−shift wavelets (n = 9) 1.6 Wavelets at level:
1.4
4 3
2
1
1.2 Scaling fn. at level 4 1 0.8 0.6 0.4 0.2 0 −8
−6
Nick Kingsbury (Univ. of Cambridge)
−4
−2 0 2 frequency (sample freq / 16) Dual-Tree Wavelets and Hilbert pairs
4
6
8 UDRC 2013 (3)
26 / 60
DT CWT in 2-D
The DT CWT in 2-D – why good directional selectivity? When the DT CWT is applied to 2-D signals (images), it has the following features: It is performed separably, using 2 trees for the rows of the image and 2 trees for the columns – yielding a Quad-Tree structure (4:1 redundancy). The 4 quad-tree components of each coefficient are combined by simple sum and difference operations to yield a pair of complex coefficients. These are part of two separate subbands in adjacent quadrants of the 2-D spectrum. This produces 6 directionally selective subbands at each level of the 2-D DT CWT. The next slide shows basis functions of these subbands at level 4, and compares them with the 3 subbands of a 2-D DWT. The DT CWT is directionally selective because the complex filters can separate positive and negative frequency components in 1-D, and hence separate adjacent quadrants of the 2-D spectrum. Real separable filters cannot do this! Nick Kingsbury (Univ. of Cambridge)
Dual-Tree Wavelets and Hilbert pairs
UDRC 2013 (3)
27 / 60
DT CWT in 2-D
HL2 (z) H0 (z2 ) G0 (z2 )
0.4 0.2
2-D Basis Functions at level 4
0 −0.2 −15 0.1
DT CWT real part
−10
−5
0
5
10
15
0
20
40
60
Level 4 Sc. funcs. 0.05 0
15
45
75
−75
−45
−15(deg)
−60 0.1
−40
−20
Level 4 Wavelets
Magnitude
0.05 0 −0.05
DT CWT imaginary part
−60
Real −40
−20
Imaginary 0
20
40
60
Time (samples)
Real DWT
e jω1 x . e jω2 y = e j(ω1 x+ω2 y ) only one wave! 90
45(?)
0
Basis functions of 2-D Q-shift complex wavelets (top), and of 2-D real wavelet filters (bottom), all illustrated at level 4 of the transforms. The complex wavelets provide 6 directionally selective filters, while real wavelets provide 3 filters, only two of which have a dominant direction. The 1-D bases, from which the 2-D complex bases are derived, are shown to the right. Nick Kingsbury (Univ. of Cambridge)
Dual-Tree Wavelets and Hilbert pairs
UDRC 2013 (3)
28 / 60
DT CWT in 2-D
Frequency Responses of 2-D Q-shift filters at levels 3 and 4 Contours of 2−D DT CWT freq responses at levels 3 & 4 4
3
2
1
Contours shown at −1 dB and −3 dB.
0
−1
−2
−3
−4 −4
Nick Kingsbury (Univ. of Cambridge)
−3
−2
−1 0 1 frequency (sample freq / 32)
Dual-Tree Wavelets and Hilbert pairs
2
3
4
UDRC 2013 (3)
29 / 60
Q-shift filter design
Q-shift Filter Design Requirements - H0 (z) - ↓ 2 - ↑ 2 - G0 (z) ?Y (z) X (z) + 6 - H1 (z) - ↓ 2 - ↑ 2 - G1 (z) Fig 2: 2-band analysis and reconstruction filter bank. H1 (z) = z −1 G0 (−z)
No aliasing:
G1 (z) = zH0 (−z) ;
Perfect reconstruction:
H0 (z)G0 (z) + H0 (−z)G0 (−z) = 2
Orthogonality:
G0 (z) = z −k H0 (z −1 )
Group delay '
1 4
sample period for H0 .
Good smoothness properties when iterated over scale.
Nick Kingsbury (Univ. of Cambridge)
Dual-Tree Wavelets and Hilbert pairs
UDRC 2013 (3)
30 / 60
Q-shift filter design
Filter Design — Delay To get 2n-tap lowpass filters, H0 (z) and G0 (z), with delays:
1 4
and
3 4
sample
Design a 4n-tap symmetric lowpass filter HL2 (z) with half the required bandwidth and a delay of 12 sample; Subsample HL2 (z) by 2:1 to get H0 (z) and G0 (z). 0.6 0.4 0.2 0 −0.2 −15
−10
−5
0
5
10
15
Fig 3: Impulse response of HL2 (z) for n = 6. The H0 and G0 filter taps0.6are shown as circles and crosses respectively. 0.4 Nick Kingsbury (Univ. of Cambridge)
Dual-Tree Wavelets and Hilbert pairs
UDRC 2013 (3)
31 / 60
Q-shift filter design
Filter Design – Perfect Reconstruction (PR) For PR and orthogonality, H0 (z) G0 (z) = H0 (z) H0 (z −1 ) must have no terms in z 2k except the term in z 0 . Therefore H0 (z 2 ) H0 (z −2 ) must have no terms in z 4k except the term in z 0 . But HL2 (z) = H0 (z 2 ) + z −1 H0 (z −2 ) and so HL2 (z) HL2 (z −1 ) = 2 H0 (z 2 ) H0 (z −2 ) + z −1 H02 (z −2 ) + z H02 (z 2 ) | {z } odd powers of z only Therefore HL2 (z) HL2 (z −1 ) must have no terms in z 4k except the term in z 0 . Hence we can include PR as a direct design constraint on HL2 (z) HL2 (z −1 ).
Nick Kingsbury (Univ. of Cambridge)
Dual-Tree Wavelets and Hilbert pairs
UDRC 2013 (3)
32 / 60
Q-shift filter design
Filter Design — Smoothness To obtain smoothness when iterated over many scales: Ensure that the stopband of H0 (z) suppresses energy at frequencies where unwanted passbands appear from subsampled filters operating at coarser scales. Consider the combined frequency response of H0 over just two scales: H0 (z) H0 (z 2 )|z=e jω = H0 (e jω ) H0 (e 2jω ) If the stopband of H0 (e jω ) covers ωs ≤ ω ≤ π, then the unwanted transition band and passband of H0 (e 2jω ) will extend from π − ω2s to π. For H0 (e jω ) to suppress the unwanted bands of H0 (e 2jω ): ωs ≤ π −
Nick Kingsbury (Univ. of Cambridge)
ωs 2
and so ωs ≤
2π 3
Dual-Tree Wavelets and Hilbert pairs
(see fig 4)
UDRC 2013 (3)
33 / 60
Q-shift filter design
Gain of filters for N = 12
i = 1, nzhpi = 1
10 T
0 −10 H0 (z) −20 −30 dB −40
HL2 (z) H0 (z)H0 (z2 )
−50 −60 −70 20
30
40
−80 0
0.2
0.4 0.6 Frequency ( ω / π )
0.8
1
Fig 4: Frequency responses of HL2 (z) (blue), H0 (z) (green), H0 (z) H0 (z 2 ) (magenta), and the gain correction matrix T (red) for n = 6 (12 taps for H0 ). Nick Kingsbury (Univ. of Cambridge)
Dual-Tree Wavelets and Hilbert pairs
UDRC 2013 (3)
34 / 60
Q-shift filter design
Optimization for MSE in the frequency domain We have now reduced the ideal design conditions for the length 4n symmetric lowpass filter HL2 to be: Zero amplitude for all the terms of HL2 (z) HL2 (z −1 ) in z 4k except the term in z 0 , which must be 1 (these are quadratic constraints on coef vector hL2 ); Zero (or near-zero) amplitude of HL2 (e jω ) for the stopband, (these are linear constraints on hL2 ).
π 3
≤ω≤π
If all constraints were linear, the LMS error solution for hL2 could be found using a matrix pseudo-inverse method. Therefore we linearise the problem and iterate. If hL2 at iteration i is hi = hi−1 + ∆hi , then hi ∗ hi = (hi−1 + ∆hi ) ∗ (hi−1 + ∆hi ) = hi−1 ∗ (hi−1 + 2∆hi ) + ∆hi ∗ ∆hi Since ∆hi becomes small as i increases, the final term can be neglected and the convolution (∗) is expressed as a linear function of ∆hi . Nick Kingsbury (Univ. of Cambridge)
Dual-Tree Wavelets and Hilbert pairs
UDRC 2013 (3)
35 / 60
Q-shift filter design
Hence we solve for ∆hi such that: Ci−1 (hi−1 + 2∆hi )
=
[0 . . . 0 1]T
F (hi−1 + ∆hi ) ' [0 . . . 0]T where Ci−1 calculates every 4th term in the convolution with hi−1 , and F evaluates the Fourier transform at M discrete frequencies ω from π3 to π (typically M ' 8n). Note that only one side of the symmetric convolution is needed in the rows of Ci−1 , and the columns of Ci−1 and F can be combined in pairs so that only the first half of the symmetric ∆hi need be solved for. To obtain high accuracy solutions to the PR constraints, we scale the equations in Ci−1 up by βi = 2i to get the following iterative LMS method for ∆hi and then hi : 2βi Ci−1 βi (c − Ci−1 hi−1 ) ∆hi = and hi = hi−1 + ∆hi F −F hi−1 where c = [0 . . . 0 1]T . Nick Kingsbury (Univ. of Cambridge)
Dual-Tree Wavelets and Hilbert pairs
UDRC 2013 (3)
36 / 60
Q-shift filter design
Two final refinements To include transition band effects, we scale rows of F by diagonal matrix Ti , the gain (at iteration i) of H0 (z 2 )/H0 (1) at frequencies corresponding to π3 ≤ ω ≤ π2 in the frequency domain of HL2 (Ti is the red curve in fig 4). To insert predefined zeros in H0 (z) or HL2 (z), we first note that a zero at z = e jπ in H0 will be produced by a pair of zeros at z = e ±jπ/2 in HL2 . We can force zeros in HL2 by forming a convolution matrix Hf such that Hf h0i = hi , where h0i is the coef vector of the filter which represents all the zeros of HL2 that are not predefined, and Hf produces convolution with the predefined zeros. Hence we now solve for ∆h0i and then hi using 2βi Ci−1 βi (c − Ci−1 hi−1 ) 0 Hf ∆hi = with hi = hi−1 + Hf ∆h0i Ti−1 F −Ti−1 F hi−1
Nick Kingsbury (Univ. of Cambridge)
Dual-Tree Wavelets and Hilbert pairs
UDRC 2013 (3)
37 / 60
Q-shift filter design
10 0 −10 −20 n=8 −30 dB −40 −50 −60 n = 16 −70 −80 0
0.2
0.4 0.6 Frequency ( ω / π )
0.8
1
Fig 5: Frequency responses of HL2 (z) for n = 8 (blue), n = 12 (green) and n = 16 (red). Each filter has one predefined zero at ω = π2 and one at ω = π. Nick Kingsbury (Univ. of Cambridge)
Dual-Tree Wavelets and Hilbert pairs
UDRC 2013 (3)
38 / 60
Q-shift filter design
Initialisation To initialise the iterative algorithm when i = 1, we must define h0 and hence C0 and T0 . This is not critical and can be achieved by a simple inverse FFT of an ‘ideal’ lowpass frequency response for HL2 (e jω ) with a root-raised-cosine transition band covering the range π π 6
