Estimating the number of components of a multicomponent nonstationary signal using the short-term time-frequency Rényi entropy

The time-frequency Rényi entropy provides a measure of complexity of a nonstationary multicomponent signal in the time-frequency plane. When the complexity of a signal corresponds to the number of its components, then this information is measured as the Rényi entropy of the time-frequency distribution (TFD) of the signal. This article presents a solution to the problem of detecting the number of components that are present in short-time interval of the signal TFD, using the short-term Rényi entropy. The method is automatic and it does not require a prior information about the signal. The algorithm is applied on both synthetic and real data, using a quadratic separable kernel TFD. The results confirm that the short-term Rényi entropy can be an effective tool for estimating the local number of components present in the signal. The key aspect of selecting a suitable TFD is also discussed.

1 Time-frequency distributions and instantaneous frequency estimation

Nonstationary signals analysis and quadratic class of time-frequency distributions
Practical signals in the various fields of engineering (telecommunications, acoustics, biomedical engineering) are nonstationary, with the instantaneous frequency (IF) being their key parameter [1].One of the fundamental information when analyzing such signals is the number of components present in the signal.When applied to a time-frequency distribution (TFD), the Rényi entropy measures the signal complexity [2,3].Signals of high complexity are composed of a large number of elementary components [2].
When dealing with highly complex signals, such as multicomponent nonstationary signals, several pieces of information are required for their characterization.Classical approaches of the time signal representation, x(t), and the frequency representation, X(f), are not best tools for obtaining those information when dealing with multicomponent signals.These representations define the signal duration, the changes of amplitude in time, as well as the entire signal frequency content.Time-frequency representations (TFRs), or TFDs, are two variable functions, C s (t, f), defined over the two-dimensional (t, f) space [4].Such a joint TFR shows how the frequency content of a signal changes in time.
One of the most popular TFDs, introduced by Wigner and extended by Ville to analytic signals [4], has been treated as a pseudo probability density function in [2,3,5] to which the Rényi entropy has been applied as a measure of signal complexity.The intuitive idea of the Wigner-Ville distribution (WVD) was to obtain a kind of instantaneous signal spectrum by performing the Fourier transform of a function related to the signal, called the kernel function K s (t, τ).The WVD of a signal s(t), denoted as W s (t, f), represents a monocomponent frequency modulated (FM) signal as a knife-edge ridge in the (t, f) plane, whose crest is the IF of the signal [4].
Let s(t) be an analytic FM signal of the form [4] s(t) = a(t)e jφ (t) , where a(t) is the instantaneous signal amplitude (assumed to be equal to one in the rest of the article), and the signal IF is defined as the time derivative of its instantaneous phase j(t) [4] The requirement for the WVD to be a knife-edge ridge can mathematically be interpreted as a series of delta functions tracking the signal IF [4]: which leads to the signal kernel definition: K s (t, τ ) = F −1 τ ←f {δ(f − f i (f ))} = e j2π f i (t)τ = e jφ (t)τ .(4) Since j'(t) is not directly available, it can be replaced by the central finite-difference approximation [4] By substituting ( 5) into (4), the kernel function and the WVD are defined as [4] Hence, the WVD can be understood as the Fourier transform of the signal kernel K s (t, f), also known as the instantaneous autocorrelation function (IAF) of s(t).
From (7), we can notice that using the IAF as the kernel function brings nonlinearity in the WVD.The effects of this nonlinearity will be most evident in the case of multicomponent signals, as explained below.Note that, in general, a component in the (t, f) domain is a ridge of energy concentration whose peaks follow the component IF law [4].
Let us consider an analytic signal of the form where are the instantaneous cross-correlation functions that will add the third term to the WVD of the two-component signal [4]: This third term in the summation in (12) is called the cross-term defined as It appears in the (t, f) plane in between the signal components, often degrading the quality of signal representation in the (t, f) plane.
The rule of interference construction in the WVD can be summarized as follows.Two points belonging to the signal will interfere to create a third point which will be located on their geometrical midpoint.The amplitude of the interference will be proportional to the double product of the amplitudes of the interfering points.In addition, the interferences oscillate perpendicularly to the line joining the two signal points, assuming both positive and negative values, with the frequency of oscillation being proportional to the distance between these two points [4,6].
It can be deduced from the general interference rule that interferences will be also present in the case of monocomponent signals with nonlinear FMs, called inner cross-terms in [4].

