This example shows how the dual-tree complex discrete wavelet transform (DT-CWT) provides advantages over the critically sampled DWT for signal, image, and volume processing. The dual-tree DWT is implemented as two separate two-channel filter banks. To gain the advantages described in this example, you cannot arbitrarily choose the scaling and wavelet filters used in the two trees. The lowpass (scaling) and highpass (wavelet) filters of one tree, $$\{{h}_{0},{h}_{1}\}$$, must generate a scaling function and wavelet that are approximate Hilbert transforms of the scaling function and wavelet generated by the lowpass and highpass filters of the other tree, $$\{{g}_{0},{g}_{1}\}$$. Therefore, the complex-valued scaling functions and wavelets formed from the two trees are approximately analytic.

As a result, the dual-tree DWT exhibits less shift variance and more directional selectivity than the critically sampled DWT with only a $${2}^{d}$$ redundancy factor for $$d$$-dimensional data. The redundancy in the dual-tree DWT is significantly less than the redundancy in the undecimated (stationary) DWT.

This example illustrates the approximate shift invariance of the dual-tree DWT, the selective orientation of the dual-tree analyzing wavelets in 2-D and 3-D, and the use of the dual-tree complex wavelet transform in image and volume denoising.

The DWT suffers from shift variance, meaning that small shifts in the input signal or image can cause significant changes in the distribution of signal/image energy across scales in the DWT coefficients. The DT-CWT is approximately shift invariant.

To demonstrate this on a test signal, construct two shifted discrete-time impulses 128 samples in length. One signal has the unit impulse at sample 60, while the other signal has the unit impulse at sample 64. Both signals clearly have unit energy ($${\ell}^{2}$$ norm).

kronDelta1 = zeros(128,1); kronDelta1(60) = 1; kronDelta2 = zeros(128,1); kronDelta2(64) = 1;

Obtain the DWT and dual-tree DWT of the two signals down to level 3 with wavelet and scaling filters of length 14. Extract the level-3 detail coefficients for comparison.

J = 3; dwt1 = dddtree('dwt',kronDelta1,J,'sym7'); dwt2 = dddtree('dwt',kronDelta2,J,'sym7'); dwt1Cfs = dwt1.cfs{J}; dwt2Cfs = dwt2.cfs{J}; dt1 = dddtree('cplxdt',kronDelta1,J,'dtf3'); dt2 = dddtree('cplxdt',kronDelta2,J,'dtf3'); dt1Cfs = dt1.cfs{J}(:,:,1)+1i*dt1.cfs{J}(:,:,2); dt2Cfs = dt2.cfs{J}(:,:,1)+1i*dt2.cfs{J}(:,:,2);

Plot the absolute values of the DWT and DT-CWT coefficients for the two signals at level 3 and compute the energy (squared $${\ell}^{2}$$ norms) of the coefficients.

figure subplot(1,2,1) stem(abs(dwt1Cfs),'markerfacecolor',[0 0 1]) title({'DWT';['Squared 2-norm = ' num2str(norm(dwt1Cfs,2)^2,3)]},... 'fontsize',10) subplot(1,2,2) stem(abs(dwt2Cfs),'markerfacecolor',[0 0 1]) title({'DWT';['Squared 2-norm = ' num2str(norm(dwt2Cfs,2)^2,3)]},... 'fontsize',10)

figure subplot(1,2,1) stem(abs(dt1Cfs),'markerfacecolor',[0 0 1]) title({'Dual-tree DWT';['Squared 2-norm = ' num2str(norm(dt1Cfs,2)^2,3)]},... 'fontsize',10) subplot(1,2,2) stem(abs(dwt2Cfs),'markerfacecolor',[0 0 1]) title({'Dual-tree DWT';['Squared 2-norm = ' num2str(norm(dt2Cfs,2)^2,3)]},... 'fontsize',10)

