A PYRAMID ALGORITHM FOR THE HAAR DISCRETE WAVELET

One area of application of the discrete wavelet transform (DWT) has been the detection and classification of physiological signals such as electroencephalography (EEG) signals. Anomalies in EEGs yield very low frequency signals which are ideal for analysis using the DWT. Anomalies in mechanical systems yield high frequency signals. The structure of the DWT makes it an un-ideal tool for the analysis of such signals. Such signals are, however, ideal for analysis using the wavelet packet transform (WPT) in which Mallat’s pyramid algorithm is applied to the multiresolution analysis (MRA) of both the approximation and detail subspaces of a signal. As a contribution to the computer-aided signal processing of non-stationary signals, this paper develops a pyramid algorithm for the discrete wavelet packet transform (DWPT) from two-scale relations for wavelet packets. The algorithm is used in the derivation of the fast Haar discrete wavelet packet transform (FHDWPT) and its inverse. It is found out that the FHDWPT is computationally as efficient as the fast Fourier transform (FFT).


INTRODUCTION
Wavelet-based digital signal processing techniques are superior to other techniques in the analysis of non-stationary signals (Mallat, 1989a;Burthiel, 2011;Mallat, 1989b). The discrete wavelet transform (DWT) is one such wavelet-based signal processing technique.
Mallat developed a pyramid algorithm which made significant improvement to the computational efficiency of the DWT (Mallat, 1989a).
The DWT has been applied to the analysis of physiological signals such as electroencephalogram (EEG) signals focusing on the detection and classification of anomalies occurring in such signals [Gell-Man, 2004;Rosso, 2001]. Anomalies in EEGs produce low frequency signals which given the structure of the DWT, are ideal for analysis using the DWT. The structure of the DWT makes it un-deal for the analysis of anomalies which produce high frequency signals. Such anomalies occur in, for example, mechanical and electrical systems [Shi, 2007]. Such anomalies are best analyzed using the wavelet packet transform (WPT) in which the wavelet decomposition yields fine gradation of both low and high frequencies in contrast to the DWT which only produces fine gradations of low frequencies [Shi, 2007]. Efficient algorithms are, therefore, needed for the computer-aided analysis of signals using the WPT.
This work develops a pyramid algorithm for the WPT starting from the two-scale relation for wavelet packets. The algorithm is applied to the development of a fast Haar discrete wavelet packet transform (FHDWPT) and its inverse. The FHDWPT is found to be computationally as efficient as the fast Fourier transform (FFT).
Multiresolution analysis (MRA)is central to the DWT (Mallat, 1989c;Coifman, et al., 1992). In MRA, a subspace V j of a multiresolution approximation is decomposed into a lower resolution approximation subspace V j-1 plus a detail subspace W j-1 . The decomposition is done by dividing the orthogonal basis 2 /$ Ф{2 ˮ − ˫{ ˫ ∈ I of Ͳ ͺ into two new orthogonal bases (1) and Ӝ2 % {2 # ˮ − ˫{ӝ ˫ ∈ I of W j-1 .
(2) which are in line with the axioms ofMRA that {Ф(t-k)} k ∈ Z is an orthonormal basis of Vo and {φ(t-k)} k ∈ Z is an orthonormal basin of W o (Mallat, 1989a). The DWT decomposition is represented by where W j-1 is the orthogonal complement of V j-1 in V j and denotes the sum of mutually orthogonal subspaces (Bulthiel, 2011).
The DWT is based on wavelet decomposition (WD). The WD splits the original signal (S), assumed to be contained in subspace V j , into two components -the approximation contained in subspace V j-1 and details contained in subspace W j-1 (Mallat, 1989a). The approximation is made up of the low frequency information in the signal and the detail contains the high frequency information.
In WD only the subspace containing the low frequency information is iteratively decomposed into lower time resolution subspaces as shown in Fig.1.

Fig. 1: The Structure of Wavelet Decomposition
Coifman, Meyer and Wickerhausen introduced wavelet packets by generalizing WD into the wavelet packet decomposition (WPD) (Coifman, et al., 1992). In the WPD both the approximation (low frequency part) and the details (high frequency part) are iteratively decomposed into lower time resolution subspaces. The WPD yields a finer gradation of frequency components in a full wavelet tree decomposition as shown in Fig. 2. The decomposition in Fig. 2 is assumed to start at some scale level N. Any node of the decomposition tree is labeled (j,n) where j is the scale level and n is the number of nodes that are on its left at the same scale level. Thus N-j≥0 is the level of the node in the tree.
Note that at each scale level of Fig.2, the first two boxes from the left correspond to the DWT decomposition (WD).
The decomposition of the wavelet subspaces in the high frequency section of Fig.2 parallels that of the DWT decomposition. In this case, to each node (j,n) of Fig.2 we associate a space W j,n which is the direct sum of two orthogonal subspaces, i.e. (Mallat, 1989a); ..
(4) Eqn (4) will form the basis of the work reported here.

Fig.2: Structure of Wavelet Packet Decomposition
The purpose of this work is to present a pyramid algorithm for the Haar Discrete Wavelet Packet Transform (HDWPT) which uses two types of matricesan orthogonal matrix whose elements are made up of the Haar scaling and wavelet coefficients and a series of permutation matrices. As a departure point we reformulate some of the ideas of Wickerhauser (Wickerhauser, 1992).

Two-Scale Relations for Wavelet Packets
The two-scale relations for scaling functions and for the mother wavelet function are given by Where {h k } k ∈Z and {g k } k ∈Z are impulse responses of low-pass and high-pass filters, respectively, known as quadrature mirror filters (QMF) related through For j=1, eqn (3) becomes From eqns (1) and (2) we note that {Ф (t-k)} k ∈Z and {φ (t-k)} k ∈Zform an orthonormal basis for V o and W o , respectively. Hence eqn (6) read together with eqn (3) shows that the V subspace can be decomposed into a direct sum of two orthogonal subspaces defined by their basis functions which are given by eqns (5a) and (5b).
This splitting algorithm for the V space leading to the DWT can also be applied to the W space representing high frequencies.
But, in order to make the development easier to understand, we first redefine the variables in eqns (5a) and (5b).
Equations (5a) and (5b) now become Unlike WD which involved the iterative splitting of the V space, WPD involves the iterative splitting of the W space. The iterative splitting of the W space is achieved by defining the following sequence of functions from W 0 (t) and W 1 (t): where now {W 2 (t-k)} k ∈Z and {W 3 (t-k)} k ∈Zare orthonormal basis functions for the two subspaces whose direct sum yields W 1,1 .
Eqns (1) and (2) are combined to form a basis for V j . Now the sequence of functions {W n (t)} n ∈Zwith their dilations and translations can be used to form an orthonormal basis for a function space.

The functions defined by
are said to form the (j,n) wavelet packet which is an orthogonal basis for the space W j,n. Eqns (11a) and (11b) thus form the two-scale relations for wavelet packets.

The Algorithm
A given signal f(t) can be decomposed using the basis functions given in eqn (12). The decomposition coefficients are given by Dilating eqns (11a) and (11b) by 2 j-1 and translating them by m and applying eqn (13) to the resulting equations may be shown to yield.
The coefficients ,$ and ,$ # can thus be successively obtained from some highest level coefficients , for some large N.

The Fast Haar Wavelet Packet Transform
Assuming that W o (t) = Ф(t) is the Haar scaling function and that W 1 (t) = φ(t) is the Haar mother wavelet, eqn 14 may be shown to reduce to Eqns (15) show that the (j-1) tn level coefficients are obtained from the j tn level coefficients through multiplication by the following orthogonal matrix We assume that the signal f(t) to be decomposed and all orthogonal functions of eqn (12) which are used in signal decomposition have support over the unit interval [0,1). For the orthogonal functions of eqn(12) this implies that the translation parameter, k, is bounded by 0≤ k< 2 j -1 (17) The signal f(t) is also assumed to be discretized to yield f k = f(t=t k ), k=0,1 ., M-1 M is set to be a power of 2, i.e. M=2 N say with the time values, t k , being restricted to the unit interval by setting t k = , k = 0,1, , M-1 (19) The nesting and density axioms of MRA ( Mallat, 1989a;Mallat, 1989c) are invoked to conclude that for large values, N, of the scaling parameter the approximation can be made to be as accurate as possible in L 2 (R) norm by proper choice of N.
For N large and with relation (17) satisfied, we can assume that the scaling functionW o (t) is so narrow that it covers only the first data point t o = 0. Thus by choosing a large N and translating the scaling function we can invoke the wavelet crime (Strang, et al., 1996) by setting the scaling function coefficients equal to the data point , =f k , k = 0, 1, ., M-1 We next express the relations ineqn (15) in matrix form. To facilitate that we need to define the orthogonal sum of matrices.
The orthogonal sum of two matrices A and B, denoted by A⊕B, is Now let H M represent the M x M matrix given by H M = ϴ⊕ϴ⊕ ϴ (23) where the operation is applied M/2 times and ϴ is the matrix defined in eqn (16). Now using eqns (15) and eqn (21)  WhereT represents the transpose operator for vectors or matrices.

Discrete Haar Series
The coefficients in eqn (27) may be used to write a Haar series for f(t). The problem with the WPT is the choice of the best basis to use in representing f(t) for a given application (Mallat, 1989a;Wickerhauser, 1992). The leaves of any admissible tree from Fig. 2 form a basis. An admissible tree ( Mallat, 1989a) is any tree where each node in the tree has 0 or 2 children nodes. This yields a large library of admissible bases which leads to an optimization problem in the choice of the best basis (Mallat, 1989a;Strang, et al., 1996).
The leaves of any admissible tree would form an adequate basis for the decomposition of functions f(t). Not all leaves of an admissible tree fall on the same level of a wavelet packet tree similar to the levels in Fig.2 (Mallat, 1989a;Daly, 2007). This is what makes it difficult to write a simple algebraic relationship to represent signal decomposition as is possible for the DWT. One of the few exceptions where we can write an expression to represent signal decomposition is the case where all leaves of an admissible tree fall on the same level.
Let us assume that all leaves of an admissible tree fall on level m (equivalent to scale parameter j=N-m). Then f(t) can be decomposed as The factor 1/ H is required to preserve the orthonormality of the wavelet packets.
Note that the normalized frequencies occupied by the wavelet packets of eqn (31)

EXAMPLE
To fix ideas we consider a simple example whose data is shown below We make two immediate conclusions from Table 2. First, if the decomposition is stopped at levels 2 or 1, the WPT decomposition compresses the signals energy more than the DWT as the latter has two coefficients whose values are 0.
Second, we can use Table 2 to find an admissible tree whose leaves are not at the same level and use this tree in signal decomposition. Using the nomenclature of Fig.2 we choose an admissible tree whose leaves are at nodes (2,0), (1,2) and (1,3).

CONCLUSION
A pyramid algorithm for the discrete wavelet packet transform has been derived by the dilation and translation of wavelet packets. The algorithm has been used to yield an 0(M log M) Fast Haar Wavelet Packet Transform and its inverse. The algorithm has been applied to a simple illustrative example.