Adaptive instantaneous frequency estimation based on time-frequency distributions with derivative approximation

Signal Processing 160 (2019) 99–105

Contents lists available at ScienceDirect

Signal Processing journal homepage: www.elsevier.com/locate/sigpro

Short communication

Adaptive instantaneous frequency estimation based on time-frequency distributions with derivative approximation Yazan Abdoush a,∗, Jose A. Garcia-Molina b, Giovanni E. Corazza a a b

Department of Electrical, Electronic and Information Engineering (DEI), University of Bologna, Bolgona, Italy European Space Agency (ESA), Noordwijk, the Netherlands

a r t i c l e

i n f o

Article history: Received 27 July 2018 Revised 26 November 2018 Accepted 16 January 2019 Available online 18 February 2019 Keywords: Instantaneous frequency estimation Wigner distribution S-method Time-frequency analysis Adaptive estimators

a b s t r a c t A new method for nonparametric and adaptive instantaneous frequency (IF) estimation of monocomponent signals based on time-frequency distributions (TFDs) is presented. This method uses an estimate of the IF second-order derivative to approximate the width of the TFD observation window associated with the estimation least mean squared error (MSE), which was previously derived in a closed-form expression. The derivative estimate is obtained in two steps. First, a preliminary TFD is computed and its local maxima are used as a rough estimate of the IF trajectory. Thence, the continuous wavelet transform (CWT) is introduced for calculations of high-order derivatives. The proposed method is evaluated and compared with state-of-the-art algorithms based on the intersection of confidence intervals (ICI) rule. Numerical results demonstrate that the proposed method achieves a significant improvement in estimation accuracy. © 2019 Elsevier B.V. All rights reserved.

1. Introduction Instantaneous frequency (IF) estimation is a key topic in timefrequency (TF) analysis of nonstationary signals whose spectra evolve with time. Its applications are vital in fields that include, but not limited to, radar and sonar, speech analysis, communications, seismology, just to name a few. The mathematical concept of the IF, IF estimators, and IF-based applications are addressed in several works; we will mention here two excellent review papers [1,2]. When accurate signal models are not available, a simple yet effective method for nonparametric IF estimation is established based on the maxima position of a TF distribution (TFD) [3–5]. Identified sources of error are: 1. bias, originating from the IF high-order derivatives; 2. small noise, the TFD local maxima remain within corresponding auto-terms; and 3. large noise, causing the local maxima to be located outside the signal components. The first two sources of error are analyzed in [6], whereas the last is considered in [7]. Analysis of the error caused by small disturbances demonstrates that the estimation bias and variance are monotonically increasing and decreasing functions of the TFD window size, respectively. For the Wigner distribution (WD), which is

∗

Corresponding author. E-mail addresses: [email protected] (Y. Abdoush), jose.antonio.garcia. [email protected] (J.A. Garcia-Molina), [email protected] (G.E. Corazza). https://doi.org/10.1016/j.sigpro.2019.01.027 0165-1684/© 2019 Elsevier B.V. All rights reserved.

the core TFD and from which all members of Cohen’s family can be derived by smoothing kernels, the window size that resolves the bias-variance trade-off in the least mean squared error (MSE) sense was derived in [8], but was deemed of impractical value for being a function of the IF second-order derivative, which is, of course, unknown, since the IF itself is to be estimated. An algorithm for adaptive IF estimation which does not require a priori knowledge of the IF derivatives was presented in [8,9] based on the intersection of confidence intervals (ICI) rule. In short, this method tests adaptively candidates for the window length from a predefined discrete set, and for each size, a confidence interval, within which the true IF exists with a certain probability, is defined using the asymptotic formula of the estimation variance. An estimate of the optimal width is determined as the largest width in the set for which two successive intervals have at least one point in common. The ICI-based method was proven to provide good estimation accuracy in working environments characterized by moderate to high signal-to-noise ratios (SNRs); its performance was analyzed in [10] and procedures to properly adjust its parameters were presented. It was subsequently applied with other TFDs, including higher-order ones, in [11–14]. A modification of this method was proposed in [15], where the amount of overlap between the current and previous confidence intervals is tested to improve the accuracy of estimation. Djurovic´ and Stankovic´ in [16] refined the method behavior in high-noise environments by excluding narrow