Note the four sample shift in the signal has caused an almost 6.5% change in the energy of the level-3 DWT wavelet coefficients. However, the dual-tree wavelet coefficients show only a 0.3% change in energy.

To demonstrate the utility of approximate shift invariance in real data, we analyze an electrocardiogram (ECG) signal. The sampling interval for the ECG signal is 1/180 seconds. The data are taken from Percival & Walden [3], p.125 (data originally provided by William Constantine and Per Reinhall, University of Washington). For convenience, we take the data to start at t=0.

load wecg dt = 1/180; t = 0:dt:(length(wecg)*dt)-dt; figure plot(t,wecg) xlabel('Seconds') ylabel('Millivolts')

The large positive peaks approximately 0.7 seconds apart are the R waves of the cardiac rhythm. First, decompose the signal using the critically sampled DWT and plot the original signal along with the level-2 and level-3 wavelet coefficients. The level-2 and level-3 coefficients were chosen because the R waves are isolated most prominently in those scales for the given sampling rate.

J = 6; dtDWT1 = dddtree('dwt',wecg,J,'farras'); details = zeros(2048,3); details(2:4:end,2) = dtDWT1.cfs{2}; details(4:8:end,3) = dtDWT1.cfs{3}; subplot(3,1,1) stem(t,details(:,2),'Marker','none','ShowBaseline','off') title('Level 2') ylabel('mV') subplot(3,1,2) stem(t,details(:,3),'Marker','none','ShowBaseline','off') title('Level 3') ylabel('mV') subplot(3,1,3) plot(t,wecg) title('Original Signal') xlabel('Seconds') ylabel('mV')

Repeat the above analysis for the dual-tree transform. In this case, just plot the real part of the dual-tree coefficients at levels 2 and 3.

dtcplx1 = dddtree('cplxdt',wecg,J,'dtf3'); details = zeros(2048,3); details(2:4:end,2) = dtcplx1.cfs{2}(:,1,1)+1i*dtcplx1.cfs{2}(:,1,2); details(4:8:end,3) = dtcplx1.cfs{3}(:,1,1)+1i*dtcplx1.cfs{3}(:,1,2); subplot(3,1,1) stem(t,real(details(:,2)),'Marker','none','ShowBaseline','off') title('Level 2') ylabel('mV') subplot(3,1,2) stem(t,real(details(:,3)),'Marker','none','ShowBaseline','off') title('Level 3') ylabel('mV') subplot(3,1,3) plot(t,wecg) title('Original Signal') xlabel('Seconds') ylabel('mV')

Both the critically sampled and dual-tree wavelet transforms localize an important feature of the ECG waveform to similar scales.

An important application of wavelets in 1-D signals is to obtain an analysis of variance by scale. It stands to reason that this analysis of variance should not be sensitive to circular shifts in the input signal. Unfortunately, this is not the case with the critically sampled DWT. To demonstrate this, we circularly shift the ECG signal by 4 samples, analyze the unshifted and shifted signals with the critically sampled DWT, and calculate the distribution of energy across scales.

wecgShift = circshift(wecg,4); dtDWT2 = dddtree('dwt',wecgShift,J,'farras'); sigenrgy = norm(wecg,2)^2; enr1 = cell2mat(cellfun(@(x)(norm(x,2)^2/sigenrgy)*100,dtDWT1.cfs,'uni',0)); enr2 = cell2mat(cellfun(@(x)(norm(x,2)^2/sigenrgy)*100,dtDWT2.cfs,'uni',0)); levels = {'D1';'D2';'D3';'D4';'D5';'D6';'A6'}; enr1 = enr1(:); enr2 = enr2(:); table(levels,enr1,enr2,'VariableNames',{'Level','enr1','enr2'})

`ans=`*7×3 table*
Level enr1 enr2
______ ______ ______
{'D1'} 4.1994 4.1994
{'D2'} 8.425 8.425
{'D3'} 13.381 10.077
{'D4'} 7.0612 10.031
{'D5'} 5.4606 5.4436
{'D6'} 3.1273 3.4584
{'A6'} 58.345 58.366

