## Dual-Tree Complex Wavelet Transform

1. Introduction

It turns out that, for some applications of the discrete wavelet transform, improvements can be obtained by using an expansive wavelet transform in place of a critically-sampled one. (An expansive transform is one that converts an N-point signal into M coefficients with M > N.) There are several kinds of expansive DWTs; here we describe the dual-tree complex discrete wavelet transform.

The dual-tree complex DWT of a signal x is implemented using two critically-sampled DWTs in parallel on the same data, as shown in the figure.

The transform is 2-times expansive because for an N-point signal it gives 2N DWT coefficients. If the filters in the upper and lower DWTs are the same, then no advantage is gained. However, if the filters are designed is a specific way, then the subband signals of the upper DWT can be interpreted as the real part of a complex wavelet transform, and subband signals of the lower DWT can be interpreted as the imaginary part. Equivalently, for specially designed sets of filters, the wavelet associated with the upper DWT can be an approximate Hilbert transform of the wavelet associated with the lower DWT. When designed in this way, the dual-tree complex DWT is nearly shift-invariant, in contrast with the critically-sampled DWT. Moreover, the dual-tree complex DWT can be used to implement 2D wavelet transforms where each wavelet is oriented, which is especially useful for image processing. (For the separable 2D DWT, recall that one of the three wavelets does not have a dominant orientation.) The dual-tree complex DWT outperforms the critically-sampled DWT for applications like image denoising and enhancement.

2. Filters for the Dual-Tree Complex DWT

For a technical reason we skip here, the filters used in the first stage of the dual-tree complex DWT should be different from the filters used in the remaining stages. One set of filters for the first stage are given by the function FSfarras.m. (FS for First Stage.) These filters are the same as those of farras.m, but they have specific alignments with respect to each other.

Table 4.1 Matlab function FSfarras.m

function [af, sf] = FSfarras

af{1} = [
0                  0
-0.08838834764832  -0.01122679215254
0.08838834764832   0.01122679215254
0.69587998903400   0.08838834764832
0.69587998903400   0.08838834764832
0.08838834764832  -0.69587998903400
-0.08838834764832   0.69587998903400
0.01122679215254  -0.08838834764832
0.01122679215254  -0.08838834764832
0                  0
];

sf{1} = af{1}(end:-1:1, :);

af{2} = [
0.01122679215254                  0
0.01122679215254                  0
-0.08838834764832  -0.08838834764832
0.08838834764832  -0.08838834764832
0.69587998903400   0.69587998903400
0.69587998903400  -0.69587998903400
0.08838834764832   0.08838834764832
-0.08838834764832   0.08838834764832
0   0.01122679215254
0  -0.01122679215254
];

sf{2} = af{2}(end:-1:1, :);

The function dualfilt1.m provides analysis and synthesis filters for the remaining stages of the dual-tree complex DWT. These filters were published in [KIN-1999].

Table 4.2 Matlab function dualfilt1.m

function [af, sf] = dualfilt1

% [af, sf] = dualfilt1
% filters for the 1D dual-tree wavelet transform.
% af{i}: analysis filters for tree i
% sf{i}: synthesis filters for tree i
%
% note: af{2} is the reverse of af{1}

af{1} = [
0.03516384000000                  0
0                  0
-0.08832942000000  -0.11430184000000
0.23389032000000                  0
0.76027237000000   0.58751830000000
0.58751830000000  -0.76027237000000
0   0.23389032000000
-0.11430184000000   0.08832942000000
0                  0
0  -0.03516384000000
];

af{2} = [
0  -0.03516384000000
0                  0
-0.11430184000000   0.08832942000000
0   0.23389032000000
0.58751830000000  -0.76027237000000
0.76027237000000   0.58751830000000
0.23389032000000                  0
-0.08832942000000  -0.11430184000000
0                  0
0.03516384000000                  0
];

sf{1} = af{1}(end:-1:1, :);