Separable kernel time-frequency distributions
A generalization of the WVD is given by a class of TFDs known as the quadratic class [4].Distributions from this class are defined as where g(t, f) is the time-frequency kernel filter, and the double asterisk denotes a double convolution in t and f.Each TFD belonging to the quadratic class can thus be defined by the double convolution of the WVD and the time-frequency kernel [4].Equation 14 can be rewritten in terms of the IAF function and the Dopplerlag kernel g(ν, τ) [4]: Equations 14 and 15 lead to the definition of the quadratic TFDs with a separable kernel, i.e. a function with, respectively, independent Doppler and lag kernel filters: Such a kernel function allows a precise and independent control of the smoothing of the WVD in both time and frequency sense.It has been shown that TFDs with separable kernel filters are easy to design, and they achieve a good compromise in components energy concentration and interfering terms suppression [4,7].If h (τ) and g(t) are some typical windows (e.g.Hamming windows), the obtained separable kernel quadratic distribution is known as the Smoothed Pseudo-WVD (SPWVD), which is mathematically defined as [4] In order to improve the representation quality of the quadratic separable kernel distribution, the distribution known as the modified B distribution (MBD) [4,8,9] has been proposed.Its kernel definition, which is lag independent, is [7] where k = Γ(2b)/(2 2b-1 Γ 2 (b)) is the normalizing factor, Γ(•) is the gamma function, and b is a real, positive number.It has been shown in [7] that the MBD outperforms other TFDs in terms of time-frequency resolution and interference suppression for a large class of non-stationary signals, such as newborn and adult EEG signals as well as heart rate variability signals.
The advantages of time-frequency signal representations over the classical analysis approaches will be illustrated on an example of a two-component signal whose components are linearly and sinusoidally frequency modulated.Figure 1 shows that the time representation and frequency representation used separately are not adequate for this multicomponent nonstationary signal analysis; from the time representation (Figure 1a) no information about the signal IF can be obtained, while the frequency representation (Figure 1b) does not provide any information about the arrival times of individual frequencies.Neither of these representations indicates the presence of multiple components.Figure 1c shows the spectrogram of the signal with the Hamming window of duration 65 s.The WVD of the same two-component signal is shown in Figure 1d, which is a representation highly corrupted by interference terms.The interference gets successfully reduced by the timefrequency smoothing performed by the SPWVD (Figure 1e) with the time and lag Hamming windows of durations 25 and 65 s, respectively.The SPWVD also achieves high energy concentration around the signal components IFs, when compared to the spectrogram.The SPWVD can be considered as a special case of separable kernel TFD.

The definition and limitations of the global Rényi entropy
Some useful properties of the TFDs belonging to the quadratic class refer to the preservation of signal energy in the (t, f) plane and the marginal conditions.The integration of the TFD over the entire (t, f) plane results into signal energy [4]: while the integration over frequency and time respectively leads to the instantaneous power and the energy spectrum: The information content and complexity of a probability density function can be measured by the entropy function.The TFDs from the quadratic class satisfy the marginal conditions given by Equations 21 and 22, and after the normalization of C s (t, f) with respect to the signal energy E s : Therefore, the instantaneous power and energy spectrum can be understood as one-dimensional densities of signal energy in time and frequency, while TFDs may be interpreted as 2-D probability density functions [5].Hence, we would expect the classical Shannon entropy [2] to be an acceptable tool for measuring the complexity and information content of a nonstationary signal in the (t, f) domain.Since C s n (t, f ) in (24) represents a probability density function, it is natural to expect that a multicomponent signal will have larger entropy when compared to a single pulse in the (t, f) plane.As explained in the first section, nonpositivity (due to the presence of interfering terms) is one of the characteristics of the quadratic class of TFDs.As a TFD can be negative in some regions of the (t, f) plane, the Shannon entropy cannot be used in practice as a signal complexity measure, due to the logarithm in (24).
A solution to this limitation was proposed in [2], introducing the generalized entropies of Rényi: As shown in [2], when the parameter a (for the Shannon entropy a 1) is an odd integer value, the oscillatory structure of interferences are annulled under integration.Let us next consider several examples illustrating the use of the Rényi entropy in the (t, f) domain.Note that all TFDs that are used in this article are normalized as per (23).

Example I: Case of two Gabor logons with same time duration and FM
When treating a TFD as a probability density function, multicomponent signals will yield a larger entropy value when compared to the entropy of a single component of the same signal.Let us consider the Rényi entropy of an ideal TFD, denoted as I s (t, f), of a compactly supported signal s(t) (in Figure 2 represented by a Gabor logon [3] in the (t, f) plane, being an elementary time-frequency building block).For the two-component signal s(t) + s(t -Δt), the ideal TFD (Figure 3) will take the form The Rényi entropy of the two-component signal TFD, H a (I s(t)+s(t-Δt) (t, f)), carries exactly one more bit of information when compared to the Rényi entropy of the one component signal TFD, H a (I s(t) (t, f)) (as long as the time shift Δt is larger than the time support of the Gabor logon) [2], i.e.
Thus, the number of components can be determined as When the third-order Rényi entropy (a = 3) is computed for the signals in Figures 2 and 3, the following results are obtained: This example confirms the accuracy of the Rényi entropy counting property when the entropy of one of the components is known in advance.

Example II: Case of two Gabor logons with same FM and different time duration
Let us now consider a signal composed of two components with different time durations.The signal, consisting of two Gabor logons, s 1 (t) and s 2 (t), with durations 96 and 32 s, respectively, having same frequency supports is shown in Figure 4.
As expected, the signal component s l (t) that occupies a larger region of the (t, f) plane exhibits a considerably larger value of the Rényi entropy when compared to the entropy of the component with the shorter-time support, s 2 (t).Consequently, the estimation of the number of components based on the difference of the Rényi entropies of the entire signal and one of its components fails regardless which of the two components is chosen as the reference signal.The implication is that the global Rényi entropy can correctly detect the number of components only when all components have same time supports, with arbitrary time/frequency shifts in the TF plane.

Example III: Case of a signal with different FM components of same time duration
Next, we considered a multicomponent signal whose components have different frequency modulations.An example of such a signal is given in Figure 5; the twocomponent signal is composed of a component with parabolic FM (denoted as s l (t)), and a component with a linearly decreasing FM (denoted as s 2 (t)).Both components have a time duration of 256 s.Then, Equation 28 does not have a unique solution, as Thus, the presented examples indicate that the counting property of the Rényi entropy is restricted to multicomponent signals composed of components with similar time and frequency supports only.In addition, having to know the entropy of one of the components in advance, makes this approach highly impractical.To  remedy this disadvantage, the next section presents a novel method for estimating the local number of components in a signal that can be applied independently of the components respective durations in time and frequency.In addition, the proposed automatic procedure does not require the prior knowledge of the complexity of one of the signal components.
3 A novel algorithm for estimating the number of components using the short-term Rényi entropy

Assumptions and constraints
As previously shown, the Rényi entropy is a good indicator of the signal complexity (the number of components present in the signal) only when 1. all components have the similar structure in the time-frequency plane, and 2. the Rényi entropy of a single component is available.
However, in most practical applications this will not be the case: in general, signals encountered in many engineering and multidisciplinary fields (telecommunications, acoustics, radar, sonar, etc.) are usually mixtures of signals with different durations and frequency bandwidths.This makes the number of components estimation based on the global Rényi entropy ineffective for real-world signals.
What is required is a solution for determining the number of components present in a signal in a shorttime interval.Such a method is presented here, based on the estimation of the short-time Rényi entropy.It exploits the fact that signals with similar time durations and frequency supports have similar Rényi entropies.Since the Rényi entropy of a signal is invariant to time and frequency shifts [2], it is expected that a signal represented in the time-frequency domain by two shortenergy impulses will have twice the energy in the (t, f) plane of a signal represented by a single energy impulse of the same duration, and its Rényi entropy should be larger by one bit.This suggests that instead of observing the entire (t, f) plane when detecting the number of signal components, we should focus on a finite time interval of the (t, f) plane, and compare the Rényi entropy of this short-time segment with the Rényi entropy of a reference test signal with the same time support.In this way, the number of components present in the chosen time interval will be automatically estimated.

The proposed algorithm
In essence, the proposed algorithm consists of annulling the TFD, C s (t, f), of the multicomponent signal outside the time interval Δt, centered around the time instant p, to obtain a TFD of the form where t 0 controls the length of the observed time interval.We next compare the short-term Rényi entropy of this TFD, denoted as C s p (t, f ) , with the Rényi entropy of the TFD of the reference cosine signal (an arbitrary chosen stationary cosine signal), C ref p (t, f ) , for the same Δt.In the proposed algorithm, the reference TFD, C ref (t, f), must be a TFD of the same type and with the same set of parameters as the signal TFD, C s (t, f).C ref (t, f) must also have same dimensions as C s (t, f).The selected synthetic reference signal is a cosine signal of arbitrary amplitude and arbitrary constant frequency.Since the FM affects the bandwidth of a component in the TF domain (signals with fast changes in the IF may present relatively larger bandwidths in the TF domain), a synthetic signal with constant or linear FM is used as a reference signal in order to maintain a constant bandwidth.The amplitude of the reference signal can also be arbitrary chosen since the Rényi entropy of a signal is amplitude invariant.
For each time instant p, a different time portion of the (t, f) plane is extracted, and thus a different value of the Rényi entropy is obtained.After comparing the obtained values of the Rényi entropy of the observed multicomponent signal and the reference one, a function of the instantaneous number of components, n(p), is obtained.
By isolating a short-time interval of the TFD of the observed signal and annulling the TFD outside this interval, the influences of different time and frequency supports of various components on the counting property of the Rényi entropy are removed.This is a result of the fact that different components now locally have same time durations and similar bandwidths (see Figure 6b), as it was the case in Example I. Hence, when only a short-time interval of the TF plane is observed, the signal reduces to the special case of Example I, for which it has been shown in [2] that the counting property of the Rényi entropy holds.In other words, even if the signals globally have different time and frequency supports, their Rényi entropies are locally comparable.
The first step of the algorithm consists in the thresholding of the TFD to remove noise and interference low energy peaks that may locally affect the entropy of the signal.The same must be done with the TFD of the reference signal.
The fundamental step in determining the local number of a signal components using the short-term Time-Frequency Rényi entropy is the choice of the TFD.As it has been shown in [2], the Rényi entropy is invariant to cross-term since they annulate under integration over the entire (t, f) plane with odd powers of a (25).Naturally, this is not the case in our algorithm, where only a short-time portion of the (t, f) plane is observed, and consequently the eventual presence of cross-terms would cause inaccurate results.Thus, the key requirement we put on the TFD is the minimization of the interference terms spreading in the sense of both the time and frequency axis.TFDs with separable kernels successfully reduce the interferences, by an independent smoothing of the WVD in time and frequency [4].They have also been shown to outperform other popular TFDs in terms of the time-frequency resolution [7].
As mentioned above, the next step of the proposed algorithm requires the definition of the time interval Δt.In this article, the time interval is chosen as a small odd value (so that it can be centered exactly around the central time instant p), Δt = 7 s, to insure temporal localization of events even when dealing with fast changing nonstationary signals [10].The sampling period for all examples presented in this article is T s = 1 s.The annulation of the TFD of the observed signal, C s (t, f), outside the time interval Δt is performed by the scalar product of C s (t, f) and a two-dimensional time-frequency masking window W p (t, f): where For example, for the signal in Figure 6a, W 130 (t, f) will isolate a time segment of duration Δt = 7 s of the signal TFD around the time instant p = 130 s.The SPWVD, C s (t, f), of the analyzed signal is shown in Figure 6a, and the time slice of the SPWVD around p = 130 s in Figure 6b.The SPWVD of the reference cosine signal, C ref (t, f), and its time slice are given in Figure 7.
The presence of a component is checked by a minimum energy criterion that requires the total energy of the time slice to be larger than a threshold value (here chosen as 0.001 * E s /(N/Δt), where E s is the signal energy, N is the signal length, and Δt is the considered time interval of C s (t, f)).Once C s p (t, f ) and C ref p (t, f )  are obtained, and the presence of the component is detected, the number of components at the observed time instant p is calculated as By varying p from p = 4 to p = max(t) -3, and repeating the previous steps of the algorithm for each p, a function of the local numbers of components is  obtained.The flowchart in Figure 8 graphically illustrates the steps of the algorithm.

