EURASIP Journal on Applied Signal Processing 2002:11, 1289–1295 c ○ 2002 Hindawi Publishing Corporation Some Results on the Wavelet Packet Decomposition of Nonstationary Processes

Wavelet/wavelet packet decomposition has become a very useful tool in describing nonstationary processes. Important examples of nonstationary processes encountered in practice are cyclostationary processes or almost-cyclostationary processes. In this paper, we study the statistical properties of the wavelet packet decomposition of a large class of nonstationary processes, including in particular cyclostationary and almost-cyclostationary processes. We first investigate in a general framework, the existence and some properties of the cumulants of wavelet packet coefficients. We then study more precisely the almost-cyclostationary case, and determine the asymptotic distributions of wavelet packet coefficients. Finally, we particularize some of our results in the cyclostationary case before providing some illustrative simulations.


INTRODUCTION
Nonstationary processes have recently received an increasing attention in the signal processing community, due to their wide applicability in modeling natural phenomena generated by physical systems with time-varying parameters. Let x(t), t ∈ R, be such a zero-mean random process. Its translated cumulants (if existing) cum n+1 x,u (t) = cum(x(u), x(t 1 + u), . . . , x(t n + u)) depend in general on the lag u. When these functions in u are uniformly almost periodic (UAP) [1], we say that the process is nth-order almost-cyclostationary. This means that the cumulants have generalized Fourier series If the dependence on u of the translated cumulants is periodic, the process is said to be nth-order cyclostationary. Cyclostationary and almost-cyclostationary processes are very useful for modeling many real signals which appear in communication, telemetry, radar, sonar, and economics [4]. In practice, the nth-order translated cumulants may also be assumed to have (n − 1)th-order finite energy, that is, to belong to L 2 (R n−1 ). Moreover, they often correspond to regular functions. This structure within the data suggests the use of a multiscale analysis, for example, a wavelet analysis, for the extraction of these higher-order statistical informations, combined with a Fourier analysis for the determination of the cyclic characteristics. In [5], a wavelet analysis for the stationary case was realized. The field of wavelet coefficients (when taken at the same resolution level) was shown to be asymptotically (i.e., for coarse resolutions) a white Gaussian noise. Similar results had already been obtained for the empirical wavelet coefficients [6], but also restricted to stationary processes. However, as already mentioned, nonstationary random processes are often encountered in many practical situations. In this paper, we are interested in higher-order statistics of the wavelet coefficients of nonstationary processes, and especially cyclostationary/almost-cyclostationary processes. Other relevant works concerning nonstationary processes may be found in [7,8,9,10,11,12].
In Section 2 of this paper, we introduce the necessary assumptions and notations. In Section 3, we investigate the existence and some properties of wavelet coefficient cumulants for a wide class of nonstationary processes. Section 4 is devoted to the almost-cyclostationary case and Section 5 provides some specific results in the cyclostationary case. In Section 6, we give some simulations in order to illustrate our results. Finally, some conclusions are drawn in Section 7.

HYPOTHESES AND NOTATIONS
Throughout the paper, the symbols N, N * , Z, and R denote, respectively, the sets of nonnegative integers, positive integers, integers, and reals. Standard notations are also used to denote the Banach space of n-dimensional (n ∈ N * ) functions and sequences L p (R n ) and p (Z n ), respectively, 1 ≤ p ≤ ∞. We consider the general case of an M-band real orthonormal wavelet packet decomposition [13,14] related to a paraunitary perfect reconstruction filter bank with impulse responses (h 0 (k)) k∈Z , . . . , (h M−1 (k)) k∈Z in 2 (Z). The associated wavelet packet functions will be denoted in the sequel by W m (t), m ∈ N. We also recall that the wavelet packets satisfy the following two-scale equation: for all m ∈ N, for all We also have Note that W 0 (t) corresponds to the scaling function associated with a multiresolution analysis (MRA) of L 2 (R) [15] while W m (t), m ∈ {1, . . . , M − 1}, represent the M-band mother wavelets. The wavelet packet coefficients of the process x(t) at the resolution level j ∈ Z and in mth frequency bin, will be denoted by Their (n + 1)th-order cumulants (if existing) at multiresolution level j = ( j 0 , j 1 , . . . , j n ) are, for k = (k 0 , . . . , k n ) ∈ Z n+1 , cum n+1 cj,m k = cum c j0,m0 k 0 , c j1,m1 k 1 , . . . , c jn,mn k n .
Finally, we make the following hypothesis.
x,u (·) has finite energy, that is, belongs to L 2 (R n ), and this energy is bounded with respect to u by a constant C n+1 .

EXISTENCE AND SOME PROPERTIES OF THE WAVELET COEFFICIENT CUMULANTS
In this part we first show that, under Hypothesis 1, the wavelet coefficient cumulants are properly defined, and using (2) and (3), we obtain a recursive algorithm for calculating these cumulants. Proposition 1. Under Hypothesis 1, the wavelet coefficient cumulants are well defined.

