Received Date: December 24, 2013; Accepted Date: January 20, 2014; Published Date: January 27, 2014
Citation: Mader M, Klatt J, Amtage F, Hellwig B, Mader W, et al. (2014) Spectral and Higher-Order-Spectral Analysis of Tremor Time Series. Clin Exp Pharmacol 4:149. doi: 10.4172/2161-1459.1000149
Copyright: © 2014 Mader M, et al. This is an open-access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.
Visit for more related articles at Journal of Clinical & Experimental Pharmacology
Objectives: The aim of this study was to investigate linear and nonlinear properties of tremor time series. To this end, we applied linear (second order) and nonlinear (higher order) spectral and cross-spectral analysis to 58 Electroencephalographic (EEG) and Electromyographic (EMG) tremor recordings of seven essential and five Parkinsonian Tremor patients.
Methods: Second and third order spectral analysis was performed on two types of data. First, data was simulated from a model mimicking the nonlinear properties of the tremor time series. Limitations of linear second order spectral analysis are illustrated in those simulations. Those limitations can be overcome by nonlinear third order spectral analysis. Second, tremor recordings from the trembling hand and contralateral motor area of the brain of the tremor patients were analyzed both by second and third order spectral analysis.
Results: Linear spectral analysis suggested nonlinearity of dynamics and interactions of processes measured by EEG and EMG. Applying bispectral analysis those nonlinearities were investigated. We found that a measure for nonlinearity based on bispectra was significant for most EMG recordings, as well as for the interaction of EEG and EMG.
Conclusions: Linear spectral analysis is a powerful tool when assessing spectral properties of time series. However, linear techniques fail to reveal nonlinear couplings. In the application of nonlinear spectral analysis to tremor time series, we showed that the dynamics of the hand muscle as well as the interaction of hand muscle and brain are governed by nonlinearities.
Higher order spectral analysis; Nonlinearity; Bispectrum; Tremor
A powerful tool for the reconstruction of process properties from measured data is spectral analysis . Univariate spectral analysis quantifies the frequency content of each process. Bivariate spectral analysis quantifies the interaction of two or more processes, revealing functional connectivity . Mostly, linear spectral analysis is applied, specifically the auto-spectrum in the univariate case and the coherence in the bivariate case [3-5]. An auto-spectral peak reveals that a process oscillates at the frequency at which the peak occurs. Peaks in the coherence spectrum reveal a linear interaction of two processes at the peak frequency. Strictly speaking, linear spectral analysis is restricted to revealing linear properties of a process. Nonlinear properties of the processes, such as quadratic phase coupling, must be investigated by nonlinear spectral analysis. As we show, higher harmonic spectral peaks cannot be distinguished from peaks which occur due to an independent sub-process by linear analysis. Thus independent sub-processes may be erroneously postulated, when truly nonlinear couplings are present.
A solution to this has been provided in the 1960s in terms of the extension of linear or second order spectral analysis to higher orders [1,6]. The latter enables quantification of higher order properties. Third order spectral analysis, e.g., quantifies non-gaussianity and quadratic phase couplings within processes [1,6]. By means of third order spectral analysis, higher harmonics can be revealed
Despite its benefits, higher order spectral analysis has remained limited in its influence onto the neurosciences [7,8]. Especially bivariate higher order spectral analyses for the quantification of nonlinear couplings of different brain areas  or coupling of the brain and limbs are rare. In this article, we apply, for the first time to our knowledge, both linear and nonlinear spectral analysis to investigate univariate and bivariate properties of the motor area in the brain and the muscles of trembling limbs of a cohort of tremor patients.
Tremor is defined as involuntary rhythmic movements of distinct body parts, predominantly the upper limbs. Oscillatory neuronal activity in the brain network generates the trembling of the limbs. Neuro-degenerative diseases, structural brain abnormalities, or toxicmetabolic conditions may lead to different forms of tremor. Common causes for limb tremor are Parkinson’s disease or essential Tremor.
An interaction of neural and muscle activity - inferred based on the Electroencephalogram (EEG) and Electromyogram (EMG) data - at the tremor frequency is expected, if neural malfunction is the cause of trembling muscles . This is reflected in an EMG-EEG coherence peak at the tremor frequency . Coherence peaks at twice the tremor frequency have been observed and interpreted - based on linear spectral analysis - as independent from the tremor-based peaks at the tremor frequency . However, such peaks are not distinguishable from higher harmonics by second order spectral analysis.
In a simulation study based on a model of coupled sinusoids we show that peaks similar to the ones occurring in tremor data can be higher harmonics. They can be detected by the third order spectral, i.e., bispectral, analysis. Applying bispectral analysis to EEG data from the motor area of the brain and EMG data from the muscles of the trembling hand of the tremor patient we uncovered that tremor is governed by nonlinear processes and interactions. We found that most coherence peaks at twice the tremor frequency are caused by higher harmonics rather than independent processes.
The article is organized as follows. In material and methods, the data of the study is introduced (in material) and methods of spectral analysis are summarized (in methods). To this end, first linear (in linear spectral analysis) and nonlinear (in nonlinear spectral analysis) spectral analysis is reviewed. The power of their combination is investigated in a simulation study (in combination of linear and nonlinear spectral analysis). In results, the results of linear (in linear spectral analysis of tremor data) and nonlinear (in nonlinear spectral analysis of tremor data) spectral analysis of the tremor data are shown. We conclude with a discussion of our findings in conclusion.
In the first part of this section we summarize the properties of the data and exclusion criteria for the analysis. In the second part linear and nonlinear spectral analysis are summarized.
Tremor data was recorded in 2-5 min segments from seven patients with Essential Tremor (ET) and five patients with Parkinsonian Tremor (PT). Each segment consisted of the rectified EMG of a patient’s trembling wrist muscle and the contralateral surface EEG of the motor area of the brain. The wrist muscles, from which the EMGs were recorded, were extensor and exor. The surface EEG-electrodes contralateral to the trembling hand were the C3 (left) and the C4 (right) electrodes, respectively. For each patient three to six such segments of simultaneous EEG and EMG were recorded. For patients with bilateral tremor, recordings from left and right wrist muscles were recorded. Segments without clear tremor in the EMG were excluded, leading to 58 segments, 30 of which were from ET, 28 from PT patients. For more details of the patients and selection criteria, see [12,13]. The data are available for download .
Spectral analysis is one of the most frequently applied tools in data analysis. It quantifies linear and nonlinear properties of a process in the frequency domain. Its application to medical data ranges back to the 1940s , when the foundation of linear spectral analysis was laid [6,15]. Linear spectral analysis is the most frequently applied type of spectral analysis. It quantifies the linear properties concerning the frequency content of the investigated processes. Nonlinear extensions have been developed in the 1960s [1,6,16]. The two types of spectral analysis are summarized in Sections (linear spectral analysis) and (nonlinear spectral analysis).
Both linear and nonlinear spectral analysis may be subdivided into univariate and bivariate analysis tools . Univariate spectral analysis quantifies the spectral properties of a single process, while bivariate spectral analysis quantifies interaction of two or more processes in the frequency domain. The subsequent summary follows , if not indicated otherwise.
Linear spectral analysis: For the sake of simplicity, we here consider two zero-mean processes x1 (t) and x2 (t). The auto-spectrum of xi (t), i=1, 2 is
It is proportional to the Fourier transform of the auto-covariance vari (τ)=E [xi (t) xi (t + τ)] of the process xi at time lags τ. The E [.] denotes the expected value.
of the two processes x1 and x2 is given by the Fourier transform of the cross-covariance . It is the bivariate generalization of the auto-spectrum. The amplitude of the cross- and auto-spectrum quantifies the covariance of the sub-processes at a given frequency ω.
Normalization of the absolute value of the cross-spectrum by the auto-spectra of the two processes, yields the coherence,
It quantifies linear interaction of two sub-processes. The coherence attains its maximum value 1 at frequency ω if x1 is a linear function of x2 at that frequency. The minimum value 0 is attained if the sub-processes are independent of each other.
In applications, the process itself is unknown. Respective spectra have to be estimated from measured data x (j), instead. Here, j refers to time when sampled at a sampling frequency . For the estimation, the data of length N is cut into m independent segments of length L=N / m, which is rounded off in case N / m is no integer. Each segment is tapered in order to reduce leakage of power into neighboring frequency bins in the frequency domain. To this end, the data is multiplied by special windowing functions (for review, see . We used a raised cosine as windowing function with taper parameter . The periodograms and cross-periodograms ,
of each tapered segment of data xa1 and xa2 are computed, squared and averaged. In case of the periodogram a1=a2=1, 2, while in case of the cross-periodogram a1=1 and a2=2, or vice versa. This yields an asymptotically consistent estimate of the respective spectrum [1,6,17].
Nonlinear spectral analysis: The idea of nonlinear, i.e., higher order, spectral analysis is to generalize the definition of the spectrum, Equation (1), by Fourier transforming k-th order momen ts
Since here we consider only two processes x1 and x2, the indexes a1,…, ak are either 1 or 2, depending on whether auto- or cross-covariances are considered. The k-th order moment, Eq. (5), of the two processes x1, x2 are considered at different time lags τj , j = 1,..., k-1 , in contrast to only one lag as for the second order covariance cov12(τ ), which in the univariate case is the auto-variance vari (τ ), i = 1, 2 . For k=2, the Fourier transform of the second order moment yields second order spectra, auto- or cross-spectrum for or , respectively, in Equations (1) and (2). For k=3, the Fourier transform of Eq. (5) results in a third order spectrum, the bispectrum,
where the indexes a1, a2, a3 denote the indexes 1, 2 of the processes x1, x2 used in the covariance obtained from Eq. (5). There are eight triple-combinations of 1 and 2 yielding a set of eight bispectra: First of all, there are two auto-bispectra B111 and B222, which are symmetric in their frequencies. Second, there are six cross-bispectra , where (.)' denotes transposition. This yields six independent auto- and cross-bispectra, since some of the cross-bispectra are the transposed others.
The bispectrum Ba1, a2, a3 quantifies the coupling of processes xa1 at frequency ω1 and xa2 at ω2 onto xa3 at ω3=ω1 + ω2. Within this article, only the absolute value of the bispectrum is considered.
Since in applications the true bispectra Ba1a2a3 are unknown, they are estimated. Averaging and taking the squared of the respective biperiodogram
of m=N/L blocks of tapered data with sub-process indexes a1, a2, a3 as listed in the eight possibilities of bispectra above, i.e., 111, 112 etc. Estimation thus is analogue to the linear case in linear spectral analysis. It yields a consistent estimate .
The bispectrum is sensitive, both to nonlinearity and non- Gaussianity. To ensure that a rejection of the null hypothesis is due to nonlinearity rather than non-Gaussianity, any data within this article was transformed to follow a Gaussian distribution prior to spectral analysis. To ensure Gaussian distribution of the data x(1),…, x(N), N Gaussian random numbers z(1), . . . , z(N) were drawn. Both x and z were sorted according to their value, yielding . Then were assigned to each other, An original data point x(j), being the i-th highest value in x, was thus substituted by the i-th highest Gaussian random number, .
Kim and Powers introduced a normalized estimated version of the bispectrum, the bicoherence 
It is estimated by averaging biperiodograms Pa1a2a3 to obtain an estimate of the bispectrum in the numerator, as well as averaging periodograms denoted by E[.] in the denominator. Similar to the coherence in the linear case, the bicoherence is normalized to values between zero and one. However, other than the coherence, which is restricted to the bivariate case, the bicoherence exists for the univariate case, as well. In the univariate case, the auto-bicoherence (a1=a2=a3, i.e., B111 or B222) refers to nonlinear self-coupling of a single process. The cross-bicoherence (B112, B212 etc.) quantifies nonlinear interaction of two processes x1 and x2 in the bivariate case.
In order to test for the statistical significance of bicoherence, block-bootstrap may be employed [19,20]. To this end, first the data is segmented into a set of m independent segments. Second, m segments are drawn randomly with replacement from this set, yielding a bootstrap realization for which the null hypothesis of independence applies. Independence refers to absent self- and cross-coupling. Third, for each block of the bootstrap realization the periodograms and biperiodograms of interest are derived. Fourth, the periodograms and the biperiodograms of all m segments are squared and averaged yielding estimates of the bicoherence of the bootstrap realization according to Eq. (8). This procedure is repeated r times, in order to generate r bootstrapped bicoherences. Fifth, empirical histograms of the bicoherences of all r bootstrap realizations at a given frequency are compiled. Finally, for a given significance level α , the 1-α quantile is considered the critical value. This refers to the α. r-th highest bootstrapped bicoherence value. If the bicoherence of the original data exceeds the critical value, the original bicoherence value exceeds zero, significantly.
Note that the null hypothesis of independence needs to be fulfilled in the bootstrap-realizations. In the case of cross-bicoherence, this is ensured by drawing random segments from the data of the considered processes at different time points. In the case of auto-bicoherence, this is achieved by ensuring segments of independent time points within the data. For all statistical tests within this article, we performed the block-bootstrap as outlined with 20 bootstrap samples and tested for a significance level of α =5%. Thus, if a bicoherence value derived from the data exceeded the highest bicoherence value obtained from the bootstrap samples, the bicoherence was considered significant with an error rate of 5%.
Combination of linear and nonlinear spectral analysis: In order to illustrate the power of bispectral analysis, consider the sinusoidal model 
with respective frequencies ω and υ, normally distributed random phases ηi (t) and ξi (t), as well as amplitudes C1 and D1for i=1; 2. Phases η3, ξ3, and frequencies ω3, υ3 represent coupling as described in the following. An auto-bispectral peak in B111 at the frequency-tuple (ω1, ω2) occurs if the sine functions of x1 are phase-locked, i.e., both ω3=ω1 + ω2 and η3 (t)=η1 (t) + η2 (t) hold. A cross-bispectral peak B122 (ω1, υ2) occurs if frequencies and phases are locked according to ω1 + υ2=υ3 and η1 (t) + η2 (t)=ξ3 (t). On the contrary, no peak occurs in the crossbispectrum B112 (ω1, ω2), if
To assess the limits and power of spectral and bispectral analysis, two settings were investigated: bispectral self-coupling, ω3=ω1 + ω2 and η3=η1 + η2, and bispectral cross-coupling, ω3=ω1+ υ2 and η3=η1+ ξ2. Within all simulations, we used m=120 blocks of length L=2500 data points for spectral estimation.
In order to investigate auto-bispectral properties, a realization of process x1(t) of Equation (9) with C1=C2=C3=1 was simulated. Gaussian white noise of variance 25 was added to simulate measurement noise, which is present in most applications. This yielded a signal to noise ratio of approximately one. The frequencies were ω1=4 Hz, ω2=9 Hz and ω3=13 Hz. The phases η1,2 were chosen randomly from a uniform distribution on the interval [0, 2π]. Random phase-shifts η were drawn for each segment, however, phase shifts remained constant within each segment for which the periodograms and biperiodograms were computed.
Auto-bispectral coupling was simulated by phase coupling
Accordingly, the estimated auto-bicoherence had a peak at the frequency tuple (4 Hz, 9 Hz), see Figure 1(a). Note that the autobispectrum is symmetric with respect to the angle bisector. Without phase-coupling but independent random phases η3, no such peaks occurred, as shown in Figure 1(b). While in case of phase-coupling (a), the auto-bispectrum is significant at the oscillation frequency, it is not in case of absent coupling (b).
Figure 1: Auto-bicoherence of a sinusoidal process as in Eq. (9) with (a) and without (b) phase coupling (Eq. (11)). In the case of present phase-coupling at frequency ω1=4 Hz and ω2=9 Hz (a), the auto-bicoherence exhibits a clear peak, which is significant. For absent phase-coupling (b) no peak occurs, the auto-bicoherence is not significant.
Analogously, the bicoherence can be assessed for two identical frequencies ω1=ω2, instead of two different ones, as above. The auto-bicoherence then resolves whether an oscillation at ω3=2ω1 is independent or occurs due to phase coupling induced by nonlinearities, leading to a higher harmonic. Similarly, the cross-bicoherence quantifies phase-coupling of two processes as modeled by the functions of sinusoids (Equations (9) and (10)).
In the second setting, we simulated such cross-coupling with identical frequencies ω1=ω2=υ1=υ2=4Hz and independent random phases η1,η2,ξ1,ξ2, respectively. The third components of both processes x1 and x2 were modeled as sines at the frequencies ω1+ υ2=ω3 and υ1+ υ2=υ3. Only process x2 coupled into the phase of process x1, such that the phase of x1 was
The phase ξ3 of process x2 was independent uniformly distributed noise. Since the signal to noise ratio does not need to be of the same order in both processes, the variance of the observational noise of x2 was reduced to one and amplitudes were set C1=C2=D1=D2=D3=1 with coupling amplitude C3=3. This yielded a signal-to-noise-ratio of approximately one for sub-process x1 and 26 for x2. The coherence exhibited clear peaks (Figure 2) at both 4 and 8 Hz. Due to the three times higher amplitude C3 of the coupling term in process x1, the coherence peak was higher at 8 than at 4 Hz. Cross-bispectral analysis correctly revealed a higher harmonic at 8 Hz due to a significant peak in the cross-bicoherence Bcoh121, Figure 3(a). None of the other crossbicoherences (Figure 3(b-d)) or auto-bicoherences was significant at the oscillation frequency.
Figure 3: Cross-bicoherence of a sinusoidal set of two processes as in Eqs. (9) and (10) with nonlinear phase coupling (Eq. (12)) from x1 and x2 onto x1, only. Only a peak in the cross-bicoherence Bcoh121 occurs at the oscillation frequency of 4 Hz (a). Other cross-bicoherences at (4 Hz, 4 Hz) are insignificant (b-d).
In order to evaluate the power of the statistical test for bicoherence based on block-bootstrap, we simulated 100 realizations of the sinusoidal model, Eqs. (9) and (10), for both auto- and cross-couplings for increasing coupling strengths. Figure 4 shows the percentage of significant phase auto-couplings (a) and cross-couplings (b) for increasing coupling strengths C3. The size of the test is correct, such that according to the significance level (black, dotted), only 5% false positive conclusions occur when no coupling, C3=0, exists. The steep increase with coupling strength for all scenarios shows that present couplings are detected reliably.
Figure 4: Power of the statistical test of auto-bicoherences in case of auto- (a) and cross-coupling (b) of increasing coupling strengths (x-axes). For absent coupling, the test obeys the significance level of 5% (black, dotted), since the percentage of rejections of the null hypothesis is approximately 5%.
To summarize, bispectral analysis is capable of revealing higher harmonics both in the univariate and bivariate case. Our simulations showed that due to such higher harmonics larger values of the coherence at twice the oscillation frequency are possible.
We investigated second (in linear spectral analysis of tremor data) and third order spectral properties (in nonlinear spectral analysis of tremor data) of the tremor data. We Gaussianized the data as described in combination of linear and nonlinear spectral analysis prior to spectral analysis. For the statistics based on block-bootstrap, the data was divided into segments of L=2500; 5000; 10000 data points, corresponding to 2.5; 5 or 10 seconds. Those block lengths were chosen such that successive segments were independent based on the correlation structure given by the autocorrelation function. Details about the block length selection are given in Appendix A (Figure 5). Results of linear and nonlinear spectral analysis are given in Tables 1 and 2.
The asterisk * refers to the one segment in which no coherence peak at the tremor frequency occurred. ‘ID’ refers to patient identity number, ‘type’ to type of tremor (ET, PT), ‘seg’ to segment number, ‘side’ to hand side (L=left, R=right), ‘muscle’ to type of muscle from which it was recorded (‘E’=extensor, ‘F’=exor), ‘TF’ to tremor frequency in Hz, ‘Group’ to which the segment was assorted based on linear cross-spectral properties (A=no coherence peak at twice the tremor frequency, B=coherence peak at twice the tremor frequency), ‘BL’ to block length for spectral analysis, i.e. for duration of decay of the autocorrelation function, ‘BE’ and ‘BM’ to auto bicoherence of the EEG and EMG, respectively, ‘XB1’ to cross-bicoherence of the EEG and EMG onto the EEG, ‘XB2’ to cross-bicoherence of the EMG and EEG onto the EMG, ‘XB3’ to cross-bicoherence of the EMG onto the EEG, ‘XB4’ to cross-bicoherence of the EEG onto the EMG.
Table 1: Results of spectral analysis in case of consistent spectral peaks at tremor frequency and twice the tremor frequency.
‘ID’ refers to patient identity number, ‘type’ to type of tremor (ET, PT), ‘seg’ to segment number, ‘side’ to hand side (L=left, R=right), ‘muscle’ to type of muscle from which it was recorded (‘E’=extensor, ‘F’=exor), ‘TF’ to tremor frequency in Hz, ‘Group’ to which the segment was assorted based on linear cross-spectral properties C=patients with inconsistent spectral peaks at the tremor and twice the tremor frequency, ‘BL’ to block length for spectral analysis, i.e. for duration of decay of the autocorrelation function, ‘BE’ and ‘BM’ to auto bicoherence of the EEG and EMG, respectively, ‘XB1’ to cross-bicoherence of the EEG and EMG onto the EEG, ‘XB2’ to cross-bicoherence of the EMG and EEG onto the EMG, ‘XB3’ to cross-bicoherence of the EMG onto the EEG, ‘XB4’ to cross-bicoherence of the EEG onto the EMG.
Table 2: Results of spectral analysis in case of inconsistent spectral peaks (group C) at tremor frequency and twice the tremor frequency.
Figure 5: Autocorrelation functions of three exemplary segments of EMG data of tremor 1 patients. The time after which the autocorrelation function is decayed determines the block lengths L=2.5 s (a, patient 2 segment 2), 5 s (b, patient 3 segment 3), or 10 s (c, patient 8 segment 2) used for spectral estimations.
Linear spectral analysis of tremor data
Linear spectral analysis was applied to all 58 segments included in this study. An exemplary set of EEG (a) and EMG (b) data as well as the respective estimated spectra (b,d) are shown in Figure 6. All segments exhibited auto-spectral peaks in the EMG. In each segment, the tremor frequency was defined as the lowest frequency in the range of 4 to 8 Hz at which the auto-spectrum of the EMG had a clear peak.
Based on auto- and cross-spectral analysis of the EEG and EMG recordings, the segments were categorized in three major groups. Group A and B were distinguished from group C as to whether auto-spectral and cross-spectral peaks at the tremor frequency and twice the tremor frequency were consistent. A segment was denoted consistent if both the frequency of the second auto-spectral peak deviated less than 0.5 Hz from twice the tremor frequency, and the first coherence peak deviated less than 0.5 Hz from the tremor frequency. Otherwise a segment was denoted inconsistent. Group C contains two patients with at least one inconsistent segment, while group A and B contained only consistent segments.
Group A and B were distinguished from each other by the existence of a cross-spectral peak at twice the tremor frequency. While group A contained segments without coherence peaks at twice the tremor frequency, group B contained segments with such peaks. Group A consisted of 16 segments of six ET patients, group B consisted of 33 segments of four ET and four PT patients. All but one segment (* in Table 1) of group B exhibited a coherence peak not only at twice the tremor frequency but also at the tremor frequency. In nine of these 33 segments, the coherence peak at twice the tremor frequency was larger than the one at the tremor frequency. As simulated by the sinusoidal model described in Section “Combination of Linear and Nonlinear Spectral Analysis”. , this higher coherence peak at twice the tremor frequency can be modeled by a strong bispectral coupling onto the target process, which however can be quantified by bispectral analysis, only. While group A only consisted of segments of ET patients, the fraction of ET segments (11 of 33) in group B was considerably lower than that of PT segments (22 of 33).
All segments of consistent PT patients were in group B. However, out of all consistent ET patients only one was exclusively assigned to group A. All segments of the other five consistent ET patients were assigned to both groups A and B.
Group C contained six segments of one PT and three segments of one ET patient. Five of those nine segments were inconsistent with respect to tremor and twice the tremor frequency, four of which from the PT patient (patient 11 segments 1-4), the other one from the ET patient (patient 12 segment 2). In two of the six PT-segments (segment 1-2), the second coherence peak occurred at 1 Hz higher than twice the tremor frequency. In two other segments (segment 3-4), the first coherence peak occurred at the tremor frequency but the second peak occurred at a more than 1 Hz higher frequency than twice the tremor frequency. From the ET patient, two segments were consistent (segment 1 and 3). The first coherence peak in the inconsistent segment (segment 2) was at a more than 0.5 Hz lower frequency than the tremor frequency.
As described in the next section, we performed nonlinear spectral and cross-spectral analysis of the EEG and EMG to investigate the interaction properties of tremor and twice the tremor frequency.
Nonlinear spectral analysis of tremor data
The EEG-EMG segments of the three groups specified by linear spectral analysis were further investigated by bicoherence analysis. A bicoherence peak was considered significant if the critical value (significance level α=5%) from 20 bootstrap realizations was exceeded. We evaluated the auto- and cross-bicoherences at the tuple (tremor frequency, tremor frequency). Results are summarized in the above mentioned Table 1 for group A and B and Table 2 for group C. Figure 7 shows the auto- (a-b) and cross-bicoherences (c-f) of one patient, exemplarily. For brevity, we denoted the cross-bicoherence significant whenever one of the cross-bicoherences as exemplarily shown in Figure 7 was significant at the respective frequencies.
Auto-bicoherences of all but two rectified EMG-segments were significant at the tremor frequency. The two insignificant autobicoherences occurred in segments of group C. They both occurred in consistent segments (patient 11 segment 6 and patient 12 segment 3). On the contrary, the auto-bicoherences of EEG segments were insignificant at the tremor frequency in all but six segments. Thus, there was hardly any evidence for the nonlinearity of the EEG, while the rectified EMG exhibited nonlinear properties.
Cross-bicoherences of the 58 EEG and EMG segments were significant in 48 cases. Only four of the segments without significant cross-bicoherences were from group A and B, one from group B. On the contrary, hardly any significant cross-bicoherences were found in group C. In both consistent segments of the ET patient, cross-bicoherence was significant, while it was insignificant in the consistent segments of the PT patient. On the contrary, cross-bicoherence was insignificant in all but one (patient 11 segment 4) inconsistent segments.
Throughout all groups investigated, we therefore conclude that a majority (83%) of the EEG and EMG segments are phase crosscoupled. Most segments (32 segments) exhibited cross-bicoherence peaks indicating a phase cross-coupling from the EMG and EEG onto the EMG, as well as from the EMG onto the EEG (27 segments). On the contrary, phase cross-couplings from the EEG and EMG onto the EEG (13 segments) as well as from the EEG to the EMG (17 segments) were seldom.
While 81% of the segments in group A had significant crossbicoherences 97% of the segments had significant cross-bicoherences in group B. The nine segments with larger coherence peaks at twice the tremor frequency than at the tremor frequency all had significant crossbicoherences, as expected from the simulations of the sinusoidal model in Section “Combination of Linear and Nonlinear Spectral Analysis”. In contrast to the groups A and B, where spectrum and coherence had consistent peak frequencies, the percentage of significant crossbicoherences in group C was considerably lower. Only 33% of the segments in group C had a significant cross-bicoherence. The low fraction of significant cross-bicoherences is in accordance with the fact that the segments in group C exhibited inconsistent linear spectral properties.
We therefore conclude that tremor is governed by nonlinear phase cross-couplings of the motor area of the brain and the trembling muscle. In almost all cases when a coherence peak at twice the tremor frequency occurred, cross-bicoherence analysis revealed it as a higher harmonic.
Based on a sinusoidal system modeling phase coupling of one or two processes we showed the limitations of second order spectral analysis and how they can be overcome by third order spectral analysis. For the statistical evaluation of third order spectral analysis, we presented a bootstrap approach. It is based on a randomization of the data recorded, such that no additional measurements are necessary. As our simulations show, the employed bootstrap method is powerful in detecting phase-coupling within one and between sub-processes. At the same time it keeps the coverage correct if phase-coupling is absent.
Based on the simulations, we analyzed the data of twelve patients of Essential Tremor (ET) and Parkinsonian Tremor (PT) with respect to their spectral properties. Both linear and nonlinear spectral properties were extracted from 58 segments of simultaneous EEG and rectified EMG recordings with clear tremor activity of approximately 2-5 min duration.
Auto- and cross-bispectral analysis of the tremor data showed that peaks in spectra and cross-spectra at double the tremor frequency are higher harmonics, not independent processes or couplings.
Linear spectral analysis revealed that more than half of the segments exhibit coherent EEG and EMG activity at twice the tremor frequency. In 90% of those cases bicoherence analysis revealed the nonlinearity of the coupling of the EEG and EMG at the tremor frequency. This means that coherence peaks at twice the tremor frequency are linked to the pathological interaction at the tremor frequency. In nine of the 58 investigated segments, the coherence was higher at twice the tremor frequency than at the tremor frequency. In all of those segments a significant cross-bicoherence peak was found at the tuple of tremor frequencies. In a sinusoidal model we showed, how such higher coherence peaks at twice the tremor frequency may be induced due to phase-couplings based on nonlinearities.
Furthermore, our analyses showed that the coupling is governed by a phase coupling from the EMG and EEG onto the EMG (32 segments) and from the EMG to EEG (27 segments) as opposed to couplings from the EEG and EMG onto the EEG (17 segments) and from the EMG to the EEG (13 segments). The tendency of fewer connections onto the EEG is explainable by the error-in-variables problem. According to this problem, connections from processes with low variance onto processes with large variance are expected to be underestimated. Besides crossbicoherence analysis, auto-bicoherence analysis was performed on EEG and EMG. We showed that the EMG is nonlinear, while nonlinearity could not be detected in the EEG in general.
The major difference of ET and PT segments is the occurrence of a coherence peak at twice the tremor frequency. All but one PT patient had a coherence peak at twice the tremor frequency in all the segments included in the study. In contrast, only five ET patients had coherence peaks at twice the tremor frequency in some of their segments. In two ET patients a coherence peak at twice the tremor frequency did not occur in any of the segments. Bicoherence properties of ET and PT patients however were comparable.
Other than in previous studies based on higher order spectra in tremor analysis [8,9] we combined the analysis of second and third order spectra of the EEG and EMG, analyzing both the properties of EEG and EMG separately and reconstructing their second order interactions. Jakubowski et al.  presented a classifier based on 30 autospectral features. They classified ET and PTs based on auto-spectral and higher order auto-spectral analysis of the accelerometer placed at the trembling hand. They aimed at classifying ET and PT patients based on 30 features. Marceglia et al.  performed cross-spectral and crossbispectral analysis of STN-local field potentials and the EMG for one ET patient, only. This patient had a significant cross-bicoherence at the tremor-frequency from the EMG and EEG onto the EMG. This is also the type of coupling we found in our study of twelve ET and PT patients most frequently. To the first time to our knowledge, first and second order auto- and cross-spectral analysis was performed on an extended set of ET and PT segments.
From the cohort we investigated we conclude, that tremor, both of ET and PT patients, is governed by nonlinear interactions of the brain and hand. We quantified nonlinear interaction by bispectral analysis both within and across EEG and EMG tremor time series. We further conclude that both linear and nonlinear spectral analyses are valuable to reconstruct spectral properties of processes. This method is not restricted to temporal EEG and EMG recordings. On the contrary, we propose to analyze both spatial and temporal data from other imaging techniques analogously.
That way, nonlinear spectral properties of the processes investigated can be revealed and interpreted.