The experiment
The algorithm presented in the previous section exploits the short-term Rényi entropy to detect the number of components present in a short-time interval of the analyzed signal.In this section, we test its performance on both synthetic and real-life signals examples.In all reported tests, the two-dimensional time-frequency window, W p (t, f), has time duration Δt = 7 s (the sampling interval is T s = 1 s).The SPWVD of both the observed and reference signals (chosen as a cosinusoidal signal with the constant frequency of 0.1 Hz) is calculated with its time and lag smoothing windows [4] chosen as Hamming windows with durations of 15 s.In general the recommended windows width durations are between N/10 and N/20, where N is the duration of the signal.With this length of the time and lag windows, the interference present in the WVD is significantly smoothed, not affecting the estimation of the short-term Rényi entropy of the signal.Also, our extensive simulations have shown that the most accurate results, in terms of correctly detecting the local number of components, are obtained for the choice of the Rényi entropy parameter a = 7, as illustrated in Figure 9b.
Figure 9a shows a three-compound signal, whose components have parabolic frequency modulations and different time durations.As seen from Figure 9b, the proposed algorithm has correctly detected the number of components present at each time instant for three different values of the parameter a. Figure 9c shows the performance of the algorithm for three different values of Δt, confirming that Δt = 7 s gives the most accurate results.
The results obtained for the signal simulating the sum of two echoes from a rotating object, are shown in Figure 10.From Figure 10b we can observe that when the components overlap in the (t, f) domain, the number of components detected by the algorithm is 1.In other words, the algorithm reports the number of different regions in the (t, f) plane that are occupied by the signal  at a given time moment, independently of its amplitude in the (t, f) plane.The fact that the Rényi entropy of a signal is invariant to its amplitude allows the use of a cosine signal of unit amplitude as a reference.
Figure 11a shows a signal corrupted by additive white Gaussian noise (SNR = 10 dB).The results in Figure 11b are the number of components for the noise-free signal (solid line), the noisy signal for SNR = 15 dB (dashed line), and the noisy signal for SNR = 10 dB (dotted line).The presence of the noise has not brought significant deteriorations in the instantaneous number of components estimation, confirming the algorithm robustness to additive white noise.

