Accurate time-frequency estimation in sαs noise via memory-dependent derivative

Abstract

This letter presents a time-frequency estimation approach based on memory-dependent derivative to obtain accurate spectrograph interpolation information. The memory correlation derivative is the convolution of a time-varying signal with a dynamic weighting function over a past time period with respect to a common derivative. Considering the described method, discrete data from previous times can be derived to estimate the signal values at the current time and to reduce the effect of noise. Fourier transforms with different scales and delay transforms are used as kernel functions to obtain energy-concentrated time-frequency curves with higher resolution and without frequency leakage. Besides, the memory-dependent derivative with adjustable scale factor is used to overcome time-frequency grid mismatches. Furthermore, differing from the phase accumulation manner of conventional time-frequency estimation, ℓ₁-norm suppresses the heavy-tailed effect from outliers, thus the robustness of estimator can be enhanced greatly. By suitably choices of scale factor, the estimator can be tuned to exhibit high resolution in targeted regions of the time-frequency spectrum.

Introduction

Extracting the time-frequency trajectory of a non-linear frequency modulated signal from small samples with noise is a foundational problem in statistical signal processing. One challenging task is that the resolution of time-frequency cell is constrained by the time bandwidth product. ¹ This indicates that the linear time-frequency analysis (TFA) cannot achieve an arbitrarily high time-frequency resolution at the same time. The other task is that TFA method designed for Gaussian noise performs poorly in the presence of outliers. The noise with heavy-tailed distribution, like impulsive noise,^2,3 will make the correlation function severely deviated from the theoretical value. ⁴ Therefore, the motivation of this letter is that time-frequency estimator designed for large samples and Gaussian noise perform poorly in the presence of outliers in a small sample.

To improve the performance of TFA methods, many endeavors have been made in the past many years, and plenty of effective methods are proposed. Generally, these methods can be divided into two strategies: time-frequency reassignment and parameterized TFA.^5–7 Time-frequency reassignment method, for example, synchrosqueezing transform (SST), belongs to a post-processing technology, and it can improve the readability of time-frequency representation obviously. ⁸ If the modulated frequency of a signal can be known in advance, the reassignment operator aims to sharpen the time-frequency representation by assigning the average of energy in a domain to the gravity center of these energy contributions. The parameterized TFA method, for example, general parameterized TFA (GPTFA), ⁹ is inspired by the fact that different analyzed window can achieve the different time-frequency resolution. The parameterized model needs to be estimated from conventional TFA, so it is needed the initialized TFA with a fair energy concentration.^10–17

Although these advanced TFA technologies can provide more precise insights into the complex structure of a signal, their inherent limitations cannot be ignored, for example, the accuracy is disturbed by time-frequency cells and outliers.^18–22 Minimum ℓ₁-norm optimization model has found extensive applications in linear parameter estimations. The phase accumulation used by conventional method is converted to ℓ₁-norm to guarantee robust frequency localization in SαS noise.^23–28

Considering that impulse noise is instantaneous, the average power within a set time is adopted to reduce the effects of instantaneous impulses partially. Memory-dependent derivative with adaptive kernel function guarantees exact localization for the time-frequency trajectory, so the SαS noise can be suppressed. Furthermore, minimizing ℓ1-norm of frequency residual within sliding window eliminate noise interference, for example, impulsive noise.^29–33

This paper is organized as following. In Section 2, the FM signal model is described. Section 3 discusses non-linear prediction model via memory-dependent derivative. In Section 4, time-frequency estimation with memory-dependent derivative is presented. Section 5 presents some processing results based on numerical simulations. The last section concludes this paper.

FM signal model

An analytic signal $x\in C^L$ with the time-varying instantaneous frequency can be given by:

$$x(t)=a{e}^{i{\int }_0^t\phi (u)\mathrm{du}}$$

(1)

where $a\in R^{+}$ and $\phi (u)$ is the amplitude and instantaneous frequency respectively, u stands for time, and $i=\sqrt{-1}$ .

The conventional linear TFA is established on the assumption of the considered signal being piece-wisely linear frequency modulation in a short time. For example, in the short time $[u-\tau ,u]$ , the instantaneous frequency $\phi (u)$ can be approximated as a constant, that is,

$${\phi (u)|}_{u\in [u-\tau ,u]}\approx {\omega }_0$$

(2)

where ${\omega }_0$ is the first-order approximation of the true frequency.