Note that the wavelet coefficients at levels 3 and 4 show approximately 3% changes in energy between the original and shifted signal. Next, we repeat this analysis using the complex dual-tree wavelet transform.

dtcplx2 = dddtree('cplxdt',wecgShift,J,'dtf3'); cfs1 = cellfun(@squeeze,dtcplx1.cfs,'uni',0); cfs2 = cellfun(@squeeze,dtcplx2.cfs,'uni',0); cfs1 = cellfun(@(x) complex(x(:,1),x(:,2)),cfs1,'uni',0); cfs2 = cellfun(@(x) complex(x(:,1),x(:,2)),cfs2,'uni',0); dtenr1 = cell2mat(cellfun(@(x)(norm(x,2)^2/sigenrgy)*100,cfs1,'uni',0)); dtenr2 = cell2mat(cellfun(@(x)(norm(x,2)^2/sigenrgy)*100,cfs2,'uni',0)); dtenr1 = dtenr1(:); dtenr2 = dtenr2(:); table(levels,dtenr1,dtenr2, 'VariableNames',{'Level','dtenr1','dtenr2'})

`ans=`*7×3 table*
Level dtenr1 dtenr2
______ ______ ______
{'D1'} 4.936 4.936
{'D2'} 6.6691 6.6691
{'D3'} 12.682 12.611
{'D4'} 8.3891 8.4808
{'D5'} 5.8868 5.8791
{'D6'} 3.053 3.0415
{'A6'} 58.384 58.382

The dual-tree transform produces a consistent analysis of variance by scale for the original signal and its circularly shifted version.

The standard implementation of the DWT in 2-D uses separable filtering of the columns and rows of the image. The LH, HL, and HH wavelets for Daubechies' least-asymmetric phase wavelet with 4 vanishing moments (sym4) are shown in the following figure.

figure J = 5; L = 3*2^(J+1); N = L/2^J; Y = zeros(L,3*L); dt = dddtree2('dwt',Y,J,'sym4'); dt.cfs{J}(N/3,N/2,1) = 1; dt.cfs{J}(N/2,N/2+N,2) = 1; dt.cfs{J}(N/2,N/2+2*N,3) = 1; dwtImage = idddtree2(dt); imagesc(dwtImage) axis xy axis off title({'Critically Sampled DWT';'2-D separable wavelets (sym4) -- LH, HL, HH'})

Note that the LH and HL wavelets have clear horizontal and vertical orientations respectively. However, the HH wavelet on the far right mixes both the +45 and -45 degree directions, producing a checkerboard artifact. This mixing of orientations is due to the use of real-valued separable filters. The HH real-valued separable filter has passbands in all four high frequency corners of the 2-D frequency plane.

The dual-tree DWT achieves directional selectivity by using wavelets that are approximately analytic, meaning that they have support on only one half of the frequency axis. In the dual-tree DWT, there are six subbands for both the real and imaginary parts. The six real parts are formed by adding the outputs of column filtering followed by row filtering of the input image in the two trees. The six imaginary parts are formed by subtracting the outputs of column filtering followed by row filtering.

The filters applied to the columns and rows may be from the same filter pair, $$\{{h}_{0},{h}_{1}\}$$ or $$\{{g}_{0},{g}_{1}\}$$, or from different filter pairs, $$\{{h}_{0},{g}_{1}\},\{{g}_{0},{h}_{1}\}$$. The following code shows the orientation of the 12 wavelets corresponding to the real and imaginary parts of the complex oriented dual-tree DWT.