Performance on real data
Finally, a real test signal (bat sound representing a natural sonar system [11]) has been analyzed to show the performance of the short-term Rényi entropy when applied to different TFDs.In many applications [12][13][14], the local number of components is an integer.In this case, the introduction of a threshold is required.The thresholding is also useful in the case when different components overlap for very short periods in time (as it is the case with some of the components in Figure 12).Thus, a threshold needs to be introduced in order to improve the sensibility of the algorithm.Our simulations have shown that the optimal threshold for the detection of the first component is 0.1, and for the ith component it is (i -1) + 0.2. Figure 12 presents the results obtained from the bat signal SPWVD (with the time and lag windows of duration 33 s), the MBD [4] (b = 0.1), and the Choi-Williams distribution (CWD) [4] (s = 2), with the solid line representing the thresholded results, and the dashed line representing the results without thresholding.In real-life signals, low energy components can be expected (as it is, for example, the minor component that appears in Figure 12a, c and  12e), thus for the information about their presence not to be neglected, a normalization of the TFD is required.Figure 12 shows that the MBD (Figure 12c) has most correctly detected the number of signal components, followed by the SPWVD (Figure 12a).The MBD shows major sensitivity to the weak spectral component.The CWD, due to the presence of pronounced cross-terms, performs poorly, showing unwanted oscillations in the local number of components.

Strengths and limitations of the algorithm
The results in this section illustrate a high accuracy of the presented algorithm in detecting the local number of signal components.Even in the case of noisy signals or low energy components, the algorithm has correctly  estimated the components number.The presented method provides the signal analyst with information on the minimum number of components present in the signal.All this can be useful in various applications that require components separation and extraction [12][13][14].
The choice of the TFD is crucial for the successful performance of the algorithm.TFDs with separable kernels are recommended in order to avoid the undesirable influence of both inner and outer artifacts of the signal TFD on the results.In this article, the MBD was shown to be a good choice for the real data local component number estimation.

Conclusion
This article proposes a method for estimating the local number of signals components.It is based on the shortterm Rényi entropy of signals in the time-frequency plane.Using the Rényi entropy of a short-term segment of a TFD of a multicomponent nonstationary signal, relative to the short-term Rényi entropy of a reference signal, the number of components present in the signal can be accurately estimated.The proposed method does not require any a priori information about the analyzed signal, nor the knowledge of the Rényi entropy of one of the signal components.The method was tested on various synthetic signals, including signals embedded in additive white Gaussian noise, and its use in practice was illustrated on a real-life signal.The method is sensitive to the selection of the TFD.The presented results indicate that the MBD, being an example of Separable kernel TFDs, is a good choice of a TFD when the proposed method is applied in practical situations [4].These results show that the proposed algorithm can be useful in many applications that require component count and component separation; and, it can be a preferred alternative to other methods such as the Empirical Mode Decomposition [15].

Figure 2 AFigure 3
Figure 2 A Gabor logon in time (top), frequency (left), and joint time-frequency domain using the SPWVD (center). n

Figure 4
Figure 4 Two Gabor logons with different time durations, shown in time, frequency, and time-frequency (SPWVD).

Figure 5
Figure 5 A two-component signal composed of a parabolic FM component and an LFM component, in time, frequency, and joint time-frequency domain (SPWVD).

Figure 6
Figure 6 (a) SPWVD of a two-component signal with linear and sinusoidal FMs; (b) 7 s time slice of the SPWVD around p = 130 s.

Figure 7
Figure 7 (a) SPWVD of the reference cosine signal; (b) 7 s time slice of the SPWVD around p = 130 s.

Figure 8
Figure 8 Flowchart of the proposed algorithm.

Figure 10
Figure 10 (a) SPWVD of a signal simulating the sum of two echoes from a rotating object, (b) number of components in time.

Figure 11
Figure 11 (a) SPWVD of a signal composed of an LFM component and two Gabor logons embedded in additive white noise with SNR = 10 dB; (b) number of components in time for the noise-free signal (solid), signal corrupted by additive white noise (SNR = 15 dB (dashed), and SNR = 10 dB (dotted)).

Figure 12 (
Figure 12 (a) SPWVD of a bat echo location signal; (b) number of components in time (thresholded results-solid line, results without thresholding-dashed line); (c) MBD of a bat echo location signal; (d) number of components in time (thresholded results-solid line, results without thresholdingdashed line); (e) CWD of a bat echo location signal; (f) number of components in time (thresholded results-solid line, results without thresholding-dashed line).