- Email: [email protected]
Dynamic spectral analysis of event-related potentials
Electroencephalography and clinical Neurophysiology 108 (1998) 251–259
Dynamic spectral analysis of event-related potentials Dmitriy Melkonian a ,*, Evian Gordon a, Christopher Rennie b, Homayoun Bahramali a a
Department of Psychological Medicine, Cognitive Neuroscience Unit, Westmead Hospital and the University of Sydney, Westmead, NSW 2145, Australia b Department of Medical Physics, Westmead Hospital and the University of Sydney, Westmead, NSW 2145, Australia Accepted for publication: 27 August 1997
Abstract This paper presents a new method for the identification of individual event related potential (ERP) components in both frequency and time domains. Using the similar basis function (SBF) algorithm the method provides a time to frequency transform, representing a frequency domain equivalent of the component waveform. Notable features of the SBF algorithm are that it allows for unevenly spaced sampled functions in both the time and frequency domains, and estimates of spectral densities are obtained by numerical computation of finite Fourier integrals. Application of this method to ERP data from 20 normal subjects demonstrated a similar shape of component amplitude frequency characteristics for traditional late component waveforms (N1, P2, N2 and P3). On this basis, a low-frequency band was found where the component amplitude frequency characteristic was described by a Gaussian function, while the component phase frequency characteristic was a linear function of frequency. These relationships are interpreted as frequency domain equivalents of the component. Transformed to the time domain, they provided an analytical description of the ERP as the sum of positive- and negative-going monopolar waves. The study points to similar mechanisms underlying these component waveforms, and analytically defines dynamic properties for the components both in the frequency and time domains. 1998 Elsevier Science Ireland Ltd. Keywords: Event-related potential; Dynamic spectral analysis; Fourier analysis
1. Introduction It is now well accepted that both time and frequency domain analysis of component structure of scalp potentials provides useful neurophysiological, neuropsychological and clinical information. The most common example is the representation of EEG in terms of different frequency bands. Most clinical studies apply spectral analysis to relatively long sections of activity, which obscures (due to the non-stationary character of EEG) rapid dynamic changes in the signal and can hide characteristic but transient spectral profiles. For these reasons dynamic spectral analysis (DSA) has become increasingly popular for various kinds of brain potential investigations (De Weerd and Kap, 1981; Norcia et al., 1986; Florian and Pfurtscheller, 1995). By calculating the spectrum over a succession of short time intervals, DSA not only allows the investigation of time courses of frequency spectra and their short-term changes, but also pro* Corresponding author. Tel.: +61 2 98456835; fax: +61 2 96357734; e-mail: [email protected]
0168-5597/98/$19.00 1998 Elsevier Science Ireland Ltd. All rights reserved PII
vides information on time courses of peak frequency and bandwidth of the dominant rhythms. Despite differences in methodology and computation, the general problems of DSA may be reduced to two basic issues: (i) the technique of local spectral decomposition of a signal (short time spectral analysis), and (ii) the choice of the short time intervals in which the signal characteristics may be considered as stationary. With respect to late component event-related EEG signals, the importance of stationarity of the data for DSA has recently been emphasized by Florian and Pfurtscheller (1995). They suggested a parametric approach that fits a sequence of autoregressive models to segments of EEG, all of the same length. Stationarity had to be assumed within each segment, although parameter variations between segments could be accommodated. With respect to the traditional late component auditory event related potential (ERP) waveform, our group have previously shown that stationary parts of the evoked potential correspond to individual ERP components (Haig et al., 1997). Consequently, to ensure the stationarity of data for DSA, the segments of data need to coincide with each com-
EEP 96182
252
D. Melkonian et al. / Electroencephalography and clinical Neurophysiology 108 (1998) 251–259
ponent of the waveform (i.e., the duration of each segment should be chosen according to the form of the ERP). Conventional tools of DSA, autoregressive models and spectral analysis based on the fast Fourier transformation (FFT) are not sufficiently adaptable in this respect. The goal of the present study was to develop and validate a new method of short-time spectral analysis of ERPs. It addresses the issue of stationarity by aligning the intervals for spectral analysis to the ERP component waveforms. With regard to the technique of spectral decomposition, the algorithm described below provides a novel means for an explicit frequency domain analysis of separate ERP component waveforms.
2. Methods In this section we first describe an efficient algorithm for evaluation of the Fourier transform of a discretely sampled time series, then we discuss how to deal with practical issues arising from application of the algorithm to ERP components, including offsets in time and voltage, and irregular sampling. Finally the application of the algorithm to inverse Fourier transforms is discussed. 2.1. Review of the similar basis function (SBF) algorithm Different problems of EEG and ERP frequency analysis (including transforms from time to frequency domain and inverse transforms from the frequency to the time domain) may be reduced to numerical computation of Fourier sine and/or cosine finite integrals:
R(u, L) =
L
0
y(x)cosuxdx, I(u, L) =
L
y(x)sinuxdx
(1)
0
where y(x) is an empirical function represented by its sampled values yn = y(xn) (n = 0,1,...,N), xn [ [0,L]. The FFT is strictly applicable only to periodic functions that are sampled with 2n evenly-spaced samples. The usual remedies of zero-padding and windowing introduce problems of their own, particularly when applied to transient, aperiodic segments of signals, and strongly limit the application of FFT to DSA of evoked brain potentials. Fourier integrals Eq. (1) are free of these complications. A widely accepted approximate method of carrying out these integrals for an aperiodic function is based on the use of a piecewise-linear approximation (Bracewell, 1965). This approach has previously been employed for spectral decomposition of a number of characteristic biopotentials, including electroretinograms and surface-negative potentials of the cerebral cortex (Melkonian, 1987; Roitbak et al., 1987). A general principle of piecewise-linear approximation is illustrated in Fig. 1A where the dotted curve represents fragment of y(x) and y0, y1, y2, y3 and y4 are the values of y(x) at points x0, x1, x2, x3 and x4. The approximating function h(x) is plotted as straight solid lines connecting
the samples y0, y1, y2, y3 and y4. Though the integrals Eq. (1) for the piecewise-linear function have well known analytical solution, its direct application requires extensive computation. Hence, for computational purposes it is highly desirable to make use of suitable specialized algorithms. The SBF algorithm (Melkonian, 1987) provides an effective means of evaluating the Fourier integrals of h(x). It involves representation of h(x) by the sum of similar basis functions of the following form N −1
h(x) = ∑ an r(x=xn + 1 )
(2)
n=0
where an are coefficients and r(x) is the triangular basis function (TBF) in Fig. 2A: r(x) = 1 − x for x [ [0, 1] and r(x) = 0 for x Ó [0, 1]. Given that x0 = 0, xn + 1 . xn and yN = 0, the coefficients an are determined by the condition hn = yn for n = 0, :::, N − 1
(3)
where hn = h(xn). A general procedure of the evaluation of the interpolation coefficients consists of substituting Eq. (3) into Eq. (2) for n = 0,...,N − 1 and finding the solution of the resulting system of N linear equations by conventional matrix methods. While this set of linear equations is easily solved by standard methods, there is a further simplification possible due to the fact that r(xm/xn) = 0 for m ≥ n. The geometric principle found useful in this connection may be seen by reference to Fig. 1. Fig. 1A portrays that approximating function may be represented in the interval [x0, x4] as a sum of triangular elements of two types: (i) the element E0 which may be represented in terms of TBS as y0r(x/x1), (ii) the elements E1, E2, E3 and E4 (the expansion may be continued in the same fashion for successive positive or negative values y5,...,yN − 1) which are the triangles of the same construction as the triangle abc in Fig. 1B. Given 0a = xn − 1, 0d = xn, 0c = xn + 1, bd = 1, ai = ae and 0g = 0f − 0h, the triangle abc may be considered as the sum of triangles 0hc, 0fd and 0ga, i.e., for x [ [xn − 1, xn + 1] as the function ln (x) = an r(x=xn + 1 ) − bn r(x=xn ) + gn r(x=xn − 1 ) where an = 0h = xn + 1 =Dxn + 1 , bn = 0f = xn
Dxn + 1 + Dxn Dxn + 1 Dxn
gn = 0g = xn − 1 =Dxn , Dxn = xn − xn − 1 In these terms, N −1
h(x) = y0 r(x=x1 ) + y1 [a1 r(x=x2 ) − b1 r(x=x1 )] + ∑ yn ln (x) n=2
Comparison with Eq. (2) shows that interpolation coefficients are readily found to be an = an yn − bn + 1 yn + 1 + gn + 2 yn + 2 (n = 0, :::, N − 1)
(4)
where a0 = 1 and it is assumed that yN = 0 and yN + 1 = 0.
D. Melkonian et al. / Electroencephalography and clinical Neurophysiology 108 (1998) 251–259
The system of simultaneous linear equations has thus been reduced to a set of independent linear equations. Now that h(x) is described in terms of a linear combination of TBFs, we turn now to the matter of finding the Fourier transform of h(x). Fourier cosine- and sine-integrals of TBF (Fig. 2B) are as follows:
R C (u) =
∞
r(x)cosuxdx =
0
R S (u) =
∞ 0
1 − cosu u2
u − sinu r(x)sinuxdx = u2
(5a)
(5b)
Given xn ≤ L, these formulas are valid for all TBFs in Eq. (2). Therefore, on the basis of the scaling theorem: N −1
R(u, L) ≈ ∑ an xn + 1 R C (uxn + 1 ), n=0
(6)
N −1
I(u, L) ≈ ∑ an xn + 1 R S (uxn + 1 ) n=0
2.2. Component Fourier transform (CFT) The problem of a dynamic spectral analysis can be considered as computation of Fourier exponential finite transform V (jq, t, t) =
t+t
n(y) e − jqy dy
(7)
t
where v(t) is the ERP voltage, [t, t + t] is an arbitrary time interval, q = 2pf, f is frequency and j = Î − 1. Let t = tC and t = tC, where tC is the time of beginning of component C (C may stand for N1, P2, N2 or P3 in the case of auditory ERPs) and tC is the duration of component C. Then, we may define the Fourier transform of the component C as
WC (jq) =
tC 0
Approximation of wC(t) is performed in the interval [0, tC] which implies that the approximating function h(tC) = 0. If wC(tC) = 0, then the interpolation condition at this point is satisfied and the SBF algorithm is immediately applicable to computing of Eq. (8). Otherwise, the intermediate function gC(t) = wC(t) − wC(tC), is introduced. Substitution of wC(t) = gC(t) + wC(tC) into Eq. (8) provides:
WC (jq) =
tC
0
gC (t) e − jqt dt +
wC (tC ) (1 − e − jqtC ) jq
The SBF algorithm may now be used for computing the Fourier transform of gC(t). Component frequency characteristic is typically calculated for a logarithmic frequency scale, i.e., for discrete frequencies qk = q0ck (k = 0,1,...,K) where q0 is the lowest value of angular frequency that is of interest, and c . 1 is a constant which defines the sampling interval in the logarithmic frequency scale. 2.3. Irregular sampling
These formulae are weighted sums of relatively simple non-periodic functions Eqs. (5a) and (5b). In summary, the SBF algorithm consists of two stages: derivation of interpolation coefficients from the discretely sampled waveform (Eq. (4)), and evaluation of the Fourier transform at any frequency u (Eq. (6)). Note that the same method may be applied to obtain a time function from a discretely sampled frequency domain signal.
253
wC (t) e − jqt dt
(8)
where wC(t) = v(t + tC) and WC(jq) is a component complex frequency characteristic. WC(jq) can be represented in terms of real functions. We define WC (jq) = R C (q) − jIC (q) = AC (q)jJC (q) where R C(q) is the component real frequency characteristic (CRFC), IC(q) is the component imaginary frequency characteristic (CIFC), AC(q) is the component amplitude frequency characteristic (CAFC) and JC(q) is the component phase frequency characteristic (CPFC).
In the most common situation the scalp potential, v(t), is sampled at evenly spaced intervals in time. Therefore, w(t) may be considered as being initially defined in [0, tC] by its evenly sampled values w(iD) (i = 0,...,I = tC /D), where D denotes the sampling interval. The choice of D is usually decided by the highest possible frequency components of the scalp potential. Consequently evenly sampled data, as a rule, contain a large volume of redundant information. Since the SBF algorithm allows both evenly and unevenly sampled data, we have an additional stage of processing, where the approximating function h(x) is determined by an automatic procedure which carries out data redundancy reduction using a first order linear interpolation algorithm (Cacers, 1965). In this case, each point xn in Eq. (2) is now considered to be a non-redundant point which may be represented as xn = MnD, where Mn is integer (M0 = 0, MN ≤ N). For these non-redundant points Eq. (3) is valid as before. Suppose that Mn + 1 − Mn = mn, where mn ≥ 1. For all m ≤ mn, i.e., for all redundant points, the procedure simply requires the condition |wC[(Mn + m)D] − h[(Mn + m)D]| , e, where e is a positive constant. Accordingly, e represents a parameter which defines the accuracy of approximation. 2.4. Inverse Fourier transform (IFT) Both R C(q) and IC(q) contain the information that makes it possible to restore wC(t) via the inverse Fourier cosine or sine transforms. We shall consider only the latter on the assumption that IC(q) is computed at points qn = q0cn (n = 0,...,N − 1, c . 1) and that the two following conditions are satisfied: (i) a piecewise- linear approximation function HC(q), which coincides with IC(q) at points qn, is a sufficiently accurate representation of IC(q) for q [ [q0, qN]; (ii) IC(q) ≈ 0 for q ≥ qN. Then the IFT is defined as
254
zC (t) =
D. Melkonian et al. / Electroencephalography and clinical Neurophysiology 108 (1998) 251–259
2 p
∞ 0
IC (q)sinqtdq ≈
2 p
qN
0
HC (q)sinqtdq
(9)
where zC(t) ≈ wC(t) for t [ [0, tC] and zC(t) ≈ 0 for t . tC. The SBF algorithm is equally applicable to numerical computation of Eq. (9) for discrete tn . 0. For the first section of the algorithm, finding interpolation coefficients, it is useful to note that due to the fact that qn follows a geometric progression, an, bn and gn in Eq. (4) do not depend on ‘n’; a = c/(c − 1), b = (c + 1)/(c − 1) and g = 1/(c − 1).Therefore, the interpolation coefficients are: an = [cIC (qn ) + (c + 1)IC (qn + 1 ) + IC (qn + 2 )]=(c − 1) (n = 0, :::, N − 1) For the second section, an effective computational scheme is provided by the choice of the following points in the time domain, ti = t0ci (t0 . 0, i = 0,1,...,I). From Eqs. (9) and (6): 2 zC (ti ) ≈ p
N −1
n=0
an qn + 1 IC (Qci + n + 1 )
(10)
where Q = q0t0. The reduction in computations provided by this formula emerges from the following rationale. To compute zC(ti) (single point) it is necessary to compute IC(q) at N points. Consequently, the general case of computation of zC(t) at I + 1 points requires evaluation of IC(q) at N(I + 1) points. However, Eq. (10) provides the conditions under which it is necessary to compute only N + I different values of IC(q).
3. Results The method described above was applied to data from 20 normal subjects, consisting of 12 males (age between 20 and 69 years) and 8 females (age between 21 and 55 years). The ERP data were averaged waveforms to 40 target stimuli in a conventional auditory oddball paradigm (Haig et al., 1995), recorded from the midline (Fz, Cz and Pz sites) with linked ears as reference. Waveforms from each electrode were analyzed in the same fashion. Each record consisted of 250 points, with 200 ms baseline prior to 800 ms data recorded after the stimulus. To identify the late components (N1, P2, N2 and P3), their amplitudes and latencies were first estimated automatically, using a method that determined the change in gradient of each component within established latency windows (Haig et al., 1995).To estimate tC and tC parameters, the tC and tC + tC instants are considered as the beginning and the end of the identified monopolar waveform (positive or negative with respect to the component under analysis) within windows, which we determined empirically from 50 normal subjects: N1 - from 40 to 180 ms; P2 - from 128 to 220 ms; N2 - from 148 to 304 ms; P3 from 208 to 488 ms. The processing techniques were first tested on a number
of typical ERPs by carrying out data redundancy reduction followed by forward then inverse transforms. Comparison of the results with the original waveforms led to the following choice of parameters: (i) for redundancy reduction with negligible loss of signal e = 0.7 (mV) (by this means the size of the records was reduced from 250 samples per second to an average of 133 samples per second); (ii) for the inverse Fourier transform a frequency band from 0.05 to 200 Hz, and c = 1.023293 (100 samples per decade) were chosen. The necessity of such a large frequency band was dictated by the required accuracy of IFT computations. If the computation of IFT is not necessary, then the frequency characteristics can be computed for a narrower frequency band. To illustrate the use of DSA, an ERP curve was chosen (Fig. 3A) whose N1, P2, N2 and P3 component waveforms were unambiguously identified by crossings with the time axis (arrows a, b, c, d and e). For the N1 component Fig. 3B and C illustrates CAFC and CIFC, respectively (the beginning and end of the epoch for analysis are marked by arrows a and b). Fig. 4A illustrates the generally accepted procedure of DSA using a succession of intervals of increasing duration running through the response. In accord with Eq. (7) the amplitude frequency characteristic V(q,0,tn) = Mod[V(jq,0,tn)] was computed for tn = nDt, where Dt = 8 ms and n = 1,...,50. A general pattern of V(q,0,tn) at fixed time, tn, is typical for a frequency domain representation of the brain evoked potentials (Stampfer and Basar, 1985). In particular, characteristic multiple resonant peaks can be recognized at t . 150 ms. Fig. 4B shows application of CFT to estimate the spectral content of the four main late components, N1, P2, N2 and P3. A striking property of component transform is that the resulting CAFC waveforms reveal a stereotyped shape (for all the ERPs analyzed), illustrated in Fig. 3B (N1 component) and Fig. 4B. It appears that CAFC has a relatively simple form, similar to the frequency response of a low-pass filter. After the initial plateau, CAFC monotonically decreases with increasing frequency, until a boundary frequency qB = 2pfB. The following parameters were chosen to quantify the CAFC: (i) amplitude coefficient, kA = AC(0) which represents geometrically the area under the component waveform; (ii) a halfpower frequency, qH = 2pfH, where the gain (relation of output power to input power) was 0.5 (−3 dB), i.e., AC(qH)/AC(0) = 0.707; (iii) a boundary frequency, fB. Evaluation of kA (mV⋅s) and fH (Hz) for different components showed the following means and standard errors (SE). For kA: N1 - 0.450 (0.034); P2 - 0.202 (0.022); N2 - 0.414 (0.058); P3 - 1.677 (0.161). For fH: N1 - 8.58 (0.34); P2 13.00 (0.79); N2 - 9.51 (0.70); P3 - 4.27 (0.28). These data illustrate the variance in CAFC parameters among separate groups of components and statistically significant intergroup differences between mean parameter values. In particular, it may be seen that the frequency composition of P3 could be attributed to a substantially lower band of frequencies than for the other component waveforms. However, it is important to provide a measure of similar-
D. Melkonian et al. / Electroencephalography and clinical Neurophysiology 108 (1998) 251–259
Fig. 1. (A) How the piecewise-linear approximating function (solid curve) may be represented as a sum of triangular elements. (B) The triangle abc can be regarded as the sum of triangles 0hc, 0fd and 0ga.
ity (or dissimilarity) of CAFC shapes. For this purpose a relative coefficient called the extension ratio has been defined as r = fB/fH. For different components, the means and SE for r are as follows: N1 - 2.80 (0.14); P2 - 2.57 (0.08); N2 - 2.78 (0.09); P3 - 2.94 (0.11).The absence of statistically significant differences between these groups points to a similar frequency composition of CAFC for all component waveforms. Considering all components together (n = 168),
Fig. 2. (A) The basis function. (B) Cosine (curve 1) and sine (curve 2) Fourier transformations of the basis function.
255
Fig. 3. (A) The solid curve plots an original ensemble averaged ERP over the time interval from 0 (stimulus presentation) to 400 ms post-stimulus. Dash-dot, dotted, dash-dot-dot and broken curves show functions zTC(t − tC) computed for N1, P2, N2 and P3, respectively. (B) The solid curve plots the function AN1(q)/AN1(0), i.e., the CAFC for N1 normalized by the amplitude. The parameters are as follows: AN1(0) = 0.594 mV⋅s, fH = 7.55 Hz; fB = 18.20 Hz. The dotted curve is a Gaussian function. (C) The solid curve plots the CIFC, IN1(q), while the dotted curve represents the approximation suggested by Eq. (14). The abscissae in B and C are in frequency f = q/2p.
the estimates for r were as follows: mean = 2.78, SE = 0.06, maximum value = 7.75, minimum value = 1.69. For the same group, estimation of AC(fB) showed a mean value = 0.049, SE = 0.003, maximum value = 0.190 and minimum value = 0.001. Therefore in the frequency band [0, fB] the CAFC is highly attenuated (−26 dB on average). The point that emerges is that CFT causes a filtering effect on the ERP component waveform, dividing it into the two bands: a low-frequency band at f , fB and a highfrequency band at f . fB. Further evidence of this is provided by best fit calculations using different approximating functions to CAFC. It appeared that independently of the type of the component, the CAFC is reliably described by Gaussian function (dotted curve in Fig. 3B) in the frequency range from zero to approximately 1.8fH. This finding is emphasized by the procedure of normalization using the kA and qH as the basis units. A normalized CAFC is defined as: AC*(n) = AC(q/qH)/kA, where n = q/qH is the relative
256
D. Melkonian et al. / Electroencephalography and clinical Neurophysiology 108 (1998) 251–259
frequency. Fig. 5A shows that after normalization the four CAFCs from Fig. 4B practically coincide in the range of relative frequencies from 0 to 1.8. Over this range it is difficult to distinguish in the figure between CAFCs and the dotted line which represents (with some constant multiplier) the Gaussian function, Q(u) = exp( − qu2 )
(11)
where q = 0.347. This value of q was derived from Q(1) = 0.707. It is important to note that the error of approximation depends on the relative frequency, n. To highlight this point consider the absolute value of deviation between the curves, e(n) = lAC*(n) − Q(n)l. Let eM1 and eM2 denote the maximum values of e(n) for the bands of relative frequencies [0, 1] and [1, 1.8], respectively. There is a high degree of consistency between the curves in the band of relative
frequencies [0, 1]: the mean value of eM1 evaluated for the whole group of 168 components was 0.00426 (0.426%) with standard deviation (SD) equal to 0.001. The corresponding estimates for the band of relative frequencies [1, 1.8] showed a very good degree of agreement: mean = 0.0483 (4.83%) and SD = 0.016. In the range of frequencies n . 1.8 a more complicated dependency is evident. Analysis of phase characteristics shows that increase of the relative frequency from n = 0 to about 2.0 is accompanied by a nearly proportional increase of the normalized CPFC, JC*(n) = JC(q/qH), from 0 to about p. For such a range the solid curves 1 and 2 in Fig. 5B represent JN1*(n) and JP3*(n), respectively. The curves only slightly deviate from linearity (to emphasize a linear character of dependencies the relative frequency is shown in linear abscissa scale). Accordingly, we considered a linear approximating function to JC*(n), LpC (n) = kJ n
Fig. 4. (A) Time course of amplitude frequency characteristic (AFC) V(q,0, tn) = Mod[V(jq,0, tn)] shows gradual changes with time (time scale corresponds to t) in the amplitudes of the ERP frequency components and reveals multiple resonant peaks that are difficult to interpret. (B) Using the same data as in A, CAFCs for each component (N1, P2, N2 and P3) show the frequency composition of the ERP component. The qualitative similarity is evident in the frequency domain images of the components. Quantitative differences in component waveform shape were found to be reflected by the fH parameter. Frequency scales in A and B are for f.
(12)
where kJ is a coefficient representing slope. The kJ coefficients have been evaluated for all component waveforms using the standard procedure for regressions through the origin (the values JC*(nj) from j = 0 to j = 157, where nj = n0cj, n0 = 0.05 and c = 1.023293, were used as initial data). No statistically significant differences have been found for the slopes among the ERP components (N1, P2, N2 and P3). For the group of component waveforms (n = 168), the mean for kJ is 1.86 (SD = 0.235). Linear dependencies (dash-dot and dotted curves) in Fig. 5B are characteristic examples of approximation. Deviation between the CPFC and approximating function typically increases with the increase of n. In the interval [0, 1.8] the mean of the absolute value of maximum deviation of CPFC from linear dependency for the whole group of component waveforms (n = 168) was 0.0916 (3.1% in relationship to the corresponding value of CPFC). We therefore found reasonable agreement for both normalized CAFCs and normalized CPFCs (ERP components N1, P2, N2, P3) with analytical dependencies Eqs. (11) and (12) in the relative frequency range from n = 0 to n about 1.8. To provide an interpretation of this result in terms of a classical additive model of evoked potentials (Lopes da Silva and Mars, 1987; Regan, 1989), let us consider wC(t) as the sum of two separate sources of activity. For t [ [0, tC] wC(t) = wTC(t) + wRC(t), where wTC(t) is a transient term linked in time and polarity to the stimulus, and wRC(t) is a random term (a statistically independent noise component usually associated with background brain activity). We can assume that relationships Eqs. (11) and (12) are relevant just to WTC(jq) = ATC(q)exp[−jJTC(q)] and consequently our study suggests that ATC (q) = kA exp( − aq2 ), JTC (q) = bq where a =
q/q2H
(13)
and b = kJ/qH. It follows from this that
ITC (q) = kA exp( − aq2 )sin(bq)
(14)
D. Melkonian et al. / Electroencephalography and clinical Neurophysiology 108 (1998) 251–259
and Fig. 3C shows that Eq. (14) provides reasonable approximation to IN1(q) in the frequency range from 0 to 1.8qH. Taking ITC(q) as a reliable frequency domain equivalent of the time domain data, we may expect that its inverse transform, zTC(t) defined by Eq. (9), approximates wTC(t) in a satisfactory manner. The existence of an analytical solution of integral Eq. (9) for function Eq. (14) (Erdelyi, 1954) provides the following formula for the transient term of a particular ERP component; for t ≥ 0 + − (t) − zTC (t) = kA {exp[ − (t − b2 =4a] zTC (t) = zTC
p
− exp[ − (t + b)2 =4a]}(2 pa) − 1
(15) The equation is presented in the form of a d’Alembert’s solution which is characteristic of wave-like phenomena (Richards and Williams, 1972). According to this representation each component waveform may be considered as a + (t) and negativesuperposition of positive-going wave zTC − going wave zTC(t). Fig. 3A illustrates that Eq. (15) provides reasonable agreement with experimental curves. The dashdot, dotted, dash-dot-dot and broken curves, represent functions zTC(t − tC) calculated from Eq. (15) for N1, P2, N2 and P3 components, respectively. The parameters kA (mV⋅s), a (s2), b (s) and tC (s) are as follows: N1 - 0.594, 0.000154, 0.0376, 0.064; P2 - 0.291, 0.000036, 0.0181, 0.14; N2 0.806, 0.000103, 0.0330, 0.184 and P3 - 1.785, 0.000227, 0.0490, 0.252.
4. Discussion The main findings reported in this paper depended critically on the methodological innovation of component Fourier transform, which was designed to provide a means for frequency domain representation of specific ERP component waveforms. The limitations of FFT for short segment spectral analysis have led previous authors to use non-traditional methods of spectral analysis, autoregressive models (Lopes da Silva and Mars, 1987; Florian and Pfurtscheller, 1995) and the wavelet transform (Bertrand et al., 1994; Schiff et al., 1994). The SBF algorithm provides a tool for DSA which does not lose the connection with the classical theory of Fourier analysis. To highlight the differences between the FFT and SBF algorithms, the FFT is based on a classical theory of Fourier series which deals with expansions of periodic functions into a sum of sines and cosines for discrete frequencies that are multiples of the fundamental harmonic. The SBF algorithm on the other hand, is based on a theory of Fourier integrals which deals with non-periodic functions. Therefore, the notion of discrete harmonics is superseded by a frequency domain representation of a signal which can be dealt with as a continuous frequency function (spectral density) defined
257
for a semi-infinite interval of frequencies from 0 to ∞. Processing empirical data of finite length disguises these differences because the computation deals with discrete values of the time and frequency domain characteristics for physically meaningful time intervals and bands of frequencies. However for the purposes of this study the SBF algorithm was of primary importance for: (i) effective DSA from evenly or unevenly sampled data on segments of arbitrary length, (ii) computation of frequency characteristics for a logarithmic scale of frequencies, and (iii) transformation from the frequency to the time domain (to control for distortions of a short-time spectral representation of the ERP waveforms). Refinements in DSA revealed a clear change in spectrotemporal patterns of ERP data provided by the CFT (illustrated in Fig. 4). In light of a general approach to the treatment of spectro-temporal structure of the brain evoked potentials (De Weerd and Kap, 1981), the deep qualitative changes in temporal evolution of the frequency domain characteristics (Fig. 4A) may be considered as an indication of the inherently non-stationary character of the ERP data. In contrast, a simple and regular form of CAFCs (Fig. 4B) provides evidence for the property of local stationarity of ERP component waveforms. The most interesting property of CFT is that it causes a specific ‘filtering effect’ on a component waveform. By this we do not mean the existence of certain fixed bands of ‘natural’ frequencies. Rather, the study shows that the half-power frequency, fH, strongly depends on the individual characteristics of the component waveform and has different estimates for groups of N1, P2, N2 and P3 components. However, this work establishes a specific form of the frequency domain patterns of CAFC and CPFC depending on the relative frequency n = f/fH. It is striking that the Gaussian function Eq. (11) and the linear function Eq. (12) provide a practically exact description of CAFC and CPFC, respectively, for the relative frequency range from 0 to approximately n = 1.8. What is of particular interest for this finding is that frequency characteristics in this band are almost free of noise contributions which are inevitable with ERPs. It may be concluded that in the relative frequency domain the frequency counterparts of noise represent components with relative frequencies higher than 1.8. Therefore, applied to the classical additive model of evoked potentials, the point that we do consider to have established, is that in the range of relative frequencies from 0 to 1.8, the CAFC and CPFC can be considered as the frequency domain description of the transient term, wTC(t).The consistent relationship of these two major characteristics provided an explicit analytical description Eq. (15) for the transient term of the component waveform in the time domain. There is a satisfactory degree of agreement of time characteristics (see typical example in Fig. 3A), since the equation parameters kA, a and b were derived entirely from the frequency domain analysis, without any additional adjustment to make them fit the original component waveforms.
258
D. Melkonian et al. / Electroencephalography and clinical Neurophysiology 108 (1998) 251–259
they tend to produce a Gaussian distribution even if the source itself is not normally distributed. Therefore, the empirically grounded appearance of the Gaussian function in the wave-like description of an ERP component, may be relevant to the ‘mass’ effect of ERP multiple sources. The CFT technique, in addition to providing reliable estimation of the frequency content of the ERP components, provides the investigator with dynamic model of the component both in the frequency and time domains. The model provides a unique set of parameters, which can be considered as the simplest frequency domain description of ERP components, and is at the same time able to reproduce the main features of transient ERP behavior. The model therefore provides a compressed representation of the ERP component wave form. The four parameters, kA, a, b and tC represent dynamic properties of ERP component waveforms, complementing conventional time and amplitude ERP component measures. Moreover, the physiological and dynamic dimensions of specific ERP components provided by our model, may help to elucidate the nature of specific ERP components and determine their interrelationship with background EEG activity in health and disease states.
Fig. 5. (A) Solid lines show the superposition of the four CAFCs from B after their normalization in both amplitude and frequency. The dotted line shows the Gaussian function described by Eq. (11). (B) Solid curves 1 and 2 plot CPFCs for N1 and P3 components, respectively. The dash-dot and broken curves show the linear approximating Eq. (12) with kJ = 1.78 and kJ = 1.92, respectively.
The time domain equation for ERP component waveforms Eq. (15) is in agreement with a previous study of our group (Haig et al., 1997) and with the study of Verleger and Wascher (1995) in that it contains the Gaussian function. However a single Gaussian function has a number of inconsistencies, including: (i) the function is non-zero at all times which is impossible for a stimulus induced process; (ii) the shape of the function is symmetric which is not always the case for ERP component waveforms. Eqn (15) is free of these limitations. Indeed, the function has a zero value at t = 0 and a non-symmetric shape with steep rise and slowly decaying tail. However, the number of free parameters remains the same. A number of authors consider brain potentials as wavelike phenomena (Freeman, 1975; Nunez, 1994). From this point of view it is interesting to note that the time domain description of the component waveform may be presented as a superposition of positive and negative going-waves. Considerable evidence exists that dynamics of scalp potentials depends upon the types of synaptic activity and spatiotemporal characteristics of multiple current sources and sinks in the soma-dendritic membrane of cortical pyramidal cells (Mitzdorf, 1985). From the fields of probability and statistics, it is well known that if the sources of the process are multiple, and if they have some degree of variability,
References Bertrand, O., Bohorquez, J. and Pernier, J. Time-frequency digital filtering based on an invertible wavelet transform: an application to evoked potentials. IEEE Trans. Biomed. Eng., 1994, 41: 77–88. Bracewell, R. The Fourier Transform and Its Applications. McGraw-Hill, New York, 1965. De Weerd, J.P.C. and Kap, J.I. Spectro-temporal representation and timevarying spectra of evoked potentials. Biol. Cybern., 1981, 41: 101– 117. Cacers, C.A. Biomedical Telemetry. Academic Press, New York, 1965. Erdelyi, A. (Ed.). Tables of Integral Transforms, Vol. 1. McGraw-Hill, New York, 1954, p. 78. Florian, G. and Pfurtscheller, G. Dynamic spectral analysis of eventrelated data. Electroenceph. clin. Neurophysiol., 1995, 95: 393–396. Freeman, W.J. Mass Action in the Nervous System. Academic Press, New York, 1975. Haig, A.R., Gordon, E., Rogers, G. and Anderson, J. Classification of single-trial ERP sub-types: application of globally optimal vector quantization using simulated annealing. Electroenceph. clin. Neurophysiol., 1995, 94: 288–297. Haig, A.R., Rennie, C. and Gordon, E. The use of Gaussian component modelling to elucidate average ERP component overlap in schizophrenia. J. Psychophysiol., 1997, 11: 173–187. Lopes da Silva, F.H. and Mars, N.J. Parametric methods in EEG analysis. In: A.S. Gevins and A. Remond (Eds.), Methods of Analysis of Brain Electrical and Magnetic Signals. EEG Handbook (Revised Series). Elsevier, Amsterdam, 1987, pp. 243–260. Melkonian, D.S. Transients in Neuronal Systems (in Russian). Publ. House of Armenian Acad. of Sciences, Yerevan, 1987. Mitzdorf, U. Current source density method and application in cat cerebral cortex: 14 investigation of evoked potentials and EEG phenomena. Physiol. Rev., 1985, 65: 37–100. Norcia, A.M., Sato, T., Shinn, P. and Mertus, J. Methods for the identification of evoked response components in the frequency and combined time/frequency domains. Electroenceph. clin. Neurophysiol., 1986, 65: 212–226.
D. Melkonian et al. / Electroencephalography and clinical Neurophysiology 108 (1998) 251–259 Nunez, P.L. Neocortical Dynamics and Human EEG Rhythms. Oxford University Press, New York, 1994. Regan, D. Human Brain Electrophysiology: Evoked Potentials and Evoked Magnetic Fields in Science and Medicine. Elsevier, New York, 1989. Richards, J.P.G. and Williams., R.P. Waves. Penguin Books, Middlesex, 1972. Roitbak, A.I., Fanardjian, V.V., Melkonian, D.S. and Melkonian, A.A. Contribution of glia and neurons to the surface-negative potentials of the cerebral cortex during its electrical stimulation. Neuroscience, 1987, 20: 1057–1067.
259
Schiff, S.J., Aldroubi, A., Unser, M. and Sato, S. Fast wavelet transformation of EEG. Electroenceph. clin. Neurophysiol., 1994, 91: 442–455. Stampfer, H.G. and Basar, E. Does frequency analysis lead to better understanding of human event related potentials. Int. J. Neurosci., 1985, 26: 181–196. Verleger, R. and Wascher, E. Fitting ex-Gauss functions to P3 waveshapes: an attempt at distinguishing between real and apparent changes of P3 latency. J. Psychophysiol., 1995, 9: 146–158.
Dynamic spectral analysis of event-related potentials Dmitriy Melkonian a ,*, Evian Gordon a, Christopher Rennie b, Homayoun Bahramali a a
Department of Psychological Medicine, Cognitive Neuroscience Unit, Westmead Hospital and the University of Sydney, Westmead, NSW 2145, Australia b Department of Medical Physics, Westmead Hospital and the University of Sydney, Westmead, NSW 2145, Australia Accepted for publication: 27 August 1997
Abstract This paper presents a new method for the identification of individual event related potential (ERP) components in both frequency and time domains. Using the similar basis function (SBF) algorithm the method provides a time to frequency transform, representing a frequency domain equivalent of the component waveform. Notable features of the SBF algorithm are that it allows for unevenly spaced sampled functions in both the time and frequency domains, and estimates of spectral densities are obtained by numerical computation of finite Fourier integrals. Application of this method to ERP data from 20 normal subjects demonstrated a similar shape of component amplitude frequency characteristics for traditional late component waveforms (N1, P2, N2 and P3). On this basis, a low-frequency band was found where the component amplitude frequency characteristic was described by a Gaussian function, while the component phase frequency characteristic was a linear function of frequency. These relationships are interpreted as frequency domain equivalents of the component. Transformed to the time domain, they provided an analytical description of the ERP as the sum of positive- and negative-going monopolar waves. The study points to similar mechanisms underlying these component waveforms, and analytically defines dynamic properties for the components both in the frequency and time domains. 1998 Elsevier Science Ireland Ltd. Keywords: Event-related potential; Dynamic spectral analysis; Fourier analysis
1. Introduction It is now well accepted that both time and frequency domain analysis of component structure of scalp potentials provides useful neurophysiological, neuropsychological and clinical information. The most common example is the representation of EEG in terms of different frequency bands. Most clinical studies apply spectral analysis to relatively long sections of activity, which obscures (due to the non-stationary character of EEG) rapid dynamic changes in the signal and can hide characteristic but transient spectral profiles. For these reasons dynamic spectral analysis (DSA) has become increasingly popular for various kinds of brain potential investigations (De Weerd and Kap, 1981; Norcia et al., 1986; Florian and Pfurtscheller, 1995). By calculating the spectrum over a succession of short time intervals, DSA not only allows the investigation of time courses of frequency spectra and their short-term changes, but also pro* Corresponding author. Tel.: +61 2 98456835; fax: +61 2 96357734; e-mail: [email protected]
0168-5597/98/$19.00 1998 Elsevier Science Ireland Ltd. All rights reserved PII
vides information on time courses of peak frequency and bandwidth of the dominant rhythms. Despite differences in methodology and computation, the general problems of DSA may be reduced to two basic issues: (i) the technique of local spectral decomposition of a signal (short time spectral analysis), and (ii) the choice of the short time intervals in which the signal characteristics may be considered as stationary. With respect to late component event-related EEG signals, the importance of stationarity of the data for DSA has recently been emphasized by Florian and Pfurtscheller (1995). They suggested a parametric approach that fits a sequence of autoregressive models to segments of EEG, all of the same length. Stationarity had to be assumed within each segment, although parameter variations between segments could be accommodated. With respect to the traditional late component auditory event related potential (ERP) waveform, our group have previously shown that stationary parts of the evoked potential correspond to individual ERP components (Haig et al., 1997). Consequently, to ensure the stationarity of data for DSA, the segments of data need to coincide with each com-
EEP 96182
252
D. Melkonian et al. / Electroencephalography and clinical Neurophysiology 108 (1998) 251–259
ponent of the waveform (i.e., the duration of each segment should be chosen according to the form of the ERP). Conventional tools of DSA, autoregressive models and spectral analysis based on the fast Fourier transformation (FFT) are not sufficiently adaptable in this respect. The goal of the present study was to develop and validate a new method of short-time spectral analysis of ERPs. It addresses the issue of stationarity by aligning the intervals for spectral analysis to the ERP component waveforms. With regard to the technique of spectral decomposition, the algorithm described below provides a novel means for an explicit frequency domain analysis of separate ERP component waveforms.
2. Methods In this section we first describe an efficient algorithm for evaluation of the Fourier transform of a discretely sampled time series, then we discuss how to deal with practical issues arising from application of the algorithm to ERP components, including offsets in time and voltage, and irregular sampling. Finally the application of the algorithm to inverse Fourier transforms is discussed. 2.1. Review of the similar basis function (SBF) algorithm Different problems of EEG and ERP frequency analysis (including transforms from time to frequency domain and inverse transforms from the frequency to the time domain) may be reduced to numerical computation of Fourier sine and/or cosine finite integrals:
R(u, L) =
L
0
y(x)cosuxdx, I(u, L) =
L
y(x)sinuxdx
(1)
0
where y(x) is an empirical function represented by its sampled values yn = y(xn) (n = 0,1,...,N), xn [ [0,L]. The FFT is strictly applicable only to periodic functions that are sampled with 2n evenly-spaced samples. The usual remedies of zero-padding and windowing introduce problems of their own, particularly when applied to transient, aperiodic segments of signals, and strongly limit the application of FFT to DSA of evoked brain potentials. Fourier integrals Eq. (1) are free of these complications. A widely accepted approximate method of carrying out these integrals for an aperiodic function is based on the use of a piecewise-linear approximation (Bracewell, 1965). This approach has previously been employed for spectral decomposition of a number of characteristic biopotentials, including electroretinograms and surface-negative potentials of the cerebral cortex (Melkonian, 1987; Roitbak et al., 1987). A general principle of piecewise-linear approximation is illustrated in Fig. 1A where the dotted curve represents fragment of y(x) and y0, y1, y2, y3 and y4 are the values of y(x) at points x0, x1, x2, x3 and x4. The approximating function h(x) is plotted as straight solid lines connecting
the samples y0, y1, y2, y3 and y4. Though the integrals Eq. (1) for the piecewise-linear function have well known analytical solution, its direct application requires extensive computation. Hence, for computational purposes it is highly desirable to make use of suitable specialized algorithms. The SBF algorithm (Melkonian, 1987) provides an effective means of evaluating the Fourier integrals of h(x). It involves representation of h(x) by the sum of similar basis functions of the following form N −1
h(x) = ∑ an r(x=xn + 1 )
(2)
n=0
where an are coefficients and r(x) is the triangular basis function (TBF) in Fig. 2A: r(x) = 1 − x for x [ [0, 1] and r(x) = 0 for x Ó [0, 1]. Given that x0 = 0, xn + 1 . xn and yN = 0, the coefficients an are determined by the condition hn = yn for n = 0, :::, N − 1
(3)
where hn = h(xn). A general procedure of the evaluation of the interpolation coefficients consists of substituting Eq. (3) into Eq. (2) for n = 0,...,N − 1 and finding the solution of the resulting system of N linear equations by conventional matrix methods. While this set of linear equations is easily solved by standard methods, there is a further simplification possible due to the fact that r(xm/xn) = 0 for m ≥ n. The geometric principle found useful in this connection may be seen by reference to Fig. 1. Fig. 1A portrays that approximating function may be represented in the interval [x0, x4] as a sum of triangular elements of two types: (i) the element E0 which may be represented in terms of TBS as y0r(x/x1), (ii) the elements E1, E2, E3 and E4 (the expansion may be continued in the same fashion for successive positive or negative values y5,...,yN − 1) which are the triangles of the same construction as the triangle abc in Fig. 1B. Given 0a = xn − 1, 0d = xn, 0c = xn + 1, bd = 1, ai = ae and 0g = 0f − 0h, the triangle abc may be considered as the sum of triangles 0hc, 0fd and 0ga, i.e., for x [ [xn − 1, xn + 1] as the function ln (x) = an r(x=xn + 1 ) − bn r(x=xn ) + gn r(x=xn − 1 ) where an = 0h = xn + 1 =Dxn + 1 , bn = 0f = xn
Dxn + 1 + Dxn Dxn + 1 Dxn
gn = 0g = xn − 1 =Dxn , Dxn = xn − xn − 1 In these terms, N −1
h(x) = y0 r(x=x1 ) + y1 [a1 r(x=x2 ) − b1 r(x=x1 )] + ∑ yn ln (x) n=2
Comparison with Eq. (2) shows that interpolation coefficients are readily found to be an = an yn − bn + 1 yn + 1 + gn + 2 yn + 2 (n = 0, :::, N − 1)
(4)
where a0 = 1 and it is assumed that yN = 0 and yN + 1 = 0.
D. Melkonian et al. / Electroencephalography and clinical Neurophysiology 108 (1998) 251–259
The system of simultaneous linear equations has thus been reduced to a set of independent linear equations. Now that h(x) is described in terms of a linear combination of TBFs, we turn now to the matter of finding the Fourier transform of h(x). Fourier cosine- and sine-integrals of TBF (Fig. 2B) are as follows:
R C (u) =
∞
r(x)cosuxdx =
0
R S (u) =
∞ 0
1 − cosu u2
u − sinu r(x)sinuxdx = u2
(5a)
(5b)
Given xn ≤ L, these formulas are valid for all TBFs in Eq. (2). Therefore, on the basis of the scaling theorem: N −1
R(u, L) ≈ ∑ an xn + 1 R C (uxn + 1 ), n=0
(6)
N −1
I(u, L) ≈ ∑ an xn + 1 R S (uxn + 1 ) n=0
2.2. Component Fourier transform (CFT) The problem of a dynamic spectral analysis can be considered as computation of Fourier exponential finite transform V (jq, t, t) =
t+t
n(y) e − jqy dy
(7)
t
where v(t) is the ERP voltage, [t, t + t] is an arbitrary time interval, q = 2pf, f is frequency and j = Î − 1. Let t = tC and t = tC, where tC is the time of beginning of component C (C may stand for N1, P2, N2 or P3 in the case of auditory ERPs) and tC is the duration of component C. Then, we may define the Fourier transform of the component C as
WC (jq) =
tC 0
Approximation of wC(t) is performed in the interval [0, tC] which implies that the approximating function h(tC) = 0. If wC(tC) = 0, then the interpolation condition at this point is satisfied and the SBF algorithm is immediately applicable to computing of Eq. (8). Otherwise, the intermediate function gC(t) = wC(t) − wC(tC), is introduced. Substitution of wC(t) = gC(t) + wC(tC) into Eq. (8) provides:
WC (jq) =
tC
0
gC (t) e − jqt dt +
wC (tC ) (1 − e − jqtC ) jq
The SBF algorithm may now be used for computing the Fourier transform of gC(t). Component frequency characteristic is typically calculated for a logarithmic frequency scale, i.e., for discrete frequencies qk = q0ck (k = 0,1,...,K) where q0 is the lowest value of angular frequency that is of interest, and c . 1 is a constant which defines the sampling interval in the logarithmic frequency scale. 2.3. Irregular sampling
These formulae are weighted sums of relatively simple non-periodic functions Eqs. (5a) and (5b). In summary, the SBF algorithm consists of two stages: derivation of interpolation coefficients from the discretely sampled waveform (Eq. (4)), and evaluation of the Fourier transform at any frequency u (Eq. (6)). Note that the same method may be applied to obtain a time function from a discretely sampled frequency domain signal.
253
wC (t) e − jqt dt
(8)
where wC(t) = v(t + tC) and WC(jq) is a component complex frequency characteristic. WC(jq) can be represented in terms of real functions. We define WC (jq) = R C (q) − jIC (q) = AC (q)jJC (q) where R C(q) is the component real frequency characteristic (CRFC), IC(q) is the component imaginary frequency characteristic (CIFC), AC(q) is the component amplitude frequency characteristic (CAFC) and JC(q) is the component phase frequency characteristic (CPFC).
In the most common situation the scalp potential, v(t), is sampled at evenly spaced intervals in time. Therefore, w(t) may be considered as being initially defined in [0, tC] by its evenly sampled values w(iD) (i = 0,...,I = tC /D), where D denotes the sampling interval. The choice of D is usually decided by the highest possible frequency components of the scalp potential. Consequently evenly sampled data, as a rule, contain a large volume of redundant information. Since the SBF algorithm allows both evenly and unevenly sampled data, we have an additional stage of processing, where the approximating function h(x) is determined by an automatic procedure which carries out data redundancy reduction using a first order linear interpolation algorithm (Cacers, 1965). In this case, each point xn in Eq. (2) is now considered to be a non-redundant point which may be represented as xn = MnD, where Mn is integer (M0 = 0, MN ≤ N). For these non-redundant points Eq. (3) is valid as before. Suppose that Mn + 1 − Mn = mn, where mn ≥ 1. For all m ≤ mn, i.e., for all redundant points, the procedure simply requires the condition |wC[(Mn + m)D] − h[(Mn + m)D]| , e, where e is a positive constant. Accordingly, e represents a parameter which defines the accuracy of approximation. 2.4. Inverse Fourier transform (IFT) Both R C(q) and IC(q) contain the information that makes it possible to restore wC(t) via the inverse Fourier cosine or sine transforms. We shall consider only the latter on the assumption that IC(q) is computed at points qn = q0cn (n = 0,...,N − 1, c . 1) and that the two following conditions are satisfied: (i) a piecewise- linear approximation function HC(q), which coincides with IC(q) at points qn, is a sufficiently accurate representation of IC(q) for q [ [q0, qN]; (ii) IC(q) ≈ 0 for q ≥ qN. Then the IFT is defined as
254
zC (t) =
D. Melkonian et al. / Electroencephalography and clinical Neurophysiology 108 (1998) 251–259
2 p
∞ 0
IC (q)sinqtdq ≈
2 p
qN
0
HC (q)sinqtdq
(9)
where zC(t) ≈ wC(t) for t [ [0, tC] and zC(t) ≈ 0 for t . tC. The SBF algorithm is equally applicable to numerical computation of Eq. (9) for discrete tn . 0. For the first section of the algorithm, finding interpolation coefficients, it is useful to note that due to the fact that qn follows a geometric progression, an, bn and gn in Eq. (4) do not depend on ‘n’; a = c/(c − 1), b = (c + 1)/(c − 1) and g = 1/(c − 1).Therefore, the interpolation coefficients are: an = [cIC (qn ) + (c + 1)IC (qn + 1 ) + IC (qn + 2 )]=(c − 1) (n = 0, :::, N − 1) For the second section, an effective computational scheme is provided by the choice of the following points in the time domain, ti = t0ci (t0 . 0, i = 0,1,...,I). From Eqs. (9) and (6): 2 zC (ti ) ≈ p
N −1
n=0
an qn + 1 IC (Qci + n + 1 )
(10)
where Q = q0t0. The reduction in computations provided by this formula emerges from the following rationale. To compute zC(ti) (single point) it is necessary to compute IC(q) at N points. Consequently, the general case of computation of zC(t) at I + 1 points requires evaluation of IC(q) at N(I + 1) points. However, Eq. (10) provides the conditions under which it is necessary to compute only N + I different values of IC(q).
3. Results The method described above was applied to data from 20 normal subjects, consisting of 12 males (age between 20 and 69 years) and 8 females (age between 21 and 55 years). The ERP data were averaged waveforms to 40 target stimuli in a conventional auditory oddball paradigm (Haig et al., 1995), recorded from the midline (Fz, Cz and Pz sites) with linked ears as reference. Waveforms from each electrode were analyzed in the same fashion. Each record consisted of 250 points, with 200 ms baseline prior to 800 ms data recorded after the stimulus. To identify the late components (N1, P2, N2 and P3), their amplitudes and latencies were first estimated automatically, using a method that determined the change in gradient of each component within established latency windows (Haig et al., 1995).To estimate tC and tC parameters, the tC and tC + tC instants are considered as the beginning and the end of the identified monopolar waveform (positive or negative with respect to the component under analysis) within windows, which we determined empirically from 50 normal subjects: N1 - from 40 to 180 ms; P2 - from 128 to 220 ms; N2 - from 148 to 304 ms; P3 from 208 to 488 ms. The processing techniques were first tested on a number
of typical ERPs by carrying out data redundancy reduction followed by forward then inverse transforms. Comparison of the results with the original waveforms led to the following choice of parameters: (i) for redundancy reduction with negligible loss of signal e = 0.7 (mV) (by this means the size of the records was reduced from 250 samples per second to an average of 133 samples per second); (ii) for the inverse Fourier transform a frequency band from 0.05 to 200 Hz, and c = 1.023293 (100 samples per decade) were chosen. The necessity of such a large frequency band was dictated by the required accuracy of IFT computations. If the computation of IFT is not necessary, then the frequency characteristics can be computed for a narrower frequency band. To illustrate the use of DSA, an ERP curve was chosen (Fig. 3A) whose N1, P2, N2 and P3 component waveforms were unambiguously identified by crossings with the time axis (arrows a, b, c, d and e). For the N1 component Fig. 3B and C illustrates CAFC and CIFC, respectively (the beginning and end of the epoch for analysis are marked by arrows a and b). Fig. 4A illustrates the generally accepted procedure of DSA using a succession of intervals of increasing duration running through the response. In accord with Eq. (7) the amplitude frequency characteristic V(q,0,tn) = Mod[V(jq,0,tn)] was computed for tn = nDt, where Dt = 8 ms and n = 1,...,50. A general pattern of V(q,0,tn) at fixed time, tn, is typical for a frequency domain representation of the brain evoked potentials (Stampfer and Basar, 1985). In particular, characteristic multiple resonant peaks can be recognized at t . 150 ms. Fig. 4B shows application of CFT to estimate the spectral content of the four main late components, N1, P2, N2 and P3. A striking property of component transform is that the resulting CAFC waveforms reveal a stereotyped shape (for all the ERPs analyzed), illustrated in Fig. 3B (N1 component) and Fig. 4B. It appears that CAFC has a relatively simple form, similar to the frequency response of a low-pass filter. After the initial plateau, CAFC monotonically decreases with increasing frequency, until a boundary frequency qB = 2pfB. The following parameters were chosen to quantify the CAFC: (i) amplitude coefficient, kA = AC(0) which represents geometrically the area under the component waveform; (ii) a halfpower frequency, qH = 2pfH, where the gain (relation of output power to input power) was 0.5 (−3 dB), i.e., AC(qH)/AC(0) = 0.707; (iii) a boundary frequency, fB. Evaluation of kA (mV⋅s) and fH (Hz) for different components showed the following means and standard errors (SE). For kA: N1 - 0.450 (0.034); P2 - 0.202 (0.022); N2 - 0.414 (0.058); P3 - 1.677 (0.161). For fH: N1 - 8.58 (0.34); P2 13.00 (0.79); N2 - 9.51 (0.70); P3 - 4.27 (0.28). These data illustrate the variance in CAFC parameters among separate groups of components and statistically significant intergroup differences between mean parameter values. In particular, it may be seen that the frequency composition of P3 could be attributed to a substantially lower band of frequencies than for the other component waveforms. However, it is important to provide a measure of similar-
D. Melkonian et al. / Electroencephalography and clinical Neurophysiology 108 (1998) 251–259
Fig. 1. (A) How the piecewise-linear approximating function (solid curve) may be represented as a sum of triangular elements. (B) The triangle abc can be regarded as the sum of triangles 0hc, 0fd and 0ga.
ity (or dissimilarity) of CAFC shapes. For this purpose a relative coefficient called the extension ratio has been defined as r = fB/fH. For different components, the means and SE for r are as follows: N1 - 2.80 (0.14); P2 - 2.57 (0.08); N2 - 2.78 (0.09); P3 - 2.94 (0.11).The absence of statistically significant differences between these groups points to a similar frequency composition of CAFC for all component waveforms. Considering all components together (n = 168),
Fig. 2. (A) The basis function. (B) Cosine (curve 1) and sine (curve 2) Fourier transformations of the basis function.
255
Fig. 3. (A) The solid curve plots an original ensemble averaged ERP over the time interval from 0 (stimulus presentation) to 400 ms post-stimulus. Dash-dot, dotted, dash-dot-dot and broken curves show functions zTC(t − tC) computed for N1, P2, N2 and P3, respectively. (B) The solid curve plots the function AN1(q)/AN1(0), i.e., the CAFC for N1 normalized by the amplitude. The parameters are as follows: AN1(0) = 0.594 mV⋅s, fH = 7.55 Hz; fB = 18.20 Hz. The dotted curve is a Gaussian function. (C) The solid curve plots the CIFC, IN1(q), while the dotted curve represents the approximation suggested by Eq. (14). The abscissae in B and C are in frequency f = q/2p.
the estimates for r were as follows: mean = 2.78, SE = 0.06, maximum value = 7.75, minimum value = 1.69. For the same group, estimation of AC(fB) showed a mean value = 0.049, SE = 0.003, maximum value = 0.190 and minimum value = 0.001. Therefore in the frequency band [0, fB] the CAFC is highly attenuated (−26 dB on average). The point that emerges is that CFT causes a filtering effect on the ERP component waveform, dividing it into the two bands: a low-frequency band at f , fB and a highfrequency band at f . fB. Further evidence of this is provided by best fit calculations using different approximating functions to CAFC. It appeared that independently of the type of the component, the CAFC is reliably described by Gaussian function (dotted curve in Fig. 3B) in the frequency range from zero to approximately 1.8fH. This finding is emphasized by the procedure of normalization using the kA and qH as the basis units. A normalized CAFC is defined as: AC*(n) = AC(q/qH)/kA, where n = q/qH is the relative
256
D. Melkonian et al. / Electroencephalography and clinical Neurophysiology 108 (1998) 251–259
frequency. Fig. 5A shows that after normalization the four CAFCs from Fig. 4B practically coincide in the range of relative frequencies from 0 to 1.8. Over this range it is difficult to distinguish in the figure between CAFCs and the dotted line which represents (with some constant multiplier) the Gaussian function, Q(u) = exp( − qu2 )
(11)
where q = 0.347. This value of q was derived from Q(1) = 0.707. It is important to note that the error of approximation depends on the relative frequency, n. To highlight this point consider the absolute value of deviation between the curves, e(n) = lAC*(n) − Q(n)l. Let eM1 and eM2 denote the maximum values of e(n) for the bands of relative frequencies [0, 1] and [1, 1.8], respectively. There is a high degree of consistency between the curves in the band of relative
frequencies [0, 1]: the mean value of eM1 evaluated for the whole group of 168 components was 0.00426 (0.426%) with standard deviation (SD) equal to 0.001. The corresponding estimates for the band of relative frequencies [1, 1.8] showed a very good degree of agreement: mean = 0.0483 (4.83%) and SD = 0.016. In the range of frequencies n . 1.8 a more complicated dependency is evident. Analysis of phase characteristics shows that increase of the relative frequency from n = 0 to about 2.0 is accompanied by a nearly proportional increase of the normalized CPFC, JC*(n) = JC(q/qH), from 0 to about p. For such a range the solid curves 1 and 2 in Fig. 5B represent JN1*(n) and JP3*(n), respectively. The curves only slightly deviate from linearity (to emphasize a linear character of dependencies the relative frequency is shown in linear abscissa scale). Accordingly, we considered a linear approximating function to JC*(n), LpC (n) = kJ n
Fig. 4. (A) Time course of amplitude frequency characteristic (AFC) V(q,0, tn) = Mod[V(jq,0, tn)] shows gradual changes with time (time scale corresponds to t) in the amplitudes of the ERP frequency components and reveals multiple resonant peaks that are difficult to interpret. (B) Using the same data as in A, CAFCs for each component (N1, P2, N2 and P3) show the frequency composition of the ERP component. The qualitative similarity is evident in the frequency domain images of the components. Quantitative differences in component waveform shape were found to be reflected by the fH parameter. Frequency scales in A and B are for f.
(12)
where kJ is a coefficient representing slope. The kJ coefficients have been evaluated for all component waveforms using the standard procedure for regressions through the origin (the values JC*(nj) from j = 0 to j = 157, where nj = n0cj, n0 = 0.05 and c = 1.023293, were used as initial data). No statistically significant differences have been found for the slopes among the ERP components (N1, P2, N2 and P3). For the group of component waveforms (n = 168), the mean for kJ is 1.86 (SD = 0.235). Linear dependencies (dash-dot and dotted curves) in Fig. 5B are characteristic examples of approximation. Deviation between the CPFC and approximating function typically increases with the increase of n. In the interval [0, 1.8] the mean of the absolute value of maximum deviation of CPFC from linear dependency for the whole group of component waveforms (n = 168) was 0.0916 (3.1% in relationship to the corresponding value of CPFC). We therefore found reasonable agreement for both normalized CAFCs and normalized CPFCs (ERP components N1, P2, N2, P3) with analytical dependencies Eqs. (11) and (12) in the relative frequency range from n = 0 to n about 1.8. To provide an interpretation of this result in terms of a classical additive model of evoked potentials (Lopes da Silva and Mars, 1987; Regan, 1989), let us consider wC(t) as the sum of two separate sources of activity. For t [ [0, tC] wC(t) = wTC(t) + wRC(t), where wTC(t) is a transient term linked in time and polarity to the stimulus, and wRC(t) is a random term (a statistically independent noise component usually associated with background brain activity). We can assume that relationships Eqs. (11) and (12) are relevant just to WTC(jq) = ATC(q)exp[−jJTC(q)] and consequently our study suggests that ATC (q) = kA exp( − aq2 ), JTC (q) = bq where a =
q/q2H
(13)
and b = kJ/qH. It follows from this that
ITC (q) = kA exp( − aq2 )sin(bq)
(14)
D. Melkonian et al. / Electroencephalography and clinical Neurophysiology 108 (1998) 251–259
and Fig. 3C shows that Eq. (14) provides reasonable approximation to IN1(q) in the frequency range from 0 to 1.8qH. Taking ITC(q) as a reliable frequency domain equivalent of the time domain data, we may expect that its inverse transform, zTC(t) defined by Eq. (9), approximates wTC(t) in a satisfactory manner. The existence of an analytical solution of integral Eq. (9) for function Eq. (14) (Erdelyi, 1954) provides the following formula for the transient term of a particular ERP component; for t ≥ 0 + − (t) − zTC (t) = kA {exp[ − (t − b2 =4a] zTC (t) = zTC
p
− exp[ − (t + b)2 =4a]}(2 pa) − 1
(15) The equation is presented in the form of a d’Alembert’s solution which is characteristic of wave-like phenomena (Richards and Williams, 1972). According to this representation each component waveform may be considered as a + (t) and negativesuperposition of positive-going wave zTC − going wave zTC(t). Fig. 3A illustrates that Eq. (15) provides reasonable agreement with experimental curves. The dashdot, dotted, dash-dot-dot and broken curves, represent functions zTC(t − tC) calculated from Eq. (15) for N1, P2, N2 and P3 components, respectively. The parameters kA (mV⋅s), a (s2), b (s) and tC (s) are as follows: N1 - 0.594, 0.000154, 0.0376, 0.064; P2 - 0.291, 0.000036, 0.0181, 0.14; N2 0.806, 0.000103, 0.0330, 0.184 and P3 - 1.785, 0.000227, 0.0490, 0.252.
4. Discussion The main findings reported in this paper depended critically on the methodological innovation of component Fourier transform, which was designed to provide a means for frequency domain representation of specific ERP component waveforms. The limitations of FFT for short segment spectral analysis have led previous authors to use non-traditional methods of spectral analysis, autoregressive models (Lopes da Silva and Mars, 1987; Florian and Pfurtscheller, 1995) and the wavelet transform (Bertrand et al., 1994; Schiff et al., 1994). The SBF algorithm provides a tool for DSA which does not lose the connection with the classical theory of Fourier analysis. To highlight the differences between the FFT and SBF algorithms, the FFT is based on a classical theory of Fourier series which deals with expansions of periodic functions into a sum of sines and cosines for discrete frequencies that are multiples of the fundamental harmonic. The SBF algorithm on the other hand, is based on a theory of Fourier integrals which deals with non-periodic functions. Therefore, the notion of discrete harmonics is superseded by a frequency domain representation of a signal which can be dealt with as a continuous frequency function (spectral density) defined
257
for a semi-infinite interval of frequencies from 0 to ∞. Processing empirical data of finite length disguises these differences because the computation deals with discrete values of the time and frequency domain characteristics for physically meaningful time intervals and bands of frequencies. However for the purposes of this study the SBF algorithm was of primary importance for: (i) effective DSA from evenly or unevenly sampled data on segments of arbitrary length, (ii) computation of frequency characteristics for a logarithmic scale of frequencies, and (iii) transformation from the frequency to the time domain (to control for distortions of a short-time spectral representation of the ERP waveforms). Refinements in DSA revealed a clear change in spectrotemporal patterns of ERP data provided by the CFT (illustrated in Fig. 4). In light of a general approach to the treatment of spectro-temporal structure of the brain evoked potentials (De Weerd and Kap, 1981), the deep qualitative changes in temporal evolution of the frequency domain characteristics (Fig. 4A) may be considered as an indication of the inherently non-stationary character of the ERP data. In contrast, a simple and regular form of CAFCs (Fig. 4B) provides evidence for the property of local stationarity of ERP component waveforms. The most interesting property of CFT is that it causes a specific ‘filtering effect’ on a component waveform. By this we do not mean the existence of certain fixed bands of ‘natural’ frequencies. Rather, the study shows that the half-power frequency, fH, strongly depends on the individual characteristics of the component waveform and has different estimates for groups of N1, P2, N2 and P3 components. However, this work establishes a specific form of the frequency domain patterns of CAFC and CPFC depending on the relative frequency n = f/fH. It is striking that the Gaussian function Eq. (11) and the linear function Eq. (12) provide a practically exact description of CAFC and CPFC, respectively, for the relative frequency range from 0 to approximately n = 1.8. What is of particular interest for this finding is that frequency characteristics in this band are almost free of noise contributions which are inevitable with ERPs. It may be concluded that in the relative frequency domain the frequency counterparts of noise represent components with relative frequencies higher than 1.8. Therefore, applied to the classical additive model of evoked potentials, the point that we do consider to have established, is that in the range of relative frequencies from 0 to 1.8, the CAFC and CPFC can be considered as the frequency domain description of the transient term, wTC(t).The consistent relationship of these two major characteristics provided an explicit analytical description Eq. (15) for the transient term of the component waveform in the time domain. There is a satisfactory degree of agreement of time characteristics (see typical example in Fig. 3A), since the equation parameters kA, a and b were derived entirely from the frequency domain analysis, without any additional adjustment to make them fit the original component waveforms.
258
D. Melkonian et al. / Electroencephalography and clinical Neurophysiology 108 (1998) 251–259
they tend to produce a Gaussian distribution even if the source itself is not normally distributed. Therefore, the empirically grounded appearance of the Gaussian function in the wave-like description of an ERP component, may be relevant to the ‘mass’ effect of ERP multiple sources. The CFT technique, in addition to providing reliable estimation of the frequency content of the ERP components, provides the investigator with dynamic model of the component both in the frequency and time domains. The model provides a unique set of parameters, which can be considered as the simplest frequency domain description of ERP components, and is at the same time able to reproduce the main features of transient ERP behavior. The model therefore provides a compressed representation of the ERP component wave form. The four parameters, kA, a, b and tC represent dynamic properties of ERP component waveforms, complementing conventional time and amplitude ERP component measures. Moreover, the physiological and dynamic dimensions of specific ERP components provided by our model, may help to elucidate the nature of specific ERP components and determine their interrelationship with background EEG activity in health and disease states.
Fig. 5. (A) Solid lines show the superposition of the four CAFCs from B after their normalization in both amplitude and frequency. The dotted line shows the Gaussian function described by Eq. (11). (B) Solid curves 1 and 2 plot CPFCs for N1 and P3 components, respectively. The dash-dot and broken curves show the linear approximating Eq. (12) with kJ = 1.78 and kJ = 1.92, respectively.
The time domain equation for ERP component waveforms Eq. (15) is in agreement with a previous study of our group (Haig et al., 1997) and with the study of Verleger and Wascher (1995) in that it contains the Gaussian function. However a single Gaussian function has a number of inconsistencies, including: (i) the function is non-zero at all times which is impossible for a stimulus induced process; (ii) the shape of the function is symmetric which is not always the case for ERP component waveforms. Eqn (15) is free of these limitations. Indeed, the function has a zero value at t = 0 and a non-symmetric shape with steep rise and slowly decaying tail. However, the number of free parameters remains the same. A number of authors consider brain potentials as wavelike phenomena (Freeman, 1975; Nunez, 1994). From this point of view it is interesting to note that the time domain description of the component waveform may be presented as a superposition of positive and negative going-waves. Considerable evidence exists that dynamics of scalp potentials depends upon the types of synaptic activity and spatiotemporal characteristics of multiple current sources and sinks in the soma-dendritic membrane of cortical pyramidal cells (Mitzdorf, 1985). From the fields of probability and statistics, it is well known that if the sources of the process are multiple, and if they have some degree of variability,
References Bertrand, O., Bohorquez, J. and Pernier, J. Time-frequency digital filtering based on an invertible wavelet transform: an application to evoked potentials. IEEE Trans. Biomed. Eng., 1994, 41: 77–88. Bracewell, R. The Fourier Transform and Its Applications. McGraw-Hill, New York, 1965. De Weerd, J.P.C. and Kap, J.I. Spectro-temporal representation and timevarying spectra of evoked potentials. Biol. Cybern., 1981, 41: 101– 117. Cacers, C.A. Biomedical Telemetry. Academic Press, New York, 1965. Erdelyi, A. (Ed.). Tables of Integral Transforms, Vol. 1. McGraw-Hill, New York, 1954, p. 78. Florian, G. and Pfurtscheller, G. Dynamic spectral analysis of eventrelated data. Electroenceph. clin. Neurophysiol., 1995, 95: 393–396. Freeman, W.J. Mass Action in the Nervous System. Academic Press, New York, 1975. Haig, A.R., Gordon, E., Rogers, G. and Anderson, J. Classification of single-trial ERP sub-types: application of globally optimal vector quantization using simulated annealing. Electroenceph. clin. Neurophysiol., 1995, 94: 288–297. Haig, A.R., Rennie, C. and Gordon, E. The use of Gaussian component modelling to elucidate average ERP component overlap in schizophrenia. J. Psychophysiol., 1997, 11: 173–187. Lopes da Silva, F.H. and Mars, N.J. Parametric methods in EEG analysis. In: A.S. Gevins and A. Remond (Eds.), Methods of Analysis of Brain Electrical and Magnetic Signals. EEG Handbook (Revised Series). Elsevier, Amsterdam, 1987, pp. 243–260. Melkonian, D.S. Transients in Neuronal Systems (in Russian). Publ. House of Armenian Acad. of Sciences, Yerevan, 1987. Mitzdorf, U. Current source density method and application in cat cerebral cortex: 14 investigation of evoked potentials and EEG phenomena. Physiol. Rev., 1985, 65: 37–100. Norcia, A.M., Sato, T., Shinn, P. and Mertus, J. Methods for the identification of evoked response components in the frequency and combined time/frequency domains. Electroenceph. clin. Neurophysiol., 1986, 65: 212–226.
D. Melkonian et al. / Electroencephalography and clinical Neurophysiology 108 (1998) 251–259 Nunez, P.L. Neocortical Dynamics and Human EEG Rhythms. Oxford University Press, New York, 1994. Regan, D. Human Brain Electrophysiology: Evoked Potentials and Evoked Magnetic Fields in Science and Medicine. Elsevier, New York, 1989. Richards, J.P.G. and Williams., R.P. Waves. Penguin Books, Middlesex, 1972. Roitbak, A.I., Fanardjian, V.V., Melkonian, D.S. and Melkonian, A.A. Contribution of glia and neurons to the surface-negative potentials of the cerebral cortex during its electrical stimulation. Neuroscience, 1987, 20: 1057–1067.
259
Schiff, S.J., Aldroubi, A., Unser, M. and Sato, S. Fast wavelet transformation of EEG. Electroenceph. clin. Neurophysiol., 1994, 91: 442–455. Stampfer, H.G. and Basar, E. Does frequency analysis lead to better understanding of human event related potentials. Int. J. Neurosci., 1985, 26: 181–196. Verleger, R. and Wascher, E. Fitting ex-Gauss functions to P3 waveshapes: an attempt at distinguishing between real and apparent changes of P3 latency. J. Psychophysiol., 1995, 9: 146–158.