Method for Automatic Estimation of Instantaneous Frequency and Group Delay in Time–Frequency Distributions with Application in EEG Seizure Signals Analysis

Instantaneous frequency (IF) is commonly used in the analysis of electroencephalogram (EEG) signals to detect oscillatory-type seizures. However, IF cannot be used to analyze seizures that appear as spikes. In this paper, we present a novel method for the automatic estimation of IF and group delay (GD) in order to detect seizures with both spike and oscillatory characteristics. Unlike previous methods that use IF alone, the proposed method utilizes information obtained from localized Rényi entropies (LREs) to generate a binary map that automatically identifies regions requiring a different estimation strategy. The method combines IF estimation algorithms for multicomponent signals with time and frequency support information to improve signal ridge estimation in the time–frequency distribution (TFD). Our experimental results indicate the superiority of the proposed combined IF and GD estimation approach over the IF estimation alone, without requiring any prior knowledge about the input signal. The LRE-based mean squared error and mean absolute error metrics showed improvements of up to 95.70% and 86.79%, respectively, for synthetic signals and up to 46.45% and 36.61% for real-life EEG seizure signals.


Introduction
Electroencephalogram (EEG) recordings are widely used for assessing brain disorders [1][2][3]. EEG is a noninvasive approach for detecting and predicting seizures [4][5][6][7], which can be difficult to identify in infants. Recurrent seizures are the hallmark of epilepsy, one of the most prevalent neurological disorders in humans. Identifying seizures in EEG recordings typically requires real-time observation by a neurologist, leading to growing interest in automated methods for seizure detection.
EEG signals are nonstationary, and therefore, time-frequency (TF) and time-scale representations are commonly used for their analysis [6,[8][9][10][11][12][13][14][15][16][17][18][19][20][21]. A time-frequency distribution (TFD) enables us to describe signal energy simultaneously in time and frequency. However, the most widely used TFDs, the quadratic TFDs (QTFDs), create highly oscillatory artifacts, known as cross-terms, for signals with several components or at least one nonlinear frequency-modulated (non-LFM) component [22][23][24]. Although the 2D low-pass filters in the ambiguity function (AF) domain are often used to suppress cross-terms, they may also suppress important components known as autoterms, resulting in a trade-off between cross-term reduction and autoterm resolution. To overcome this trade-off, a variety of filtering methods in the TF domain have been developed [22,25] as an alternative to conventional filtering methods utilizing nonstationary signals in the time domain [26][27][28]. One of these advanced filtering methods that outperforms others is the adaptive directional TFD (ADTFD), which achieves high resolution for multicomponent signals having multiple directions of energy distribution in the TF domain, such as EEG seizure signals [7,[29][30][31][32].
Instantaneous frequency (IF) provides vital information on the time-varying spectral changes in nonstationary signals. In the TF signal analysis, a number of IF estimation methods have been developed [29,[33][34][35][36][37][38][39][40][41][42][43][44]. However, IF is not as useful to characterize signals composed of spikes, where an infinite amount of frequencies is present [45][46][47]. Such behavior is present in EEG seizure signals, for which the group delay (GD) has been proven to be a better alternative to IF [46,47]. However, the limitation of the IF estimation method proposed in [48], and hence, the approach in [46,47], is that the number of components is required as input from a user, which limits its usage when no a priori information is available and when the number of components varies over time.
In this paper, we present a novel method for generating a binary map that automatically identifies distinct regions within the signal's TFD, where either the IF or GD measure is particularly suitable for analysis. Our method utilizes information about the local number of signal components obtained from the localized Rényi entropy (LRE) information, namely the short-term Rényi entropy (STRE) [33,34] and the narrow-band Rényi entropy (NBRE) [35,36], to detect components and segment the TFD accordingly. Unlike other methods that require a priori knowledge of the number of components, our approach is applicable to signals with unknown and time-varying numbers of components.
The resulting binary map is used in a joint IF and GD estimation method, which was applied to two commonly used IF estimation algorithms in EEG signal analysis, namely image-based [22,37,38] and blind-source separation (BSS) [15,39,40]. These methods were selected for their suitability in an automatic environment, as they effectively use the component time support information from the STRE. Our approach further incorporates the frequency support information from the NBRE to estimate the GD and reduce the dependency of the IF algorithms on the STRE accuracy, which can be achieved for the considered signal examples.
Furthermore, we demonstrate that the shrinkage operator proposed in [35,36] for sparse TFD reconstruction can be effectively used for IF and GD estimation, with competitive performance compare with the image-based and BSS algorithms. We evaluate the performance of our method using mean squared error (MSE) between the local number of signal components before and after estimation, and we show the superior performance of our combined IF and GD estimation approach over IF estimation alone on both synthetic signals with additive noise and real-life EEG seizure signals.
The rest of this paper is organized as follows. Background theory, the proposed method and the EEG dataset used in this study are described in Section 2. The obtained results are thoroughly presented and discussed in Section 3. Finally, the paper's conclusions are summarized in Section 4.

Time-Frequency Signal Analysis
A multicomponent nonstationary signal, denoted as z(t), is defined as the analytic associate of a real signal s(t) as where NC is the number of components, while a i (t) and ϕ i (t) denote the instantaneous amplitude and instantaneous phase of the signals i-th component, respectively. The ideal TFD,ρ(t, f ), is a unit delta function following the crests of the ridges which represent the IF, f 0 i (t), of the i-th component: indicating the dominant frequency of the signal's i-th component at a given time. A dual (or inverse) of the IF, namely the GD, τ d i ( f ), indicates the dominant time of the signal's i-th component at a given frequency: where Z( f ) is the Fourier transform of z(t). The definitions of IF and GD are closely related, involving interchanging the time and frequency variables, with an extra minus sign in Equation (4). In practice, the IF is obtained by determining the signal's component ridge across time slices of its TFD, while the GD is determined by determining the signal's component ridge across frequency slices of the TFD. In the vast majority of instances, the ideal TFD is not achievable, because practical TFDs are not precisely localized and may be affected by cross-terms [22]. The Wigner-Ville Distribution (WVD) is widely used as the most fundamental TFD, defined as [22]: and it provides an estimate of the IF/GD for a signal with a single LFM component in the TF plane that is almost perfect. Yet, the cross-term vulnerability (when dealing with multicomponent signals) necessitates proper cross-terms suppression. Using the AF, A(ν, τ), calculated as the highly oscillatory cross-terms can be suppressed with a 2D low-pass filter, defining a QTFD class of TFDs, ρ(t, f ): where g(ν, τ) is the low-pass filter kernel in the AF. The conventional approaches to kernel design typically entail a compromise between the concentration of autoterms and the suppression of cross-terms [22].

Adaptive Directional TFD
To circumvent the above problem of the conventional TFDs, the adaptive directional TFD (ADTFD) that adjusts the direction of the smoothing kernel at each TF point in the TF plane is introduced [7,29], and it is mathematically expressed as follows: where γ θ (t, f ) is the smoothing kernel whose direction is controlled by θ, while the double asterisk denotes double convolution in t and f . In this work, we used the Extended Modified B Distribution (EMBD) as an underlying QTFD with its kernel: where α E = β E = 0.25 are the time and frequency smoothing parameters [7,29,31,32]. As γ θ (t, f ), we selected the double-derivative directional Gaussian filter (DGF) as in [7,29,31,32]: where t θ = t cos(θ) + f sin(θ) and f θ = −t sin(θ) + f cos(θ), while parameters a and b control the extent of smoothing along the time and frequency axes. The DGF has low-pass characteristics along the time axis (e −a 2 t 2 θ ), while it performs second-order differentiation along the frequency axis ab . The direction angle of γ θ (t, f ) is adapted locally for each point in the TF domain by maximizing the correlation between the γ θ (t, f ) and TF ridges as where −π/2 ≤ θ ≤ π/2. The implementation of directional smoothing results in the suppression of cross-terms and the enhancement of autoterms. The optimization of smoothing kernel parameters and shape is necessary for achieving optimal performance, as they are dependent on the signal being analyzed. Previous studies [7,29,31,32] have indicated that assigning a small value to parameter a results in intensive smoothing along the major axis, whereas a larger value for parameter b prevents the merging of close components.
To be more precise, a ∈ [2, 3], while b ∈ [5,30]. In addition to shape parameters a and b, the window length, W L, of the γ θ (t, f ) affects the performance of the ADTFD. A filter with a small W L value fails to resolve close components and eliminate cross-terms, but it preserves the energy of short-duration components. Conversely, a larger W L value achieves the opposite effect. The computational demand of the exhaustive search involving all possible combinations of (a, b, W L) is significant. Therefore, we used a method for the automatic parameter optimization of ADTFD, namely the locally adaptive-ADTFD (LO-ADTFD) proposed in [31], where the final LO-ADTFD is obtained by choosing TF points with the minimum value from a given set of ADTFDs {ρ (ad) 1 (t, f ), ρ (ad) 2 (t, f ), . . . } and their respective parameters {(a 1 , b 1 , W L 1 ), (a 2 , b 2 , W L 2 ), . . . }: where ρ (ad) k (t, f ) is the k-th ADTFD in the defined set. That way, the LO-ADTFD preserves the energy of short-duration signals while achieving high-resolution TFD with resolved close components and suppressed cross-terms. In this work, we selected the parameter (a, b) values from the following set {(3, 6), (3,8), (2,20), (2, 30)}, while W L was optimized for each (a, b) pair using the concentration measure proposed in [41]:

The Localized Rényi Entropy
The Rényi entropy, denoted by R(ρ(t, f )), is a comprehensive metric for signal complexity in the TF plane [42][43][44], defined as where for the odd integer parameter α R > 2 the cross-terms get integrated out from the QTFD, ρ(t, f ), which is normalized with respect to its total energy [34,42].
To refine the global approach mentioned earlier, the counting attribute of the Rényi entropy has been utilized to extract the local number of signal components from ρ(t, f ) using STRE [34]. To achieve this, the Rényi entropy of the extracted signal's TFD is compared with that of a reference TFD with a known number of component as follows: where t denotes localization using time slices, and t 0 is the observed time slice, while ρ(t, f ) and ρ ref (t, f ) denote the considered and reference TFD, respectively. The time-localization operator Γ t 0 sets all TFD samples to zero, except those in the vicinity of t 0 : where Θ t is the parameter controlling the time-window length. The reference signal is a cosine signal with an amplitude of 1 and a constant normalized frequency of 0.1, providing a reference energy of a single component in each time slice [34].
In [35,36], more research reveals the weaknesses of STRE for specific signal types and introduces NBRE to counteract them. Using NBRE, one may determine the local number of signal components per frequency slice, f 0 , by substituting the frequency-localization operator for the time-localization parameter in Equation (16) with the frequency-localization operator: where Θ f is the frequency window length. The reference signal is a delta function centered at t = 15 [36]. Note that ρ(t, f ) and ρ ref (t, f ) have to be obtained with the same TFD with reduced interference in order for the comparison to be valid. The comparison between local numbers of components obtained from STRE and NBRE are shown in Figure 1 for the synthetic signal z 4LFM , with N t = 256 samples composed of four LFM components with different amplitudes embedded in additive Gaussian noise (AWGN) with a signal-to-noise ratio (SNR) SNR = 3 dB. Figure 1b,c show the reasoning behind introducing NBRE in [35,36]-an inaccurate increase in the local number of signal components, NC t (t), is evident for signals whose components are more aligned with the frequency axis (i.e., deviate from the method's reference component), such as the signal z 4LFM , whose LO-ADTFD is shown in Figure 1a. On the other hand, when signals have components more aligned with the time axis, STRE presents more accurate estimations of the local number of signal components [35,36].

Image-Based Method
The first method for estimating the IF is an image processing method [37], which is divided into two sequential steps. In the first step, local peaks of the TFD are detected using its first and second derivatives as which generates a binary (t, f ) image consisting of ones on all peak locations and zeros on all other points. This step usually provides peaks that do not belong exclusively to the autoterms. Hence, in the second step, the IFs of signal components are extracted by using the m-connectivity criterion derived from image processing. In this work, we used m = 10, as in [7,22,37], which defines a 10-neighborhood set for a detected peak at location (x, y) as {(x − 1, y), (x − 1, y − 1), (x − 1, y + 1), (x − 1, y + 2), (x − 1, y − 2), (x + 1, y), (x + 1, y − 1), (x + 1, y + 1), (x + 1, y + 2), (x + 1, y − 2)}. According to this criterion, the points above and below the IF curves are not included, meaning that there can be only one frequency at any time instant for any given signal component. Finally, a threshold must be specified for the minimum time duration of a valid signal component so that the final TFD contains only components that meet the threshold. In our example, where no prior knowledge of the input signal is available, the threshold is determined by the least component time support information provided by the STRE. This method demonstrated good computing efficiency and performance for real-world signals without requiring previous knowledge of the components' IF laws and amplitudes [7,22,37,38].

Blind-Source Separation Method
The second IF estimation method applied in this study is the blind-source separation (BSS) method, which is an efficient method for the localization and extraction of components from multicomponent signals in the TF domain [39]. The term "blind" refers to the provision of a mixture of statistically independent components without prior knowledge of its structure or number of components. The STRE's information on the time supports of components is included in the version of the approach utilized in this study, namely the BSS-STRE [40], thus removing the need for several thresholds required by the original method [39]. The steps of the BSS component extraction method are summarized below.
First, the TFD of a signal is computed, ρ(t, f ), and the corresponding NC t (t) from the STRE is obtained. The algorithm then locates the largest TFD peak at (t 0 , f 0 ) and calculates the adaptive neighboring component frequency band, ∆B = B L + B R , for the time slice ρ(t 0 , f 0 − B L : f 0 + B R ). Next, the component is extracted at time-slice t 0 , and NC t (t 0 ) is reduced by 1. Following that, the previous steps of the method are repeated in both directions around t 0 as t 0 ← t 0 − 1 and t 0 ← t 0 + 1 until the component edges are detected by the first derivation of |NC t (t)|, |NC t (t)| = 0. If there is at least one component remaining in ρ(t, f ), the above steps are repeated for the succeeding component. The full pseudocode of this method may be found in [15,40].
The BSS method results in several TFDs, each containing a single extracted component from which the component IF is estimated in a separate vector. In order to be comparable with the B (IM) (t, f ), all estimated IFs will be displayed within a single binary TFD, denoted by B (BSS) (t, f ).
Note that the ADTFD, or more precisely, its automatically optimized version LO-ADTFD, is used as the underlying TFD in both estimation methods for this work, as its superior performance has been demonstrated in numerous IF estimation applications studies [38,45,46].

Limitations of the Considered IF Estimation Algorithms Based on STRE
The accuracy of IF estimation methods is highly dependent on the accuracy of STRE. The findings in [35,36] indicate that STRE is unsuitable for signals whose components deviate from the time axis, resulting in artificially increased NC t (t) and reduced estimation accuracy. Since this issue may occur in real-life signals, the question arises as to which IF estimation approach should be used for such signals. An incorrect NC t (t) in the imagebased IF estimation method can lead to either too low a threshold, resulting in interference being classified as a true signal component, or too high a threshold, causing some true signal components to be rejected. In the case of the BSS method, an incorrect NC t (t) can result in incomplete extraction of signal components, and higher NC t (t) values may cause the estimation of interference IFs.
In this paper, we aim to demonstrate the negative consequences of using an inappropriate localization approach in estimating the IFs of signal components. Specifically, we show that such an approach can cause estimated IFs to be discontinuous and shifted away from the true component ridge. Although polynomial functions can approximate discontinuous estimation samples, the approximation error increases as more estimated samples do not belong to the autoterms. This phenomenon is particularly problematic for the image-based method, which requires a larger m-connectivity criterion to connect true signal components. Otherwise, discontinuous signal components are often classified as interference and fail to meet the threshold criteria. Increasing the m-connectivity criterion is not recommended, as it can lead to the linking of interference terms. Our paper aims to address these limitations by proposing the use of the frequency localization approach and the estimation of the GD with information from the NBRE method for certain signals. By doing so, we can minimize the negative effects of an inappropriate localization approach on IF estimation.

The Proposed Rényi-Entropy-Based Method for Component Alignment Detection towards Time or Frequency Axis
In this section, we present a novel method for the automatic detection of TFD regions that require a time or frequency localization strategy. Such a strategy involves using TFD analysis with time or frequency slices, linked to an estimate of IFs or GDs. To detect such regions, we implemented the STRE and NBRE methods, which are sensitive to estimate errors when dealing with signal components that vary from the methods' respective reference components. In this method, this sensitivity is turned into an advantage, as a considerable increase in the local number of signal components can indicate the need for an alternate localization strategy in a certain TFD region.
Next, the proposed method validates the detected increase in the local number of signal components by comparing the quality of IF and GD estimations in the identified TFD region. We observed that an incorrect localization approach can lead to discontinuous and inaccurate IF or GD estimations.
To estimate the IF and GD trajectories, we employed an approach based on the shrinkage operator proposed in [35,36] for sparse TFD reconstruction. The shrinkage operator, denoted with shrink t, f , operates independently for each time and frequency slice and removes samples that do not belong to the autoterms. The autoterms are locally associated with the NC t (t) or NC f ( f ) largest areas, where an area is calculated as a sum of samples between the minima to the left and right of the detected local maxima. This operator involves parameters δ t and δ f , which control the number of samples around local maxima classified as autoterms [35,36]. Hence, by applying the shrinkage operator on desired TFD with parameters δ t = δ f = 1, we may extract only local maxima belonging to signal autoterms, which basically represent the IF and GD estimations of signal components: where the t, f notation denotes shrinkage performed over time or frequency slices, while ρ t (t, f ) and ρ f (t, f ) denote the signal's estimated IFs and GDs (or ridges), respectively. To demonstrate the proposed method's steps, we created a synthetic signal, z mix (t), with N t = 256 samples, consisting of two constant FM and four non-LFM components with various different directions and time/frequency supports. Figure 2 illustrates the LO-ADTFD of the signal, as well as its estimated IFs, ρ t (lo) (t, f ), and GDs, ρ Figure 2b displays four discontinuous non-LFM components in ρ t (lo) (t, f ) whose alignment deviates from the time axis. Alternatively, Figure 2c depicts the same occurrence for ρ f (lo) (t, f ) but with inverse results. Now, the identical four non-LFM components exhibit stronger connectivity than two constant FM components that are aligned with the time axis.
To evaluate the quality of estimated IFs and GDs, we propose a metric based on the number of continuously connected regions of TFD samples denoted by N r . Specifically, a TFD sample at location (x, y) is considered to be part of a region if it is connected to at least one sample in its 8-neighborhood set: The resulting components are then counted to obtain N r . This way, the proposed metric detects and penalizes discontinuities in the estimated IF or GD trajectories. Higher values of N r indicate lower consistency among the estimated components (i.e., autoterms), reflecting a lower quality of the estimated IF or GD trajectories. In our proposed method, the result of the N r metric serves as a decision-making factor in each TFD region under consideration.
The objective of our proposed method is to generate a binary component alignment map, denoted by BM(t, f ), that distinguishes TFD regions with components suitable for time or frequency localization using ones and zeros, respectively. We provide a comprehensive description of the proposed method's steps below: 1.
The first step in the proposed approach involves calculating the TFD for a given signal in the time domain, followed by estimation of the IFs and GDs using the shrinkage operator, resulting in ρ t (t, f ) and ρ f (t, f ), as shown in Figure 2 , the proposed method assumes that signal components are primarily aligned with the time axis and generates BM(t, f ) using ones. Additionally, the STRE method is employed to calculate NC t (t), which is used to investigate local component behavior.
, the proposed method generates BM(t, f ) using zeros and uses the NBRE with NC f ( f ) as the initial localization approach. 3.
The proposed algorithm examines the input NC t (t) or NC f ( f ) for pronounced local maxima, which may indicate an inadequate local component for the current STRE or NBRE approach. To identify such maxima, we first locate all local maxima within NC t (t) (or NC f ( f )), followed by the calculation of the difference in the local number of signal components, denoted as ∆NC, between the observed maximum and the minima to the left and right. We consider all ∆NC ≥ 1.50 as "suspicious" intervals that require further analysis. The chosen threshold value of 1.50 is based on the desire to detect components that deviate from the respective time or frequency axis more than the LFM component, with a starting and stopping normalized frequency at 0 and 0.5, respectively. This component showed marginal accuracy for both approaches, with mean numbers of local components obtained as 1 48. If all ∆NC values are less than 1.50, the algorithm outputs a BM(t, f ) consisting of only ones or zeros and terminates. Otherwise, the algorithm proceeds to the next step for further analysis. 4.
Next, the algorithm identifies a segment of time (or frequency) slices from NC t (t) (or NC f ( f )), where the edges are defined by local minima satisfying ∆NC ≥ 1.50. An example of NC f ( f ) for the signal z mix (t) is shown in Figure 3a, where the first segment is indicated by red dashed lines at frequency bins f 1 and f 2 . The same segment is then extracted from a TFD, within which an unsuitable signal component for the current localization approach may be present. Figure 3b illustrates an example of ρ (lo) (t, f ) with the segment constrained by the previously detected frequency bins f 1 and f 2 , where a constant FM component needs to be further detected in the subsequent algorithm steps as unsuitable for the frequency localization approach. 5.
At this point, additional localization is performed within the segmented TFD by computing the LRE in the opposite direction of the previous step. In particular, . We follow the same procedure as in the previous two steps, identifying all local maxima and minima with ∆NC ≥ 1.50. Then, we detect all segments with borders satisfying ∆NC ≥ 1.50, which define a 2D TF region or block within ρ t (t, f ) and ρ t (t, f ) that are evaluated in reference to the N r value. If N r is lower for a TFD block within ρ t (t, f ) than in ρ f (t, f ), it implies that a time localization technique is more suited, and BM(t, f ) in positions specified with a TFD block is set to 1; otherwise, it is set to 0. All the remaining estimated signal components that are saturated inside minima with ∆NC < 1.50 belong to the current localization approach (which differs from the localization approach in the previous step), and the corresponding TFD block within the final BM(t, f ) changes its values from 1 ↔ 0. Finally, the N r value is compared in ρ t (t, f ) and ρ f (t, f ) for time samples or frequency bins where no component was detected (i.e., NC f ( f ) = 0 or NC t (t) = 0), and the BM(t, f ) is set to 0 or 1 based on the lower N r value. Figure 4 illustrates the procedure described in this section. Firstly, an LO-ADTFD segment bounded in the range [ f 1 , f 2 ] is extracted and subjected to an opposite LRE approach, such as STRE, to obtain the NC t (t), as shown in Figure 4a,b. The red dashed lines in Figure 4b denote the segment detected using the NC t (t), which defines a TFD block (t 1 : t 2 , f 1 : f 2 ) for evaluation in ρ t (lo) (t, f ) and ρ f (lo) (t, f ), shown in Figure 4c,d, respectively. The estimation of GDs within TFD blocks in red dashed lines exhibits superior connectivity than IF estimation, leading to the BM(t 1 : t 2 , f 1 : f 2 ) to remain unchanged from the initial values of zero. Since the constant FM component detected inside the green dashed lines does not produce inaccuracies in the NC t (t), as shown in Figure 4b, the BM(t 3 : t 4 , f 1 : f 2 ) is changed to one. 6.
Steps 4 and 5 are repeated for the remaining detected segments in the NC t (t) (or NC f ( f )) input, i.e., until all ∆NC ≥ 1.50 are examined.

Component Extraction Using the Component Alignment Map
The BM(t, f ) map obtained provides the means to extract signal components in the TFD. Extraction of these components is not achieved individually but through two sets of components: one suitable for localization through time slices and the other suitable for localization through frequency bins. The extraction process is simplified by multiplying the BM(t, f ) map with the TFD using the operator κ{·}, which we defined in two ways. To derive the signal components corresponding to the localized approach through time slices, the operator κ t {·} is used: where, by multiplying the BM(t, f ) map with TFD, only the regions defined by the units in BM(t, f ) are retained in TFD. Likewise, the operator κ f {·} is employed to obtain the components corresponding to the localized approach through frequency bins: by means of which only the TFD regions that are defined by zeros in BM(t, f ) are kept in TFD. Consequently, when max{BM(t, f )} = 1 and min{BM(t, f )} = 0, the input TFD can be split into two TFDs, κ t {ρ(t, f )} and κ f {ρ(t, f )}, enabling the local number of signal components to be computed using the STRE and NBRE methods, respectively. Since both κ t {ρ(t, f )} and κ f {ρ(t, f )} are expected to contain components corresponding to the chosen localization approach, more precise estimates of the local number of signal components can be obtained compared with those obtained from the original TFD. Figure 6 illustrates the results of applying the operators κ t and κ f to the signal z mix (t). The obtained κ t {ρ (lo) (t, f )} and κ f {ρ (lo) (t, f )}, shown in Figure 6a,b, respectively, contain components suitable for the chosen localization approach. This is further corroborated by the NC t (t) and NC f ( f ) estimates in κ t {ρ (lo) (t, f )} and κ f {ρ (lo) (t, f )}, respectively, which lack significant inaccurate local maxima. Notably, the NC t (t) and NC f ( f ) estimates obtained from the split TFDs are considerably more precise than those obtained from the original TFD, as demonstrated in Figure 6c,d. However, it is important to emphasize that the estimates of the local number of signal components should be interpreted in conjunction with BM(t, f ), since they do not represent estimates of the entire TFD but rather of the TFD regions specified within BM(t, f ).

Method for an Automatic Estimation of IF and GD
Upon completion of the necessary prerequisites, we propose a new method that can automatically estimate both the IF and GD of signal components in a TFD. This method leverages the use of the binary map BM(t, f ) to identify the TFD regions that require IF or GD estimation. Meanwhile, the behavior of the signal components, or autoterms, is defined by the local number of signal components obtained through the STRE and NBRE methods applied on TFDs with extracted components, κ t {ρ(t, f )} and κ f {ρ(t, f )}, respectively.
The proposed method is composed of the following steps: 1. First, the input signal's TFD is processed using the κ t, f {·} operator to obtain two TFDs: one TFD, κ t {ρ(t, f )}, containing signal components suitable for IF estimation using a time-slice approach, and another TFD, κ f {ρ(t, f )}, containing signal components suitable for GD estimation using a frequency-slice approach. 2.
Next, the IFs are estimated from κ t {ρ(t, f )} using any IF estimation algorithm, the result of which is denoted as B t (t, f ).

3.
A matrix transpose is applied to the discrete version of κ f {ρ(t, f )}, which interchanges the time and frequency axes and enables GD estimation in frequency slices using the same IF estimation algorithm approach, resulting in B f (t, f ).

4.
Finally, the IF and GD estimations are combined within a resulting binary TFD, B(t, f ), which is an output of the proposed method summing B t (t, f ) and B f (t, f ). The proposed method's steps are visualized in Figure 7. We implemented the aboveproposed method in the considered IF estimation methods, which may be now considered as IF/GD estimation algorithms utilizing both LRE methods, namely the image-based STRE-NBRE and BSS-STRE-NBRE.
Apply IF estimation algorithm: It is worth noting that the combined estimates of the IF and GD may be extracted from the estimates generated by the operator shrink t, f {·}, previously used in the proposed method for BM(t, f ) as In the following section, the performance of the B (shrink) (t, f ) estimation is compared with that of the image-based and BSS methods.

EEG Dataset Description
Seizure signals in EEG recordings are often modeled as multicomponent piecewise FM signals: However , this model does not take into consideration the spikes or short-duration transients that are regularly seen in EEG readings. In order to account for these spikes, the updated signal model that we use in this study is mathematically given as follows [7,30]: where T i is a time-varying shift. Analyzing signals that contain both rhythmic and spike features using traditional TF techniques is challenging, as they have energy distributed along both the time and frequency axes. Smoothing along the frequency axis can eliminate cross-terms created by spikes, but it also reduces the resolution of the sinusoidal signal components. In addition, such signals require the use of a combined IF and GD estimation method, as estimating only the IFs fails to recover spike features [47,48].
We used a database of 200 EEG seizure segments, which were previously uploaded as supplementary material in [38] and have been used in [7,30,32,47], from which an illustrative example was chosen, as used in additional studies [22,31]. The data and relevant code are publicly available at https://github.com/nabeelalikhan1/EEG-Classification-IFand-GD-features (accessed on 1 October 2022). The EEG recordings were obtained from newborns at the NICU of the Royal Brisbane and Women's Hospital, Brisbane, Australia, using the MEDELEC Profile System. Twelve electrodes were placed according to the international 10-20 standard, which were used to construct a 20-channel bipolar montage. The recordings underwent prefiltering using an analog bandpass filter with a bandwidth of 0.5 to 70 Hz. The signal was then sampled to 256 Hz before being digitally resampled to 32 Hz, as the majority of the signal energy is typically found below 12 Hz. The resulting signal segment is 8 s in duration and acquired at a sampling rate of 32 Hz, resulting in a total of N t = 256 samples [7,30,32,38,47]. Previous studies have shown that a differentiator filter can be used to whiten the EEG background and enhance the signature of spikes in EEG signals [30,47,49,50]. Hence, the proposed IF/GD estimation method's performance was tested on EEG seizure signal without, denoted as z EEG (t), and with a differentiatior filter, denoted as z EEG filt (t).

Results and Discussion
We compared the performance of the combined IF and GD estimation approach with that of the IF estimation approach alone for both synthetic signals, z LFM (t) and z mix (t), and real-life EEG seizure signals, z EEG (t) and z EEG filt (t). It is important to note that for IF estimation only, algorithms use the STRE method applied to an input LO-ADTFD. However, for IF and GD estimation, the algorithms utilize the proposed method with BM(t, f ), along with the STRE and NBRE methods applied to the extracted components from an input LO-ADTFD using the proposed operators κ t and κ f . We calculated the STRE and NBRE using the parameter α R = 3, with Θ t = Θ f = 11 for the synthetic signals z LFM (t) and z mix (t) and Θ t = Θ f = 5 for the signals z EEG (t) and z EEG filt (t), in order to capture spike features more precisely, which have been shown to be stable in [33,35,36,51].
To supplement the visual inspection of the results, it is important to quantify the impact of missing estimated IFs and GDs on the connectivity of components in the signal. Since the loss of components can occur at any point in the TFD, it is necessary to use a performance indicator that can monitor local components. To achieve this, we utilized an LRE-based indicator that has been shown to be effective in detecting reconstructed TFDs with discontinuous autoterms in prior work [36]. Specifically, we employed two performance indicators that measure the error between the local number of signal components in the original LO-ADTFD with fully preserved auto terms (ρ (lo) (t, f )) and the TFD with estimated IFs/GDs (B(t, f )), using mean squared error (MSE) given as and mean absolute error (MAE) given as Higher values of MSE t, f and MAE t, f indicate a greater amount of missing estimated IFs and/or GDs in B(t, f ), suggesting that an inappropriate estimation strategy has been em-ployed. Normalizing the local component count enabled a fair MSE and MAE comparison across signals.

Results for Synthetic Signals
The efficacy of the proposed BM(t, f ) was evaluated on the synthetic signal example z LFM (t), with the results presented in Figure 8. Since all four LFM components of the signal z LFM (t) deviate from the time axis, the estimated GDs, shown in Figure 8b, offer a better connection than the estimated IFs, shown in Figure 8a. Consequently, the proposed method generates the BM(t, f ) map, correctly highlighting the TFD regions along the signal's components for the frequency localization approach, as shown in Figure 8c,d. In the case of the signal z LFM (t), Figure 9a demonstrates that the image-based STRE method was incapable of linking TF peaks that deviate from the time axis, resulting in the rejection of all four LFM autoterms and the estimation of interference terms that were almost parallel to the time axis. However, when utilizing the proposed combined IF and GD estimation, the autoterms were fully estimated with high connectivity (as shown in Figure 9b). Nevertheless, some interference and noise terms were also estimated due to the obtained spike in NC f ( f ) estimation shown in Figure 1c, which provided a frequency support threshold that is too small. While the BSS-STRE method did not experience difficulties with missing components, the estimated IFs were highly discontinuous, and not all of them belonged to the autoterms, as shown in Figure 9c. Conversely, Figure 9d demonstrates that the combined IF and GD estimation resulted in significantly improved component connectivity, with nearly all samples belonging to the autoterms. Similar results were obtained for the signal z mix (t). The image-based STRE method was unable to connect non-LFM components (or parts of components) that deviated from the time axis, as illustrated in Figure 10a. Meanwhile, the BSS-STRE method produced IF estimates that were discontinuous and dislocated from the true component ridge, as seen in Figure 10c. However, the proposed combined IF and GD estimation significantly improved the performance of both methods, as shown in Figure 10b,d. The performance of the B (shrink) (t, f ) obtained by the shrinkage operator is depicted in Figure 11. For both signal examples, the estimated component ridges in the TFD showed superior performance compared with those obtained using the image-based method, while being very similar to the BSS method, demonstrating high component connectivity and belonging to the autoterms of the signal.
The results presented in Table 1 demonstrate a significant reduction in MSE t, f and MAE t, f when using the proposed method for combined IF and GD estimation with STRE and NBRE information. For the signals z LFM (t) and z mix (t), the image-based method's estimation improved by 71.25% and 81.11% in terms of MSE t, f and 50.09% and 62.70% in terms of MAE t, f , respectively. The BSS method's estimation also improved by 92.95% and 83.17% in terms of MSE t, f and 82.37% and 65.07% in terms of MAE t, f , respectively.
Moreover, the results demonstrate that the B (shrink) (t, f ) obtained using the shrinkageoperator-based method for IF and GD estimation outperformed the image-based STRE and BSS-STRE methods-MSE t, f and MAE t, f values improved by up to 95.62% and 86.52%, respectively, for the signal z LFM (t) and 90.01% and 72.79%, respectively, for the signal z mix (t). Moreover, the results show that the obtained B (shrink) (t, f ) is competitive with the B (BSS) (t, f ) obtained using the BSS-STRE-NBRE algorithm, showing a reduction in MSE t, f for z LFM (t) by 3.20%, while the rest of the indicators are slightly in favor of the BSS-STRE-NBRE method (by up to 2.01%). (a) (b) Figure 11. Estimated IFs and GDs obtained in B (shrink) (t, f ) using the shrinkage operator for the signals: (a) z LFM (t); (b) z mix (t).

Sensitivity to Noise
In this study, the proposed shrinkage-operator-based, BSS-STRE-NBRE and imagebased STRE-NBRE methods for estimating IF and GD were evaluated for their robustness to noise. Synthetic signals, including z LFM (t) and z mix (t), were embedded in AWGN with an SNR that varied between 0 and 10 dB in 1000 independent simulations. The F1 score metric, which combines precision and recall, was used to evaluate the algorithm's performance, given as where true positives (TP) and true negatives (TN) indicate the number of samples that were correctly estimated as a signal or noise/interference component, respectively, while false positives (FP) and false negatives (FN) refer to the number of noise/interference or signal samples that were incorrectly identified as signal or noise, respectively. Note that F1 values range from 0 to 1, with higher values indicating better performance.
To further validate the F1 score metric, the 2D MSE between the noise-free and noisy TFDs was calculated using the scaled and squared Frobenius norm as follows: The F-norm value was defined as the squared norm of the difference between the two TFDs, divided by the total number of time-frequency bins, N t × N f . A lower F-norm value indicates better performance.
We evaluated the accuracy of the local number of signal component estimates by computing the MSE between NC t (t) and NC f ( f ) for both noise-free and noisy signals using Equation (3). The results, depicted in Figure 12, show that the LRE methods produce stable estimates for SNR values above 1 dB. Moreover, the noise sensitivity of all IF/GD estimation algorithms considered is determined by the sensitivity of LRE methods, as indicated by the F1 and F-norm values shown in Figures 13 and 14. Notably, the imagebased STRE-NBRE method is more sensitive to noise than the shrinkage-operator-based and BSS-STRE-NBRE methods, as even minor inaccuracies in the estimated local number of signal components can result in poor threshold values, as seen in Figure 9b.  Figure 15a presents the time-domain representation of the signal z EEG (t), revealing several spikes. These spikes are effectively captured in the TF domain using LO-ADTFD, as demonstrated in Figure 15b,c for the original and filtered signals, z EEG (t) and z EEGfilt (t), respectively, which also reveal an additional single sinusoidal component. Notably, the differentiator filter significantly reduced the background noise and enhanced the desired spike components, leading to a cleaner signal, z EEG filt (t).  Figure 16 demonstrates that the proposed BM(t, f ) effectively identified and separated the spike and sinusoidal components in both signals, z EEG (t) and z EEG filt (t). As a result, Figure 17 displays the extracted sinusoidal and spike components in separate TFDs. The local number of components, NC t (t) and NC f ( f ), estimated from these TFDs, shows a significant reduction in inaccurate local maxima compared with those obtained from the input LO-ADTFD, as illustrated in Figure 18.

Results for Real-Life EEG Seizure Signals
(a) (b) (c) (d) Figure 18. Comparison between the local number of signal components obtained by STRE and NBRE in starting TFD (dashed red line) and from extracted components using the proposed operators κ t and κ f (solid blue line) for the signals: (a,b) z EEG (t); (c,d) z EEG filt (t). Figure 19 displays the estimated IFs and GDs of the existing approach in [47] for the signals z EEG (t) and z EEGfilt (t). The IFs are estimated from the signal in the time domain, while the GDs are obtained from the Fourier transform of a signal using the duality property, which transposes a signal in the TF domain [47]. The results indicate that the approach is effective in estimating the IFs and GDs for the filtered EEG signal z EEGfilt (t) used in the original study [47]. However, the approach proved unsuitable for estimating the GDs of the unfiltered signal z EEG (t) due to the presence of background noise and the inconsistent number of components over the TF domain, leading to inaccurate GD estimates, as illustrated in Figure 19d. Furthermore, it is worth noting that the user needs to provide the global number of components to obtain these estimates, which may present a practical limitation when acquiring an unknown signal. Figures 20 and 21 compare the performance of the IF estimation methods against the proposed mutual IFs and GDs estimations. As the spike and sinusoidal components are perpendicular to each other, their estimated IFs exhibit completely opposite performances. Specifically, all IF estimation methods successfully estimated the IFs of the sinusoidal component, while the IF estimates of the spike components were highly discontinuous and, especially for the unfiltered signal z EEG (t), indistinguishable from the background samples. However, estimating the GDs for the spike components significantly improved their connectivity and overall preservation, as shown in Figures 20 and 21, using all considered IF and GD estimation methods.   The numerical results of the mutual IF and GD estimations, presented on illustrative z EEG (t) and averaged on a dataset containing 200 examples, are summarized in Table 2. The results confirm that the combined IF and GD estimation approach outperforms the IF estimation alone. Specifically, for the EEG dataset of unfiltered z EEG (t) signals, the imagebased STRE-NBRE and BSS-STRE-NBRE algorithms improved the MSE t, f and MAE t, f indicators by up to 42.23% and 30.08%, respectively, compared with the image-based STRE and BSS-STRE algorithms. Improvements were obtained for the dataset of filtered EEG signals z EEGfilt (t) also, where the MSE t, f and MAE t, f indicators were reduced by up to 34.96% and 33.41%, respectively, considering the same algorithms, namely the imagebased STRE-NBRE and BSS-STRE-NBRE. Furthermore, the obtained results show the superiority of the shrinkage operator to the image-based STRE and BSS-STRE algorithms, with improvements of up to 46.45% and 36.61% for the dataset of z EEG (t) signals and up to 31.71% and 30.35% for the dataset of z EEGfilt (t) signals in terms of MSE t, f and MAE t, f , respectively. Again, the shrinkage operator approach was shown to be competitive with the BSS-STRE-NBRE algorithm.

Interpretation of Obtained Results
The experimental results show that the proposed method, which generates a binary map BM(t, f ), successfully identifies regions in the TFD that contain signal components requiring different time or frequency localization approaches for all considered synthetic and real-life EEG signals. Furthermore, it is shown that estimating the local number of signal components on extracted TFDs significantly reduces local estimation inaccuracies compared with the original estimate on the TFD with all components present.
The results obtained when comparing mutual IF and GD estimation with IF estimation alone show a significant improvement in estimated component connectivity and preservation for all considered synthetic and real-life EEG signals. The image-based STRE-NBRE and BSS-STRE-NBRE methods effectively used frequency support information from NBRE, providing IF and GD estimates without requiring prior knowledge about the signal. The results also show that IF and GD estimates can be obtained using the shrinkage operator derived from sparse reconstruction, which has competitive performance compared with the BSS algorithm, while both outperforming the image-based algorithm. Additionally, the considered algorithms' robustness to noise was connected with the STRE and NBRE methods, where the BSS-STRE-NBRE and shrinkage-operator-based estimations outperformed the image-based STRE-NBRE algorithm for all considered SNRs. The image-based algorithm's dependence on LRE accuracy is the reason behind this, where even the smallest error can cause a threshold that is too small or too large for the method.
The advantage of the BSS algorithm over the shrinkage operator is that it extracts components one by one using the double-directional approach. This implies that the estimated samples follow a line, which is evident when comparing Figures 9d and 11a. However, this can be a disadvantage if interference is falsely chosen as a component due to its higher maxima than the autoterm's, as the BSS algorithm will force the estimation of an interference. In this case, tracking the largest local surfaces instead of only local maxima in the shrinkage operator can avoid some interference samples, as shown in Figures 20d and 21a.
The results show that the proposed method is feasible to be used for estimating IFs and GDs of EEG seizure signals, z EEG (t) and z EEG filt (t). In the case of the filtered signal z EEG filt (t), the estimated IFs and GDs using the proposed method are competitive with the approach in [47], with the significant advantage that the proposed method does not require the number of components to be set in advance. However, for the unfiltered signal z EEG (t), the proposed method outperforms the approach in [47] and shows feasibility for signals whose number of components changes over time.
It should be noted that when extracting intersecting components that require different localization approaches, a small portion of components near the intersection point may be extracted in a different TFD. This phenomenon is evident in Figure 17a,b for the unfiltered EEG signal z EEG (t), where we observe that small parts of spike components have been extracted alongside the sinusoidal component. This behavior can be attributed to the practical calculation of STRE and NBRE, where the sliding window with size Θ t or Θ f within the stable range defined in [33,51] detects a component's time or frequency support as a few samples or bins more than the ideal. Consequently, the TFD region borders in the proposed BM(t, f ) are slightly wider than the actual component's time or frequency support to accommodate for this behavior.

Conclusions
The analysis of signals that exhibit both rhythmic and spike features, such as EEG seizure signals, presents a significant challenge when utilizing conventional TF methods. In order to extract valuable components that are distributed across both the time and frequency axes, a comprehensive method is necessary. In this paper, we introduced a novel method for automatically estimating the IF and GD of a signal in the TF domain. In order to define TFD regions requiring a different time or frequency localization strategy, we proposed a method for generating a binary map BM(t, f ) based on the information from LRE methods. An increase in the local number of signal components obtained using the LRE methods was indicative of the presence of a component that may require a different localization approach than what was observed, whereas measuring IF and GD estimates with the proposed measure N r was effective for identifying discontinuous estimates.
Through the implementation of the suggested BM(t, f ), we successfully extracted components that necessitate either a time or frequency localization approach, thereby yielding more accurate evaluations of the numbers of local components using the STRE and NBRE methods. The STRE method's reduced accuracy for certain signals prompted modifications to image-based and BSS IF estimation algorithms, enabling them to efficiently incorporate the NBRE method and decrease their dependence on the STRE method. The proposed method yielded a notable enhancement in performance and facilitated the simultaneous estimation of IF and GD.
The results obtained demonstrate that the proposed method's combined IF and GD estimation outperforms the IF estimation alone. This was demonstrated through the analysis of noisy synthetic and real-life EEG seizure signals with characteristic rhythmic and spike features. In contrast to current methodologies, the proposed approach does not necessitate an a priori understanding of an input signal and is applicable to signals whose number of components varies with time or frequency.
The following research efforts will focus on the advancement of characteristics that can distinguish and categorize EEG signals from the surrounding environment, utilizing the IF and GD evaluations derived from this study. Furthermore, a research area of interest involves the creation of a concentration measure for TFDs utilizing the estimated IFs and GDs. The primary objective of this measure will be to impose a penalty for the lack of autoterms that occurs in signal processing using advanced TF techniques. Funding: This research was funded by the University of Rijeka under the projects "Analysis and Classification of Non-Stationary Signals Using Advanced Deep Learning Methods" (uniri-mladitehnic-22-16), "Computer-Aided Digital Analysis And Classification of Signals" (uniri-tehnic- , and the ZIP UNIRI project "Development of computational algorithms based on digital signal processing and artificial intelligence methods with application in modern transport systems".
Institutional Review Board Statement: Not applicable.

Informed Consent Statement: Not applicable.
Data Availability Statement: Not applicable.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript: