Rule-Based EEG Classifier Utilizing Local Entropy of Time–Frequency Distributions

Electroencephalogram (EEG) signals are known to contain signatures of stimuli that induce brain activities. However, detecting these signatures to classify captured EEG waveforms is one of the most challenging tasks of EEG analysis. This paper proposes a novel time–frequency-based method for EEG analysis and characterization implemented in a computer-aided decision-support system that can be used to assist medical experts in interpreting EEG patterns. The computerized method utilizes EEG spectral non-stationarity, which is clearly revealed in the time–frequency distributions (TFDs) of multicomponent signals. The proposed algorithm, which is based on the modification of the Rényi entropy, called local or short-term Rényi entropy (STRE), was upgraded with a blind component separation procedure and instantaneous frequency (IF) estimation. The method was applied to EEGs of both forward and backward movements of the left and right hands, as well as to EEGs of imagined hand movements, which were captured by a 19-channel EEG recording system. The obtained results show that in a given virtual instrument, the proposed methods efficiently distinguish between real and imagined limb movements by considering their signatures in terms of the dominant EEG component’s IFs at the specified subset of EEG channels (namely, F3, F4, F7, F8, T3, and T4). Furthermore, computing the number of EEG signal components, their extraction, and IF estimation provide important information that shows potential to enhance existing clinical diagnostic techniques for detecting the intensity, location, and type of brain function abnormalities in patients with neurological motor control disorders.