sf{2} = af{2}(end:-1:1, :);

3. Matlab Implementation

The Matlab function, dualtree.m, computes the J-scale dual-tree complex DWT w of a signal x. This function is very similar to the function dwt.m. It repeatedly calls the function afb.m to implement the analysis filter bank. The wavelet coefficients are also stored as a cell array w. For j = 1..J, w{j}{1} is the high frequency subband signal produced at stage j in the upper DWT, w{j}{2} is the high frequency subband signal produced at stage j in the lower DWT. The function idualtree.m implements the inverse transform.

Table 4.3 Matlab function dualtree.m

function w = dualtree(x, J, Faf, af)

% Dual-tree Complex Discrete Wavelet Transform
%
% w = dualtree1D(x, J, Faf, af)
%
% OUTPUT
%   w{j}: scale j (j = 1..J)
%   w{j}{i}: tree i (i = 1,2)
%
% INPUT
%   x   : 1-D signal
%   J   : number of stages
%   Faf : first stage filter
%   af  : filter for remaining stages

% normalization
x = x/sqrt(2);

% Tree 1
[x1 w{1}{1}] = afb(x, Faf{1});
for k = 2:J
[x1 w{k}{1}] = afb(x1, af{1});
end
w{J+1}{1} = x1;

% Tree 2
[x2 w{1}{2}] = afb(x, Faf{2});
for k = 2:J
[x2 w{k}{2}] = afb(x2, af{2});
end
w{J+1}{2} = x2;

The perfect reconstruction property of the dual-tree complex DWT is verified in the following example. We create a random input signal x of length 512, apply the dual-tree complex DWT and its inverse. We then show that the reconstructed signal y is equal to the original signal x by computing the maximum value of |x-y|.

#### EXAMPLE (VERIFY PERFECT RECONSTRUNCTION)

>> x = rand(1, 512);      % Test signal
>> J = 4;                 % number of stages
>> [Faf, Fsf] = FSfarras; % 1st stage anal. & synth. filters
>> [af, sf] = dualfilt1;
>> w = dualtree(x, J, Faf, af);
>> y = idualtree(w, J, Fsf, sf);
>>  err = x - y;
>> max(abs(err))

ans =

2.1651e-008


The complex wavelet associated with the dual-tree complex DWT can be computed using the following Matlab code fragment. To compute the real part of the complex wavelet, we set all coefficients to zero, except for one coefficient in the upper DWT, and then compute the inverse transform. To compute the imaginary part, the exception is a single coefficient in the lower DWT.

This figure was obtained using the following code fragment. Different sets of filters will produce different wavelets.

Table 4.4 Matlab function dualtree_eg1.m

J = 5;  % J: number of stages

% get filters
[Faf, Fsf] = FSfarras;
[af, sf] = dualfilt1;

x = zeros(1,256);  % zero signal

% Compute dual-tree complex DWT of zero signal
w = dualtree(x, J, Faf, af);
% Set a single (real) coefficient to 1
w{5}{1}(4) = 1;
% Compute the inverse transform
y1 = idualtree(w, J, Fsf, sf);

% Compute dual-tree complex DWT of zero signal
w = dualtree(x, J, Faf, af);
% Set a single (imaginary) coefficient to 1
w{5}{2}(4) = 1;
% Compute the inverse transform
y2 = idualtree(w, J, Fsf, sf);

% Display real and imaginary parts and magnitude
n = [1:256]/256;
plot(n,y1,n,y2,n,sqrt(y1.^2+y2.^2))
title('Complex 1-D wavelet')
xlabel('t');
ylabel('\psi(t)');


Summary of programs:

1. FSfarras.m - Compute filters for stage 1
2. dualfilt1.m - Compute filters for remaining stages
3. dualtree.m - forward 1-D dual-tree DWT
4. idualtree.m - inverse 1-D dual-tree DWT

Introduction
Separable DWT

Dual-tree DWT

Denoising

References
People