J = 4; L = 3*2^(J+1); N = L/2^J; Y = zeros(2*L,6*L); wt = dddtree2('cplxdt',Y,J,'dtf3'); wt.cfs{J}(N/2,N/2+0*N,2,2,1) = 1; wt.cfs{J}(N/2,N/2+1*N,3,1,1) = 1; wt.cfs{J}(N/2,N/2+2*N,1,2,1) = 1; wt.cfs{J}(N/2,N/2+3*N,1,1,1) = 1; wt.cfs{J}(N/2,N/2+4*N,3,2,1) = 1; wt.cfs{J}(N/2,N/2+5*N,2,1,1) = 1; wt.cfs{J}(N/2+N,N/2+0*N,2,2,2) = 1; wt.cfs{J}(N/2+N,N/2+1*N,3,1,2) = 1; wt.cfs{J}(N/2+N,N/2+2*N,1,2,2) = 1; wt.cfs{J}(N/2+N,N/2+3*N,1,1,2) = 1; wt.cfs{J}(N/2+N,N/2+4*N,3,2,2) = 1; wt.cfs{J}(N/2+N,N/2+5*N,2,1,2) = 1; waveIm = idddtree2(wt); imagesc(waveIm) axis off title('Complex Dual-Tree 2-D Wavelets')

The top row of the preceding figure shows the six directional wavelets of the real oriented dual-tree wavelet transform. The second row shows the imaginary parts. Together the real and imaginary parts form the complex oriented dual-tree wavelet transform. The real and imaginary parts are oriented in the same direction. You can use `dddtree2`

with the `'realdt'`

option to obtain the real oriented dual-tree DWT, which uses only the real parts. Using the real oriented dual-tree wavelet transform, you can achieve directional selectivity, but you do not gain the full benefit of using analytic wavelets such as approximate shift invariance.

The approximate analyticity and selective orientation of the complex dual-tree wavelets provide superior performance over the standard 2-D DWT in the representation of edges in images. To illustrate this, we analyze test images with edges consisting of line and curve singularities in multiple directions using the critically sampled 2-D DWT and the 2-D complex oriented dual-tree transform. First, analyze an image of an octagon, which consists of line singularities.

load woctagon figure imagesc(woctagon) colormap gray title('Original Image') axis equal axis off

Decompose the image down to level 4 and reconstruct an image approximation based on the level-4 detail coefficients.

dtcplx = dddtree2('cplxdt',woctagon,4,'dtf3'); dtDWT = dddtree2('dwt',woctagon,4,'farras'); dtcplx.cfs{1} = zeros(size(dtcplx.cfs{1})); dtcplx.cfs{2} = zeros(size(dtcplx.cfs{2})); dtcplx.cfs{3} = zeros(size(dtcplx.cfs{3})); dtcplx.cfs{5} = zeros(size(dtcplx.cfs{5})); dtDWT.cfs{1} = zeros(size(dtDWT.cfs{1})); dtDWT.cfs{2} = zeros(size(dtDWT.cfs{2})); dtDWT.cfs{3} = zeros(size(dtDWT.cfs{3})); dtDWT.cfs{5} = zeros(size(dtDWT.cfs{5})); dtImage = idddtree2(dtcplx); dwtImage = idddtree2(dtDWT); subplot(1,2,1) imagesc(dtImage) axis equal axis off colormap gray title('Complex Oriented Dual-Tree') subplot(1,2,2) imagesc(dwtImage) axis equal axis off colormap gray title('DWT')

Next, analyze an octagon with hyperbolic edges. The edges in the hyperbolic octagon are curve singularities.

load woctagonHyperbolic figure imagesc(woctagonHyperbolic) colormap gray title('Octagon with Hyperbolic Edges') axis equal axis off

Again, decompose the image down to level 4 and reconstruct an image approximation based on the level-4 detail coefficients for both the standard 2-D DWT and the complex oriented dual-tree DWT.

