# Why does MATLAB's dualtree3 function not return the lowpass subband of the imaginary tree?

MATLAB has a function named dualtree3 which computes 3-D dual-tree complex wavelet transform. This function only retains the last low-pass coefficients/subband of the real tree. And, it does not save the low-pass subband of the imaginary tree. How can you justify this?

• In [1], the authors say that when the input signal is real, there is a symmetry between imaginary and real parts. However, you should note that this symmetry exists for detail subbands (or highpass coefficients). This MATLAB's function seems to remove the computations for the lowpass coefficients. ---- [1] Selesnick, Ivan W., Richard G. Baraniuk, and Nick C. Kingsbury. "The dual-tree complex wavelet transform." IEEE signal processing magazine 22.6 (2005): 123-151. Commented Sep 10, 2020 at 4:16
• I don't clearly see the link between the title (complexity reduction) and the low-pass vs Hilbert behavior Commented Sep 14, 2020 at 17:45
• You are right @Laurent. At first, I thought that when the function does not return one of the two low-pass parts, it actually reduces the amount of computations. Now, I changed the title to better show my question. Commented Sep 15, 2020 at 8:06

TL,DR: the low-pass component (approximation coefficients, a) has a size bigger than expected ($$2^3$$ times). So I guess that the 8 avatars of the approximation subbands are gathered into one.

First, I did not take enough time to check the codes, so this may be a partial answer. While performing a dualtree3 decomposition:

zr = rand(64,64,64);
[a,d] = dualtree3(zr,2);


and its inversion:

xrec = idualtree3(a,d);


the $$\ell_\infty$$ norm is very low, typically of the order of $$10^{-15}$$. So, numerically speaking, the reconstruction is "perfect". As such, I won't go for the ideas that the low-pass approximations are "very similar". This would yield higher relative differences. It may not matter a lot for practical volume processing

Now, let us look at the subbands. For d, th details subbands, one gets:

  1×2 cell array
{32×32×32×28 double}    {16×16×16×28 double}


so you recover the $$28$$ directional details, that's okay. Their size are respectively $$2^3$$ and $$4^3$$ those of the original data. Now, look at the approximation:

a         32x32x32            262144  double


Oddly, the size does corresponds to what we could expect. Each of the low-pass approximations (8 in 3D) ought to be of size 16×16×16, similarly to the 2nd level of details. Instead, it looks like the 8 low-pass approximations are packed into a $$2^3$$ bigger array.

Finally, if you consider that the high-pass subbands are complex (hence have twice the array size), the global size of the 3D dual-tree, divided by the size of the orignal cube, is

(2*prod(size(d{1}))+2*prod(size(d{2}))+prod(size(a)))/prod(size(zr))


which yields $$8$$, as expected. I obtained the same results with three level of decomposition. To me, all the information is here, with the expected redundancy. No savings, unfortunately. Gathering approximation coefficients is not uncommon. It is rare in 2D (because one gets a shape factor). This is possible in 3D because the $$8$$ instances of approximations can nicely fit into a cube. I have no mention of that in the literature in mind.

Prof. Nick Kingsbury kindly provide an answer to my question!.

In 1-D, the lowpass basis functions (scaling functions) from the two trees (a) and (b) of the dual tree WT tend to look very similar to each other, apart from a shift of half the output sample period between them. Hence it usually makes most sense to regard the two sets of lowpass samples as interleaved samples from a single 2x over-sampled real basis function.

On the other hand, the highpass basis functions (wavelets) from (a) and (b) look approximately like a Hilbert pair of functions. so they are most logically interpreted as the real and imaginary parts of an analytic basis function.

Hence the over-sampling by 2 in the lowpass channel replaces the spitting into real + imaginary channels for the bandpass channels of a multi-scale wavelet transform. Similar things apply when 4 trees are used for the 2-D dual-tree WT, or when 8 trees are used for the 3-D version.

I hope this helps. Do refer to our 2005 paper in IEEE Signal Processing Magazine for more details.

[UPDATE at 9/18/2020]

Since I quoted Prof. Kingsubury's answer above, I should also note one more thing. He mentioned in his email that he is not familiar with dualtree3 function in MATLAB. However, based on my question, he gave the above general answer.

• Not fully convinced by now because the inversion result is close to perfect. I have updated my contribution Commented Sep 18, 2020 at 13:07
• I added an update. Still I am working on your answer. Thanks. Commented Sep 18, 2020 at 16:25

The @LaurentDuval is right.

The MATLAB's dualtree3 function returns all detail and approximation subbands which are expected when performing 3D dual-tree complex wavelet transform (cplxdual3D). However, instead of returning separate approximation subbands, dualtree3 function groups them all in a single cube.

In the following, I explain it by using both MATLAB's built-in dualtree3 function and cplxdual3D function from Prof. Ivan Selesnick's Wavelet toolbox.

## Example 1: One stage decomposition

We know that real 3D discrete wavelet transform (DWT) returns $$1$$ approximate subband and $$7$$ detail subbands.

In practice, cplxdual3D is implemented using eight real 3D discrete wavelet transforms. Therefore, after performing one stage of cplxdual3D, we expect to have $$1\times8$$ approximate subbands and $$7\times8$$ detail subbands.

J = 1; % number of decomposition stage
[Faf,Fsf] = FSfarras; % filterbanks used in the first stage
[af,sf] = dualfilt1; % filterbanks used in the remaining stages

zr = rand(64,64,64); % Input volume

w = cplxdual3D(zr,J,Faf,af);


This function returns all subbands in a cell array named w, in which the approximate subbands can be accessed by w{J+1}{m}{n}{p}, where $$m,n,p\in\{1,2\}$$.

For the above example:

m = 1; n = 1; p = 1;
size(w{J+1}{m}{n}{p})

ans =
32    32    32


So, we have $$32 \times 32 \times 32 \times 8$$ entries as approximate coefficients.

NOTE: I do not know why but MATLAB's dualtree3 function could not support one stage decomposition.

## Example 2: Two stage decomposition

J = 2
w = cplxdual3D(zr,J,Faf,af);

size(w{J+1}{1}{1}{1})

ans =
16    16    16


As mentioned in the above example, here we also have $$8$$ approximate subbands. So we have $$16 \times 16 \times 16 \times 8$$ entries as approximate coefficients.

If we run a similar experiment with MATLAB's dualtree3 function, we see that dualtree3 groups all approximate coefficients in a cube of size $$32 \times 32 \times 32$$:

J = 2;
[a,d] = dualtree3(zr,J);
size(a)

ans =
32    32    32


Additional Note: At this time, I do not know in which order dualtree3 puts the approximate coefficients into the cube. In order to find it, you should perform cplxdual3D and dualtree3 with the same filter banks and compare the results of decompositions for a small and simple example.

• Yes, this one level impossibility is weird. I should have thought it was from Selesnick word (the Farras filter). Commented Sep 19, 2020 at 8:21