In fact, the window function needs to be short enough for strongly modulated frequency signal, otherwise, the instantaneous frequency is still time-varying. However, small samples within a short windows not only leads to a poor time-frequency resolution, but also sensitive to noise at a low signal-to-noise ratio (SNR) level. Specifically, the impulsive noise whose variance or energy is infinite causes large error for the signal within short window. Thus, this letter proposes an accuracy and robust time-frequency estimation method for non-linear FM signal.

Non-linear prediction model via memory-dependent derivative

In this section, the memory-dependent derivative as a predict model of FM signal is derived from the integral of Taylor series at expansion points over an interval. Linear prediction is a mathematical operation where the current value of a linear frequency signal is estimated as a weighted function of previous samples. The most common representation is given by

$$\overline{x}(t)={\int }_{t-\tau }^{t}h(t-\tau )\overline{x}(v)\mathrm{dv}$$

(3)

where v stands for time, $\overline{x}(t)$ is the predicted signal value, $\overline{x}(v)$ is the previous observed values, and $h(t-v)$ is the kernel function.

Most prediction method can model linear frequency signal, but fails for FM signal. One of the simplest and effective specifications is to make $x(t)$ piecewise linear. The Taylor series of FM signal converges for all $t$ in the interval $(v-\sigma ,v+\sigma )$ ( $\sigma >0$ is infinitesimal) and the sum of the series is equal to $x(t)$ .

In this case, the first-order Taylor series expansion as a linear function approximates $x(t)$ within interval $(v-\sigma ,v+\sigma )$ , that is,

$$x(t)\approx x(v)+(t-v)x\prime (v)$$

(4)

Additionally, the second and higher order are ignored because the derivative amplifies the noise inevitably and the resulting inaccuracy worsens for higher derivatives. ¹⁰

Next, as expansion point v slides in the short time interval $v\in [t-\tau ,t]$ , we can obtain a series of linear functions converged to $x(t)$ . Liking the linear prediction method as shown in equation (4), the current value $x(t)$ can be estimated by using the weighted sum of first-order Taylor series expanded in the previous period $v\in [t-\tau ,t]$ . The prediction value is an average of linear superposition of a series of local linear convergence values as follows

$$x(t)\approx \frac{1}{\tau }{\int }_{t-\tau }^{t}x(v)+(t-v)x\prime (v)\mathrm{dv}$$

(5)

Inserting equation (1) and $x^{\prime }(v)\approx i{\omega }_{0}x(v)$ into equation (5), it can be rewritten as

$$\begin{array}{c}\multicolumn{1}{c}{x(t)\approx \frac{1}{\tau }{\int }_{t-\tau }^t[t-v-\frac{i}{{\omega }_0}]x\prime (v)\mathrm{dv}}\\ =\frac{1}{\tau }{\int }_{t-\tau }^t[t-v-\frac{i}{a{\omega }_0}]\mathrm{ia}{\omega }_0{e}^{i{\omega }_{0}v}\mathrm{dv}\hfill \end{array}$$

(6)

which indicates that $x(t)$ can be estimated by its previous values within $[t-\tau ,t]$ .

In fact, ${\omega }_0$ in equation (6) is unknown for the problem of instantaneous frequency estimation. Thus, insteading of $t-v-\frac{i}{a{\omega }_0}$ , a kernel function h should be chosen properly to present the time-frequency characteristics of x. Enlightened by these, we have introduced the concept of memory-dependent derivative to predict the non-linear FM signal in a distinct way. Generalizing the kernel function of equation (6) to the general form, the memory-dependent derivative $D_{\tau }[·]$ ¹¹ can be defined as follows

$$D_{\tau }[x(t)]=\frac{1}{\tau }{\int }_{t-\tau }^{t}h(t-v)x\prime (v)\mathrm{dv}$$

(7)

which is an integral form of $x^{\prime }(v)$ with a kernel function $h(t-v)$ on a slipping interval $[t-\tau ,t]$ . A memory-dependent derivative is characterized by its kernel function which can be chosen according to the frequency function of signal.

Time-frequency estimation with memory-dependent derivative

The kernel function determines the accuracy of the prediction model for FM signal, so a memory-dependent derivative is characterized by its kernel function. A fundamental issue is to design an appropriate kernel function which improve accuracy of time-frequency estimation result.

Given a rectangular window $g(t-v)$ to truncate the signal to instead of the integral interval $[t-\tau ,t]$ , the spectrum expression of $D_{\tau }[x(t)]$ can be obtained through Fourier delay transformation as follows

