From Wavelets to Wavelet Packets

In the orthogonal wavelet decomposition procedure, the generic step splits the approximation coefficients into two parts. After splitting we obtain a vector of approximation coefficients and a vector of detail coefficients, both at a coarser scale. The information lost between two successive approximations is captured in the detail coefficients. Then the next step consists of splitting the new approximation coefficient vector; successive details are never reanalyzed.

In the corresponding wavelet packet situation, each detail coefficient vector is also decomposed into two parts using the same approach as in approximation vector splitting. This offers the richest analysis: the complete binary tree is produced as shown in the following figure.

The idea of this decomposition is to start from a scale-oriented decomposition, and then to analyze the obtained signals on frequency subbands.

Wavelet Packets in Action: An Introduction

This example illustrates certain differences between wavelet analysis and wavelet packet analysis.

Fs = 1e3;
frq = 128;
t = 0:1/Fs:2;
sig = sin(2*pi*frq*t);

figure
level = 6;
wpt = wpdec(sig,level,"sym8");
[spec,tm,frq] = wpspectrum(wpt,Fs,"plot");

windowsize = 128;
window = hanning(windowsize);
noverlap = windowsize-1;
stft(sig,Fs,Window=window,OverlapLength=noverlap, ...
    FrequencyRange="onesided")
colorbar off
colormap("default")

Fs = 1e3;
frq1 = 64;
frq2 = 128;
t = 0:1/Fs:2;
sig = sin(2*pi*frq1*t) + ...
    sin(2*pi*frq2*t);

clf
wpt = wpdec(sig,level,"sym8");
[Spec,Time,Freq] = wpspectrum(wpt,Fs,"plot");

stft(sig,Fs,Window=window,OverlapLength=noverlap, ...
    FrequencyRange="onesided")
colorbar off
colormap("default")

Fs = 500;
frq1 = 16;
frq2 = 64;
t = 0:1/Fs:4;
sig = sin(2*pi*frq1*t).*(t<2) + ...
    sin(2*pi*frq2*t).*(t>=2);

clf
wpt = wpdec(sig,level,"sym8");
[Spec,Time,Freq] = wpspectrum(wpt,Fs,"plot");

stft(sig,Fs,Window=window,OverlapLength=noverlap, ...
    FrequencyRange="onesided")
colorbar off
colormap("default")

Fs = 1e3;
frq = 128;
t = 0:1/Fs:2;
sig = sin(2*pi*frq*t.^2);

clf
wpt = wpdec(sig,level,"sym8");
[Spec,Time,Freq] = wpspectrum(wpt,Fs,"plot");

stft(sig,Fs,Window=window,OverlapLength=noverlap, ...
    FrequencyRange="onesided")
colorbar off
colormap("default")

sig = wnoise("quadchirp",10);
len = length(sig);
t = linspace(0,5,len);
Fs = 1/t(2);

clf
wpt = wpdec(sig,level,"sym8");
[Spec,Time,Freq] = wpspectrum(wpt,Fs,"plot");

stft(sig,Fs,Window=window,OverlapLength=noverlap, ...
    FrequencyRange="onesided")
colorbar off
colormap("default")

Building Wavelet Packets

The computational scheme for wavelet packets generation is easy when using an orthogonal wavelet. We start with the two filters of length $$2N$$, $$h(n)$$ and $$g(n)$$, that correspond to the wavelet.

Now by induction, let us define the sequence of functions $$\{W_n(x),n=0,1,2,\dots \}$$ by:

where $$W_0(x)=\varphi (x)$$ is the scaling function and $$W_1(x)=\psi (x)$$ is the wavelet function.

For example, for the Haar wavelet we have

and

The equations become

and

$$W_0(x)=\varphi (x)$$ is the Haar scaling function, and $$W_1(x)=\psi (x)$$ is the Haar wavelet. Both functions are supported in the closed interval [-1,1]. Then we can obtain $$W_{2n}$$ by adding half-scaled versions of $$W_n$$ with distinct supports [0,1/2] and [1/2,1].

Use the wpfun function to plot the $$W$$-functions associated with the Haar wavelet on a 1/2^5 grid.

[wfun,xgrid] = wpfun("db1",7,5);
t = tiledlayout(2,4);
for k=0:7
    nexttile
    plot(xgrid,wfun(k+1,:),LineWidth=2)
    title("W"+num2str(k))
    ylim([-2 2])
