A Wavelet Packet-Based Algorithm for the Extraction of Neural Rhythms

Abstract

Neural rhythms are associated with different brain functions and pathological conditions. These rhythms are often clinically relevant for purposes of diagnosis or treatment, though their complex, time-varying features make them difficult to isolate. The wavelet packet transform has proven itself to be versatile and effective with respect to resolving signal features in both time and frequency. We propose a signal analysis technique, called neural rhythm extraction (NRE) that incorporates wavelet packet analysis along with a threshold-based scheme for separating rhythmic neural features from non-rhythmic ones. We applied NRE to rat in vitro intracellular recordings and human scalp electroencephalogram (EEG) signals, and were able to isolate and classify individual neural rhythms in signals containing large amplitude spikes and other artifacts. NRE is capable of discriminating signal features sharing similar time or frequency localization, as well as extracting low-amplitude, low-power rhythms otherwise masked by spectrally dominant signal components. The algorithm allows for independent retention and reconstruction of rhythmic features, which may serve to enhance other analysis techniques such as independent component analysis (ICA), and aid in application-specific tasks such as detection, classification or tracking.

Appendix

The wavelet transform10 is the correlation of a signal, x(t), with a mother wavelet function, ψ(t), that is shifted and scaled, yielding wavelet coefficients c _σ,τ :

$$ c_{{\sigma ,\tau }} = \frac{1} {{{\sqrt \sigma }}}{\int\limits_{ - \infty }^\infty {x(t)\psi {\left( {\frac{{t - \tau }} {\sigma }} \right)}dt} }. $$

(A.1)

If discrete scaling and translation are applied such that σ = 2^j and τ = 2^j k, where k and j are integers, then a basis of orthonormal wavelet functions for L ²(R) can be obtained30:

$$ \psi _{{j,k}} = {\sqrt {2^{{ - j}} } }\psi (2^{{ - j}} t - k). $$

(A.2)

Restating the coefficients in terms of j and k, the signal can be represented analytically as a summation of wavelet coefficients and basis functions:

$$ x(t) = {\sum\limits_{j = - \infty }^\infty {{\sum\limits_{k = - \infty }^\infty {c_{{j,k}} \psi _{{j,k}} } }} }. $$

(A.3)

Practical implementation of the discrete wavelet transform (DWT) is accomplished using a quadrature mirror filter (QMF) cascade.10,30 This requires the basis functions have compact support and the signal have finite energy and bandwidth. Let ϕ(t) be a lowpass scaling function associated with the wavelet function ψ(t), such that the following expressions hold:

$$ \phi (t) = {\sqrt 2 }{\sum\limits_k {h_{k} \phi (2t - k)} } $$

(A.4)

$$ \psi (t) = {\sqrt 2 }{\sum\limits_k {g_{k} \phi (2t - k)} } $$

(A.5)

where h _k and g _k are lowpass and highpass QMF coefficients, respectively. The filter coefficients can be calculated directly by using Eqs. (A.4) and (A.5) and the scale and shift orthogonality of ϕ(t) and ψ(t) to form the inner products

$$ h_{k} ={\left\langle {\phi (t),\phi_{- 1,k}} \right\rangle } = {\sqrt 2 }{\int\limits_{ - \infty }^\infty {\phi (t) \cdot \phi (2t - k)dt} } $$

(A.6)

$$ g_{k} ={\left\langle {\psi (t),\phi_{-1, k}} \right\rangle } = {\sqrt 2 }{\int\limits_{ - \infty }^\infty {\psi (t) \cdot \phi (2t - k)dt} } $$

(A.7)

for all finite k ∈ Z. The QMF pair has the special relationship

$$ g_{k} = ( - 1)^{k} h_{{1 - k}} $$

(A.8)

in addition to satisfying the following orthogonality and normalization properties:

$$ {\sum\limits_k {h_{{k - 2m}} h_{{k - 2{m}^{\prime} }} } } = \delta _{{m{m}^{\prime} }} $$

(A.9)

$$ {\sum\limits_k {g_{{k - 2m}} g_{{k - 2{m}^{\prime} }} } } = \delta _{{m{m}^{\prime} }} $$

(A.10)

$$ {\sum\limits_k {g_{{k - 2m}} h_{{k - 2{m}^{\prime} }} } } = 0 $$

(A.11)

$$ {\sum\limits_k {h_{k} } } = {\sqrt 2 } $$

(A.12)

$$ {\sum\limits_k {g_{k} } } = 0 $$

(A.13)

for all $m,{m}^{\prime} \in Z,$ where $\delta _{{m{m}^{\prime} }}$ is the Kronecker delta.10 Alternatively, when the filter coefficients are known, the wavelet and scaling functions can be produced from the coefficients.