$$\begin{array}{c}\multicolumn{1}{c}{F\left\{D_{\tau }\left[x\left(t\right)\right]\right\}=\frac{1}{\tau }{\int }_{-\infty }^{+\infty }[{\int }_{-\infty }^{+\infty }g(t-v)h(t-v)x\prime \left(v\right)\mathrm{dv}]e^{-i\omega t}\mathrm{dt}}\\ =\frac{1}{\tau }{\int }_{-\infty }^{+\infty }x\prime \left(v\right)\left[{\int }_{-\infty }^{+\infty }g(t-v)h(t-v)e^{-i\omega t}\mathrm{dt}\right]\mathrm{dv}\hfill \\ =\frac{1}{\tau }{\int }_{-\infty }^{+\infty }x\prime (v)F\{g(t-v)h(t-v)\}\mathrm{dv}\hfill \\ =\frac{1}{\tau }{\int }_{-\infty }^{+\infty }x\prime (v)F\{g(t-v)\}*F\{h(t-v)\}\mathrm{dv}\hfill \end{array}$$

(8)

where $F\{·\}$ and * denotes the Fourier transform and convolution operator, respectively. According to the definition of Fourier transform, we define these Fourier transform pairs as follows

$$\begin{array}{c}F\{x\prime (v)\}=\overline{X}(\tilde{\overline{\omega }}),F\{g(t-v)\}=e^{-i\omega v}G(\omega ),\hfill \\ F\{h(t-v)\}=e^{-i\omega v}H(\omega )\hfill \end{array}$$

(9)

Substituting equation (9) into equation (8), it can be rewritten as

$$\begin{array}{c}\multicolumn{1}{c}{\mathrm{Equation}\left(8\right)=\frac{1}{\tau }{\int }_{-\infty }^{+\infty }x\prime \left(v\right)\{e^{-i\omega v}G\left(\omega \right)*e^{-i\omega v}H\left(\omega \right)\}\mathrm{dv}}\\ =\frac{1}{\tau }{\int }_{-\infty }^{+\infty }x\prime \left(v\right){\int }_{-\infty }^{+\infty }{e}^{-i\kappa v}G\left(\kappa \right)e^{-\mathrm{iv}(\omega -\kappa )}H(\omega -\kappa )d\kappa \mathrm{dv}\hfill \\ =\frac{1}{\tau }{\int }_{-\infty }^{+\infty }G\left(\kappa \right)H(\omega -\kappa ){\int }_{-\infty }^{+\infty }x\prime \left(v\right)e^{-\mathrm{iv}\omega }\mathrm{dvd}\kappa \hfill \\ =\frac{1}{\tau }{\int }_{-\infty }^{+\infty }G\left(\kappa \right)H(\omega -\kappa )\overline{X}\left(\overline{\omega }\right)d\kappa \hfill \end{array}$$

(10)

Traditionally, the TFA can discretize the continuous space to a finite number of time-frequency cells presented the time-frequency trajectory. This discretized model is simple and easy to handle analytically, but it inevitably brings the drawback that the restriction of uncertainty principle causes the energy dispersion. Especially, in the case of short window size, the conventional TF representation is not sufficient to generate a more energy concentrated results.

Utilizing the Fourier scale transformation, an kernel function with tunable grids is designed to generate the more energy-concentrated result by introducing the scale coefficient $\rho >0$ . The kernel function is given by

$$h(t)=\frac{\sqrt{2}}{2\hat{\omega }}{e}^{-i\frac{\rho t-4\hat{\omega }}{4\rho }t}$$

(11)

where $\hat{\omega }$ is estimated frequency. Its Fourier transform can be expressed as

$$F\{h(t)\}={\int }_{-\infty }^{+\infty }\frac{\sqrt{2}}{2\hat{\omega }}{e}^{-i\frac{\rho t-4\hat{\omega }}{4\rho }t}{e}^{-i\omega t}\mathrm{dt}=\frac{1}{\hat{\omega }}{e}^{-i{(\frac{\hat{\omega }}{\rho }-\omega )}^2+i\frac{\pi }{4}}$$

(12)

Obviously, the frequency domain grid can be refined for $\rho >1$ . By using the scale transform, the basis mismatch problem can be resolved in exact location of instantaneous frequency. For signals with well-separated frequencies, we show the $\rho$ is roughly proportional to the number of frequencies, up to polylogarithmic factors. Instantaneous frequency estimation is possible even though kernel function is not inherently continuous at all and does not satisfy any sort of restricted isometry conditions.

Inserting equation (12) into equation (10), we can obtain as follows