100

Y. Abdoush, J.A. Garcia-Molina and G.E. Corazza / Signal Processing 160 (2019) 99–105

windows that produce large errors and applying a median filter to the estimates returned by middle-size windows to remover possible outliers. Despite the fact that the ICI-based method is of central importance in adaptive IF estimation of monocomponent signals, as accentuated in [2], the following limitations remain. First, since the method nature is based on statistical hypothesis testing, the IF estimation accuracy depends, obviously, on the cardinality of the set of windows used by the ICI rule. Second, if we make use of a larger set of window widths, we may improve on accuracy but at the expense of higher computational complexity, compromising the method practicality. Tackling the previous shortcomings, in this communication, we take a different approach to adaptive IF estimation. Essentially, we argue that the closed-form formula which relates the optimal width (in the MSE sense) to the IF secondorder derivative can be used directly as a criterion for width selection once a reasonable estimate of the second-order derivative is available. To this end, we first construct a preliminary TFD (PTFD), which is a WD or one based on the S-method (SM) [17], and use the peaks of its dominant ridges as a rough estimate of the IF trajectory. Thence, the continuous wavelet transform (CWT) is introduced to accommodate an estimate of the IF second-order derivative. Although numerical differentiation of noisy process is known to be a classical ill-posed problem, we show that, by a proper selection of a wavelet mother, the CWT combines smoothing with high-order differentiation, providing a good approximation of derivatives. Results on simulated data validate that the proposed method can outperform its counterparts built based on the ICI rule. 2. Background theory 2.1. IF estimation model

ω

y(n) which is given by [8]: h/2−1



wh (k )y(n + k )y∗ (n − k ) exp (− j2ωk )

(1)

k=−h/2

in which wh (k) ≡ wh (kT) is a real-valued symmetric window of length h. The estimation error, defined as ω ˆ (n ) = ω ˆ ( n ) − ω ( n ), is analyzed in [8]; therein, it is demonstrated that the error has the following bias and variance, respectively:

σ2 1 Var(h ) = σ 2 (h ) ∼ = c2 2 3 , A h

Bias(n, h ) ∼ = c1 h2 ω (n ),

(2)

where c1 and c2 are window-dependent parameters. The contradictory behavior between the bias and the variance with respect to the length h is made clear by (2), and the fact that the estimation bias is dependent on the IF itself establishes that any reasonable choice of the length h should be data-driven. Optimizing h by minimizing the estimation MSE, defined as the sum of the squared bias and the variance, results in [9]:

 hopt (n ) =

3c2 σ2 /A2

(2c1 ω (n ) )2

1 Pe (h ) = 1 − √ 2πσWD



× exp

1 / 7 .

(3)



+∞





z 1 − 0.5 erfc √ 2σWD

−∞

−(z − AWD )2 2 2σWD

h−1



dz,

(4)

2 = E σ 2 (2A2 + σ 2 ) with E and E dewhere AWD = E0 A2 and σWD 1  0 1  h/2−1  2−1 2 (k ), respectively. noting k=−h/2 wh (k ) and h/ w k=−h/2 h

3. Proposed method We present here an alternative algorithm for adaptive IF estimation. 3.1. Outline and motivation Essentially, we try to obtain an estimate of the IF second-order derivative, which, once available, can be used directly to adjust the window size by a direct application of (3). The main argument in support of this alternative approach comes from observing that, by letting a reasonable estimate of the IF second-order derivative control the width selection according to (3), variations of the MSE around its minimal value become significantly slower compared with the case in which the MSE changes directly as a function of h. To clarify this point, we begin by writing the MSE in two forms. In the first one, the bias and the variance are given by the set of Eq. (2), resulting in:



Let y(n ) = x(n ) +  (n ) be a noisy signal, where x(n) ≡ x(nT) is a sampled version of the continuous-time and noise-free signal x(t ) = A exp( jφ (t )) with T being the sampling interval.  (n) ≡  (nT) is a complex-valued additive white Gaussian noise (AWGN) of variance σ2 . The signal phase is φ (t) and A is a constant real-valued amplitude. By definition, the signal IF is the first-order derivative of its phase [1–3]: ω (n )  φ  (t )|t=nT . Here, the IF is estimated based on the maxima position of the WD according to: ωˆ (n ) = argmax WDy (n, ω ), where WDy (n, ω) denotes the WD of

