The MATLAB function, doubledualtree_f1D.m, computes the double-density complex DWT, w, of an input signal x. This function is similar to the forward double-density DWT, double_f1D.m, in that it repeatedly calls the afb3.m function to implement the analysis filter bank. MATLAB's cell array data structure is used to store each of the subband signals such that for j = 1 - J, w{j}{1} is the high frequency subband signals produced in the upper DWT at stage j and w{j}{2} is the high frequency subband signals produced in the lower DWT at stage j.

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

% Forward 1-D Double-Density Dual-Tree Complex DWT
%
% USAGE:
%    w = doubledualtree_f1D(x, J, Faf, af)
% INPUT:
%   x - 1-D signal
%   J - number of stages
%   Faf - filters for the first stage 
%   af - filters for the remaining stages
% OUTPUT:
%   w - DWT coefficients
%      w{j}{1}, j = 1..J - real part 
%      w{j}{2}, j = 1..J - imaginary part 
%      w{J+1}{d} - lowpass coefficients, d = 1,2

% normalization
x = x/sqrt(2);

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

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

The inverse double-density complex DWT is computed with the doubledualtree_i1D.m MATLAB function. To verify the perfect reconstruction condition, we write another small piece of MATLAB code as shown below.

>> x = rand(1, 512);	      		% create generic 1-D test signal
>> J = 4;	                 	% number of stages (4)
>> [Faf, Fsf] = FSdoubledualfilt; 	% filters for the first stage
>> [af, sf] = doubledualfilt;		% filters for the remaining stages
>> w = doubledualtree_f1D(x, J, Faf, af);% analysis filter banks
>> y = doubledualtree_i1D(w, J, Fsf, sf);% synthesis filter banks
>> err = x - y;				% compute error signal
>> max(abs(err))			% verify that error is zero

ans =

  3.0898e-013

First, we generate a random 1-D input signal x of length 512. The analysis and synthesis filter bank is then computed first from the FSdoubledualfilt.m function and then from the doubledualfilt.m function. The 4-scale forward double-density complex DWT is computed with the doubledualtree_f1D.m function to generate the forward transform w. We then apply the inverse double-density complex DWT to w using the doubledualtree_i1D.m MATLAB function. To ensure that the perfect reconstruction condition has been met, we then find the difference between the input signal x and the output signal y. As you can see, the differences found here are again significantly close to zero, and can actually be considered as zero.

The complex wavelets associated with the double-density complex DWT can be computed using the following MATLAB code fragment. To compute the real part of the complex wavelet, we set all the 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.

>> J = 5;  				% number of stages
>> [Faf, Fsf] = FSdoubledualfilt;	% filters for first stage
>> [af, sf] = doubledualfilt;		% filters for remaining stages
>> x = zeros(1,256);  			% create 1-D zero signal
>> w = doubledualtree_f1D(x, J, Faf, af); % compute forward transform
>> w{6}{1}(4) = 1;			% set one real coefficient to 1
>> y1 = doubledualtree_i1D(w, J, Fsf, sf);% compute the inverse transform 
>> w = doubledualtree_f1D(x, J, Faf, af); % compute forward transform
>> w{6}{2}(4) = 1;			% set one imaginary coefficient to 1
>> y2 = doubledualtree_i1D(w, J, Fsf, sf);% compute the inverse transform 
>> n = [1:256]/256;			% display Re{w},Im{w}, & mag.
>> figure(1)				% display the lowpass filter
>> plot(n,y1,n,y2,n,sqrt(y1.^2+y2.^2))
>> title('Complex 1-D Scaling Function') 
>> xlabel('t');
>> ylabel('\phi(t)');

>> J = 5;  				% number of stages
>> [Faf, Fsf] = FSdoubledualfilt;	% filters for the first stage
>> [af, sf] = doubledualfilt;		% filters for the remaining stages
>> x = zeros(1,256);  			% create a 1-D zero signal
>> w = doubledualtree_f1D(x, J, Faf, af); % compute forward transform
>> w{5}{1}{1}(4) = 1;			% set one real coefficient to 1
>> y1 = doubledualtree_i1D(w, J, Fsf, sf);% compute the inverse transform 
>> w = doubledualtree_f1D(x, J, Faf, af); % compute forward transform
>> w{5}{2}{1}(4) = 1;			% set one imaginary coefficient to 1
>> y2 = doubledualtree_i1D(w, J, Fsf, sf);% compute the inverse transform 
>> n = [1:256]/256;			% display Re{w},Im{w}, & mag.
>> figure(2)				% display the 1st highpass filter
>> plot(n,y1,n,y2,n,sqrt(y1.^2+y2.^2))
>> title('Complex 1-D Wavelet') 
>> xlabel('t');
>> ylabel('\psi(t)');

>> J = 5;  				% number of stages
>> [Faf, Fsf] = FSdoubledualfilt;	% filters for the first stage
>> [af, sf] = doubledualfilt;		% filters for the remaining stages
>> x = zeros(1,256);  			% create a 1-D zero signal
>> w = doubledualtree_f1D(x, J, Faf, af); % compute forward transform
>> w{5}{1}{2}(4) = 1;			% set one real coefficient to 1
>> y1 = doubledualtree_i1D(w, J, Fsf, sf);% compute the inverse transform 
>> w = doubledualtree_f1D(x, J, Faf, af); % compute forward transform
>> w{5}{2}{2}(4) = 1;			% set one imaginary coefficient to 1
>> y2 = doubledualtree_i1D(w, J, Fsf, sf);% compute the inverse transform 
>> n = [1:256]/256;			% display Re{w},Im{w}, & mag.
>> figure(3)				% display the 2nd highpass filter
>> plot(n,y1,n,y2,n,sqrt(y1.^2+y2.^2))	
>> title('Complex 1-D Wavelet') 
>> xlabel('t');
>> ylabel('\psi(t)');

The following figure illustrates the eight wavelets generated from the above MATLAB code fragment.

Figure 10.   1-D Double-Density Complex Wavelets.

Summary of Programs:

1. FSdoubledualfilt.m - Set of filters for the first stage
2. doubledualfilt.m - Set of filters for the remaining stages
3. afb3.m - Analysis filter bank
4. sfb3.m - Synthesis filter bank
5. doubledualtree_f1D.m - Forward 1-D double-density complex wavelet transform
6. doubledualtree_i1D.m - Inverse 1-D double-density complex wavelet transform
7. cshift.m - Circular shift