Proposition 2. Under Hypothesis 1, the time-shifted cumulant cum n+1
cq,j 0 ,m 0 ; j,m (·) defines a sequence that belongs to 2 (Z n ). By using the Plancherel formula, the cumulant is equal to Note that this integral (on R n+1 ) is (as the first one) absolutely convergent. So, from Fubini theorem, it is equal to and, consequently [16], Using Cauchy inequality and the orthonormality of the wavelet packets written in frequency domain, as well as the Parseval equality, we obtain thatΓ n+1 cq,j 0 ,m 0 ; j,m (ν) belongs to L 2 ([−1/2, 1/2] n ) and its norm is upper bounded by C n+1 M j0/2 W m0 L 1 (R) . Consequently the Fourier transform of cum n+1 cq,j 0 ,m 0 ; j,m (·) isΓ n+1 cq,j 0 ,m 0 ; j,m (·) and cum n+1 cq,j 0 ,m 0 ; j,m (·) belongs to 2 (Z n ). This result will be used in Section 4 to determine the cyclospectra of the wavelet packet coefficient field in the almost-cyclostationary case.
From the two-scale equations (2) and (3), a scale recursive algorithm can be proposed to compute the higher-order statistics of the wavelet packet coefficients.
The cumulants of the wavelet packet coefficients at multiresolution level j + 1 = ( j 0 + 1, . . . , j n + 1) are obtained from those of the coefficients at level ( j r , j r + 1) = ( j 0 , . . . , j r−1 , j r + 1, . . . , j n + 1) through the multidimensional filter bank shown in Figure 1 involving filters with impulse responses and r-dimensional decimators by a factor M.
Proof. We have So the result is obtained from (2), applied on the r first wavelet packets, if assuming a finite length for the original filters h i , i = 0, . . . , M − 1. Otherwise, from the fact that (2) holds in L 1 (R) sense, as well as L 2 (R) sense, and using some technical arguments, we obtain the same result.
This result could be useful in providing alternatives to statistical methods based on prior modelling of the wavelet coefficients. This formulation provides both analysis and reconstruction equations. The analysis equation allows to express the cumulants at a multiresolution level from those at the lower levels, and the reconstruction equation allows to do the reverse operation. If we suppose that we have (by means of asymptotic estimation techniques or other methods) the cumulant field at a multiresolution level ( j * 0 , j * 1 , . . . , j * n ), these equations allow to calculate recursively wavelet coefficient cumulant field corresponding to different scales ( j 0 , j 1 , . . . , j n ).
Similar equations have been obtained in [5] but they are restricted to the stationary case and they can only be used to calculate recursively intra-scale higher-order statistics.
We can also derive an energy-conservation property through resolution levels for shifted cumulants.
In the following, we focus on the almost-cyclostationary case. We give frequency characteristics, and we obtain important asymptotic results for the wavelet coefficient cumulants.

THE ALMOST-CYCLOSTATIONARY CASE
We assume that the polyspectrum of the signal Γ n+1 x,u (ν) has a Fourier expansion with cyclic frequencies {η s } s∈Z [4]. This is the case, for example, if the cumulants have Fourier expansions in L 2 (R n ) sense (with respect to t) (In order to simplify the notations, we do not recall the dependence on n in A s .) So the polyspectrum of the signal is expanded, in L 2 (R n ) sense (with respect to ν), as whereÂ s (ν) is the n-dimensional Fourier transform of A s (t). This expansion allows to obtain a simple expression for the shifted polyspectrum of the wavelet packet coefficients. We denote this polyspectrum by Γ n+1 cq,j 0 ,m 0 ; j,m (ν), which is the n-dimensional Fourier transform of the shifted cumulant cum n+1 cq,j 0 ,m 0 ; j,m (k) with respect to k. This polyspectrum is well defined thanks to Proposition 2. Proposition 4. Under Hypothesis 1, and assuming that for almost all ν, s∈Z |Â s (ν)| < ∞, the translated polyspectrum is given by Proof. In the proof of Proposition 2, we have seen that the polyspectrum of the decomposition coefficients is So using Hypothesis 1(2) and the fact that s∈Z |Â s (ν)| < ∞, the desired expression of the polyspectrum is obtained.
From this expression of the polyspectrum, we note that, if x is (n + 1)th-order almost-cyclostationary with spectrum S x = {η s } s∈Z , then cum n+1 c(j,...,j),m (q, k + q) is almost-periodic (in q) [1], and its spectrum is {M j · η s } s∈Z . The coefficients of the expansion with respect to s ∈ Z are the high-order cyclospectra of the wavelet coefficients. For coefficients defined at different resolutions, this property of cyclostationarity no longer holds. Now, consider the asymptotic behaviour of the wavelet coefficients. First, we obtain the following bounds on the shifted polyspectrum, and then we deduce bounds on the translated cumulants of the wavelet packet coefficients.
In the following, we will make use of the 2nd-order wavelet packets [5] as defined below Obviously, we have We establish the following lemma by iterating the analysis equation (2).
Note that these results may be useful in computing the 2nd-order wavelet packets I 2 ( j0, j1),(m0,m1) (k 0 , k 1 ) which are closely related to the asymptotic correlations of the wavelet coefficients as shown below. Proposition 6. With the same notations as in the previous propositions, let C j,m (k) = (c j0,m0 (k 0 ), . . . , c jp,mp (k p )) be an array of wavelet packet coefficients of a zero-mean almostcyclostationary process with fixed values for the differences j i − j 0 , i = 1, . . . , p. Then C j,m (k) converges in law to a Gaussian vector as j = min i=0,...,p ( j i ) → ∞, provided that Cj,m(k) [l, l ] = I 2 ( jl, j l ),(ml,m l ) k l , k l .
We denote by R j,m,k (u) the remainder of order 3 of the Taylor-series expansion of ψ Cj,m(k) (u). From Proposition 4 we deduce that incum c j0,m0 k 0 n0 , . . . , c jp,mp K p np n 0 ! · · · n p ! u n0 Consequently, it is sufficient to find the asymptotic statistics of order 2. From the expression of the polyspectrum in Proposition 4, we easily obtain the correlation matrix by further using that, as W m belongs to L 1 (R), Note 2. The convergence could be more directly proved by invoking the theorem of moments, but the interest of our proof is to provide the convergence rate. Some examples of signals satisfying the assumptions made in this section are: • AM signals: where f 0 and f 1 are incommensurable, a(t) and b(t) are nth-order independent stationary signals such that cum n+1 a (t) and cum n+1 b (t) belong to L 2 (R n ) and their polyspectra are bounded by (n + 1)!K n+1 . These conditions are generally satisfied when a(t) and b(t) are linear processes obtained by filtering white noises with bounded-frequency response filters. • PM signals: x(t) = cos(2π f 0 t +φ(t)), with φ(t) almostcyclostationary. • Real AM signals: is almost-cyclostationary and a(t) is stationary and independent of φ(t).