$$T{F}_x(v,\hat{\omega }|\rho )=\frac{{\omega }_0}{\tau \hat{\omega }}{e}^{-i{(\frac{\hat{\omega }}{\rho }-\omega )}^2+i\frac{\pi }{4}}\overline{X}(\overline{\omega }){\int }_{-\infty }^{+\infty }{e}^{-i\kappa v}G(\kappa )d\kappa$$

(13)

Then, equation (13) is given by

$$T{F}_x(v,\hat{\omega }|\rho )=\frac{\widehat{\tilde{\overline{\omega }}}}{\tau \hat{\omega }}{e}^{-i{(\frac{\hat{\omega }}{\rho }-\overline{\omega })}^2+i\frac{\pi }{4}}{\int }_{-\infty }^{+\infty }{e}^{-i\kappa v}G(\kappa )d\kappa$$

(14)

The energy of result in the time-frequency space divided into many small grid is well-concentrated, and it is not restricted by the Heisenberg uncertainty principle because it is based on an inner product operator with the refined kernel function.

Furthermore, because a peak at the true frequency at $\omega =\overline{\omega }$ and some null for the side-lobes is presented at non-targeted cells, the estimated frequency $\overline{\omega }$ within window function can be calculated the derivative of $T{F}_x(v,\hat{\omega }|\rho )$ with respect to $\hat{\omega }$ as follows

$$\frac{\partial T{F}_x(v,\hat{\omega }|\rho )}{\partial \hat{\omega }}=-\frac{i2(\frac{\hat{\omega }}{\rho }-\overline{\omega })}{\rho }T{F}_x(v,\hat{\omega }|\rho )=0$$

(15)

As the window function $g(t-v)$ moves, residuals of time-frequency trajectories in the whole time-frequency spectrum can be readily formulated as

$$||\frac{\hat{\omega }}{\rho }-\overline{\omega }|{|}_1=||i\frac{\rho }{2}{\mathrm{TF}}_x^{-1}(v,\hat{\omega }|\rho )\frac{\partial T{F}_x(v,\hat{\omega }|\rho )}{\partial \hat{\omega }}|{|}_1$$

(16)

The ℓ₀ or ℓ₂ norm is sensitive to outlier observations. This sensitivity is the principal reason for exploring alternative norms. Procedures for the proposed method that involve the norm have been developed to increase robustness. The ℓ₁ norm has been applied in numerous variations of frequency estimation. ℓ₁ norm is an attractive alternative to ℓ₀ norm or ℓ₂ norm because it provides globally optimal solutions in polynomial time and can impart robustness in the presence of outliers and is indicated for models where standard Gaussian assumptions about the noise may not apply.^34–38 Therefore, this equation (16) proposes an ℓ₁ norm procedure based on the efficient calculation of the optimal solution of the ℓ₁ norm best-fit hyperplane problem.

Additionally, a proper scale factor can ensure a good concentration of energy in the time-frequency presentation. Thus, the optimal scale factor ${\rho }_{\mathrm{opt}}$ can be determined by calculating the derivative of $T{F}_x(v,\hat{\omega }|\rho )$ with respect to $\rho$ as follows

$$\frac{\partial T{F}_x(v,\hat{\omega }|\rho )}{\partial \rho }=-\frac{i2\hat{\omega }(\frac{\hat{\omega }}{\rho }-\overline{\omega })}{{\rho }^2}T{F}_x(v,\hat{\omega }|\rho )=0$$

(17)

Due to both of $T{F}_x(v,\hat{\omega }|\rho )$ and $\rho$ are positive, the optimal scale factor ${\rho }_{\mathrm{opt}}$ is obtained by solving equation (17) as follows

$${\rho }_{\mathrm{opt}}=\frac{\overline{\omega }}{\hat{\omega }}$$

(18)

Optimal scaling factor ${\rho }_{\mathrm{opt}}$ is determined by the frequency deviation. Specifically, the estimated frequency is equal to the true, that is, $\hat{\omega }=\overline{\omega }$ , at ${\rho }_{\mathrm{opt}}\to +\infty$ .

Considered an FM signal x(t) and a short and sliding rectangular window $g(t-v)$ , the estimated instantaneous frequency $\hat{\omega }$ at time v is obtained by minimizing equation (16) as follows

$$\underset{\hat{\omega },\rho }{\min }||\hat{\omega }-\overline{\omega }|{|}_1=||{\mathrm{iTF}}_x^{-1}(v,\hat{\omega }|{\rho }_{\mathrm{opt}})\frac{\partial T{F}_x(v,\hat{\omega }|{\rho }_{\mathrm{opt}})}{\partial \hat{\omega }}+\frac{i}{\hat{\omega }}|{|}_1$$

(19)