dtcplx = dddtree2('cplxdt',woctagonHyperbolic,4,'dtf3'); dtDWT = dddtree2('dwt',woctagonHyperbolic,4,'farras'); dtcplx.cfs{1} = zeros(size(dtcplx.cfs{1})); dtcplx.cfs{2} = zeros(size(dtcplx.cfs{2})); dtcplx.cfs{3} = zeros(size(dtcplx.cfs{3})); dtcplx.cfs{5} = zeros(size(dtcplx.cfs{5})); dtDWT.cfs{1} = zeros(size(dtDWT.cfs{1})); dtDWT.cfs{2} = zeros(size(dtDWT.cfs{2})); dtDWT.cfs{3} = zeros(size(dtDWT.cfs{3})); dtDWT.cfs{5} = zeros(size(dtDWT.cfs{5})); dtImage = idddtree2(dtcplx); dwtImage = idddtree2(dtDWT); subplot(1,2,1) imagesc(dtImage) axis equal axis off colormap gray title('DT-CWT') subplot(1,2,2) imagesc(dwtImage) axis equal axis off colormap gray title('DWT')

Note that the ringing artifacts evident in the 2-D critically sampled DWT are absent in the 2-D DT-CWT of both images. The DT-CWT more faithfully reproduces line and curve singularities.

Because of the ability to isolate distinct orientations in separate subbands, the dual-tree DWT is often able to outperform the standard separable DWT in applications like image denoising. To demonstrate this, use the helper function `helperCompare2DDenoising`

. The helper function loads an image and adds zero-mean white Gaussian noise with $$\sigma =25$$. For a user-supplied range of thresholds, the function compares denoising using soft thresholding for the critically sampled DWT, the real oriented dual-tree DWT, and the complex oriented dual-tree DWT. For each threshold value, the root-mean-square (RMS) error and peak signal-to-noise ratio (PSNR) are displayed.

```
numex = 3;
helperCompare2DDenoising(numex,0:2:100,'PlotMetrics');
```

Both the real oriented and the complex oriented dual-tree DWTs outperform the standard DWT in RMS error and PSNR.

Next, obtain the denoised images for a threshold value of 25, which is equal to the standard deviation of the additive noise.

```
numex = 3;
helperCompare2DDenoising(numex,25,'PlotImage');
```

With a threshold value equal to the standard deviation of the additive noise, the complex oriented dual-tree transform provides a PSNR almost 4 dB higher than the standard 2-D DWT.

The ringing artifacts observed with the separable DWT in two dimensions is exacerbated when extending wavelet analysis to higher dimensions. The DT-CWT enables you to maintain directional selectivity in 3-D with minimal redundancy. In 3-D, there are 28 wavelet subbands in the dual-tree transform.

To demonstrate the directional selectivity of the 3-D dual-tree wavelet transform, visualize example 3-D isosurfaces of both 3-D dual-tree and separable DWT wavelets. First, visualize the real and imaginary parts separately of two dual-tree subbands.

helperVisualize3D('Dual-Tree',28,'separate');

helperVisualize3D('Dual-Tree',25,'separate');

The red portion of the isosurface plot indicates the positive excursion of the wavelet from zero, while blue denotes the negative excursion. You can clearly see the directional selectivity in space of the real and imaginary parts of the dual-tree wavelets. Now visualize one of the dual-tree subbands with the real and imaginary plots plotted together as one isosurface.

helperVisualize3D('Dual-Tree',25,'real-imaginary');

The preceding plot demonstrates that the real and imaginary parts are shifted versions of each other in space. This reflects the fact that the imaginary part of the complex wavelet is the approximate Hilbert transform of the real part. Next, visualize the isosurface of a real orthogonal wavelet in 3-D for comparison.

`helperVisualize3D('DWT',7);`

The mixing of orientations observed in the 2-D DWT is even more pronounced in 3-D. Just as in the 2-D case, the mixing of orientations in 3-D leads to significant ringing, or blocking artifacts. To demonstrate this, examine the 3-D DWT and DT-CWT wavelet details of a spherical volume. The sphere is 64-by-64-by-64.