WDy (n, ω ) =

If the signal is highly contaminated by noise, the previous analysis loses its validity. In such a case, it was demonstrated in [7] that the estimation error exhibits an impulsive nature, and the WD peaks in any location outside the auto-term position with the following probability:

MSE1 (n, h ) = c1 ω (n )h2

2

+

c2 σ2 /A2 . h3

(5)

While in the second formulation, h in (2) is substituted by  (n ))2 ]1/7 with ω  (n ) denoting an estimate of [(3c2 σ2 /A2 )/(2c1 ω the IF second-order derivative:

 

c σ 2 /A2 4 1/7

2   (n ) = 3 ω (n ) + 4 ω  (n ) MSE2 ω , (6)

 (n ) 8 ω







2



2

where c = c13 c22 /108. At a given time instant, the minimum MSE, denoted by MSEmin (n ), results from (5) by the substitu (n ) with ω  (n). tion h = hopt (n ), or from (6) by replacing ω Now, to understand how the functions in (5) and (6) varies  (n ) deviating around MSEmin (n ) as a result of h(n) and ω from their respective optimal values, we define the following

two quantities: L1 (ρ ) = MSE1 n, ρ hopt (n ) /MSEmin (n ) and L2 (ρ ) =





MSE2 ρω (n ) /MSEmin (n ), where the deviations of h(n) and  (n ) from hopt (n) and ω  (n), respectively, are described by a ω real parameter ρ . It can be shown that the ratios L1 (ρ ) and L2 (ρ ) change as functions of ρ according to: L1 (ρ ) = (4 + 3ρ 7 )/(7ρ 3 ) and L2 (ρ ) = (3 + 4ρ 2 )/(7ρ 8/7 ). These expressions, in turn, make clear that, as ρ shifts away from 1, L1 (n) and L2 (n) increase above their stationary point (at ρ = 1) with a rate being significantly larger for L1 (n) than it is for L2 (n), as illustrated in Fig. 1. We conclude that, the IF estimation approach in which an estimate of the IF second-order derivative controls the size of the TFD observation window using (3) can tolerate error in the derivative estimation larger than the error an approach based on adaptive length selection from a discrete scheme (like the ICI method) can endure in its search for the most suitable length. The evident flatness of L2 (ρ )  (n ) does around its stationary point indicates that the estimate ω not have to be highly accurate in order to achieve a low MSE, but

Y. Abdoush, J.A. Garcia-Molina and G.E. Corazza / Signal Processing 160 (2019) 99–105

101

30

3.2. Derivative approximation by the CWT

25

Despite being of fundamental importance in various fields of study, numerical differentiation (ND) is known to be an ill-posed problem, meaning that small perturbations lead to large error in the approximate solution. In this context, several works have reported the potential offered by the CWT to derivative calculation [18–20], which originates from its intrinsic capability to combine smoothing with differentiation, as briefly explained in what follows. It should be mentioned first that the CWT itself, and a closely related transform known as the S-transform (ST), are well-studied IF estimators [21,22], but here, the CWT will be used for a completely different purpose, which is derivative approximation. The CWT of a continuous-time signal f(t) is defined by

20 15 10 5 0

0

0.5

1

1.5

2

2.5

3

rather reasonable. To obtain such an estimate, the method proposed in this communication constructs first a PTFD and uses its dominant ridges as a rough estimate of the IF trajectory. Then, the IF second-order derivative is estimated by the CWT. This estimate is finally inserted into (3) to adjust the window size of an improved TFD, whose local peaks provide a more accurate estimation of the IF.

where the analyzing wavelet (assumed real in the previous definition) results from a mother wavelet ψ (u) by dilation and translation through s and t, respectively. The transform in (7) can be written equivalently using the convolution operator according to: √ CWT f (t, s; ψ ) = 1/ s f (t ) ∗ ψ (−t/s ). For a wavelet characterized by m vanishing moments, there exists a function θ (t) whose mthorder derivative relates to this wavelet by [23]:

Spectral Response Amplitude

CWT f (t, s; ψ ) =

Second-order Derivative

Fig. 1. Variations of L1 (ρ ) and L2 (ρ ) (in dB). Function L2 (ρ ) increases above the stationary point at ρ = 1 much slower compared with L1 (ρ ). Note, for example, L1 (ρ = 2 ) is more than 5 times larger than L2 (ρ = 2 ).

0.2 0.1 0 -0.1 -0.2 0

100

200

300

400



+∞ −∞



1 u−t f (u ) √ ψ s s

 du,

Gaussian Derivative Savitzky-Golay

0.03

0.02

0.01

0

-3

-2

-1

0

1

2

(a)

(b) Second-order Derivative

Second-order Derivative

0.01 0.01

Theoretical

0.005 0 -0.005 -0.01 100

200

300

Time Sample

(c)

3

Frequency (rad/s)

Time Sample

0

(7)

400

500

CWT Theoretical

0.005 0 -0.005 -0.01

0

100

200

300

400

500

Time Sample

(d)

Fig. 2. Example of computing the second-order derivative of noisy sine waveform (SNR = 15 dB). (a) Derivative computed by conventional finite-difference method. (b) Spectral response amplitude of filter s (t) given by second-order derivative of a Gaussian and SG second-order differential filter (note the high side lobes of SG filter). Parameters of filters are set such that they have the same cutoff frequency. Derivative computed by: SG differential filter in (c), and the CWT in (d). Most accurate result is provided by the CWT.



Y. Abdoush, J.A. Garcia-Molina and G.E. Corazza / Signal Processing 160 (2019) 99–105 +∞

−∞

t i ψ (t ) dt = 0, ⇐⇒

i = 0, . . . , m − 1

ψ (t ) = (−1 )m

Estimate A2 and σ2

m

d θ (t ). dt m

(8)

It was also proven that θ (t) satisfies the conditions of a smoothing function by having a fast decay and nonzero integral, which establishes that θ (t) behaves as a lowpass filter in the spectral domain. By the theory of convolution, the CWT of f(t) can be written now as:

CWT f (t, s; ψ ) = sm+1/2 Rearranging (9) gives:





1 dm −t f (t ) ∗ θ dt m s s





CWT f (t, s; ψ ) dm 1 −t = f (t ) ∗ m θ dt s s sm+1/2



Choose a suitable hint satisfying the accuracy requirement on Pe (hint )



.

= f (t ) ∗ s (t ).

Construct a PTFD using hint and obtain a rough estimate of the IF trajectory

(9)

(10)

Estimate the IF-second order derivative using the CWT

Eqs. (9) and (10) clearly demonstrate that the CWT through a mother wavelet with m vanishing moments is intrinsically equivalent to the mth-order derivative of a lowpass-filtered version of the analyzed signal (up to a constant factor). We remark that the lowpass behavior exhibited by s (t) is certainly desirable if the signal to be differentiated is noisy. This behavior is controlled by the dilation parameter, which sets the filter cutoff frequency. With increasing s, the spectral response of s (t) becomes narrower, suppressing a larger amount of noise power. On the other hand, excessive smoothing leads to loss of important derivative details. Therefore, a trade-off should be adopted when selecting a suitable value for s, but, in general, the sole constraint on the chosen mother wavelet is given in (8). Note also that the mth-order derivative can be computed by one transform procedure. For the purpose of this work, we selected the second-order derivative of a Gaussian as √ mother wavelet, namely: θ (t ) = 1/ π exp(−t 2 ). Example of computing the second-order derivative of a noisy sine waveform with a comparison using Savitzky-Golay (SG) filter [24] is illustrated in Fig. 2.

The PTFD is aimed to provide a rough estimate of the IF trajectory, which, once available, is used by the CWT for derivative calculations. We stress again that this preliminary estimation does not have to be highly accurate, since, as explained before, an IF estimate with small MSE can be obtained even with an error in the derivative estimate, provided that this error is not very large. Indeed, it turns out that only large errors, which dislocate the PTFD peaks far from the corresponding auto-terms, impair the quality of the final IF estimates. To avoid the occurrence of such errors, based on an estimate of A and σ  obtained following [8, (23, 24)], the PTFD window size, denoted henceforth by hint , is determined as the smallest in a predefined discrete set of window widths such that the probability of large error is smaller than a threshold: Pe (hint ) < Pth . By satisfying this requirement on Pe (hint ), large errors occur with a small probability so as to allow for the smoothing offered by the CWT to effectively alleviate their impact on derivative calculations. We remark that the reasons behind selecting the narrowest window among those that satisfy the previous accuracy requirement are: 1. to avoid IF estimation with large bias, because large biasedness in the initial IF estimate is more detrimental to derivative calculations than the variance, whose impact can be better mitigated by smoothing; 2. abrupt changes in the IF trajectory cannot be revealed by a wide observation window, leading to serious errors in derivative estimation. In the numerical realization, Pth is set to 0.1% at high values of SNR (SNR ≥ 5 dB), while at lower SNRs (SNR < 5 dB), to avoid selection of very large windows, which are prone to inner-artifacts, Pth is increased to 1/3 (similar probability thresholds were used in [16]). For signals with highly

ˆ opt (n) Compute an adaptive TFD using h

Estimate the IF using the adaptive TFD Fig. 3. Flowchart of proposed method for adaptive IF estimation.

nonlinear IF laws in high-noise environments, the WD is contaminated by aggravating interfering inner-terms, which constitute in this case the dominant source of estimation errors. As a possible remedy to inner-terms, we may apply the proposed algorithm with the SM, which is a reduced-interference TFD cantoning the appealing features of the WD and the short-time Fourier transform (STFT) (see [17] for more details and implementation). We note that the SM was also used for adaptive IF estimation in [16]. The proposed method is summarized by the flowchart in Fig. 3. 4. Numerical evaluation In this section, the proposed IF estimation method is evaluated and compared with that based on the ICI rule [8] and its modified version described in [16] with the extra condition of Lerga and Sucic [15]. The performance is quantified by means of the normalized -10 NMSE (dB)

3.3. Preliminary TFD

Use the derivative estiˆ opt (n) mates to calculate h

-10

Original ICI (WD) Modified ICI (WD) Proposed (WD) Original ICI (SM) Modified ICI (SM) Proposed (SM)

-20 -30

NMSE (dB)

102

-40 -50 -60 -5

0

5

10

SNR (dB) (a): Signal A

15

20

25

Original ICI (WD) Modified ICI (WD) Proposed (WD) Original ICI (SM) Modified ICI (SM) Proposed (SM)

-20 -30 -40 -50 -5

0

5

10

15

20

25

SNR (dB) (b): Signal B

Fig. 4. Comparison between IF estimation methods in terms of NMSE using two test signals. All considered methods are implemented using WD and SM.

Y. Abdoush, J.A. Garcia-Molina and G.E. Corazza / Signal Processing 160 (2019) 99–105

103

0 0

100 200 300 400 500

1 0

Time (S)

0

100

ICI 10

2

0 0

100 200 300 400 500

100 200 300 400 500

Time (S)

Time (S)

Time (S)

Proposed

ICI

Proposed

-4

0.01

0.04

Time (S)

ICI

Poposed

0.02 0 0

100 200 300 400 500

Time (S)

200

0.005 0 0

400

Time (S)

200

400

Time (S)

ICI

Proposed

100 50 0 0

100 200 300 400 500

100 200 300 400 500

Time (S)

Time (S)

ICI

400

Width

Width

Width

500

0 0

2

0 0

500

1

0 0

100 200 300 400 500

400

Width

0 0

300

Error

5

Error

Error

10

-4

200

IF (rad/S)

2

Error

2

IF (rad/S)

IF (rad/S)

IF (rad/S)

3

200 0 0

100 200 300 400 500

Time (S)

Proposed

ICI

100 50 0 0

100 200 300 400 500

Time (S) Proposed

Fig. 5. IF estimation for Signal A. Proposed method is compared with the modified ICI-based algorithm. First and second columns from left: SNR = 11 dB and WD is used for TFD. Third and fourth columns from left: SNR = 0 dB and SM is used for TFD. First row from top: estimated IF (solid line) and true IF (red dashed line). Second row: absolute error. Third row from top: window width. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

MSE (NMSE) of the IF estimates over a large number of MonteCarlo simulations (here is set to 100). In each iteration, AWGN is added to a noiseless signal at SNR (SNR = A2 /σ2 ) ranging from 0 to 25 dB with a 2-dB step. For numerical evaluation, we use two signals1 , Signal A and Signal B, with IF laws, ωA (n) and ωB (n), given according to:

ωA ( n ) =

⎧ 2.58π ⎪ , ⎪ ⎪ 3 ⎪ ⎪ ⎪ ⎪ 1.5π π ⎪ ⎪ + cos (4π (n1 − 1 ) ), ⎨

3 3 1 . 5 π π ⎪ ⎪ + cos (2π (n1 − 2 ) ), ⎪ ⎪ 3 3 ⎪ ⎪ ⎪   ⎪ ⎪ ⎩ 2.5π + 0.5π − 2.5π ( n1 − 3 ), 3 3

for n ∈ {0, 127} for n ∈ {128, 255} for n ∈ {256, 383} for n ∈ {384, 511} (11)

1 ωB (n ) = 0.5π + 0.4π cos(4π n2 ) − cos(12π n2 ) +

1 cos(20π n2 ) , 3

3

n ∈ {0, 255},

(12)

where n1 = n/128, n2 = (n − 128 )/256. The amplitude of both signals is set to unity.

1 Similar signals are used in [2], and in [25] to evaluate TFDs with complex-time argument, which improve the energy concentration along nonlinear IF laws.

4.1. Simulation results In all experiments, the proposed method and the ICI-based algorithms are implemented using the WD and the SM as well. Results of the NMSE (in decibels) of IF estimation are shown in Fig. 4 as functions of SNR. The following conclusions can be drawn from the experimental results. 1. The proposed method using the WD outperforms the ICI-based algorithms at SNRs larger than 5 dB, whereas at lower SNRs, the modified ICI method provides more accurate IF estimates. The performance shortcoming in high-noise environments is due to emphatic inner-interference presented in the WD, which impairs the validity of (3). 2. When the TFDs are implemented using the SM instead of the WD for inner-terms suppression, the proposed method regains its effectiveness, improving on the ICI-based algorithms in terms of the accuracy of the IF estimates. Note that the SM has also lead to significant improvements in the results of the ICI-based methods at SNRs below 5 dB. 3. At SNRs larger than 5 dB, where the inner-terms are less prominent, the accuracy of the proposed method estimates is slightly degraded when the SM is used instead of the WD (this degradation is more evident with Signal B). Nonetheless, these estimates are still more accurate than those of the ICI-based algorithms. Therefore, the suggestion is to use the SM only at SNRs less than 5 dB; otherwise, the WD is preferred. We illustrate in Fig. 5 a single trial of IF estimation using Signal A at two values of SNR (0 and 11 dB). Fig. 6 replicates the results in Fig. 5 for Signal B. The PTFDs used in the proposed method to obtain the results in Figs. 5 and 6 are illustrated in Fig. 7. Further details on the numerical realization of the methods are provided in a supplementary document.

Y. Abdoush, J.A. Garcia-Molina and G.E. Corazza / Signal Processing 160 (2019) 99–105

50

0 0

100 150 200 250

50

ICI 10

50

100

150

IF (rad/S)

ICI

Proposed

0.1

100

150

Time (S)

ICI

Poposed

0.05 0 0

200

50

0.05

0 0

100 150 200

50

Proposed 50

100 0 0

50

Time (S)

ICI

100 150 200

Time (S)

ICI

100 150 200 250

Time (S)

50

Time (S)

200

0 0

100 150 200 250

Error

Error

50

Width

Width

Proposed -4

50

50

100 150 200 250

Time (S)

Time (S)

100

50

Time (S)

2 0 0

200

200

0 0

0 0

100 150 200 250

4

1 0 0

50

2

Time (S)

Width

10

Error

Error

2

-3

0 0

100 150 200 250

Time (S)

2

Width

0 0

2

IF (rad/S)

2

IF (rad/S)

IF (rad/S)

104

0 0

100 150 200 250

50

Time (S)

Proposed

100 150 200 250

Time (S)

ICI

Proposed

Fig. 6. IF estimation for Signal B. Proposed method is compared with the modified ICI-based algorithm. First and second columns from left: SNR = 11 dB and WD is used for TFD. Third and fourh columns from left: SNR = 0 dB and SM is used for TFD. First row from top: estimated IF (solid line) and true IF (red dashed line). Second row: absolute error. Third row from top: window width. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

3 IF (rad/S)

IF (rad/S)

3 2 1 0

2 1 0

0

100

200

300

400

500

0

100

(a): Signal A, SNR = 11 dB.

400

500

3 IF (rad/S)

IF (rad/S)

300

(b): Signal A, SNR = 0 dB.

3 2 1 0

200

Time (S)

Time (S)

0

50

100

150

200

Time (S) (c): Signal B, SNR = 11 dB.

250

2 1 0

0

50

100

150

200

250

Time (S) (d): Signal B, SNR = 0 dB.

Fig. 7. PTFDs for Test Signal A and Test Signal B at low and high SNRs. PTFDs are used in the proposed method to obtain rough IF estimate. First column from left: TFDs constructed using WD with hint = 16. Second column from left: TFDs constructed using SM with hint = 32. These PTFDs are used to obtain the results in Figs. 5 and 6. All PTFDs are displayed with the same colormap.

Y. Abdoush, J.A. Garcia-Molina and G.E. Corazza / Signal Processing 160 (2019) 99–105 Table 1 Comparison between adaptive IF estimation methods in terms of computation time. SNR is set to 11 dB, and WD is used as TFD. Results are given in seconds, and experiments were run on Windows-7 Intel Core i-5 platform with 8 GB RAM. Proposed method is the most computationally efficient. Signal

Original ICI

Modified ICI

Proposed method

Signal A (N = 512) Signal B (N = 256)

2.09 0.81

4.59 1.94

0.55 0.34

Another advantage of the proposed method is reduced computational complexity, since only two TFDs are implemented (PTFD and final TFD), each of which requires computational complexity of order O(N2 log2 N), besides the CWT for derivative calculation, which can be realized with O(N) computational complexity. With the ICI-based algorithm, on the other hand, various TFDs should be constructed, depending on the size of the used set of widths. A comparison between the developed method and the ICI-based algorithms in terms of computation time is provided in Table 1. Results demonstrate the computational efficiency of the proposed method with respect to its counterparts. 5. Conclusion A new method for nonparametric and adaptive IF estimation based on TFDs is presented. This method uses an estimate of the IF second-order derivative, provided by the CWT, to approximate the optimal window width of a WD. At moderate to high SNRs, numerical results demonstrated the method superiority over its counterparts based on the ICI rule. Furthermore, at lower SNRs and with application of the SM instead of the WD for suppression of inner-terms, the proposed method improved the accuracy of IF estimation. Moreover, the computational complexity of the developed method is significantly less than that of the ICI-based algorithms. Future research could involve 1. applying the presented method with high-order TFDs, 2. extending it for IF estimation of multicomponent signals using an algorithm for component extraction and tracking (as in [26]), and comparing it with the quasimaximum-likelihood (QML) [27] estimator of polynomial phase signals. Declaration of interest None. Supplementary material Supplementary material associated with this article can be found, in the online version, at doi:10.1016/j.sigpro.2019.01.027. References [1] B. Boashash, Estimating and interpreting the instantaneous frequency of a signal. II. Algorithms and applications, Proc. IEEE 80 (4) (1992) 540–568.

105

[2] L. Stankovic´ , I. Djurovic´ , S. Stankovic´ , M. Simeunovic´ , S. Djukanovic´ , M. Dakovic´ , Instantaneous frequency in time–frequency analysis: enhanced concepts and performance of estimation algorithms, Digit. Signal Process. 35 (2014) 1–13. [3] B. Boashash, Time Frequency Signal Analysis and Processing: A Comprehensive Reference, second, Elsevier, London, UK, 2003. [4] E. Sejdic´ , I. Djurovic´ , J. Jiang, Time-frequency feature representation using energy concentration: an overview of recent advances, Digit. Signal Process. 19 (1) (2009) 153–183. [5] L. Cohen, Time-frequency distributions-a review, Proc. IEEE 77 (7) (1989) 941–981. [6] V.N. Ivanovic, M. Dakovic, L. Stankovic, Performance of quadratic time-frequency distributions as instantaneous frequency estimators, IEEE Trans. Signal Process. 51 (1) (2003) 77–89. [7] I. Djurovic, L. Stankoic, Influence of high noise on the instantaneous frequency estimation using quadratic time-frequency distributions, IEEE Signal Process. Lett. 7 (11) (20 0 0) 317–319. [8] V. Katkovnik, L. Stankovic, Instantaneous frequency estimation using the Wigner distribution with varying and data-driven window length, IEEE Trans. Signal Process. 46 (9) (1998) 2315–2325. [9] L.J. Stankovic, V. Katkovnik, Algorithm for the instantaneous frequency estimation using time-frequency distributions with adaptive window width, IEEE Signal Process. Lett. 5 (9) (1998) 224–227. [10] L. Stankovic, Performance analysis of the adaptive algorithm for bias-to-variance tradeoff, IEEE Trans. Signal Process. 52 (5) (2004) 1228–1234. [11] V. Katkovnik, L. Stankovic´ , Periodogram with varying and data-driven window length, Signal Process. 67 (3) (1998) 345–358. [12] L. Stankovic, V. Katkovnik, Instantaneous frequency estimation using higher order L-Wigner distributions with data-driven order and window length, IEEE Trans. Inf. Theory 46 (1) (20 0 0) 302–311. [13] B. Barkat, B. Boashash, L. Stankovic, Adaptive window in the PWVD for the IF estimation of FM signals in additive Gaussian noise, in: Acoustics, Speech, and Signal Processing, 1999. Proceedings., 1999 IEEE International Conference on, 3, IEEE, 1999, pp. 1317–1320. [14] Z.M. Hussain, B. Boashash, Adaptive instantaneous frequency estimation of multicomponent FM signals using quadratic time-frequency distributions, IEEE Trans. Signal Process. 50 (8) (2002) 1866–1876. [15] J. Lerga, V. Sucic, Nonlinear IF estimation based on the pseudo WVD adapted using the improved sliding pairwise ICI rule, IEEE Signal Process. Lett. 16 (11) (2009) 953–956. [16] I. Djurovic, L. Stankovic, Modification of the ICI rule-based IF estimator for high noise environments, IEEE Trans. Signal Process. 52 (9) (2004) 2655–2661. [17] L. Stankovic, A method for time-frequency analysis, IEEE Trans. Signal Process. 42 (1) (1994) 225–229. [18] C.-L. Fu, X.-L. Feng, Z. Qian, Wavelets and high order numerical differentiation, Appl. Math. Modell. 34 (10) (2010) 3008–3021. [19] A. Gentile, A. Messina, On the continuous wavelet transforms applied to discrete vibrational data for detecting open cracks in damaged beams, Int. J. Solids Struct. 40 (2) (2003) 295–315. [20] X. Shao, C. Ma, A general approach to derivative calculation using wavelet transform, Chemom. Intell. Lab. Syst. 69 (1–2) (2003) 157–165. [21] E. Sejdic, I. Djurovic, et al., Quantitative performance analysis of scalogram as instantaneous frequency estimator, IEEE Trans. Signal Process. 56 (8) (2008) 3837–3845. [22] E. Sejdic, L. Stankovic, M. Dakovic, J. Jiang, Instantaneous frequency estimation using the S-transform, IEEE Signal Process. Lett. 15 (2008) 309–312. [23] S. Mallat, A wavelet tour of signal processing: the sparse way, Academic Press, 2008. [24] R.W. Schafer, What is a Savitzky-Golay filter? [Lecture notes], IEEE Signal Process. Mag. 28 (4) (2011) 111–117. [25] L. Stankovic, Time-frequency distributions with complex argument, IEEE Trans. Signal Process. 50 (3) (2002) 475–486. [26] J. Lerga, V. Sucic, B. Boashash, An efficient algorithm for instantaneous frequency estimation of nonstationary multicomponent signals in low SNR, EURASIP J. Adv. Signal Process. 2011 (1) (2011) 725189. [27] I. Djurovic, L. Stankovic, Quasi-maximum-likelihood estimator of polynomial phase signals, IET Signal Proc. 8 (4) (2014) 347–359.