Home | About | Journals | Submit | Contact Us | Français |

**|**HHS Author Manuscripts**|**PMC2993774

Formats

Article sections

- Abstract
- 1. Introduction
- 2. Theory
- 3. Illustrative examples
- 4. Discussion
- 5. Experimental methods
- References

Authors

Related links

J Neurosci Methods. Author manuscript; available in PMC 2011 December 15.

Published in final edited form as:

Published online 2010 October 7. doi: 10.1016/j.jneumeth.2010.09.019

PMCID: PMC2993774

NIHMSID: NIHMS243596

David Hsu, MD, PhD,^{1} Murielle Hsu, PhD,^{1} Heidi L. Grabenstatter, PhD,^{2} Gregory A. Worrell, MD, PhD,^{3} and Thomas P. Sutula, MD, PhD^{1}

Corresponding author: David Hsu, Department of Neurology, Rm 7178, 1685 Highland Av, Madison WI 53705-2281 USA, Email: hsu/at/neurology.wisc.edu, Phone 608-263-8851, Fax 608-263-0412

The publisher's final edited version of this article is available at J Neurosci Methods

See other articles in PMC that cite the published article.

The damped-oscillator pseudo-wavelet is presented as a method of time-frequency analysis along with a new spectral density measure, the *data power*. An instantaneous phase can be defined for this pseudo-wavelet, and it is easily inverted. The data power measure is tested on both computer generated data and *in vivo* intrahippocampal electrophysiological recordings from a rat. The data power spectral density is found to give better time and frequency resolution than the more conventional total energy measure, and additionally shows intricate time-frequency structure in the rat that is altered in association with the emergence of epilepsy. With epileptogenesis, the baseline theta oscillation is severely degraded and is absorbed into a broader gamma band. There are also broad 600 Hz and 2000 Hz bands which localize to hippocampal layers that are distinct from those of the theta and gamma bands. The 600 Hz band decreases in prominence with epileptogenesis while the 2000 Hz band increases in prominence. The origins of these high frequency bands await further study. In general, we find that the damped-oscillator pseudo-wavelet is easy to use and is particularly suitable for problems where a wide range of oscillator frequencies is expected.

Time-frequency analysis involves monitoring the changes in the frequency spectrum of a system over time. It is of importance in nearly every field of science and engineering. The old work-horse of frequency analysis, the fast Fourier transform (FFT), can be applied to this problem by segmenting data into shorter time windows. By monitoring changes in the power spectrum in each time window, one can determine the presence or absence of activity at certain frequencies within that time window. The segmentation, however, introduces artifact and precludes finer time resolution within time windows. There have been many discussions of these issues and many refinements of the FFT approach (Thomson 1982; Mitra and Pesaran 1999; Groechenig 2001). These refinements are all somewhat labor intensive, requiring judgments to be made about window sizes and windowing functions, which may be different for different frequency ranges. Mathematical generalization of time windowed FFT power spectra to a wide range of time-frequency distribution functions has been accomplished (Cohen 1989; Oonincx 2000). At least one of these time-frequency distribution functions has found application to the study of electroencephalograms (EEG) (Zaveri, Williams et al. 1992). Each of these time-frequency distribution functions has its own strengths and weaknesses.

Over the last 20 years, wavelet transforms have also become an important new way of doing time-frequency analysis (Daubechies 1990; Farge 1992; Graps 1995; Cohen and Kovacevic 1996). A wavelet consists of an oscillatory waveform that has a fairly well-defined frequency and which exists only for a brief period of time. By convolving time series data with a suitably chosen wavelet, one can determine whether an oscillation of a certain frequency is present at a certain interval in time, in a way that is more convenient and less susceptible to artifact than time windowed FFT. There are many kinds of wavelets, and many applications of wavelet transforms to EEG analysis (Le Van Quyen, Foucher et al. 2001; Senhadji and Wendling 2002; Le Van Quyen and Bragin 2007; Allen and MacKinnon 2009; Storti, Formaggio et al. 2009).

A defining property of wavelets is the *admissibility criterion*, a consequence of which is that the mean of the wavelet when averaged over all time must equal zero. This criterion insures that a stable inverse transform exists (Farge 1992; Cohen and Kovacevic 1996). A stable inverse transform is critical for reliable signal transmission and reconstruction. However, in many fields of science and engineering, one is not interested in signal reconstruction. Rather, one may only wish to *detect* whether oscillations of certain frequencies appear, at what times and for what duration of time. For example, in the brain, it is known that oscillations in the theta (4–8 Hz) and gamma ranges (30–100 Hz) are associated with cognitive activity. Oscillations from these two frequency ranges are sometimes phase-coupled, such that the faster gamma frequencies ride entirely on the crests (or troughs) of slower theta rhythms (Lisman 2005; Canolty, Edwards et al. 2006; Axmacher, Henseler et al. 2010). Of increasing clinical interest are high frequency oscillations (HFO’s) in the range of 200–500 Hz. These oscillations tend to occur more frequently in brain regions that are epileptogenic, and so may be useful as a marker for tissue that should be surgically resected in people with refractory epilepsy (Worrell, Gardner et al. 2008; Zijlmans, Jacobs et al. 2009; Bragin, Engel et al. 2010). For these applications, one might consider relaxing the admissibility criterion.

Relaxing the admissibility criterion results in a wide class of waveforms called pseudo-wavelets (or quasi-wavelets), which have found application in the study of turbulence and other complex phenomena (Qiu, Pawu et al. 1995; Wilson, Ostashev et al. 2009). Here we describe a particularly simple pseudo-wavelet that has a uniquely simple physical interpretation. Our pseudo-wavelet corresponds to using a mathematical model of a frictionally damped harmonic oscillator to detect data oscillations of the same frequency. We have previously referred to this pseudo-wavelet as the damped-oscillator oscillator detector (DOOD) (Hsu, Hsu et al. 2007). Here we develop the damped oscillator pseudo-wavelet in more detail and apply it to time series analysis of the EEG. We illustrate our approach with computer generated oscillations as well as data from *in vivo* rat hippocampal recordings. The rat recordings show a richness in time-frequency dynamics which invites further exploration.

The motivation for the DOOD pseudo-wavelet comes from considering how one might construct a mathematical oscillator detector based on analogy to mechanical oscillators. One first constructs a set of mathematical harmonic oscillators each of which has its own natural frequency of oscillation. Next, consider the time series data as a driving force that acts upon each mathematical harmonic oscillator. Start with the mathematical oscillators all at rest, and then let the data drive activity in the mathematical oscillators. If there is an oscillation in the data (a “data oscillator”) that is nearly resonant with one of the mathematical oscillators, then that mathematical oscillator will steadily gain in energy. That mathematical oscillator will continue to gain in energy as long as the data oscillator remains on. If the data oscillator is turned off, the energy of the mathematical oscillator will plateau. If the data oscillator is turned back on again, the mathematical oscillator will either gain or lose energy, depending on whether the data oscillator is in-phase or out-of-phase with the mathematical oscillator. Sudden increases or decreases in energy are both markers of a resonant oscillation in the data. Therefore, we will choose to monitor the time rate of change of the energy of each mathematical oscillator, i.e., its power.

Time resolution can be improved in this scheme by allowing friction to act on the mathematical oscillators. Friction induces the mathematical oscillators to approach a steady state. The steady state will be different depending on whether or not there is resonant driving by data oscillators. The higher the friction, the faster the mathematical oscillator approaches steady state, and the better the time resolution of when data oscillators turn on and off. The cost of better time resolution, of course, is poorer frequency resolution. This counterbalance between time and frequency resolution is a manifestation of the uncertainty principle.

To begin, let *x _{data}* (

Next let *x*(*n*, *t*) represent the coordinate of the *n ^{th}* mathematical oscillator, (

Eq (1)

The simplest non-transient solution for Eq (1) can be written in terms of a complex wavefunction *ψ*(*n*, *t*):

Eq (2)

with
. Taking *ψ _{R}*(

Eq (3)

Eq (4)

This solution can be verified by taking the first time derivative of Eq (2) and expressing *ψ*(*n*, *t*) in terms of *ψ _{R}* (

The wavefunction *ψ*(*n*, *t*) can be calculated easily using a recursion relationship. To derive this relationship, take one incremental time step *δt*:

Eq (5)

For very small time steps *δt*, Eq (5) can be rewritten

Eq (6)

Equation (6) allows us to calculate *ψ*(*n*, *t*+*δt*) rapidly given *h*(*t*+*δt*) and *ψ*(*n*, *t*) without ever having to recalculate *ψ*(*n*, *t*) or to keep track of *h*(*t*) from prior time steps. Note that the exponential factor only has to be calculated once, before taking the first time step. Once we know *ψ*(*n*, *t*+*δt*), we can calculate *x*(*n*, *t*+*δt*) and (*n*, *t*+*δt*) using Eqs (3) and (4).

The instantaneous energy of the *n ^{th}* mathematical oscillator has terms that involve the friction and terms that are independent of the friction. Since we are interested only in changes in the energy due to coupling to

Eq (7)

The first term is the kinetic energy term, and the second term is the potential energy term. Taking the time derivative of Eq (7) and inserting Eq (1), we find the time rate of change of the instantaneous energy, i.e., the instantaneous power:

Eq (8)

The term in Eq (8) containing *g*(*n*) is a term due to friction; we shall ignore it. The term in *h*(*t*) is the power due to the data driving force. We take this term to be our measure of spectral density. To emphasize that this measure is that part of the total power due only to the data driving force, we refer to it as the *data power*. The data power is then written as:

Eq (9)

In comparison, the more traditional wavelet and pseudo-wavelet approaches identify the total energy of the wavefunction *ψ*(*n*, *t*) as the spectral density, which is given as:

Eq (10)

In the absence of friction, the total energy *E*(*n*, *t*) plateaus when there are no active data oscillators. In the presence of friction, the total energy will decay back to baseline when data oscillators are turned off with a rate given by the friction constant *g*(*n*). One might also imagine looking at the total power, given by *Ė*(*n*, *t*). This quantity drops back to baseline faster than the total energy, but it contains a term that depends only on the friction and not on the data oscillator, i.e., the second term in Eq (8). Therefore, neither *E*(*n*, *t*) nor *Ė*(*n*, *t*) is as sensitive to state changes in the data oscillator as is *S*(*n*, *t*). In the first example below, we will see that *S*(*n*, *t*) can show a sudden drop almost immediately after a data oscillator is turned off, within one period of the resonant frequency. Although *S*(*n*, *t*) is still subject to the uncertainty principle, time resolution is typically much better with *S*(*n*, *t*) than with *E*(*n*, *t*).

In summary, the DOOD algorithm consists of the following steps. First, start with all the mathematical oscillators at rest, *ψ*(*n*, *t*)= 0 at *t*= 0, with the frequency of each mathematical oscillator being *f*(*n*)=*ω*(*n*)/2*π*. Second, use Eq (6) to calculate *ψ*(*n*, *t*) at the next instant in time. Third, calculate (*n*, *t*) using Eq (4). Fourth, calculate the data power spectral density using Eq (9). Fifth, one may not wish to save the data power spectral density at every instant in time. For instance, one may wish to know the data power spectral density only at time intervals given by *t _{dump}*. In this case, one can average

The friction constant *g*(*n*) can be used to sharpen the time resolution of the data power spectral density *S*(*n*, *t*). If a data oscillator turns on, drives up the spectral power to a new level, and then turns off, the friction will bring the data power spectral density back down to baseline on a time scale of 1/*g*(*n*). Larger friction will bring down the spectral power faster, thus giving better time resolution as to when the data oscillator turned off. Larger friction will also cause broadening of the spectral peak, however. This manifestation of the uncertainty principle is expected and unavoidable. It can be demonstrated by inspecting Eq (2). If one inserts a pure sinusoid for *h*(*t*) in Eq (2), one finds a Lorentzian line shape with a peak at the frequency of the sinusoid of *h*(*t*) and a half-width-at-half-maximum given by *g*(*n*). The larger is *g*(*n*), the wider is this Lorentzian peak. That is, the spectral bandwidth of the *n ^{th}* mathematical oscillator is

There is no single right answer to this question. In choosing the *g*(*n*)’s, we have to decide whether we wish to have better time resolution or better frequency resolution. One way to choose the *g*(*n*)’s is to consider how one distributes the frequencies of the mathematical oscillators, i.e., how one chooses the *f*(*n*)’s. If one chooses the *f*(*n*)’s on a conventional linear grid, with *f*(*n* +1)= *f*(*n*)+ *δf*, then it is most reasonable to choose *g*(*n*)/2*π*=*δ f*, so that the spectral bandwidths of neighboring mathematical oscillators overlap, leaving no gap between. Here the absolute frequency resolution is given by *δ f*, while the *relative* frequency resolution is given by *g*(*n*)/*ω*(*n*)=*δ f*/*f*(*n*). Note that in this case, the relative frequency resolution goes like ~1/*n*, i.e., the relative spectral bandwidth of each mathematical oscillator *decreases* for higher frequency oscillators.

One might reasonably argue that, rather than setting the absolute frequency resolution constant for every mathematical oscillator, one should set the relative frequency resolution constant, e.g., by taking [*f*(*n* +1) − *f*(*n*)]/*f*(*n*)= *g*_{0}, where *g*_{0} is a small positive number, the same for every mathematical oscillator. Such a choice then implies that the frequency grid is *geometric*, *f*(*n* +1) = *f*(*n*)[1+ *g*_{0}].

When the frequency range that needs to be spanned for a given problem is small, say, of 1 to 2 orders of magnitude, then it does not matter much whether one chooses a linear or geometric frequency grid. However, when the frequency range that needs to be spanned is 3–4 orders of magnitude, or even larger, then computational economy suggests that one should adopt a geometric frequency grid. In this way, one can span the greatest frequency range with the smallest number of mathematical oscillators, while yet maintaining the same relative spectral resolution of each mathematical oscillator.

Here we consider computational cost and compare with the cost of traditional time-windowed fast Fourier transform (FFT). Let *N _{tot}* be the total number of time steps in the entire time series, and let

Eq (11)

The DOOD algorithm Eq (6) yields a power spectrum at every single time step. Its cost therefore goes like

Eq (12)

Comparing Eqs (11) and (12), one sees that FFT is much faster than DOOD for large *N _{f}*. Of course, it is not entirely fair to compare

Eq (11′)

Comparing Eq (12) to Eq (11′), one sees that the DOOD algorithm is computationally *less expensive* than FFT, when time resolution at the level of the sampling rate is desired.

Another situation where the DOOD algorithm will be faster than FFT is where one wishes to monitor a single well defined frequency *f _{data}* with spectral width

It is also possible to implement the DOOD algorithm using Eq (2) directly. In this approach, one segments the total data time series into time windows, as one would for conventional FFT. Within each time window, one can calculate Eq (2) rapidly by taking the FFT of both sides of Eq (2), evaluating Eq (2) in Fourier space, then back transforming to recover Eq (2) in the time domain. The computational cost in this case is of the same order as that of straight FFT. In practice, we find DOOD fast enough using Eq (6) that we do not resort to this hybrid DOOD-FFT algorithm. It is conceivable that for situations where time resolution is not a high priority, that one may adopt a hybrid DOOD-FFT approach.

The sampling rate *f _{S}* of the raw data places an upper limit on the highest frequency oscillator that can be calculated reliably, because if there is an oscillation in the data that is of the same frequency or higher, then it will be artifactually “folded” into the dynamics of oscillators of somewhat lower frequency. This well-known phenomenon applies equally to FFT and to the DOOD pseudo-wavelet. It is advisable to take mathematical oscillators of frequency no higher than

A wavelet transform has the form:

Eq (13)

Here * _{w}* (

Eq (14)

Here * _{w}*(

However, violation of the admissibility criterion does not necessarily preclude an inverse; it merely means that the general inversion formula discovered for admissible wavelets is not applicable. We will come back to this point shortly.

Comparing Eqs (14) with (2), one sees that the DOOD transform is reminiscent of a wavelet transform, if one identifies the pseudo-wavelet kernel as

Eq (15)

Because * _{pw}* (

Nonetheless, for every DOOD pseudo-wavelet *n*, an inverse does exist. For formal completeness, we show how to calculate this inverse. We shall first find *h*(*n*, *t*), which represents the reconstruction of the total signal *h*(*t*) using only the *n ^{th}* damped oscillator. The best reconstruction of the total signal will then be

To proceed, define the one-sided Fourier transform of *ψ*(*n*, *t*) as usual:

Eq (16)

Taking the one-sided Fourier transform of both sides of Eq (2) then yields:

Eq (17)

Eq (18)

where (*n*, *s*) * _{R}* (

Eq (19)

Since the spectral density of interest is the data power, we construct a variety of correlation functions from the data power spectral density. The pair correlation function is defined:

Eq (23)

If we wish to extract phase information, we have to be careful not to introduce phase shifts through the sign of the data power. This kind of artifact is most easily prevented by squaring the data power before taking the time average:

Eq (24)

Eq (25)

Eq (26)

Eq (27)

Equations (24)–(27) can be used to check whether the phase of oscillator *n* modifies the data power amplitude of oscillator *m*. If oscillator *m* tends to be more activated at a particular phase of oscillator *n*, then *σ*_{0}(*m*, *n*) will be relatively larger for that pair of (*m,n*) and *θ*_{0}(*m*, *n*) will give the relative phase of *n* during which the square magnitude of the data power of *m* is largest.

Recall that *x _{data}* (

In general, as we shall see in the examples, one should perform both X-DOOD and V-DOOD when performing time-frequency analysis, in order to explore both low frequency and high frequency structure. However, because V-DOOD requires taking a numerical derivative of the data, there is the potential problem that this numerical derivative may magnify the effects of noise. We shall see that, in practice, the noise magnification introduced by V-DOOD is quite tolerable. In fact, in our examples, V-DOOD does a better job of detecting low frequency oscillations than X-DOOD does of detecting high frequency oscillations, and so V-DOOD may be preferable to X-DOOD as an initial screening tool.

Finally, on the subject of taking time derivatives and magnifying noise, it is noteworthy that the velocity of the *n ^{th}* mathematical oscillator, (

Here we explore the DOOD algorithm with both computer generated data and real data from rats. We will use both X-DOOD and V-DOOD, anticipating that X-DOOD will be more sensitive to low frequency activity and V-DOOD will be more sensitive to high frequency activity.

Consider computer generated data given by the equation:

Eq (28)

Equation (28) describes the dynamics of two sinusoidal oscillators with frequencies of *f*(1)= 7 Hz and *f*(2)= 60 Hz oscillators. The frequencies are chosen to model a single theta (7 Hz) and a single gamma (60 Hz) oscillator. These oscillators are coupled so that the 60 Hz oscillator is more activated on the crests of the 7 Hz oscillator, and less activated in the troughs. Here *R*(*t*) represents a Gaussian random number with root mean square deviation of one. Both the 7 Hz and 60 Hz data oscillators are turned on at time *t*=12 s and off at time *t*=14 s. The sampling rate is taken to be *f _{S}*= 400 Hz which corresponds to time steps of

The X-DOOD and V-DOOD oscillators are chosen with frequencies given by *f*(*n*)= *n* Hz, with *n*=1 to 100. Recall that taking *h*(*t*) = *x _{data}* (

In Fig. 1(a), we plot *h*(*t*) = *x _{data}* (

Damped-oscillator pseudo-wavelet analysis of computer generated data using X-DOOD. In these graphs, the computer generated data driving force *h*(*t*), data power *S*(*f*, *t*) and total energy measure *E*(*f*, *t*) are all rescaled so that the maximum value of each **...**

In Fig. 1(b) we show the data power spectral density of the 7 Hz mathematical oscillator, *S*(7 Hz, *t*), and the corresponding total energy spectral density, *E*(7 Hz, *t*). The friction is taken to be zero, *g*(*n*)= 0 for all *n*. Also shown for reference is *h*(*t*) = *x _{data}* (

Also note the “beating effect” in the amplitude of *S*(7 Hz, *t*) in Fig. 1(b). The higher amplitude beats occur on the crests of the data time series, *h*(*t*) = *x _{data}* (

As an aside, we may also predict that periodic spike trains in the data may be interpreted as a periodic phenomenon that only “pushes” but does not “pull”. In this case, if the spikes occur at a frequency of *f*_{0}, then *S*(*f*_{0}, *t*) will exhibit beats *at the fundamental*, i.e., at *f*_{0}.

To return to Fig. 1(b), note the fine structure on the crests of the data time series, *h*(*t*). This fine structure corresponds to the 60 Hz data oscillator, which as described by the exponential factor in Eq. 28 is amplified on the crests of the 7 Hz data oscillator and damped in amplitude in the troughs. The 60 Hz oscillation shows up clearly on the FFT power spectrum of *S*(7 Hz, *t*) shown in Fig. 1(c). In Fig. 1(c), we see not only a 60 Hz peak, but also satellite peaks at 60 +/− 7 Hz, 60 +/− 14 Hz, etc. The satellite peaks are due to coupling of the 60 Hz oscillation to the 7 Hz oscillation, and are expected. The coupling arises from the exponential factor in Eq. 28. Expansion of this exponential results in sinusoids with frequencies at 60 +/− 7 Hz, 60 +/− 14 Hz, etc.

In Fig. 1(d), we show the data power of the 60 Hz oscillator, *S*(60 Hz, *t*), and the corresponding total energy spectral density, *E*(60 Hz, *t*). The friction is still zero, *g*(*n*)= 0 for all *n*. This figure is similar to Fig. 1(b) except a slow 7 Hz modulation of amplitudes is evident in *S*(60 Hz, *t*), and the fast beats are much more prominent. Since the resonant fundamental frequency for this mathematical oscillator is 60 Hz, we expect to find 120 Hz beats in its power spectrum. The FFT power spectrum shown in Fig. 1(e) confirms that there are indeed 120 Hz beats, as well as associated satellites. It also shows the expected peaks at 7 Hz, its first harmonic at 14 Hz, and a doublet at 60 +/− 7 Hz. Figure 1(f) shows that turning on a bit of friction, *g*(*n*)/2*π*=1 Hz for all *n*, helps damp out the energy trail after *E*(*f*, *t*) at 60 Hz. Time resolution is still much better with *S*(*f*, *t*).

Recall that the X-DOOD algorithm corresponds to taking *h*(*t*) = *x _{data}* (

Time-frequency and phase-amplitude analysis applied to computer generated data. The friction is *g*(*n*)=1 Hz for all *n*, except for (e) and (f) where *g*(*n*)= 0 for all *n*. In these graphs, *S*(*f*, *t*) and *E*(*f*, *t*) are rescaled so that the maximum value is equal to **...**

In Fig. 2(e), we show the phase-amplitude correlation function, *σ*_{0} (*f*_{2}, *f*_{1}) from Eq. (27). If the 60 Hz oscillation tends to be activated at a certain phase of the 7 Hz oscillation, then *σ*_{0}(*f*_{2}, *f*_{1}) will show a peak at *f*_{1}= 7 Hz and *f*_{2}= 60 Hz. The phase frequency is *f*_{1} while the amplitude frequency is *f*_{2}. It is seen from Fig. 2(e) that there is indeed a peak at *f*_{1}= 7 Hz and *f*_{2}= 60 Hz, along with a couple of satellite peaks. There is a stronger peak at *f*_{1} = *f*_{2} = 7 Hz which is the correlation of the 7 Hz oscillation with itself. The phase-amplitude correlation function shown in Fig. 2(e) is calculated with the X-DOOD algorithm with *g*(*n*)= 0 for all *n*. If we were to use the V-DOOD algorithm, the peak at *f*_{1} = 7 Hz and *f*_{2} = 60 Hz would be larger than the self-correlation peak at *f*_{1} = *f*_{2} = 7 Hz, because V-DOOD favors higher frequency peaks. The phase of the 7 Hz oscillator at which the amplitude of the 60 Hz oscillator is largest can be read off the graph for *θ*_{0}(*f*_{2}, *f*_{1}), as defined in Eq. (26). In Fig. 2(f) we see that *θ*_{0}(60 Hz, 7 Hz) ≈ 0, which is the correct phase relationship as shown in Fig. 1(a).

In summary, what we have learned from testing the DOOD algorithm on computer generated data are the following:

- When there are oscillators of multiple frequencies present in a given data set, V-DOOD is better for detecting the higher frequency oscillators than X-DOOD, and using X-DOOD alone may cause one to miss the higher frequency oscillations. In fact, V-DOOD does a fairly good job of detecting low frequency oscillations, and may be preferable to X-DOOD as an initial screening tool.
- The data power measure
*S*(*f*,*t*) gives better time resolution, especially of when an oscillator turns off, than the more conventional total energy measure*E*(*f*,*t*). - The data power measure
*S*(*f*,*t*) exhibits a beating effect when a data oscillator is present. This beating effect may help us distinguish an oscillatory signal from non-oscillatory driving forces when dealing with real biological data. - As a point for future study, it may also be possible to analyze the fine structure of these beats in
*S*(*f*,*t*), in order to extract further information, such as whether there is coupling between data oscillators of different frequencies. - The instantaneous phase
*θ*(*f*,*t*) is a stable, reliable and easily calculated quantity. - The phase-amplitude correlation function,
*σ*_{0}(*f*_{2},*f*_{1}), may be a useful function for detecting phase-modulated coupling of a low frequency oscillation to a higher frequency oscillation. This kind of coupling is thought to be an important mechanism for coding higher order information in the brain (Lisman 2005).

That V-DOOD may be preferable to X-DOOD as an initial screening tool is significant because transforming the raw data represented by the time series *x _{data}* (

A 16-channel silicon probe was inserted into the right hippocampus of two adult Sprague-Dawley rats. Please see Experimental Methods section below for details. The electrode channels span area CA1, the dorsal blade of the dentate gyrus and the hilus of the dentate gyrus. Two weeks after electrode implantation, status epilepticus was induced with repeated small (5 mg/kg) intraperitoneal doses of kainic acid. Field potentials were sampled at a frequency of *f _{S}*=12,207.03 Hz and filtered online between 2 Hz and 6000 Hz. Recordings were saved just prior to injection of kainic acid, immediately after injection of kainic acid, and at weekly intervals out to week 6 after kainic acid injection. By six weeks out, the rats were clearly epileptic, with frequent clinical seizures noted. Each recording was at least 2 minutes in duration. Sample interictal traces from one of these rats at baseline (prior to kainic acid injection) and six weeks after kainic acid injection are shown in Figs. 3(a) and 3(b).

Sample local field potentials recorded from rat hippocampus with an electrode with 16 contacts. The first contact is in the alveus (bottom tracing) and the last contact is in the hilus of the dentate gyrus (top tracing). See Experimental Methods for a **...**

The choice of this experimental system to study was motivated by the focus of our research groups on changes in brain dynamics with epileptogenesis. In this paper, we will not explore all of these changes in detail, but rather our main intent here is to demonstrate the use of the DOOD algorithm. The results of this initial analysis are presented here for one of these rats.

By inspection, it appears that the electrode outputs from channels 2, 4, 6, 12, 14, and 16 are very similar to one another. We shall refer to these as group I channels. We suspect these channels reflect activity in neuronal layers, although they do not match up exactly with later histology (see the Experimental Methods section below). Similarly, the outputs from the remaining channels are more or less similar to one another. We suspect these channels are in the non-neuronal layers, and we will refer to them as group II channels. Detailed analysis of all 16 channels through the six weeks after kainic acid injection will be presented in a future paper. For now, we show results from analyzing channels 16 and 1 only. Based on histology, we believe channel 1 to be in the alveus while channel 16 is in the hilus of the dentate gyrus.

For this data, we wish to cover a wide frequency range, and so we use a geometric frequency grid, as discussed earlier in Section 2.1. The frequency of the *n ^{th}* mathematical oscillator is given by

Shown in Fig. 4(a) and 4(b) are *S*(*f*, *t*) using V-DOOD with *t _{dump}*=1 sec for channels 16 and 1 at baseline, just prior to kainic acid injection. The theta oscillation at 7 Hz is the most prominent finding. There is also a broad gamma frequency band centered at 60 Hz, and intermittently there are bursts of higher frequency activity centered at 600 Hz. A yet higher frequency band centered at 2000 Hz may also be discerned. The theta and gamma bands are most prominent in channel 16 while the 2000 Hz band is better seen in channel 1. The 600 Hz band is seen clearly in both channels 16 and 1. The 600 Hz activity appears to be intermittent while the 7 Hz, 60 Hz and 2000 Hz activity appear to be nearly continuous.

Time-frequency V-DOOD spectra from intrahippocampal rat recordings at baseline, prior to kainic acid injection. In these graphs, *S*(*f*, *t*) is rescaled by subtracting out the mean and dividing by the standard deviation. The mean and standard deviation are **...**

To appreciate the finer temporal structure, we show in Figs. 4(c) and 4(d) the square of the data power, *S*^{2}(*f*, *t*), averaged over shorter time windows of duration *t _{dump}*=5 msec. Squaring the data power before time averaging prevents large cancellations between positive and negative variations in mathematical oscillator energy. Seen in Figs. 4(c) and 4(d) is a strong 7 Hz oscillation which now exhibits the expected “beating” effect. The gamma band can now also be seen to be composed of many distinct gamma frequency oscillations which come and go. The time-frequency patterns created by the theta and gamma frequencies are quite intricate. Figure 4(d) is a magnification of the burst of 600 Hz activity at around time

Are the 600 Hz and 2000 Hz bands “real” or could they be due to external artifact? Although Fig. 4(d) shows that the burst of 600 Hz activity is time correlated with increased activity at both 7 Hz and in the 60 Hz band, note the absence of a “beating” effect in the 600 Hz band, as compared with the 7 Hz and 60 Hz activity. This finding suggests that the activity in the 600 Hz band is not oscillatory. The 600 Hz band may represent simply fast, non-oscillatory activity.

Figure 4(e) shows a magnification of the local field potential 600 Hz burst near time t = 56.4 sec, depicted in the time domain, which shows only a high frequency “buzz” in all 16 channels at the time of the 600 Hz burst. No clear evidence of external artifact is present around this time. Figure 4(f) shows the data power averaged over the first 60 secs and plotted as a function of frequency and channel number. In this figure, it is clear that the theta and gamma oscillations co-localize to the same hippocampal layers (which we suspect to be the neuronal layers), while the 600 Hz and 2000 Hz activity co-localize to a different set of layers. The only exception is channel 6, which shows activity from all four frequency bands. If the 600 Hz or 2000 Hz band activity were from external noise, one would not expect localization of these bands to distinct hippocampal layers, but rather all the layers should be involved to a similar extent, or there should be a unidirectional gradient in amplitude. The 600 Hz band may be related to the high frequency oscillations (HFO’s) observed by others (Worrell, Gardner et al. 2008; Zijlmans, Jacobs et al. 2009; Bragin, Engel et al. 2010). We do not know what the 2000 Hz band could possibly correspond to. It is much too fast to be related to neuron-to-neuron synaptic communication. It might be related to fast local ion flows, possibly in the dendritic trees, neuronal cell membranes or in the distal axonal arborizations. Future work is necessary to sort this out.

In Figs. 5(a) and 5(b) we show results of finite non-overlapping time-windowed FFT time-frequency analysis applied to the same data, with time windows chosen to be 0.67 sec in duration. Prior to application of the fast Fourier transform, the data is first converted to its first time derivative, i.e., FFT is performed on *h*(*t*) = * _{data}* (

V-FFT time-frequency power spectra with *h*(*t*) = _{data} (*t*) from channel 16 from the rat prior to kainic acid injection. Non-overlapping time windows of 0.67 sec were used, with a total time interval of 60 secs and frequency resolution of 1.49 Hz. **...**

In Fig. 6(a) we show *S*(*f*, *t*) averaged over 60 seconds from the same rat 6 weeks after kainic acid injection. At this point, the rat exhibits frequent spontaneous seizures. Figure 6(a) shows that the 7 Hz theta oscillation has been severely degraded as to be virtually undetectable. The gamma band has expanded to absorb the theta oscillation and there is no longer a clear separation between theta and gamma bands. The 600 Hz band has also been degraded, while the 2000 Hz band has grown in prominence, but note that the 2000 Hz band is still localized to distinct layers of the hippocampus, being more prominent in layers in which the gamma oscillation is less prominent.

Time-frequency V-DOOD spectra with *h*(*t*) = _{data} (*t*) from intrahippocampal rat recordings at six weeks after kainic acid injection. In these graphs, *S*(*f*, *t*) is rescaled by subtracting out the mean and dividing by the standard deviation. (a) *S*( **...**

Figure 6(b) shows *S*(*f*, *t*) for channel 16 averaged over time windows of duration *t _{dump}*= 250 msec. We see that activity is now dominated by gamma frequency discharges rather than theta. Figure 6(c) shows

To look more closely at the finer time-frequency structure, we plot the square of the data power, *S*^{2}(*f*, *t*), averaged over time windows of duration *t _{dump}*= 4 msec, shown in Figs. 6(d) and 6(e). Again one sees beautifully intricate time-frequency patterns spanning theta through gamma frequencies, but now there is not a dominant theta rhythm, and there is no clear boundary between the theta and gamma bands. It is as though gamma oscillations have invaded into theta territory and disrupted the normally dominant theta oscillation. The loss of clear separation between theta and gamma bands and the degradation of the dominant theta rhythm with epileptogenesis may impair the ability of the epileptic brain to couple gamma oscillations to the theta rhythm. If such theta-gamma coupling is important for cognitive performance (Lisman 2005), then the disruption of theta-gamma coupling may contribute to cognitive impairment in epileptic rats.

It may also appear that the gamma oscillations are not as well formed in Fig. 6(d) as in Fig. 4(c). In Fig. 4(c), each gamma oscillation is a clean vertical line, but in Fig. 6(d), the beginnings of each oscillation is a triangular clump, and each oscillation peters out quickly. We are not certain if this finding will hold up with more detailed studies, but we will return to this point in future studies.

We have also noted that the squared data power measure, *S*^{2}(*f*, *t*), appears to be more sensitive to low frequency structure and not as sensitive to high frequency structure. If one were to plot the Figs. 6(d) and 6(e) on a logarithmic scale, it is possible to see the high frequency band at 2000 Hz, but this band appears only as a small hump. It seems that *S*(*f*, *t*) is more sensitive to high frequency activity than *S*^{2}(*f*, *t*).

In Figs. 7(a) through 7(f), we show non-overlapping time windowed V-FFT results at 6 weeks after status epilepticus using time windows of 0.67 sec. These figures are in agreement with those using *S*(*f*, *t*), showing loss of the theta rhythm and expansion of the gamma band down to theta frequencies. The 2000 Hz band has grown in prominence, but it is most prominent in those channels where the gamma rhythm is weakest. However, Figs. 7(c) and 7(e) shows that V-FFT is rather poor at showing the fine time-frequency structure that is so striking when using V-DOOD *S*^{2}(*f*, *t*).

In this paper, our main goal is to introduce the DOOD algorithm. More detailed analysis of the time-frequency spectra of *in vivo* rat (and human) data will be presented in future publications, including more detailed study of the 600 Hz and 2000 Hz bands. Nonetheless, the results shown here suggest a striking richness of time-frequency structure in the normal rat, which is altered after the rat undergoes epileptogenesis. If we believe that oscillations underlie the function of the brain (Buzsaki 2006), then the fine details of the behaviors of these oscillations – when they arise, where they arise, for how long and in what relation to other oscillations – must contain interesting information about how the brain works and what happens to it in disease. It would also be of interest to know whether the 2000 Hz activity is reproducible in other experiments, and if so, to determine from where it arises. Is it generated in the dendritic tree, in cell membranes, or from somewhere else? We were initially inclined to dismiss it as artifact, but the fact that the 2000 Hz activity is most prominent in specific hippocampal layers suggests that it may not be due to artifact, or at least not due to artifact external to the hippocampus. We do have data from a second rat, prepared the same way as described here, for which 2000 Hz activity is also seen and which also grows in prominence after status epilepticus.

Others have described the same pseudo-wavelet presented here but without identifying its relationship to damped harmonic oscillators (Szu, Telfer et al. 1992; Vazquez, Mazilu et al. 2005). The spectral density tested by those authors was the total energy measure, *E*(*f*, *t*), not the data power measure, *S*(*f*, *t*), nor the squared power measure, *S*^{2}(*f*, *t*), introduced here. We expect the data power and squared power measures to be capable of better time and frequency resolution. In application to an example consisting of computer generated data, we showed that the data power measure yields better time and frequency resolution than the total energy measure. As shown in Fig. 2, the time resolution even at zero friction is remarkable. The better time resolution allows us to use lower levels of friction which in principle should yield an improved frequency resolution as well.

Does the data power measure of the damped-oscillator pseudo-wavelet give better time and frequency resolution than other wavelets or pseudo-wavelets? We have not compared the damped-oscillator pseudo-wavelet directly with other wavelets. However, unless a wavelet or pseudo-wavelet is especially chosen to match a particular application, there are not usually dramatic differences in time and frequency resolution between them. Since all other wavelet and pseudo-wavelet methods define spectral density in terms of the equivalent of the total energy measure, it is possible that the data power measure introduced here may be superior to them as well.

The velocity form of the damped-oscillator pseudo-wavelet (V-DOOD) is more sensitive to high frequency activity while the coordinate form (X-DOOD) is more sensitive to low frequency activity. This property is common to all methods of time-frequency analysis; that is, all time-frequency methods have analogous coordinate and velocity forms with the same relative sensitivities to low and high frequency activity. It is generally good practice to scan data with both forms so as not to miss either low or high frequency activity. However, our experience has been that V-DOOD is fairly good at detecting even low frequency activity while X-DOOD is relatively insensitive to higher frequency activity. Thus we favor V-DOOD over X-DOOD, at least as an initial screening tool.

In application to *in vivo* rat data, we found that the squared data power measure shows intricate time-frequency structure at low frequencies. The square of the data power brings out this structure somewhat better than the data power itself.

In summary, we have developed a novel pseudo-wavelet approach for time-frequency analysis that is intuitively appealing, easy to use, and provides striking time-frequency resolution. The damped-oscillator pseudo-wavelet is particularly useful when multiple time scales are involved spanning several orders of magnitude or more, and when one wishes to have maximal flexibility in deciding which frequency ranges to monitor. Our preliminary results suggest a richness of time-frequency structure in local field potentials in a rat brain undergoing epileptogenesis, that invites further study.

Adult, Spraque-Dawley rats were anesthetized with 2% isoflurane after pretreatment with atropine (2 mg/kg) and placed in a stereotaxic apparatus. The area of incision was injected with 0.5 ml of 0.5% bupivacaine for prolonged local anesthesia/analgesia and burr holes for the recording probe and ground screw were placed by conventional surgical techniques. The depth-recording electrode, a 16-channel silicon probe (NeuroNexus, Inc., Ann Arbor, MI), was implanted in the right hippocampus (3.0 mm posterior, 2.6 mm lateral, 3.1 mm ventral), spanning CA1, dorsal dentate gyrus, and the hilus at 100 micron intervals. A skull screw was placed in the left frontal bone to serve as ground. Dental acrylic was used to secure the electrodes according to standard chronic methods. Postoperative analgesia was administered once daily (buprenorphrine 0.01 mg/kg, subcutaneous). At two weeks following electrode implantation, kainic acid (Tocris Bioscience, Ellisville, MO) was administered in repeated, low doses (5 mg/kg, intraperitoneal) hourly until each rat experienced convulsive status epilepticus for >3 h (Hellier, Patrylo et al. 1998).

Simultaneous *in vivo* field potential recordings from 16-channel silicon probes were obtained in freely-moving rats during each recording session. Field potentials were recorded at a sampling frequency of 12,207.03 Hz and filtered on-line between 2 Hz and 6000 Hz using a Tucker Davis Technologies recording system (Alachua, FL). Off-line conversion of the electrographic data to f32 binary files was accomplished using custom-written software developed in MATLAB® (The MathWorks, Inc., Natick, MA).

After completion of long-term recordings epileptic rats were perfused with 0.01M phosphate buffered saline followed by an aqueous solution of aldehyde fixatives (4% paraformaldehyde-0.5% glutaraldehyde). The brains were removed, post-fixed overnight in same fixative solution, and sectioned on a Pelco Vibratome 1500 sectioning system into 100-μm coronal sections. The sections were wet-mounted in 0.9% saline and imaged with a digital camera (Spot II; Diagnostic Instruments, Sterling Heights, MI) on a Nikon E600 Eclipse epifluorescent microscope with ×1–4 planapochromatic objectives in order to locate and measure the recording tracts. Brightfield images were acquired at an initial 36-bit tone scale and saved as 16-bit files. The final images were prepared in Adobe Photoshop 7.0. Adjustments in tone scale, gamma, contrast, and hue and subsequent sharpening with the unsharp mask algorithm were applied to the entire image. The 16 recording sites along the probe length were reconstructed and associated with relevant hippocampal structures using histological measurements and the known inter-site intervals. Electrodes 1–3 were in stratum oriens of CA1, with electrode 1 probably in alveus. Electrodes 4–5 were in the pyramidal cell layer. Electrodes 6–7 were in stratum radiatum. Electrodes 8–9 were in the lacunosum moleculare layer. Electrodes 10–11 were in the molecular layer of the dentate gyrus. Electrodes 12–13 were in the dentate granule cell layer. Electrodes 14–16 were in the hilus of the dentate gyrus.

Institutional regulations for the care of all animals were strictly adhered to. Animals that are persistently ill after induction of seizures or that suffer frequent, prolonged seizures are euthanized.

- Time-frequency analysis using mathematical oscillators to detect data oscillators
- High time-frequency resolution
- Fast and convenient to use
- Tested on computer generated and
*in vivo electrophysiological*data

DH was supported by the NIH Loan Repayment Program. TPS was supported by NIH R01-25020. GAW was supported by NIH R01-NS63039-01.

- DOOD
- damped-oscillator oscillator detector
- FFT
- fast Fourier transform

**Publisher's Disclaimer: **This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final citable form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

- Allen DP, MacKinnon CD. Time-frequency analysis of movement-related spectral power in EEG during repetitive movements: A comparison of methods. J Neurosci Methods. 2009;186(1):107–115. [PMC free article] [PubMed]
- Axmacher N, Henseler MM, et al. Cross-frequency coupling supports multi-item working memory in the human hippocampus. Proceedings of the National Academy of Sciences. 2010;107(7):3228–3233. [PubMed]
- Berne BJ, Pecora R. Dynamic Light Scattering. New York: Wiley; 1976.
- Bragin A, Engel JJ, et al. High-frequency oscillations in epileptic brain. Current Opinion in Neurology. 2010;23(2):151–156. [PubMed]
- Buzsaki G. Rhythms of the Brain. New York: Oxford University Press; 2006.
- Canolty RT, Edwards E, et al. High gamma power is phase-locked to theta oscillations in human neocortex. Science. 2006;313(5793):1626–8. [PMC free article] [PubMed]
- Cohen A, Kovacevic J. Wavelets: The mathematical background. Proceedings of the IEEE. 1996;84(4):514–522.
- Cohen L. Time-frequency distributions - a review. Proceedings of the IEEE. 1989;77(7):941–981.
- Daubechies I. The wavelet transform, time-frequency localization and signal analysis. IEEE Transactions on Information Theory. 1990;36(5):961–1005.
- Farge M. Wavelet transforms and their application to turbulence. Annual Review of Fluid Mechanics. 1992;24:395–457.
- Graps A. An introduction to wavelets. IEEE Computational Science & Engineering. 1995;2(2):50–61.
- Groechenig K. Foundations of Time-Frequency Analysis. Boston: Birkhauser; 2001.
- Hellier JL, Patrylo PR, et al. Recurrent spontaneous motor seizures after repeated low-dose systemic treatment with kainate: assessment of a rat model of temporal lobe epilepsy. Epilepsy Res. 1998;31(1):73–84. [PubMed]
- Hsu D, Hsu M, et al. An algorithm for detecting oscillatory behavior in discretized data: the damped-oscillator oscillator detector. 2007. Aug 9, from arxiv.org/abs/0708.1341.
- Le Van Quyen M, Bragin A. Analysis of dynamic brain oscillations: methodological advances. Trends Neurosci. 2007;30(7):365–73. [PubMed]
- Le Van Quyen M, Foucher J, et al. Comparison of Hilbert transform and wavelet methods for the analysis of neuronal synchrony. J Neurosci Methods. 2001;111(2):83–98. [PubMed]
- Lisman J. The theta/gamma discrete phase code occuring during the hippocampal phase precession may be a more general brain coding scheme. Hippocampus. 2005;15(7):913–22. [PubMed]
- Mitra PP, Pesaran B. Analysis of dynamic brain imaging data. Biophys J. 1999;76(2):691–708. [PubMed]
- Oonincx PJ. Mathematical Signal Analysis: Wavelets, Wigner Distribution and a Seismic Application. Amsterdam: Centrum voor Wiskunde en Informatica; 2000.
- Qiu J, Pawu UK, et al. Pseudo-wavelet analysis of turbulence patterns in 3 vegetation layers. Boundary-Layer Meteorology. 1995;72(1–2):177–204.
- Senhadji L, Wendling F. Epileptic transient detection: wavelets and time-frequency approaches. Neurophysiologie Clinique/Clinical Neurophysiology. 2002;32(3):175–192. [PubMed]
- Storti S, Formaggio E, et al. Wavelet Analysis as a Tool for Investigating Movement-Related Cortical Oscillations in EEG-fMRI Coregistration. Brain Topography. 2009;23(1):46–57. [PubMed]
- Szu HH, Telfer B, et al. Causal analytical wavelet transform. Optical Engineering. 1992;31:1825–9.
- Thomson DJ. Spectrum estimation and harmonic analysis. Proceedings of the IEEE. 1982;70(9):1055–1096.
- Vazquez JM, Mazilu M, et al. Wavelet transforms for optical pulse analysis. Journal of Optical Society of America A. 2005;22:2890–9. [PubMed]
- Wilson DK, V, Ostashev E, et al. Quasi-wavelet formulations of turbulence and other random fields with correlated properties. Probabilistic Engineering Mechanics. 2009;24(3):343–357.
- Worrell GA, Gardner AB, et al. High-frequency oscillations in human temporal lobe: simultaneous microwire and clinical macroelectrode recordings. Brain. 2008;131(Pt 4):928–37. [PMC free article] [PubMed]
- Zaveri HP, Williams WJ, et al. Time-frequency representation of electrocorticograms in temporal lobe epilepsy. IEEE Trans Biomed Eng. 1992;39(5):502–9. [PubMed]
- Zijlmans M, Jacobs J, et al. High frequency oscillations and seizure frequency in patients with focal epilepsy. Epilepsy Res. 2009;85(2–3):287–92. [PubMed]

PubMed Central Canada is a service of the Canadian Institutes of Health Research (CIHR) working in partnership with the National Research Council's Canada Institute for Scientific and Technical Information in cooperation with the National Center for Biotechnology Information at the U.S. National Library of Medicine(NCBI/NLM). It includes content provided to the PubMed Central International archive by participating publishers. |