end
title(t,"Haar Wavelet Packets")

Starting from more regular original wavelets and using a similar construction, we obtain smoothed versions of this system of W-functions, all with support in the interval [0, 2N–1].

Use the wpfun function to plot the $$W$$-functions associated with the db2 wavelet.

[wfun,xgrid] = wpfun("db2",7,5);
t = tiledlayout(2,4);
for k=0:7
    nexttile
    plot(xgrid,wfun(k+1,:),LineWidth=2)
    title("W"+num2str(k))
    ylim([-3 3])
end
title(t,"db2 Wavelet Packets")

Wavelet Packet Atoms

Starting from the functions $$\{W_n(x),n\in \mathbb{N}\}$$ and following the same line leading to orthogonal wavelets, we consider the three-indexed family of analyzing functions (the waveforms):

where $$n\in \mathbb{N}$$ and $$(j,k)\in {\mathbb{Z}}^2$$.

As in the wavelet framework, k can be interpreted as a time-localization parameter and j as a scale parameter. So what is the interpretation of n?

The basic idea of wavelet packets is that for fixed values of j and k, W_j,n,k analyzes the fluctuations of the signal roughly around the position 2^j· k, at the scale 2^j, and at various frequencies for the different admissible values of the last parameter n.

In fact, examining carefully the wavelet packets displayed in Haar Wavelet Packets and db2 Wavelet Packets, the naturally ordered W_n for n = 0, 1, ..., 7, does not match exactly the order defined by the number of oscillations. More precisely, counting the number of zero crossings (up-crossings and down-crossings) for the db1 wavelet packets, we have the following.

So, to restore the property that the main frequency increases monotonically with the order, it is convenient to define the frequency order obtained from the natural one recursively.

As can be seen in the previous figures, W_r(n)(x) “oscillates” approximately n times.

To analyze a signal (the Linear Chirp, for instance), it is better to plot the wavelet packet coefficients following the frequency order from the low frequencies at the bottom to the high frequencies at the top, rather than naturally ordered coefficients.

To plot the coefficients in a wavelet packet tree, use the wpviewcf function.

Organizing the Wavelet Packets

The set of functions $$W_{j,n}=\{W_{j,n,k}(x),k\in \mathbb{Z}\}$$ is the (j,n) wavelet packet. For positive values of integers j and n, wavelet packets are organized in trees. The tree in the figure Wavelet Packets Organized in a Tree; Scale j Defines Depth and Frequency n Defines Position in the Tree is created to give a maximum level decomposition equal to 3. For each scale j, the possible values of parameter n are 0, 1, ..., 2^j–1.

The notation W_j,n, where j denotes the scale parameter and n the frequency parameter, is consistent with the usual depth-position tree labeling.

We have $$W_{0,0}=\{\varphi (x-k),k\in \mathbb{Z}\}$$, and $$W_{1,1}=\left\{\psi \left(\frac{x}{2}-k\right),k\in \mathbb{Z}\right\}$$.

It turns out that the library of wavelet packet bases contains the wavelet basis and several other bases. Let us have a look at some of those bases. More precisely, let V₀ denote the space (spanned by the family W_0,0) in which the signal to be analyzed lies; then {W_d,1; d ≥ 1} is an orthogonal basis of V₀.

For every strictly positive integer D, {W_D,0, {W_d,1; 1 ≤ d ≤ D}} is an orthogonal basis of V₀.

We also know that the family of functions {(W_j+1,2n), (W_j+1,2n+1)} is an orthogonal basis of the space spanned by W_j,n, which is split into two subspaces: W_j+1,2n spans the first subspace, and W_j+1,2n+1 the second one.

This last property gives a precise interpretation of splitting in the wavelet packet organization tree, because all the developed nodes are of the form shown in the figure Wavelet Packet Tree: Split and Merge.

It follows that the leaves of every connected binary subtree of the complete tree correspond to an orthogonal basis of the initial space.

For a finite energy signal belonging to V₀, any wavelet packet basis will provide exact reconstruction and offer a specific way of coding the signal, using information allocation in frequency scale subbands.

Choosing Optimal Decomposition