Introduction
Activities of neurons in the cerebral cortex cause differences in electrical potentials, which may be mapped in time, frequency, and space when captured by electrodes placed on standard positions on a patient's head. The recorded signals, referred to as an electroencephalogram (EEG), vary from a few microvolts to several hundred microvolts, providing important information on various neurological brain disorders [1][2][3][4]. Due to the high temporal resolution and low recording cost, as well as its noninvasiveness, EEG has attracted vast scientific research, which has provided plenty of physical, psychological, and pathological information on the neurological activity of the brain [5][6][7][8]. This has resulted in better understanding of both clinical diagnoses and treatment of some illnesses, as well as in progress in cognitive science [5].
Neurophysiologists often correlate functions of the central nervous system to various EEG patterns based on empirical visual EEG inspections [9,10]. However, visual pattern recognition may become tedious and time-consuming, especially in the case of the long EEG datasets that are inspected by medical experts searching for abnormalities in EEG waveforms [11]. Moreover, it is subject to misinterpretations, such as in cases of cerebral metabolic disturbances and convulsive disorders [9,12]. For example, an isoelectric EEG may be the result of selective pCO 2 increases when the brain is sufficiently supplied with O 2 [9]. On the other hand, normal pCO 2 levels may occur in the case of cerebral oxygen deficiency [9]. However, the two cases have significantly different outcome prognoses [9]. In addition, traditional visual inspections do not allow adequate EEG systematization [13]. To overcome these limitations, quantitative EEG analysis has introduced computer-aided EEG processing techniques for measuring brain dynamics objectively, such as with the methods proposed in [13][14][15][16][17].
It is well known that EEG signals are highly non-stationary stochastic processes, especially when they are induced by external stimuli [18]. Cognitive tasks are reflected in both EEG spectral changes and activities located across different brain regions [18]. Localizing transient brain activities, both in space and time or frequency, in order to detect signatures of stimuli in EEGs has been the subject of extensive research attention in the last few decades, and is still one of the most challenging tasks in EEG analysis [18].
For example, the analysis of event-related potentials, which measure the brain's response to external stimuli, is often done using time-domain processing techniques, since non-stationary EEG signals and intensive background noise make the Fourier analysis inadequate for this purpose [18,19]. Furthermore, the Fourier transform ensures only the information on the signal's frequency content, without any time resolution [19,20]. In order to overcome this shortcoming, a windowed Fourier transform with a time-localized window function was introduced [19]. However, its performance was limited by the fixed window size, which affected the time and frequency resolution of the signal's timefrequency representation [20][21][22][23][24][25][26][27]. The trade-off between time and frequency resolution motivated the development of other time-frequency signal representations, which allowed their applications in various technical and biomedical fields [20].
One way to quantify the degree of signal complexity is to consider spectral entropy from the Fourier power spectrum, which is not defined as a function of time, and is thus suitable for stationary signals only [13]. In order to overcome this limitation, Powell and Percival defined a time-evolving entropy from the short-time Fourier transform (STFT) using the Hanning window [42]. An alternative entropy-based approach aiming to analyze the dynamics of time series using a discrete probability distribution based on pointwise wavelet leaders for evaluating the Shannon entropy was introduced by Rosenblatt et al. [43,44]. A comparison of global and pointwise information-theory-based quantifiers applied to EEG signal characterization was given in [13], where the normalized Shannon wavelet entropy and the wavelet statistical complexity were chosen as global, and the wavelet leader Shannon entropy and the wavelet leader statistical complexity were chosen as the pointwise quantifiers [13]. The first was shown to compute changes in the signal's frequencies, while the latter described the signal's morphological changes by considering its regularity [13]. Bhattacharyya et al. proposed measuring the complexity of epilepsy EEG signals with a multi-scale entropy [45]. More specifically, EEG signals were decomposed into the number of sub-bands using the tunable-Q wavelet transform (TQWT), and then the K-nearest neighbor entropies were cumulatively estimated from various sub-bands [45]. The proposed entropy measure was used to classify seizure, seizure-free, and normal EEG signals [45]. The TQWT-based method upgraded with multivariate fuzzy entropy was further applied to analysis and classification of multivariate sub-band focal and non-focal EEG signals [46]. Acharya et al. proposed classification of epilepsy EEG signals into normal, interictal, and ictal classes using a continuous wavelet transform, higher-order spectra, and textures [47]. Sharma et al. proposed an automatic approach to detecting epileptic seizures utilizing the analytic time-frequency flexible wavelet transform and the fractal dimension. The method also divides EEG signals into sub-bands with the fractal dimension calculated for each sub-band [30].
In the last decade, through the development of assistive technology, analysis of EEG features has shown potential to enhance techniques for detecting types of brain functions in patients with neurological motor control disorders, together with the possibility of using brain activity independently of muscles and peripheral nerves to develop brain-computer interfaces (BCIs). In [48,49], real (executed) left-right limb movement classification from EEG data was presented. The wavelet coefficients and power spectral density (PSD) of the alpha and central beta band and the average power of the respective bands were employed as features for classification after pre-processing of the captured signals. In [50], again, the mean power of the signal evaluated in the frequency domain (PSD) in eight EEG frequency bands and a linear-discriminant-analysis-based classifier were used to distinguish between real (executed) upper limb movements. Recently, several studies [51][52][53][54][55] have achieved classification of the same limb movements with appreciable performance using EEG data, although many of these studies classified only two movements. In [56], several classifiers (LS-SVM, QP-SVM, SMO, k-NN, SVM, etc.) were employed to detect single limb movement intentions.
The comprehensive method proposed in this paper introduces an EEG information complexity quantifier called local or short-term Rényi entropy (STRE), which is applied to limb movement EEG signal analysis in the time-frequency domain. Next, the STRE was upgraded with a blind component separation for EEG component localization and extraction, followed by the IF estimation of the dominant component. Thus, the method detects, localizes, and extracts spectral changes in EEG rhythms with the purpose of detecting and distinguishing signatures of various limb movements (both imagined and real movements used as stimuli) found in EEG time-frequency equivalents. The high values of the dominant EEG component IFs at specific electrodes were found to be signatures of both imagined and real limb movements. The ability to detect limb movement signatures in the time-frequency domain of the analyzed EEGs with high precision is the main contribution of this paper. Furthermore, the proposed novel algorithm for EEG signal characterization opens the door for building an extensive computer-aided decision-support system for numerous clinical applications with the purpose of supplying neurophysiologists with additional information on limb-movement-related brain activities and associated motor control abnormalities. We hope this could lead to the improvement of their medical treatment and diagnostics from noninvasive EEG records.
The paper is structured as follows. The proposed method is described in Section 2. In particular, Section 2.1 gives a brief overview of the time-frequency representations used. Section 2.2 introduces the method for measuring signal information content with the modification of the Rényi entropy in the time-frequency domain. The IF estimation and the component extraction algorithm are described in Sections 2.3 and 2.4, while the experimental setup and data are described in Section 2.5. The experimental results are thoroughly presented in Section 3, and are discussed in detail in Section 4. Conclusions are found in Section 5.

