Wavelet Packets
The wavelet packet method is a generalization of wavelet decomposition that offers a richer signal analysis.
Wavelet packet atoms are waveforms indexed by three naturally interpreted parameters: position, scale (as in wavelet decomposition), and frequency.
For a given orthogonal wavelet function, we generate a library of bases called wavelet packet bases. Each of these bases offers a particular way of coding signals, preserving global energy, and reconstructing exact features. The wavelet packets can be used for numerous expansions of a given signal. We then select the most suitable decomposition of a given signal with respect to an entropy-based criterion.
There exist simple and efficient algorithms for both wavelet packet decomposition and optimal decomposition selection. We can then produce adaptive filtering algorithms with direct applications in optimal signal coding and data compression.
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.
Wavelet Packet Decomposition Tree at Level 3
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
The following simple examples illustrate certain differences between wavelet analysis and wavelet packet analysis.
Wavelet Packet Spectrum
The spectral analysis of wide-sense stationary signals using the Fourier transform is well-established. For nonstationary signals, there exist local Fourier methods such as the short-time Fourier transform (STFT). See Short-Time Fourier Transform for a brief description.
Because wavelets are localized in time and frequency, it is
possible to use wavelet-based counterparts to the STFT for the time-frequency
analysis of nonstationary signals. For example, it is possible to
construct the scalogram (wscalogram
)
based on the continuous wavelet transform (CWT). However, a potential
drawback of using the CWT is that it is computationally expensive.
The discrete wavelet transform (DWT) permits a time-frequency decomposition of the input signal, but the degree of frequency resolution in the DWT is typically considered too coarse for practical time-frequency analysis.
As a compromise between the DWT- and CWT-based techniques, wavelet
packets provide a computationally-efficient alternative with sufficient
frequency resolution. You can use wpspectrum
to
perform a time-frequency analysis of your signal using wavelet packets.
The following examples illustrate the use of wavelet packets
to perform a local spectral analysis. The following examples also
use spectrogram
(Signal Processing Toolbox) from the Signal Processing Toolbox™ software
as a benchmark to compare against the wavelet packet spectrum. If
you do not have the Signal Processing Toolbox software, you can
simply run the wavelet packet spectrum examples.
Wavelet packet spectrum of a sine wave.
fs = 1000; % sampling rate t = 0:1/fs:2; % 2 secs at 1kHz sample rate y = sin(256*pi*t); % sine of period 128 level = 6; wpt = wpdec(y,level,'sym8'); [Spec,Time,Freq] = wpspectrum(wpt,fs,'plot');
If you have the Signal Processing Toolbox software, you can compute the short-time Fourier transform.
figure; windowsize = 128; window = hanning(windowsize); nfft = windowsize; noverlap = windowsize-1; [S,F,T] = spectrogram(y,window,noverlap,nfft,fs); imagesc(T,F,log10(abs(S))) set(gca,'YDir','Normal') xlabel('Time (secs)') ylabel('Freq (Hz)') title('Short-time Fourier Transform spectrum')
Sum of two sine waves with frequencies of 64 and 128 hertz.
fs = 1000; t = 0:1/fs:2; y = sin(128*pi*t) + sin(256*pi*t); % sine of periods 64 and 128. level = 6; wpt = wpdec(y,level,'sym8'); [Spec,Time,Freq] = wpspectrum(wpt,fs,'plot');
If you have the Signal Processing Toolbox software, you can compute the short-time Fourier transform.
figure; windowsize = 128; window = hanning(windowsize); nfft = windowsize; noverlap = windowsize-1; [S,F,T] = spectrogram(y,window,noverlap,nfft,fs); imagesc(T,F,log10(abs(S))) set(gca,'YDir','Normal') xlabel('Time (secs)') ylabel('Freq (Hz)') title('Short-time Fourier Transform spectrum')
Signal with an abrupt change in frequency from 16 to 64 hertz at two seconds.
fs = 500; t = 0:1/fs:4; y = sin(32*pi*t).*(t<2) + sin(128*pi*t).*(t>=2); level = 6; wpt = wpdec(y,level,'sym8'); [Spec,Time,Freq] = wpspectrum(wpt,fs,'plot');
If you have the Signal Processing Toolbox software, you can compute the short-time Fourier transform.
figure; windowsize = 128; window = hanning(windowsize); nfft = windowsize; noverlap = windowsize-1; [S,F,T] = spectrogram(y,window,noverlap,nfft,fs); imagesc(T,F,log10(abs(S))) set(gca,'YDir','Normal') xlabel('Time (secs)') ylabel('Freq (Hz)') title('Short-time Fourier Transform spectrum')
Wavelet packet spectrum of a linear chirp.
fs = 1000; t = 0:1/fs:2; y = sin(256*pi*t.^2); level = 6; wpt = wpdec(y,level,'sym8'); [Spec,Time,Freq] = wpspectrum(wpt,fs,'plot');
If you have the Signal Processing Toolbox software, you can compute the short-time Fourier transform.
figure; windowsize = 128; window = hanning(windowsize); nfft = windowsize; noverlap = windowsize-1; [S,F,T] = spectrogram(y,window,noverlap,nfft,fs); imagesc(T,F,log10(abs(S))) set(gca,'YDir','Normal') xlabel('Time (secs)') ylabel('Freq (Hz)') title('Short-time Fourier Transform spectrum')
Wavelet packet spectrum of quadratic chirp.
y = wnoise('quadchirp',10); len = length(y); t = linspace(0,5,len); fs = 1/t(2); level = 6; wpt = wpdec(y,level,'sym8'); [Spec,Time,Freq] = wpspectrum(wpt,fs,'plot');
If you have the Signal Processing Toolbox software, you can compute the short-time Fourier transform.
windowsize = 128; window = hanning(windowsize); nfft = windowsize; noverlap = windowsize-1; imagesc(T,F,log10(abs(S))) set(gca,'YDir','Normal') xlabel('Time (secs)') ylabel('Freq (Hz)') title('Short-time Fourier Transform spectrum')
Building Wavelet Packets
The computation scheme for wavelet packets generation is easy
when using an orthogonal wavelet. We start with the two filters of
length 2N, where h(n)
and g(n)
,
correspond to the wavelet.
Now by induction let us define the following sequence of functions:
(Wn(x), n = 0, 1, 2, ...)
by
where W0(x)
=
φ(x)
is the scaling function
and W1(x)
=
ψ(x)
is the wavelet function.
For example for the Haar wavelet we have
and
The equations become
and
W0(x)
=
φ(x)
is the Haar scaling
function and W1(x)
=
ψ(x)
is the Haar wavelet,
both supported in [0, 1]. Then we can obtain W2n by
adding two 1/2-scaled versions of Wn with
distinct supports [0,1/2] and [1/2,1] and obtain W2n+1 by
subtracting the same versions of Wn.
For n = 0 to 7, we have the W-functions shown in the figure Haar Wavelet Packets.
Haar Wavelet Packets
This can be obtained using the following command:
[wfun,xgrid] = wpfun('db1',7,5);
which returns in wfun
the approximate values
of Wn for n =
0 to 7, computed on a 1/25 grid of the
support xgrid
.
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].
The figure db2 Wavelet Packets presents the
system of W-functions for the original db2
wavelet.
db2 Wavelet Packets
Wavelet Packet Atoms
Starting from the functions and following the same line leading to orthogonal wavelets, we consider the three-indexed family of analyzing functions (the waveforms):
where n∊N and (j,k)∊Z2.
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 the wavelet packets is that for fixed values of j and k, Wj,n,k analyzes the fluctuations of the signal roughly around the position 2j· k, at the scale 2j 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 Wn 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, Wr(n)(x)
“oscillates”
approximately n times.
To analyze a signal (the chirp of Example 2 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.
These options are also available from command-line mode when
using the wpviewcf
function.
Organizing the Wavelet Packets
The set of functions Wj,n = (Wj,n,k(x),k∊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, ..., 2j–1.
Wavelet Packets Organized in a Tree; Scale j Defines Depth and Frequency n Defines Position in the Tree
The notation Wj,n, where j denotes scale parameter and n the frequency parameter, is consistent with the usual depth-position tree labeling.
We have , and .
It turns out that the library of wavelet packet bases contains the wavelet basis and also several other bases. Let us have a look at some of those bases. More precisely, let V0 denote the space (spanned by the family W0,0) in which the signal to be analyzed lies; then (Wd,1; d ≥ 1) is an orthogonal basis of V0.
For every strictly positive integer D, (WD,0, (Wd,1; 1 ≤ d ≤ D)) is an orthogonal basis of V0.
We also know that the family of functions {(Wj+1,2n), (Wj+1,2n+1)} is an orthogonal basis of the space spanned by Wj,n, which is split into two subspaces: Wj+1,2n spans the first subspace, and Wj+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 V0, 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 the 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 = 2L can be expanded in α different ways, where α is the number of binary subtrees of a complete binary tree of depth L. As a result, (see [Mal98] page 323).
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 verifying 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 [CoiW92]); many others are available
and can be easily integrated (type help
wentropy
).
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 such that E(0) = 0 and
The (nonnormalized) Shannon entropy
so
with the convention 0log(0) = 0.
The concentration in lp norm with 1 ℜ ≤ p
so
The logarithm of the “energy” entropy
so
with the convention log(0) = 0.
The threshold entropy
if and 0 elsewhere, so E4(s) = # {i such that } is the number of time instants when the signal is greater than a threshold ε.
These entropy functions are available using the wentropy
file.
Example 1: Compute Various Entropies
Generate a signal of energy equal to 1.
s = ones(1,16)*0.25;
Compute the Shannon entropy of s.
e1 = wentropy(s,'shannon') e1 = 2.7726
Compute the l1.5 entropy of s, equivalent to
norm(s,1.5)
1.5.e2 = wentropy(s,'norm',1.5) e2 = 2
Compute the “log energy” entropy of s.
e3 = wentropy(s,'log energy') e3 = -44.3614
Compute the threshold entropy of s using a threshold value of 0.24.
e4 = wentropy(s,'threshold', 0.24) e4 = 16
Example 2: Minimum-Entropy Decomposition
This simple example illustrates the use of entropy to determine whether a new splitting is of interest to obtain a minimum-entropy decomposition.
We start with a constant original signal. Two pieces of information are sufficient to define and to recover the signal (i.e., length and constant value).
w00 = ones(1,16)*0.25;
Compute entropy of original signal.
e00 = wentropy(w00,'shannon') e00 = 2.7726
Then split
w00
using the haar wavelet.[w10,w11] = dwt(w00,'db1');
Compute entropy of approximation at level 1.
e10 = wentropy(w10,'shannon') e10 = 2.0794
The detail of level 1,
w11
, is zero; the entropye11
is zero. Due to the additivity property the entropy of decomposition is given bye10+e11=2.0794
. This has to be compared to the initial entropye00=2.7726
. We havee10 + e11 < e00
, so the splitting is interesting.Now split
w10
(notw11
because the splitting of a null vector is without interest since the entropy is zero).[w20,w21] = dwt(w10,'db1');
We have
w20=0.5*ones(1,4)
andw21
is zero. The entropy of the approximation level 2 ise20 = wentropy(w20,'shannon') e20 = 1.3863
Again we have
e20 + 0 < e10
, so splitting makes the entropy decrease.Then
[w30,w31] = dwt(w20,'db1'); e30 = wentropy(w30,'shannon') e30 = 0.6931 [w40,w41] = dwt(w30,'db1') w40 = 1.0000 w41 = 0 e40 = wentropy(w40,'shannon') e40 = 0
In the last splitting operation we find that only one piece of information is needed to reconstruct the original signal. The wavelet basis at level 4 is a best basis according to Shannon entropy (with null optimal entropy since
e40+e41+e31+e21+e11 = 0
).Perform wavelet packets decomposition of the signal s defined in example 1.
t = wpdec(s,4,'haar','shannon');
The wavelet packet tree in Entropy Values shows the nodes labeled with original entropy numbers.
Entropy Values
Compute the best tree.
bt = besttree(t);
The best tree is shown in the following figure. In this case, the best tree corresponds to the wavelet tree. The nodes are labeled with optimal entropy.
Optimal Entropy Values
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.
Decomposition Tree | Subtree Such That the Set of Leaves Is a Basis |
---|---|
Wavelet packets decomposition tree | Complete binary tree: WPT of depth D |
Wavelet packets optimal decomposition tree | Binary subtree of WPT |
Wavelet packets best-level tree | Complete binary subtree of WPT |
Wavelet decomposition tree | Left unilateral binary subtree of WPT of depth D |
Wavelet best-basis tree | Left unilateral binary subtree of WPT |
We deduce the following definitions of optimal decompositions, with respect to an entropy criterion E.
Decompositions | Optimal Decomposition | |
---|---|---|
Wavelet packet decompositions | Search among 2D trees | Search among D trees |
Wavelet decompositions | Search among D trees | Search among D trees |
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).
Entropy Condition | Action on Tree and on Entropy Labeling |
---|---|
| If (node≠root), merge and set Eopt(node) = E(node) |
Split and set |
with the natural initial condition on the reference tree, Eopt(t)
= E(t)
for
each terminal node t.
Reconstructing a Signal Approximation from a Node
You can use the function wprcoef
to
reconstruct an approximation to your signal from any node in the wavelet
packet tree. This is true irrespective of whether you are working
with a full wavelet packet tree, or a subtree determined by an optimality
criterion. Use wpcoef
if you
want to extract the wavelet packet coefficients from a node without
reconstructing an approximation to the signal.
Load the noisy Doppler signal.
load noisdopp
Compute the wavelet packet decomposition down to level 5 using
the sym4
wavelet. Use the periodization mode.
dwtmode('per'); T = wpdec(noisdopp,5,'sym4'); plot(T)
Plot the binary wavelet packet tree and click on the (4,1) doublet (node 16).
Extract the wavelet packet coefficients from node 16.
wpc = wpcoef(T,16); % wpc is length 64
Obtain an approximation to the signal from node 16.
rwpc = wprcoef(T,16); % rwpc is length 1024 plot(noisdopp,'k'); hold on; plot(rwpc,'b','linewidth',2); axis tight;
Determine the optimum binary wavelet packet tree.
Topt = besttree(T); % plot the best tree plot(Topt)
Reconstruct an approximation to the signal from the (3,0) doublet (node 7).
rsig = wprcoef(Topt,7); % rsig is length 1024 plot(noisdopp,'k'); hold on; plot(rsig,'b','linewidth',2); axis tight;
If you know which doublet in the binary wavelet packet tree
you want to extract, you can determine the node corresponding to that
doublet with depo2ind
.
For example, to determine the node corresponding to the doublet (3,0), enter:
Node = depo2ind(2,[3 0]);
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.
Quaternary Tree of 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
with
an entropy function.