From a statistical perspective, ℓ1-norm is robust to impulsive noise. ¹² Equation (19) not only effectively mitigates the effect of artifacts caused by the noise, but also produces a time-frequency estimation for the nonlinear FM signal with an excellent concentration.

Numerical simulations

We compare the performance of the memory-dependent derivative algorithm with several classical algorithms, and they are assessed by tests on generalized sinusoidal FM (GSFM) signal. The signal is given as

$$\phi (t)=350+200t+150t^2$$

(20)

The sample time is T = 1 S and the sample frequency is fs = 5000 Hz. In the simulation, symmetric $\alpha$ -stable (S $\alpha$ S) distributions is used to model the real-world noise since it represents a different noise statistical character by changing parameters.^13–16 Additionally, the average relative error of the estimated instantaneous frequency $\hat{\phi }[m]$ with respect to the true instantaneous frequency is defined as follows

$$\xi (T)=\frac{1}{T}{\int }_0^T||[\hat{\phi }(t)-\phi (t)]/\phi (t)||\mathrm{dt}$$

(21)

Figure 1 displays the performance of SST, ⁸ GPTFA ⁹ and the proposed algorithm at $\alpha =1.6$ and GSNR = −5 dB (GSNR means generalized signal-to-noise ratio). The window size is 0.02 S. The frequency modulation rate of the FM signal can be estimated by minimizing equation (21). As shown in Figure 1(b), the energy of the signal would be scattered over a relatively wide frequency band, so their estimated IF trajectories cannot show themselves off very clearly. These estimated IF trajectories fluctuate around the true IF, even deviate from the true IF significantly at some moments in Figure 1(c) and (d). However, from the results, we can see that the proposed algorithm clearly outperform GPTFA and SST, especially at the moment of impulses appearance. Time- frequency estimation generated by the proposed algorithm is more excellent concentration because its time-frequency cell is adjustable adaptively. The relative error of GPTFA, SST and the proposed algorithm are 1.63%, 10.9% and 0.54%, respectively. This shows that the proposed method is suitable for heavy tailed and low GSNR environment.

Figure 1.

(a) FM signal at SNR = −5 dB and $\alpha =1.6$ , (b) result by GPTFA, (c) result by SST, and (d) result by the proposed method.

Next, Figure 2 illustrates the average relative error with the whole sample time for SST, GPTFA, and the proposed algorithm for 200 Monte-Carlo simulations. As shown in Figure 2(a), $\xi (T)$ by the proposed algorithm is very small and is insensitive to the S $\alpha$ S noise with various degrees of heavy tailed effect, that is, the characteristic parameter $\alpha \in [0.2,2]$ every 0.2 at the fixed GSNR = 5 dB. Although $\xi (T)$ by GPTFA is insensitive as well to $\alpha$ , its error is much higher than the proposed algorithm. Besides, the estimated results of SST have larger bias as the characteristic parameter $\alpha$ increases. For non-Gaussian noise with high GSNR, results of the estimated IF by the proposed method and SST are similar for the whole data, but the proposed method is more accurate under $\alpha$  > 1.2. For the different GSNR, that is, GSNR ∈ [−10, 10] dB every 2 dB at the fixed $\alpha =1.6$ , the Figure 2(b) reveals that the proposed algorithm has a better estimation performance than the other two methods whether GSNR is high or low. The proposed method is superior to other two algorithms, which are sensitive to the low GSNR SαS noise. And the SST performance is seriously degraded under the condition of low GSNR. There is no lower bound for the estimation performance shown in Figure 2 due to the complex probability density of the alpha-stable noise. It is difficult to determine the algorithm’s lower limit for frequency estimation. However, the error pattern is shown in Figure 2.

Figure 2.

Results for the FM signal at SNR = 5 dB and $\alpha =1.6$ . (a) $\alpha \in [0.2,1.8]$ every 0.2 at GSNR = 5 dB, (b) $\mathrm{GSNR}\in [-10,10]$  dB every 2 dB at $\alpha =1.6$ .

Conclusion

This letter presents the time-frequency estimation algorithm based on memory-dependent derivative for the exact and robust instantaneous frequency estimation of nonlinear FM signal. Unlike previous work in time-frequency estimation, the instantaneous frequency is not assumed to lie on a grid, but can assume any value in the frequency domain. Furthermore, the phase accumulation used by conventional method is converted to ℓ1-norm to guarantee robust frequency localization in S $\alpha$ S noise. Simulation results demonstrate that the algorithm performs well at various levels of the heavy-tailed noise, and its estimation accuracy is high at low GSNR levels.