Signal Analysis in the Time-Frequency Domain
In order to obtain an exhaustive representation of multicomponent non-stationary signals, time-frequency distributions (TFDs) were introduced. TFDs include time dependency in the signal frequency spectrum, which is thus a two-variable function, C z (t, f ), defined over the two-dimensional (t, f ) domain (where t and f denote time and frequency, and z(t) is the EEG time series) [20,57,58]. Over the past few decades, numerous TFDs have been formulated in order to optimize the representation quality of different classes of signals. Most applications were focused on the following [20]: • Identification of signal characteristics (such as time and frequency variations and the number of signal components); • Extracting components from their mixtures and background noise; • Ability to synthesize extracted components in the time domain; • Analysis of features (such as instantaneous amplitude, frequency, and bandwidth of each component).
Until now, none of the existing TFDs have guaranteed ideal time-frequency representation (cross-term suppression and maximal time and frequency resolution for all classes of signals). In light of the above, in this paper, three TFDs are considered (namely, the spectrogram (SP), Wigner-Ville distribution (WVD), and Rihaczek distribution (RD), which are defined in the following), and the results obtained with the proposed method are compared in Section 3.

Spectrogram
The spectrogram, S z (t, f ), is a computationally simple TFD obtained by squaring the magnitude of the STFT [20,[57][58][59]: The STFT is a linear time-frequency transformation that introduces the frequency dimension in the signal representation by performing the Fourier transform only on the portion of the signal included inside the analyzing window w(t). Since the signal and the analyzing window cannot have arbitrarily small supports in time and in frequency simultaneously, the STFT and, hence, the spectrogram suffer from a limited time-frequency resolution [57,60].
It can be concluded that the spectrogram has difficulty in concentrating components' energy around the respective IFs. On the other hand, it does not produce cross-terms unless components overlap in the (t, f ) plane [59].

Wigner-Ville Distribution (WVD)
One of the most popular TFDs was introduced by Wigner [61] and was afterward extended to analytic signals by Ville [20,[62][63][64][65][66]. The WVD of an analytic signal z(t) = a(t)e jφ(t) (where a(t) is the instantaneous amplitude and φ(t) is the instantaneous phase), denoted as WVD z (t, f ), represents a mono-component signal with a(t) = 1 and quadratic phase as a knife-edge ridge in the (t, f ) plane (an elongated region of concentrated energy, the crest of which corresponds to the signal IF [20,67]) where the IF is defined as This leads to the signal kernel defined as Since φ (t) is not directly available, it can be replaced by the central finite-difference ap- , leading to the WVD given in the form [20,68,69] From Equation (6), it can be noticed that using the instantaneous autocorrelation function (IAF) as the kernel function brings nonlinearity in the WVD. The effects of this nonlinearity are most evident in the case of components with nonlinear IFs or multicomponent signals, producing undesirable interferences (often referred to as cross-terms) [20,59,68].

The Quadratic Class of Time-Frequency Distributions
To maintain the high concentration exhibited by the WVD, but to avoid or reduce unwanted interferences caused by its quadratic nature at the same time, a smoothed version of the IAF was introduced: where G(t, τ) is the time-lag kernel and * t denotes convolution in the time domain. Thus, the class of quadratic TFDs, defined as C z (t, f ) = F τ→ f {R z (t, τ)}, can be calculated as [20]: The time-lag kernel determines the performance of the TFD, allowing a trade-off between the time-frequency supports of the signal components and suppression of interferences.
In this paper, as an exponent of the quadratic class of TFDs, the RD was used. It is defined as [20] where Z( f ) is the Fourier transform of the signal z(t). The RD can also be expressed in the form: Once the signal time-frequency representation is calculated, its information content in terms of the number of signal components should be obtained. One of the ways of finding the number of frequency-modulated components based on the spectrogram is the method described in [70]. It exploits Kolmogorov complexity to model the information in the spectrogram, which is then converted into a binary map through automatic thresholding based on the minimum description length, and mode counting is performed through two-dimensional run-length encoding. However, the proposed approach is not limited to spectrograms only (as the algorithm in [70] is), and it also works for other quadratic time-frequency representations. Hence, another approach to counting the number of signal components was developed and utilized in this paper (based on the modification of the Rényi entropy, as described in the following).

Measuring TFDs' Information Content Using the Global Rényi Entropy
The entropy measure originates from physics as a measure of the disorder of thermodynamic systems. It has been widely applied in information theory as an information uncertainty estimator of probability density functions [71].
In signal processing, the time-frequency Rényi entropy has been adopted as an estimator of signal complexity. The idea of measuring signal complexity using the Rényi entropy is based on structural similarities of probability density functions and TFDs. Namely, the energy preservation property of the TFD can be expressed as [20] while the marginal conditions, similarly to the probability density functions, are calculated as [20]: Thus, the TFD properties in Equations (11)-(13) justify treating TFDs as probability density functions when quantifying their information content. Thus, the global Rényi entropy (estimated over the entire TFD plane) of the normalized TFD takes the form [71]: where α is the order of the Rényi entropy [71,72]. A multicomponent signal z M (t), embedded in additive white Gaussian noise ν(t), may be modeled as where z m (t) is the m-th signal component. When applied to TFDs, the Rényi entropy allows detection of the total number of signal components M using its counting property. The counting property of the Rényi entropy can be illustrated as follows: Consider a compactly supported signal z 0 (t) and its copy z c (t) (obtained by shifting z 0 (t) in the (t, f ) plane). The Rényi entropy of the TFD of the two-component signal z 0 (t) + z c (t) (H α,z 0 +z c ) carries exactly one bit of information more than the Rényi entropy of the one-component signal TFD (H α,z 0 ) [71].
Thus, the number of components of the multicomponent signal z M (t) can be determined as [71]: where z 1 (t) is the one-component reference signal.
Although the Rényi entropy was introduced for estimating the complexity of TFDs that satisfy marginal conditions and the energy preservation property, simulations have shown that TFDs that do not meet these conditions (such as spectrograms or the absolute value of a complex-valued RD) perform similarly to the probability density function in the case of the estimation by means of the Rényi entropy [73,74].

The Local or Short-Term Rényi Entropy (STRE)
The main shortcoming of the counting property of the global Rényi entropy estimation method is its applicability to only multicomponent signals constructed from shifted copies of one basic component. In fact, the global Rényi entropy estimation method fails in the case of multicomponent signals with components that have different time/frequency supports and frequency modulations.
In order to solve this issue, Saulig et al. proposed limiting the TFD of the analyzed signal with a short-time moving window (with duration ∆t) before evaluating the Rényi entropy [73,74]. Due to the use of a short time interval ∆t, the proposed complexity measure was named local or short-term Rényi entropy (STRE). The application of the local Rényi entropy to the signal TFD enables one to obtain the continuous time function of the number of components inside the moving time interval as where Even though, over short-term estimation intervals, time and frequency marginals are not preserved, a simplified model of a time slice of an idealized TFD was introduced in [73] to show that the counting property of the RE holds under the assumption of a positive TFD and cross-term suppression.
Based on Equation (18), an iterative algorithm for estimating the number of components in signals whose components have different amplitudes was proposed [73].
An example of estimation of the local number of signal components using the proposed iterative STRE-based method in comparison to the non-iterative algorithm presented in [74] is given in Figure 1. The signal's time series and TFD are shown in Figure 1a,b, followed by the instantaneous number of signal components obtained using the non-iterative ( Figure 1c) and iterative STRE0-based methods (Figure 1d,e for different α values). As can be seen, the iterative method unambiguously detects the number of signal components for each time instant, as well as their time supports. Furthermore, it can be observed that the proposed method also performs well in the case of components with different amplitudes, which are typical for the EEG time-frequency structures [74].
Information on the number of signal components, M p (t), obtained using the Rényi entropy allows localization and extraction of components from multicomponent signal TFDs, and then the estimation of their IFs, as described in the next subsection.

Instantaneous Frequency Estimation
Signal IF, defined as the derivative of the signal phase with respect to time, is one of the key signal parameters for providing important information concerning the timevarying spectral changes in non-stationary signals. The concept of IF found its usage in various technical fields and applications, such as seismology, machine condition monitoring, radar, sonar, communication (where it is referred to as the frequency modulation law), biomedicine, etc. [20]. Here, we apply IF estimation to biomedical signal characterization, i.e., to multichannel EEG signal analysis from their TFDs.
In the case of signal being mono-component, its IF is calculated as f 1 (t) = φ 1 (t)/2π [20]. Furthermore, the IF corresponds to an elongated energy region in the time-frequency representation of the signal and may be estimated by tracking the crest of the time-frequency "ridge", i.e., [20], A multicomponent signal z M (t) (as most of multichannel EEG signals are) can be modeled as a sum of M mono-component signals (the IF of each being f m (t)): where φ m (t) is the m-th signal component's instantaneous phase, and a m (t) is its instantaneous amplitude. The IFs can be calculated as (c) Local number of signal components M p (t) obtained using the non-iterative method [74]. (d) Local number of signal components M p (t) obtained using the iterative method with α = 3 [73]. (e) Local number of signal components M p (t) obtained using the iterative method with α = 5 (dotted) and α = 7 (dashed) [73].
In order to obtain the IFs from a multicomponent EEG signal, a component separation procedure should precede the IF estimation [20]. The procedure for localization and extraction of EEG signal components from its time-frequency representation is given in the following subsection.
Extracting signal components from a mixture of two or more statistically independent components is often referred to as blind source separation (BSS) [20,22]. The term "blind" signifies that neither the mixture structure nor source signals are known a priori. Instead, only their mixture is required (as is the case with recorded multichannel EEG signals) [20,22]. In addition, the proposed method is not limited to signals with simultaneously existing components, as is the case with the method in [75].
The challenging task of detecting and separating EEG signal components was performed in the time-frequency domain with the method described in the flowchart in Figure 2. As can be seen, the first step is the localization of the TFD maxima coordinates (t 0 , f 0 ) in the time-frequency plane, followed by the adaptive frequency bandwidth size calculation for the time slice C z (t 0 , f ).

Component Extraction Procedure
After all components are extracted from the signal's TFD, their IFs can be obtained from the TFD maxima, as reported in Section 2.3.
The described procedure was applied to all analyzed EEG signals, and the experimental results are presented in the next section.

Real data
The proposed method, which was based on the STRE and upgraded with the EEG component extraction and IF estimation, was implemented in Matlab (the Time Frequency Toolbox was used) and tested on 27,360 EEG signals recorded using 19 electrodes (placed at standard locations according to . The method was applied to reallife EEG signals from an open-access EEG database made available by the Institute for Neural Computation, University of California San Diego (Swartz Center for Computational Neuroscience) at https://sccn.ucsd.edu/~arno/fam2data/publicly_available_EEG_data. html (accessed on 5 January 2021) [76,77]. The analyzed EEG motor movement/imagined movement dataset consisted of EEG data captured for forward and backward movements of the left and right hand, as well as imagined forward and backward movements of the left and right hand, which were used as stimuli. Each movement was recorded independently (not as part of a sequence of movements) with an equal number of trails (with eyes closed and without controlling breathing or swallowing). The room was without electromagnetic shielding. EEG time series were captured by a Neurofax EEG System with the sampling frequency of 500 Hz and daisy-chain montage (records were exported using the Eemagine EEG) [76,77].
The noise was suppressed by applying hard amplitude thresholding (all TFD coefficients smaller than 5% of the TFD global maxima were set to zero in all TFDs). This is standard preprocessing that is found in many other time-frequency-based methods (for example [74,78]) based on the well-known fact that the white noise has a flat spectrum in the time-frequency domain [79,80]. The purpose of this step in our method was to enhance the algorithm's precision by reducing its sensitivity to noise and low-energy peaks in the time-frequency domain. In addition, the DC component and cross-terms caused by it were removed from the time-frequency representation prior to calculation of the number of EEG components and their IF estimations. Component localization, counting, extraction, and IF estimation were limited to the bandwidth of interest. Namely, the sub-bandwidth of 0.01-0.2 was considered instead of the entire normalized frequency range of 0-0.5, since it contains significant information on EEG spectral changes (knowing that most of the important EEG features are found in the range of 1-20 Hz [81]). EEG signals that were 5.12 s long were down-sampled to 512 samples in order to reduce computation burden (downsampling reduced sampling frequency from 500 to 100 Hz). However, down-sampling did not lead to loss of significant information, since the spectral bandwidth containing important EEG content remained unaltered. The Rényi entropy of order α = 3 was evaluated, since it was shown in [71,82] that odd entropy orders annul the effects of oscillatory interferences between the components. Selecting larger α values increases computational complexity without significant effects on the obtained results, as shown in Figure 1.

Results
Examples of the proposed method's performance when applied to the EEG spectrograms, WVDs, and RDs are given in Figures 3, 4 Figures 3i, 4i, and 5i and Figures 3m, 4m,  and 5m, respectively. Figures 3q, 4q, and 5q and Figures 3u, 4u, and 5u show the EEG time series for left-hand backward and imagined left-hand backward movement, respectively, while right-hand backward and imagined right-hand backward movement EEG time series are represented in Figures 3y, 4ac, and 5ac and Figures 3y, 4ac, and 5ac, respectively. The corresponding signal spectrograms, WVDs, and RDs are given in Figures 3b,f,j,n,r,v,z,ad, 4b,f,j,n,r,v,z,ad, and 5b,f,j,n,r,v,z,ad, respectively. As can be seen, EEG signals often have multiple components, resulting in multiple "ridges" in their time-frequency representations. The number of their components, as well as their locations in the time-frequency plane, may be used to detect limb movement signatures in recorded EEG signals, as shown in the following. Prior to their extraction, the number of signal components M p (t) as a function of time was obtained using the STRE (as described in Section 2.2) and is shown in Figure 3c,g,k,o,s,w,aa,ae for the spectrogram, in Figure 4c,g,k,o,s,w,aa,ae for the WVD, and in Figure 5c,g,k,o,s,w,aa,ae for the RD. Finally, the IFs of each EEG signal and each signal component were estimated from the extracted components using the approach described in Section 2.3. The dominant component IF, f 1 (t), was singled out for the purpose of the EEG signal characterization in the time-frequency domain (the dominant component was defined as the one with the highest amplitude in the (t, f ) plane). The estimated IFs are given in Figure 3d,h,l,p,t,x,ab,af for the spectrogram, in Figure 4c,g,k,o,s,w,aa,ae for the WVD, and in Figure 5c,g,k,o,s,w,aa,ae for the RD, respectively.
Tables 1, 2, and 3 show the results for the mean value of the dominant component IF f 1 for all analyzed EEG signals (captured by all electrodes) for the spectrogram, WVD, and RD, respectively. The first column in Tables 1-3 gives the electrode scalp locations, followed by two columns that provide the average of the dominant component IFf 1 estimated from the spectrogram, the WVD, and the RD for the left-hand forward and imagined left-hand forward movements, respectively. The following two columns in Tables 1-3 providef 1 results for the right-hand forward movement and imagined right-hand forward movement calculated from the spectrogram, the WVD, and the RD, respectively. The next two columns in Tables 1-3 give the results estimated from the spectrogram, the WVD, and the RD for the left-hand backward movement and imagined left-hand backward movement (the mean value of the dominant component IF), respectively. The last two columns in Tables 1-3 give the mean of the dominant IF for the right-hand backward and imagined right-hand backward movement, respectively, for the spectrogram, the WVD, and the RD. Simple statistical analyses, in terms of the mean, median, and standard deviation of the average of the dominant IF for all electrodes, as well as for all tested EEG signals and calculated TFDs, are given in the last three rows of Tables 1-3, respectively. In addition, a comparison of the number of components M andf 1 for the analyzed real (executed) limb movements' and imagined movements' spectrograms, WVDs, and RDs is presented in Figures 6 and 7, respectively.

EEG Analyzer Implemented in a Virtual Computer Instrument
The proposed method based on the STRE and IF estimation for limb movement EEG analysis was implemented in a virtual computer instrument (shown in Figure 8) as a groundwork for an extensive computer-aided decision-support clinical system. The instrument emulates the decision-making that is often performed by medical experts based on a visual inspection and interpretation of patterns found in EEG waveforms. This user-friendly virtual instrument for multi-channel EEG analysis implements a spectrogram (the time/frequency resolution of which is controlled by a corresponding window width), WVD, and RD as time-frequency representations. The functionalities of the proposed instrument are listed in the following.

•
Load multi-channel EEG signals: The user imports multi-channel EEG records as a MAT file (for real and imagined limb movements), which are then shown in the first row of the proposed instrument. Selecting TFD and windows' sizes and types: The user is allowed to choose a TFD from the provided list of TFDs, and depending on the chosen TFD, they are allowed to set the analyzing window width and type (Hamming, Hanning, rectangular, triangular, Gauss, and Kaiser). • STRE and TFD parameters: In addition, prior to running EEG analysis, the user is allowed to select values of the STRE sensitivity parameter, the TFD threshold for reducing noise and low-energy cross-terms, and the component extraction threshold. • TFD display options: The proposed instrument allows display of the TFDs as imagesc (TFD is displayed as an image), contour (TFD values are treated as heights above a plan), mesh (TFD is treated as colored parametric mesh), and surf (plots colored parametric surface). • References: By clicking on this button, related papers and previous works of the authors that led to the development of the proposed instrument are given.
The next section presents a detailed elaboration of the obtained results, specifying the proposed method's contributions and advantages in EEG signal characterization.

Discussion
The reported results show that EEG signals often have multiple components, with comparable complexity in terms of the number of signal components. Furthermore, the obtained results show that limb movements (both imagined and real movements) leave their signatures in EEG records in the time-frequency domain. The number of EEG components and the dominant component IF were used to capture these signatures.
More specifically, for the spectrogram, the number of EEG components in the (t, f ) plane varied from 2 to 21 for hand movements and from 1 to 20 for imagined hand movements. Extreme values were obtained for the left-hand forward and the right-hand backward movements at C3 when the maximal number of components was reached, and for the imagined right-hand backward movement at CZ when the minimal number of components was reached. Furthermore, important information was obtained from the mean of the normalized IF of the dominant EEG componentf 1 , varying from 0.05 to 0.19 (shown in Table 1). Namely, for hand movements,f 1 varied from 0.05 to 0.19, and for imagined hand movements, it varied from 0.05 to 0.18. This variation off 1 at specific electrodes was used to detect signatures of the analyzed limb movements in EEG records, as discussed in the following.
In terms of computational cost, the STRE was significantly more expensive than both component extraction and IF estimation. Namely, component extraction and IF estimation required an average time of only 0.41 s (in Matlab 2016b run on a computer with Intel i7-6700HQ CPU at 2.60 GHz with four cores and 16 GB RAM). On the other hand, calculating both the local and total number of components with the STRE required an average time of 86.99 s. In addition, calculating EEG time-frequency representations required an average time of 0.017 s. Hence, the overall computational cost of the proposed method is a limiting factor for its real-time application.
Comparable results in terms of the number of EEG components, and especially in terms of the dominant component IF, were also obtained from the WVD and RD. This proves that the proposed method is not limited only to the spectrogram, but it can also be used with other TFDs. The number of EEG components found in the (t, f ) domain for the tested limb movements varied from 6 to 39 in the case of the WVD and from 1 to 107 in the case of the RD. The reason for the increased number of components in the WVD and RD when compared to the spectrogram may be found in the quadratic nature of these time-frequency representations.
The proposed STRE-based approach was compared to a method used to extract useful information content (in terms of the number of energy clusters) from signals in the timefrequency domain [83]. This method applied the K-means algorithm for amplitude-based clustering of the TFD, with the optimal number of classes for the segmentation procedure being computed using the Davies-Bouldin criterion [84]. Following the TFD segmentation, the class containing the TFD elements with the largest amplitudes (i.e., useful information content) was used as the input data for the Hoshen-Kopelman algorithm [86]. The Hoshen-Kopelman algorithm has been widely used to estimate the number of separate energy clusters in the time-frequency domain, which, in this case, provides a quantitative measure of EEG complexity (shown in Table 4). However, since this method does not provide information on the local EEG complexity, at each time instant, as required for component extraction and IF estimation procedures, the STRE-based approach was favored in the proposed method.
The results (in terms of the number of EEG components) obtained using the STREbased method describe the variable complexity of EEG signals with respect to particular channels. Namely, distinguishing signatures of the analyzed limb movements using only the number of EEG components was shown to be insufficiently effective, as can be seen in Figure 6. Hence, we proposed an upgrade of the method with the IF estimation. In order to estimate EEG components' IFs, component extraction is required beforehand. As shown in Section 2.4, the described component extraction approach requires the information on the number of EEG components locally, at each time instant, which was previously ensured by the STRE-based method (also used to detect the total number of EEG components).
The IFs of the dominant EEG components estimated from the spectrogram, WVD, and RD provide crucial information that is necessary in order to detect limb movement signatures in EEG signals (as shown in Tables 2 and 3). Namely,f 1 varied from 0.05 to 0.16 for the WVD and from 0.05 to 0.17 for the RD. These results are similar to the results obtained using the spectrogram (whenf 1 varied from 0.05 to 0. 19), proving that the dominant EEG component IF is highly robust to the chosen TFD.
Next, the dominant EEG component IF was used to detect and distinguish between various limb movement signatures found in EEG time-frequency distributions. Moreover, understanding and distinguishing the signatures of real and imagined limb movements in EEG signals has great usability in computer-assisted decision-support systems for analyzing neurological motor control disorders. As can be seen in Figure 7, high frequencies in the dominant EEG component IF were detected at F3, F4, F7, F8, T3, and T4. Furthermore, considering the dominant EEG component IF at F7, F8, and T4, we were able to distinguish the analyzed hand movement signatures from imagined movements.
In fact, the signature of the imagined left-hand forward movement was detected in its EEG record as a high frequency of the IF of the dominant EEG component at F7 for all three analyzed TFDs. Furthermore, the imagined left-hand forward movement resulted in a largerf 1 value at F7 when compared to the corresponding executed hand movement for all three time-frequency representations calculated.
Furthermore, the signatures of the right-hand forward and imagined right-hand forward movements were found in their EEG signals as high values off 1 at F8 and F7, respectively. The dominant EEG component IF for the imagined right-hand forward movement was larger thanf 1 for the respective executed hand movement at F7, while the executed right-hand forward movement increasedf 1 when compared to the respective imagined movement at F8 (for all three TFDs).
Next, the signature of the left-hand backward movement was detected as a high frequencyf 1 at T4 (with an increased dominant EEG component IF when compared to the imagined left-hand backward movement for the spectrogram, WVD, and RD).
Finally, the signatures of the right-hand backward and imagined right-hand backward movements (for all three TFDs) were detected in the EEG as high frequencies of the dominant EEG component IF at F8 and T4 (with dominant EEG component IFs for the moved hand being larger at T4 than for the respective imagined hand movements), as well as at F7 (with largerf 1 for the moved hand than for the respective imagined hand movement), respectively.
Note that the obtained results were stable for all three analyzed TFDs. Moreover, the proposed method was shown to efficiently distinguish between real movements from the corresponding imaginary movements, as well as between backward and forward hand movements for both the left and right hand. More specifically, the provided study suggests that limb movement signals leave their signatures in EEG records, which may be detected by analyzing the dominant EEG component IF in the TFD at specific channels (namely, F3, F4, and T3, but especially F7, F8, and T4) in the time-frequency domain.
Moreover, the presented results suggest that the proposed time-frequency-based approach using the STRE and IF estimation has potential to detect signatures of various neurological phenomena, and possibly signatures of some neurological abnormalities. However, a specific neurological disorder would require extensive study, which should include experts from both the engineering and medical fields (the clinical applications exceed the scope of this paper and are planned for our future research).
Hence, the results of this study in terms of detecting signatures of various motor control EEGs in the time-frequency domain encourage hopes that the proposed method may be used in the future to also detect signatures of numerous neurological disorders in EEG time-frequency representations.

Conclusions
Detecting the signatures of various stimuli that induce brain activities in recorded EEG signals is one of the most challenging tasks of EEG analysis. For that purpose, this paper proposes a novel time-frequency-based method for multichannel EEG signal analysis as groundwork for building an extensive computer-aided decision-support system to correlate various EEG patterns with functions and dysfunctions of the central nervous system. The method introduces the modification of the Rényi entropy estimation method, upgraded with a blind component separation procedure (used for EEG component extraction) and component IF estimation (used for EEG characterization). As shown in the paper, the numbers of EEG components and, especially, their dominant IFs provide useful information that may help to detect signatures of various limb movements in EEG time-frequency representations. By considering the EEG channels F3, F4, F7, F8, T3, and T4 (with an emphasis on F7, F8, and T4), the method was able to detect signatures and distinguish between forward and backward movements of the left and right hand, as well as imagined hand movements. Specifically, the imagined left-hand forward movement was detected with the high frequency of the dominant EEG component IF at F7. High values of the dominant EEG components' IFs were found to be signatures of the imagined right-hand forward and right-hand forward movements at F7 and F8, respectively. The left-hand backward movement was detected with the high dominant EEG component IF at T4. The signatures of the right-hand backward and imagined right-hand backward movements were detected with the high frequencies of the dominant EEG components' IFs at F8 and T4 and at F7, respectively.
Furthermore, the method's performance in terms of accuracy, precision, recall, and F1score in distinguishing real from imagined hand movements was reported for all analyzed TFDs in the case of low-noise measurements. The achieved results show potential for further enhancement of clinical diagnostic methods and medical treatment of numerous neurological disorders if the signatures of these abnormalities, as expected, are detected by the proposed method in recorded EEG data, which is planned to be investigated in future work.