load sphr [A,D] = dualtree3(sphr,2,'excludeL1'); A = zeros(size(A)); sphrDTCWT = idualtree3(A,D); montage(reshape(sphrDTCWT,[64 64 1 64]),'DisplayRange',[]) title('DT-CWT Level 2 Details')

Compare the preceding plot against the second-level details based on the separable DWT.

sphrDEC = wavedec3(sphr,2,'sym4','mode','per'); sphrDEC.dec{1} = zeros(size(sphrDEC.dec{1})); for kk = 2:8 sphrDEC.dec{kk} = zeros(size(sphrDEC.dec{kk})); end sphrrecDWT = waverec3(sphrDEC); figure montage(reshape(sphrrecDWT,[64 64 1 64]),'DisplayRange',[]) title('DWT Level 2 Details')

Zoom in on the images in both the DT-CWT and DWT montages and you will see how prominent the blocking artifacts in the DWT details are compared to the DT-CWT.

Similar to the 2-D case, the directional selectivity of the 3-D DT-CWT often leads to improvements in volume denoising.

To demonstrate this, consider an MRI dataset consisting of 16 slices. Gaussian noise with a standard deviation of 10 has been added to the original dataset. Display the noisy dataset.

load MRI3D montage(reshape(noisyMRI,[128 128 1 16]),'DisplayRange',[])

Note that the original SNR prior to denoising is approximately 11 dB.

20*log10(norm(origMRI(:),2)/norm(origMRI(:)-noisyMRI(:),2))

ans = 11.2997

Denoise the MRI dataset down to level 4 using both the DT-CWT and the DWT. Similar wavelet filter lengths are used in both cases. Plot the resulting SNR as a function of the threshold. Display the denoised results for both the DT-CWT and DWT obtained at the best SNR.

[imrecDTCWT,imrecDWT] = helperCompare3DDenoising(origMRI,noisyMRI);

figure montage(reshape(imrecDTCWT,[128 128 1 16]),'DisplayRange',[]) title('DT-CWT Denoised Volume')

figure montage(reshape(imrecDWT,[128 128 1 16]),'DisplayRange',[]) title('DWT Denoised Volume')

We have shown that the dual-tree DWT possesses the desirable properties of near shift invariance and directional selectivity not achievable with the critically sampled DWT. We have demonstrated how these properties can result in improved performance in signal analysis, the representation of edges in images and volumes, and image and volume denoising.

[1] Chen, H., and N. G.
Kingsbury. "Efficient registration of nonrigid 3-D bodies." *IEEE
Transactions on Image Processing*. Vol. 21, Number 1, January
2012, pp. 262–272.

[2] Kingsbury, N. G.
"Complex Wavelets for Shift Invariant Analysis and Filtering of Signals."
*Journal of Applied and Computational Harmonic
Analysis*. Vol. 10, Number 3, May 2001, pp.
234–253.

[3] Percival, D. B., and
A. T. Walden. *Wavelet Methods for Time Series Analysis*.
Cambridge, UK: Cambridge University Press, 2000.

[4] Selesnick, I.,
Baraniuk, R. G., and N. G. Kingsbury. "The Dual-Tree Complex Wavelet Transform."
*IEEE Signal Processing Magazine*. Vol. 22, Number 6,
November 2005, pp. 123–151.

[5] Selesnick, I. "The
Double Density DWT." *Wavelets in Signal and Image Analysis: From
Theory to Practice* (A. A. Petrosian, and F. G. Meyer, eds.).
Norwell, MA: Kluwer Academic Publishers, 2001, pp. 39–66.

[6] Selesnick, I. "The
Double-Density Dual-Tree Wavelet Transform." *IEEE Transactions on
Signal Processing*. Vol. 52, Number 5, May 2004, pp.
1304–1314.

`dddtree`

| `dddtree2`

| `dualtree3`

| `idddtree`

| `idddtree2`

| `idualtree3`