Based on the organization of the wavelet packet library, it is natural to count the decompositions issued from a given orthogonal wavelet.

A signal of length $$N=2^L$$ can be expanded in $$\alpha$$ different ways, where $$\alpha$$ is the number of binary subtrees of a complete binary tree of depth L. As a result, $$\alpha \ge 2^{N/2}$$ (see [1]).

As this number may be very large, and since explicit enumeration is generally unmanageable, it is interesting to find an optimal decomposition with respect to a convenient criterion, computable by an efficient algorithm. We are looking for a minimum of the criterion.

Functions satisfying an additivity-type property are well suited for efficient searching of binary-tree structures and the fundamental splitting. Classical entropy-based criteria match these conditions and describe information-related properties for an accurate representation of a given signal. Entropy is a common concept in many fields, mainly in signal processing. Let us list four different entropy criteria (see [2]); many others are available and can be easily integrated (see wpdec). In the following expressions s is the signal and (si) are the coefficients of s in an orthonormal basis.

The entropy $$E$$ must be an additive cost function with that $$E(0)=0$$ and $$E(s)=\sum _{i}E(s_i)$$.

s = ones(1,16)*0.25;

e1 = -sum((s.^2).*log(s.^2))

e1 = 
2.7726

e2 = norm(s,1.5)^1.5

e2 = 
2.0000

e3 = sum(log(s.^2))

e3 = 
-44.3614

e4 = sum(abs(s)>0.24)

e4 = 
16

shannon = @(x)(-sum((x.^2).*log(x.^2)));

w00 = ones(1,16)*0.25;

e00 = shannon(w00)

e00 = 
2.7726

[w10,w11] = dwt(w00,"db1");

e10 = shannon(w10)

e10 = 
2.0794

[w20,w21] = dwt(w10,"db1");

e20 = shannon(w20)

e20 = 
1.3863

[w30,w31] = dwt(w20,"db1");
e30 = shannon(w30)

e30 = 
0.6931

[w40,w41] = dwt(w30,"db1")

w40 = 
1.0000

w41 = 
0

e40 = shannon(w40)

e40 = 
-4.4409e-16

t = wpdec(w00,4,"haar","shannon");

helperPlotTreeWithLabels([e00,e10,e20,e30])

bt = besttree(t);
plot(bt)

Some Interesting Subtrees

Using wavelet packets requires tree-related actions and labeling. The implementation of the user interface is built around this consideration. For more information on the technical details, see the reference pages.

The complete binary tree of depth D corresponding to a wavelet packet decomposition tree developed at level D is denoted by WPT.

We have the following interesting subtrees.

We deduce the following definitions of optimal decompositions, with respect to an entropy criterion E.

For any nonterminal node, we use the following basic step to find the optimal subtree with respect to a given entropy criterion E (where Eopt denotes the optimal entropy value).

with the natural initial condition on the reference tree, Eopt(t) = E(t) for each terminal node t.

load noisdopp
origMode = dwtmode("status","nodisp");

dwtmode("per","nodisp")
T = wpdec(noisdopp,5,"sym4");
plot(T)

wpc = wpcoef(T,16);
length(wpc)

ans = 
64

rwpc = wprcoef(T,16);
isequal(length(rwpc),length(noisdopp))

ans = logical
   1

plot(noisdopp,"k")
hold on
plot(rwpc,"b",LineWidth=2)
axis tight
legend("Original","Approximation")
hold off

Topt = besttree(T);
plot(Topt)

rsig = wprcoef(Topt,7);
plot(noisdopp,"k")
hold on
plot(rsig,"b",LineWidth=2)
axis tight
legend("Original","Approximation")
hold off

Node = depo2ind(2,[3 0])

Node = 
7

dwtmode(origMode,"nodisp")

Wavelet Packets 2-D Decomposition Structure

Exactly as in the wavelet decomposition case, the preceding 1-D framework can be extended to image analysis. Minor direct modifications lead to quaternary tree-related definitions. An example is shown the following figure for depth 2.

Wavelet Packets for Compression and Denoising

In the wavelet packet framework, compression and denoising ideas are identical to those developed in the wavelet framework. The only new feature is a more complete analysis that provides increased flexibility. A single decomposition using wavelet packets generates a large number of bases. You can then look for the best representation with respect to a design objective, using the besttree function with an entropy function.