SPECIFIC RESULTS IN THE CYCLOSTATIONARY CASE
In this section, we show how the general results established in the almost-cyclostationary case simplify in the cyclostationary case.
Under Hypothesis 1, we can easily obtain the Fourier expansion of the cumulants of a T-cyclostationary process x(t), in L 2 ([0, T] × R n ) sense, cum n+1 x,u t = s∈Z A s t e 2π(su/T) .
So the polyspectrum of the signal is expanded, in L 2 ([0, T] × R n ) sense, as whereÂ s (ν) is the n-dimensional Fourier transform of A s (t). This expansion is obtained, in the cyclostationary case, only by using Hypothesis 1. Note that the spectrum of the T-cyclostationary process is (1/T)Z.
The spectrum of the decomposition coefficients at the same resolution level j is (M j /T)Z. This implies that the decomposition coefficients field can be stationary only for Note that the condition |η s − η s | ≥ γ > 0 required in Proposition 5 is automatically satisfied by cyclostationary processes.
Finally the asymptotic correlations of the decomposition coefficients are the same as those obtained in Proposition 6 when replacing (lim D→+∞ (1/D) An important example of cyclostationary processes is given by PAM signals: x(t) = +∞ m=−∞ a m p(t − mT), where the support of the pulse function p(t) is included in the interval (−T/2, T/2), and (a m ) m is an IID sequence. Then x(t) is cyclostationary with cyclic period T.

SIMULATIONS
In this section, we present some simulation results illustrating some of the theoretical claims made in this paper. More precisely, we verify the Gaussianity of the decomposition coefficients at coarse resolution levels.
We generated a random time series X(t) corresponding to an amplitude-modulation of a χ 2 (independent stationary) signal by an almost-periodic function where ω 1 and ω 2 are incommensurable (ω 1 = 0.01 rd/s, ω 2 = πω 1 ). We generated 10000 samples of size D = 2 14 according to (36). For each realisation, we realized a wavelet decomposition on seven resolution levels. To do this, we used the Symmlet four basis. We focused on the resolutions levels j = 7 (the coarsest one), j = 5, j = 3, and j = 1.
In Figures 2a, 2b, 2c, and 2d we show the histogram of one particular coefficient at each level. We also show its approximations by generalised Gaussian distributions, with the We see that at coarse resolution levels the coefficients are well approximated by Gaussian distributions, β is close to 2. Such an approximation is no longer valid in higher levels, but generalized Gaussian distributions still provide good approximations for the wavelet decomposition.

CONCLUSIONS
In this paper, we have established fundamental results on the wavelet packet decomposition of almost cyclostationary processes. One of the main results is that the wavelet packet coefficients at a given resolution level also are almost cyclostationary. Due to the importance of this particular class of nonstationary processes in modeling communication signals as well as many other real world signals, our results should be useful in many applications.
In our future work, we plan to use the derived properties in denoising problems involving cyclostationary or almost cyclostationary noise.