DualTree Complex Wavelet Transform1. 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 criticallysampled
one. (An expansive transform is one that converts an Npoint
signal into M coefficients with M > N.)
There are several kinds of expansive DWTs;
here we describe the dualtree complex discrete wavelet
transform.
The transform is 2times expansive because for an Npoint 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 dualtree complex DWT is nearly shiftinvariant, in contrast with the criticallysampled DWT. Moreover, the dualtree 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 dualtree complex DWT outperforms the criticallysampled DWT for applications like image denoising and enhancement. 2. Filters for the DualTree Complex DWT For a technical reason we skip here, the filters used in the first stage of the dualtree 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 dualtree complex DWT. These filters were published in [KIN1999].
Table 4.2 Matlab function dualfilt1.m
function [af, sf] = dualfilt1 % [af, sf] = dualfilt1 % filters for the 1D dualtree 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 Jscale dualtree 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) % Dualtree 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 : 1D 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 dualtree complex DWT is verified in the following example. We create a random input signal x of length 512, apply the dualtree 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 xy.
The complex wavelet associated with the dualtree 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 dualtree 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 dualtree 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 1D wavelet') xlabel('t'); ylabel('\psi(t)'); Summary of programs:
