Hilbert transform for brain waves

• Dr. Walter J. Freeman, University of California, Berkeley, California

State variables from brain potentials are time series that are either recorded and digitized or derived from recordings by use of the Fourier and Hilbert transforms; they provide the primary raw materials by which models of brain dynamics are constructed and evaluated.

Introduction

Digitized records of continuously fluctuating potentials of the brain are signs of brain functions. From depth electrodes they are local field potentials (LFP). From cortical surface electrodes they are electrocorticograms (ECoG). From scalp electrodes they are electroencephalograms (EEG). For convenience "EEG" refers to all three. They provide useful mesoscopic correlates of intentional behaviors, of the rates and intervals of microscopic pulse trains, of "MEG" referring to magnetic potentials and magnetoencephalography, and of macroscopic images like those from fMRI. In the context of modeling brain dynamics they help to define states and state variables of brains.

Common practice is to assume linearity and stationarity in the dynamics of their sources, and to decompose the signals by use of PCA, ICA, ARMA, etc. into derivative state variables. Prominent among these transforms is the fast Fourier transform (FFT), which decomposes each raw signal into a family of state variables having fixed frequencies and amplitudes (limit cycles). Closely related are wavelets that are used for linear decomposition of signals varying in amplitude. The temporal resolution of wavelets and the FFT is bounded by the Nyquist criterion: the digitizing rate must be at least twice and preferably three times the highest component frequency. The duration of segments for decomposition must exceed at least one cycle of the lowest component frequency.

Brain dynamics revealed by state variables derived from EEG is noisy, nonstationary, nonlinear, and rife with temporal discontinuities from state transitions. Though the Hilbert transform (HT) like the FFT is a linear operator, it is useful for analyzing nonstationary signals by expressing frequency as a rate of change in phase, so that the frequency can vary with time. Typically multiple time-varying frequencies coexist in raw recordings. Empirical mode decomposition (EMD), the Hilbert-Huang transform (Huang and Shen, 2005), gives high spectral resolution of arbitrary frequencies. More useful for EEG is 'clinical mode decomposition' (CMD) by band pass filtering to decompose raw signals into components corresponding to the divisions of the clinical spectrum. During epochs between discontinuities in these frequencies the oscillations briefly approach limit cycles with modulation of center frequencies in amplitude (AM), frequency (FM) and phase (PM). These state variables provide sufficient temporal resolution to locate state transitions and document the relative time-invariance of recurring brief epochs, during which state variables may approach linearity and stationarity (Freeman, 2004a,b).

The HT and FFT give the same results when both transforms are applied to signals having the relatively long durations needed for the FFT and wavelets (Le Van Quyen et al., 2001; Quiroga et al., 2002). Otherwise they are complementary. The FFT and EMD give high frequency resolution. The HT and CMD give high temporal resolution of rapid changes in analytic state variables for frequency, phase and amplitude. Decomposition by either FFT or HT gives a family of state variables for each selected pass band. Each state variable defines an axis in the state space of a dynamic model, which is a finite projection from infinite brain state space. The ranges of possible observed values on the axes define the limits of the state space.

The further advantage of the HT is the sensitive access it gives to AM patterns of analytic amplitude that are correlated with intentional behaviors (Freeman, 2005). The greatest difficulty in using the HT is to distinguish physiological state transitions from spurious discontinuities in the analytic phase known as phase slip (Pikovsky, Rosenblum and Kurths, 2001). Phase slip occurs by interference between multiple overlapping signals with differing frequencies and AM patterns. EEG signals commonly have power-law distributions (see 1/f noise) in power spectral densities (PSD) that are not well handled by the FFT. Application of the HT to unfiltered EEG gives analytic phase values resembling a random walk. Band pass filtering using the FFT is necessary to get time series from the HT that are interpretable as state variables. The passbands of clinical EEG from 80 years of experience are delta 1-3 Hz, theta 3-7 Hz, alpha 7-12 Hz, beta 12-30 Hz, gamma 30-80 Hz, epsilon 80-250 Hz.

Calculation of the analytic signal from EEG with the Hilbert transform

Figure 1: Application of the HT to 64-channel EEG in an 800 ms epoch of visual cortical EEG centered at delivery of a CS at 0 ms. From Freeman [2004a], Fig. A1.03.

The analytic signal, \(V(t)\ ,\) is calculated from a filtered EEG by applying the Hilbert transform (Barlow, 1993; Pikovsky, Rosenblum and Kurths, 2001; Freeman, Burke and Holmes, 2003). The EEG of each channel, \(j\ ,\) denoted \(v_j(t)\ ,\) from a recording through a complete behavioral trial is transformed to a time series of complex numbers, \(V_j(t)\ ,\) with a real part, \(v_j(t)\ ,\) and an imaginary part, \(u_j(t)\ ,\) \[\tag{1} V_j(t) = v_j(t) + \imath u_j(t) \ ,\]

\(j = 1, 64\ ,\)

where the real part is the same as the filtered EEG ( Figure 1, A, blue curve), and the imaginary part is from the Hilbert transform of \(v_j(t)\ ,\) which is \(u_j(t)\) (red curve) given by equation (2) corresponding to the negative rate of change of \(v_j(t)\ ,\) the quadrature:

\[\tag{2} u_j(t) = \frac{1}{\pi} \mbox{ PV } \int_{-\infty}^{+\infty} V_j(t')/(t - t')\, dt' \ ,\]

where PV signifies the Cauchy Principal Value. At each digitizing step the EEG yields a point in a polar plot in the complex plane ( Figure 1, B). Details on the computation using the FFT can be found in MATLAB (see "help hilbert").

Sequences of digitized values give a trajectory of the tip of a vector rotating counterclockwise in the complex plane with elapsed time. The vector length at each digitizing step, \(t\ ,\) is the state variable for analytic amplitude (the blue curve in C), \[\tag{3} A_j(t) = \sqrt{v_j^2(t) + u_j^2(t)} \ ,\]

and the arctangent of the angle of the vector with respect to the real axis is the state variable for analytic phase (blue saw tooth in D), \[ P_j(t) = \mbox{atan}\, \frac{u_j(t)}{v_j(t)} \ ,\] \(j = 1, \ldots, 64\ .\)

The arctangent is calculated using either the two-quadrant inverse tangent function (atan in MATLAB) or the four-quadrant inverse tangent function (atan2). The atan function gives phase values that range from \(-\pi/2\) to \(\pi/2\) and on reaching \(\pi/2\) fall to \(-\pi/2\) twice in each cycle, while the atan2 function of the same data ranged from \(-\pi\) to \(\pi\) and on reaching \(\pi\) fall to \(-\pi\) once each cycle. Both give a sawtooth curve with small increments of phase lag along a diagonal from the lower bound to the upper bound, and a downward fall in one step to the lower bound on reaching the upper bound.

In order to track \(P_j(t)\) over arbitrarily long time intervals the disjoint phase sequences are straightened by adding \(\pi\) to the atan function or \(2\pi\) to the atan2 function at each reset (branch point) to get the unwrapped analytic phase, \(p_j(t)\ .\) The slope of the ramp in rad/s gives the mean frequency over the duration of the time interval fitted with a line segment. The center frequency is the peak frequency in the FFT of the filtered signal that is determined mainly by the pass band of the filter. The upward and downward deviations from the mean demarcate either phase slip or state transitions. The challenge is to distinguish them.

Figure 2: The temporal power spectral density (PSDT) conforms to \(1/f^\alpha\) except for excess power in the theta (3-7 Hz) and beta-gamma (20-80 Hz) ranges. From Freeman [2004a].

Optimization of the temporal band pass filter required for use of the clinical HT

Owing to the origin of background EEG in self-organized activity among neurons forming random nets, the canonical form of the PSD of EEG is \(1/f^\alpha\) in log-log coordinates [Freeman, 2006], where the slope of the power-law relation in subjects at rest and disengaged from the environment ranges between -2 ("brown noise", exponent \(\alpha = 2\)) and -3 ("black noise", exponent \(\alpha = 3\) [Schroeder, 1991]). Behavioral arousal in tasks brings peaks above the linear slope ( Figure 2), which indicate the emergence of order from noise, and which are seen in the empirical ranges long established in clinical EEG (here theta, 3-7 Hz, beta (13-30 Hz) and gamma (30-80 Hz). The simplest method to optimize the pass band is to bracket peaks in the PSDT with a band pass filter (20-80 Hz in Figure 2 and also for the curves in Figure 1, A). Each pass band gives a different set of state variables.

A computationally intensive method for optimization of a pass band relies on multiple simultaneously recorded signals from an array of closely spaced electrodes. The successive differences, \(\Delta p_j(t)\ ,\) divided by the digitizing interval, \(\Delta t\ ,\) approximate the analytic frequency, \(\omega_j(t)\) (a state variable!), in rad/s (Hz when divided also by \(2\pi\)rad/cycle): \[ \omega_j(t) = [p_j(t)-p_j(t-1)]/\Delta t = \Delta p_j(t)/\Delta t \ .\] A raster plot (Figure 3) of \(\Delta p_j(t)\) from 64 electrodes in an 8x8 array (5.6x5.6 mm) shows that \(\omega_j(t)\) alternates between relatively flat epochs (nearly constant frequency) and short periods of high or low (even negative) phase differences. The jumps are coordinated across the array in timing but not in direction. The instantaneous spatial standard deviation, \(SD_X(t)\ ,\) of \(\omega_j(t)\) or \(\Delta p_j(t)\) gives a statistic that is useful for locating times of onset of nearly stationary epochs for further study of analytic amplitudes (see Figure 5).

The reason for the intermittently high values of \(SD_X(t)\)is that each spatial AM pattern is separated from those before and after by state transitions. Each AM pattern is accompanied by a spatial PM (phase modulation) pattern in the form of a cone. This phase cone reveals a property common to distributed systems, that a state transition is not everywhere instantaneous but begins at a site of nucleation and spreads concentrically, like the formation of a snowflake around a grain of dust. The apex that marks the site in cortex is a random variable both in sign (lead or lag) and location. The spatial gradient in rad/mm combined with the carrier frequency in rad/s gives the phase velocity in m/s, which is comparable to the propagation velocities of cortical axons mediating the state transition (Freeman, 2004b). The oscillation at each recording site shows a discontinuity as the phase is re-set, but the direction of re-setting is different as one cone is replaced by its successor (Freeman, 2005), usually with a differing carrier frequency in the beta or gamma band (Freeman, 2006). Of course the Hilbert transform of a continuous signal always gives increasing phase values (no negative frequencies), but the discontinuities from phase re-setting allow for momentary negative phase differences.

Figure 3: The raster plot shows the successive differences of the unwrapped analytic phase, \(\Delta p_j(t)\ ,\) changing with time (left abscissa) and channel (right abscissa) from Figure 1, A. The 8 columns of 8 rows are aligned to show near-coincidence of sudden jumps or dips by fast-forward or backward rotation of the vector in Figure 1, B). When the jumps or dips are aligned with the right abscissa, they have a distribution that is centered at zero with a range of \(\pm\pi/8\) (45°), the half-power level from phase dispersion. From Freeman [2004a], Fig. A1.04.

The peaks in \(SD_X(t)\) recur at rates in the alpha (8-12 Hz) and theta ranges. This is shown by cross correlating the unfiltered EEG with \(SD_X(t)\ ,\) calculating the cospectrum, finding the frequency at which cospectral power in the alpha band is maximal, and measuring that power. Optimization of the pass band is calculating \(SD_X(t)\) repeatedly as the high cut-off filter is applied at frequencies varied in steps from 80 Hz down to 12 Hz to construct a tuning curve of maximal alpha power as a function of cut-off frequency (Freeman, Burke and Holmes, 2003). The cut-off frequency giving the highest power in the alpha range is fixed at the optimal high cut-off filter setting, and the low cut-off filter is applied in steps from 4 to 20 Hz to construct a second tuning curve to locate maximal cospectral power and specify the optimal low cut-off filter setting.

The criterion most relevant to intentional behaviors, but computationally most intensive for optimizing the upper and lower cut-off frequencies, is the level of correct classification of spatial AM patterns of brief epochs of beta or gamma oscillations with conditioned stimuli (CS) (Freeman, 2005, 2006).

Additional constraints on applying the HT to EEG are illustrated in Figure 4. In A the phase of a 20 Hz cosine in 2 ms steps is reset at 3.5 rad (28/50 ms = 0.56 cycle) to simulate a discontinuity before completion of the first cycle. In B the unwrapped phase, \(\Delta p(t)\ ,\) shows that the change to a new ramp is distributed over ±5 digitizing steps. In C the phase differences, \(\Delta p(t)\ ,\) mark the time of reset, but with distortion over several time steps that strongly resembles the end effects. Similar distortions occur in the analytic amplitude, A(t), which should be constant. This simulation shows that \(\Delta p(t)\) suffices to locate discontinuities in analytic phase, but that reliable measurements of analytic amplitude and phase cannot be made near phase slips, nor at or close to the ends of digitized segments.

Figure 4: Constraints are illustrated on the use of HT at and near phase discontinuities and ends of segments. From Freeman [2004a], Fig. A1.11.

In D the unwrapped analytic phase values, \(p_j(t)\ ,\) of 64 EEG signals after band pass filtering (20-80 Hz) are distributed about a mean phase at time zero. With elapsed time the \(p_j(t)\) diverge, forming parallel sloping lines separated by \(\pm\pi\) rad. The diverging intercepts are either due to residual low frequency oscillations having a different spatial distribution from that of the center frequencies, so that on some cycles of \(v_j(t)\) the peak of oscillation fails to reach the zero from above or below, which prevents an expected branch point, or they are due to residual high frequency oscillations that throw in extra zero crossings and branch points. The divergence is exacerbated when analytic amplitude approaches zero, as it does at state transitions shown between 60 and 80 ms. The main use of \(\Delta p_j(t)\) phase differences is to locate the onsets of state transitions by \(SD_X(t)\ .\)

Comparisons of unwrapped analytic phase across multiple signals simultaneously recorded from arrays after band pass filtering are contraindicated over windows exceeding the cycle duration of the center frequency minus the phase distribution.

Time-varying state variables derived from EEG with the HT to evaluate cortical states

Figure 5: Examples of 4 of 8 new parameters from the HT. From Freeman (2004a, Fig. A1.09).

The waveforms of the EEG from a high-density array have varying degrees of similarity, depending on the degree of spatial coherence among them. A simple measure of global coherence that avoids time-lagged correlation and derivation of phase distributions at multiple frequencies is given by the ratio, \(R_e(t)\ ,\) which is given by \(SD_T(t)\) of the mean waveform divided by the mean \(\underline{SD}_T(t)\) of the several waveforms (Figure 5, A, red curve). \[ R_e(t) = SD_T(t)/\underline{SD}_T(t) \ ,\] where the moving window used to calculate the \(SD_T\) has twice the duration of the wavelength of the center frequency of the pass band for the HT, and time \(t\) is at the center. If the waveforms all have the same shape and phase, the ratio is one, even though the amplitudes differ. If the signals are all independent, the ratio is \(1/n^{0.5}\ ,\) where \(n\) is the number of time points in the moving window.

The reciprocal, \(1/R_e(t)\ ,\) is compared with another measure of coherence (Figure 5, A, black spikes), \(SD_X(t)\) of the spatial array of \(\Delta p_j(t)\) shown in Figure 3 which, divided by the digitizing step in \(s\ ,\) approximates the analytic frequency, \(\omega(t)\ .\) The mean, \(\omega(t)\ ,\) and SD_X(t) are highly correlated. The two measures of coherence are related but independent, as shown by the divergence marked by the * symbols.

The analytic amplitude, \(A_j(t)\ ,\) for the j-th channel from equation (3) approximates the envelope of the filtered EEG (Figure 5, B). The mean squared value, \(\underline{A}^2(t)\ ,\) gives an estimate of the power dissipated by the flow of dendritic current across the fixed extracellular tissue impedance. The 64 scalar values of analytic amplitude, \(A_j^2(t)\) form a vector at each time step, \({\mathbf A}^2(t)\ .\) The vector is normalized by dividing every value by the mean, \(\underline{A}^2(t)\ ,\) to fix a point in 64-space that represents the spatial AM pattern at each time point. Successive points give a trajectory in 64-space that is projected into 2-space for visualization. The rate of pattern evolution is proportional to the Euclidean distance, \(D_e(t)\ ,\) between successive time points: \[ D_e(t) = |{\mathbf A}^2(t)| - |{\mathbf A}^2(t-1)| \ .\] \(D_e(t)\) is always ≥0 and should not be confused with the derivative of mean amplitude. An example of \(D_e(t)\) is shown as the green curve in Figure 5, C. The comparison with \(1/R_e(t)\) shows close relation in peaks and plateaus, indicating that relatively stationary epochs with stable frequency tend to have relatively stable AM patterns. Comparison of \(\underline{A}(t)\) with \(D_e(t)\) (Figure 5, D) shows that increases in \(\underline{A}^2(t)\) follow stabilization of frequency and AM pattern. The red bars in D highlight the times of maximal \(\underline{A}(t)\) in epochs of duration > 3 cycles at the center frequency used to determine the length of the window for calculating \(R_e(t)\) by equation (7). The orange bar shows an example for duration \(\leq\)3 cycles. A visual CS was given at \(t = 0\) ms; the red bars near 200 ms and 400 ms show the locations of values of \(\underline{A}^2(t)\) specifying AM patterns that were classified with respect to the CS.

In summary, the HT gives analytic state variables that reveal cinematographic frames of cortical activity, with abrupt onset by phase resetting in the beta or gamma carrier wave \(SD_X(t)\ ,\) resynchronization \(R_e(t)\) of the new carrier wave \(\omega(t)\ ,\) emergence and stabilization \(D_e(t)\) of its spatial AM pattern \({\mathbf A}^2(t)\ ,\) then massive increase in power output \(\underline{A}^2(t)\ .\) A cluster of points in \(n\)-space represents a brain state; a trajectory represents a state transition. The vector, \({\mathbf A}^2(t)\ ,\) serves as both a state variable and an order parameter (Freeman and Vitiello, 2006).

References

• Barlow JS. (1993) The Electroencephalogram: Its Patterns and Origins. Cambridge MA: MIT Press.
• Freeman WJ. (1975) Mass Action in the Nervous System. Academic Press, New York. Reprinted 2004:
• Freeman WJ. (2004a) Origin, structure, and role of background EEG activity. Part 1. Analytic amplitude. Clin. Neurophysiol. 115: 2077-2088.
• Freeman WJ. (2004b) Origin, structure, and role of background EEG activity. Part 2. Analytic phase. Clin. Neurophysiol. 115: 2089-2107. http://repositories.cdlib.org/postprints/987
• Freeman WJ. (2005) Origin, structure, and role of background EEG activity. Part 3. Neural frame classification. Clin. Neurophysiol. 116 (5): 1118-1129. http://authors.elsevier.com/sd/article/S1388245705000064
• Freeman WJ. (2006) Origin, structure, and role of background EEG activity. Part 4. Neural frame simulation. Clin. Neurophysiol. 117(3): 572-589.
• Freeman WJ, Burke BC, Holmes MD. (2003) Aperiodic phase re-setting in scalp EEG of beta-gamma oscillations by state transitions at alpha-theta rates. Hum. Brain Mapp. 19: 248-272.
• Freeman WJ, Vitiello G [2006] Nonlinear brain dynamics as macroscopic manifestation of underlying many-body field dynamics. Physics of Life Reviews 3: 93-118. http://dx.doi.org/10.1016/j.plrev.2006.02.001, http://repositories.cdlib.org/postprints/1515
• Huang NE, Shen Z, Long SR, Wu MC, Shih HH, Zheng Q, Yen N-C, Tung CC, Liu HH. (1998) The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis. Proc. R. Soc. Lond. 454: 903-995.
• Le Van Quyen M, Foucher J, Lachaux J-P, Rodriguez E, Lutz A, Martinerie J, Varela F. (2001) Comparison of Hilbert transform and wavelet methods for the analysis of neuronal synchrony. J. Neurosci. Meth. 111: 83-98.
• Pikovsky A, Rosenblum M, Kurths J. (2001) Synchronization — A Universal Concept in Non-linear Sciences. Cambridge UK: Cambridge U.P.
• Quiroga RQ, Kraskov A, Kreuz T, Grassberger P. (2002) Performance of different synchronization measures in real data: A case study on electroencephalographic signals. Physical Rev E, 6504:U645-U6 58 - art. no. 041903.
• Schroeder M. (1991) Fractals, Chaos, Power Laws: Minutes from an Infinite Paradise. San Francisco: WH Freeman.

Internal references

• Valentino Braitenberg (2007) Brain. Scholarpedia, 2(11):2918.
• James Meiss (2007) Dynamical systems. Scholarpedia, 2(2):1629.
• Paul L. Nunez and Ramesh Srinivasan (2007) Electroencephalogram. Scholarpedia, 2(2):1348.
• Roger D. Traub (2006) Fast oscillations. Scholarpedia, 1(12):1764.
• William D. Penny and Karl J. Friston (2007) Functional imaging. Scholarpedia, 2(5):1478.
• Seiji Ogawa and Yul-Wan Sung (2007) Functional magnetic resonance imaging. Scholarpedia, 2(10):3105.
• Rob Schreiber (2007) MATLAB. Scholarpedia, 2(7):2929.
• Jeff Moehlis, Kresimir Josic, Eric T. Shea-Brown (2006) Periodic orbit. Scholarpedia, 1(7):1358.
• Philip Holmes and Eric T. Shea-Brown (2006) Stability. Scholarpedia, 1(10):1838.
• Arkady Pikovsky and Michael Rosenblum (2007) Synchronization. Scholarpedia, 2(12):1459.

External links

• Walter J. Freeman's website

See Also

1/f Noise, Electroencephalogram, Functional Imaging, Hilbert Transform, Local Field Potential, Magnetoencephalogram, Oscillators, Phase Model, Synchronization, Synchronization of Chaotic Oscillators

Sponsored by: Eugene M. Izhikevich, Editor-in-Chief of Scholarpedia, the peer-reviewed open-access encyclopedia

Reviewed by: Dr. Giuseppe Vitiello, Dipartimento di Matematica e Informatica and Istituto Nazionale di Fisica Nucleare, University of Salerno, Italy

Reviewed by: Dr. Yanqing Chen, The Neurosciences Institute, San Diego, California

Accepted on: 2007-01-15 18:13:27 GMT