In general, for an N level multiresolution analysis using the DWT, the signal is represented by the expression:

$$ \begin{aligned}{} x(t) &= A_{N} (t) + D_{N} (t) + \cdots + D_{1} (t) \\ & = {\sum\limits_k {c^{{A,n}}_{k} \phi _{{N,k}} + {\sum\limits_n {{\sum\limits_k {c^{{D,n}}_{k} \psi _{{n,k}} } }} }} } \\ \end{aligned} $$

(A.14)

where A(t) and D(t) are approximation and detail signals, and $c^{{A,n}}_{k}$ and $c^{{D,n}}_{k}$ are the nth level approximation and detail coefficients, respectively. The (n + 1) level coefficients are computed by discrete convolution of the previous nth level approximation coefficients with the time-reversed QMF impulse responses, h and g, along with downsampling by 230:

$$ c^{{A,n + 1}}_{k} = [h_{{ - k}} * c^{{A,n}}_{k} ]( \downarrow 2) = {\sum\limits_m {h_{{m - 2k}} c^{{A,n}}_{m} } } $$

(A.15)

$$ c^{{D,n + 1}}_{k} = [g_{{ - k}} * c^{{A,n}}_{k} ]( \downarrow 2) = {\sum\limits_m {g_{{m - 2k}} c^{{A,n}}_{m} } }. $$

(A.16)

The downsampling operation removes redundant information in the coefficients, since the doubling of scale at each level halves the temporal resolution. By downsampling, the total number of coefficients remains roughly the same throughout the analysis. Time reversal of the filters for decomposition vs. reconstruction prevents aliasing. For signal reconstruction using the QMF cascade, the coefficients are upsampled by 2 before being passed through their respective filters and added together.

To obtain a finer resolution analysis than the DWT, and for better discrimination of higher frequency features, the discrete wavelet packet transform (DWPT) is used.8,18 Implementation of the DWPT is similar to the DWT, except that the scaling and wavelet functions of the DWT are extended into a set of orthonormal wavelet packet functions, W _m (t), which are defined by the recursion relations18:

$$ W_{{2r}} (t) = {\sqrt 2 }{\sum\limits_k {h_{k} W_{r} (2t - k)} } $$

(A.17)

$$ W_{{2r + 1}} (t) = {\sqrt 2 }{\sum\limits_k {g_{k} W_{r} (2t - k)} } $$

(A.18)

for integers r ≥ 0, where W ₀(t) = ϕ(t) and W ₁(t) = ψ(t). Multiresolution analysis gives rise to a binary tree structure, with coefficient nodes generated from a dictionary of wavelet packet atoms having discrete scaling, translation and sub-band localization:

$$ W_{{n,r,k}} = {\sqrt {2^{{ - n}} } }W_{r} (2^{{ - n}} t - k) $$

(A.19)

where the scale exponent n = 0, 1, 2, …, corresponds directly to the level in the tree, and index r = 0, 1, 2, …, 2ⁿ − 1 corresponds to the horizontal node position. If tree nodes are re-ordered in terms of increasing frequency, then r becomes a frequency index relating to the sub-band analyzed by atoms W _n,r,k at scale 2ⁿ and time 2ⁿ k. For this paper, wavelet packet atoms were derived from Meyer scaling and wavelet functions, which are orthogonal, symmetric and compactly supported in the frequency domain.24 Compact frequency support provides good frequency localization for purposes of separating rhythmic components, and symmetric filters ensure the phase behavior is linear. For purposes of computing the DWPT, the impulse responses of the Meyer functions were approximated to obtain compact support in time.

Signal decomposition via the DWPT is accomplished by sequentially splitting coefficient nodes into daughter node pairs, while synthesis involves the reverse operation. The leaves (terminal nodes) of the connected tree constitute an orthonormal basis from which the original signal can be represented by the following summation:

$$ x(t) = {\sum\limits_{{n}^{\prime} } {{\sum\limits_{{r}^{\prime} } {x_{{{n}^{\prime},{r}^{\prime} }} (t)} }} } $$

(A.20)

where the reconstructed nodal time series

$$ x_{{{n}^{\prime},{r}^{\prime} }} (t) = {\sum\limits_k {c_{{{n}^{\prime},{r}^{\prime},k}} W_{{{n}^{\prime},{r}^{\prime},k}} } } $$

(A.21)

are defined for (n′,r′) belonging to the set denoting the leaves of the tree (with coefficients $c_{{{n}^{\prime},{r}^{\prime},k}}$).