Second-order Time-Reassigned Synchrosqueezing Transform: Application to Draupner Wave Analysis

TL;DR

This paper introduces a second-order time-reassigned synchrosqueezing transform that improves time-frequency analysis of impulsive and modulated signals, demonstrated on synthetic data and pioneering analysis of Draupner wave records.

Contribution

It presents a novel second-order enhancement of the synchrosqueezing transform for better joint time-frequency representation of complex signals.

Findings

Enhanced time-frequency resolution for impulsive signals

Reversible and sharpened signal representations achieved

First detailed analysis of Draupner wave in time-frequency domain

Abstract

This paper addresses the problem of efficiently jointly representing a non-stationary multicomponent signal in time and frequency. We introduce a novel enhancement of the time-reassigned synchrosqueezing method designed to compute sharpened and reversible representations of impulsive or strongly modulated signals. After establishing theoretical relations of the new proposed method with our previous results, we illustrate in numerical experiments the improvement brought by our proposal when applied on both synthetic and real-world signals. Our experiments deal with an analysis of the Draupner wave record for which we provide pioneered time-frequency analysis results.

Tables

Table 1. Table 1: Signal reconstruction quality obtained for the reversible TFR s presented in Fig.  1 .

Method	RQF (dB)
STFT	269.27
classical synchrosqueezing	35.89
second-order vertical synchrosqueezing	23.80
time-reassigned synchrosqueezing	116.67
second-order time-reassigned synchrosqueezing	116.67

Equations

$$F_x^h(t,\omega )$$

$${\int }_{\mathbb{R}}{F}_x^h(t,\omega )dt$$

$$=\int {\int }_{{\mathbb{R}}^2}h(u{)}^{\ast }x(\tau ){\mathbf{e}}^{-j\omega \tau }dud\tau$$

$$={\int }_{\mathbb{R}}h(u{)}^{\ast }du{\int }_{\mathbb{R}}x(\tau ){\mathbf{e}}^{-j\omega \tau }d\tau$$

$$=F_h(0{)}^{\ast }{F}_x(\omega )$$

$$F_x(\omega )=\frac{1}{F_h(0{)}^{\ast }}{\int }_{\mathbb{R}}{F}_x^h(t,\omega )dt$$

$$x(t)=\frac{1}{2\pi F_h(0{)}^{\ast }}\int {\int }_{{\mathbb{R}}^2}{F}_x^h(\tau ,\omega ){\mathbf{e}}^{j\omega t}d\tau d\omega .$$

$${\hat{t}}_x(t,\omega )$$

$${\hat{\omega }}_x(t,\omega )$$

$$\int {\int }_{{\mathbb{R}}^2}|F_x^h(\tau ,\Omega ){|}^2\delta \left(t-{\hat{t}}_x(\tau ,\Omega )\right)\delta \left(\omega -{\hat{\omega }}_x(\tau ,\Omega )\right)\mathrm{d}\tau \mathrm{d}\Omega .$$

$$S_x^h(t,\omega )={\int }_{\mathbb{R}}{F}_x^h(\tau ,\omega )\delta \left(t-{\hat{t}}_x(\tau ,\omega )\right)d\tau$$

$${\int }_{\mathbb{R}}{S}_x^h(t,\omega )dt$$

$$={\int }_{\mathbb{R}}{F}_x^h(\tau ,\omega )d\tau =F_h(0{)}^{\ast }{F}_x(\omega ).$$

$$x(t)=\frac{1}{2\pi F_h(0{)}^{\ast }}\int {\int }_{{\mathbb{R}}^2}{S}_x^h(\tau ,\omega ){\mathbf{e}}^{j\omega t}d\tau d\omega .$$

$$x(t)$$

$$\text{with}{\lambda }_x(t)$$

$$\text{and}{\varphi }_x(t)$$

$${\omega }_x=\text{Im}({\tilde{\omega }}_x(t,\omega )-q_x{\tilde{t}}_x(t,\omega ))={\hat{\omega }}_x(t,\omega )-\text{Im}(q_x{\tilde{t}}_x(t,\omega ))$$

$$t_x^{(2)}=\frac{\omega -{\omega }_x}{{\alpha }_x}={\hat{t}}_x(t,\omega )+\frac{\omega -{\hat{\omega }}_x(t,\omega )}{{\alpha }_x}+\frac{{\nu }_x}{{\alpha }_x}\text{Im}({\tilde{t}}_x(t,\omega ))$$

$${\hat{t}}_x^{(2)}(t,\omega )=\{\begin{array}{ll}\frac{\omega -{\hat{\omega }}_x(t,\omega )+\text{Im}({\hat{q}}_x(t,\omega ){\tilde{t}}_x(t,\omega ))}{{\hat{\alpha }}_x(t,\omega )}& \mathrm{if}{\hat{\alpha }}_x(t,\omega )\ne 0\\ {\hat{t}}_x(t,\omega )& \mathrm{otherwise}\end{array}$$

$${\hat{t}}_x^{(2b)}(t,\omega )=\{\begin{array}{ll}{\hat{t}}_x(t,\omega )+\frac{\omega -{\hat{\omega }}_x(t,\omega )}{{\hat{\alpha }}_x(t,\omega )}& \mathrm{if}{\hat{\alpha }}_x(t,\omega )\ne 0\\ {\hat{t}}_x(t,\omega )& \mathrm{otherwise}\end{array}.$$

$${\hat{q}}_x^{(tn)}(t,\omega )$$

$${\hat{q}}_x^{(\omega n)}(t,\omega )$$

$$G(t)={\left({\int }_{\Omega }|S_x^h(t,\omega ){|}^{2}d\omega \right)}^{\frac{1}{2}}.$$

$$\hat{S}(t,\omega )=\{\begin{array}{ll}{S}_x^h(t,\omega )& \text{ if }G(t)>\Gamma \\ 0& \text{otherwise}\end{array}.$$

Full text

1 Introduction

Time-frequency and time-scale analysis [1, 2, 3, 4] aim at developing efficient and innovative methods to deal with non-stationary multicomponent signals. Among the common approaches, the Short-Time Fourier Transform (STFT) and the Continuous Wavelet Transform (CWT) [5] are the simplest linear transforms which have been intensively applied in various applications such as audio [6], biomedical [7], seismic or radar.

Unfortunately, these tools are limited by the Heisenberg-Gabor uncertainty principle. As a consequence, the resulting representations are blurred with a poor energy concentration and require a trade-off between the accuracy of the time or frequency localization. Another approach, the reassignment method [8, 9] was introduced as a mathematically elegant and efficient solution to improve the readability of a time-frequency representation (TFR). The inconvenience is that reassignment provides non-invertible TFRs which limits its interest to analysis or modeling applications.

More recently, synchrosqueezing [10, 11] was introduced as a variant of the reassignment technique due to its capability to provide sharpen and reversible TFRs. This reconstruction capability make this method continuously gaining interest since it paves the way of an infinite number of synchrosqueezing-based applications such as noise removal [12], signal components extraction or separation [10, 13, 14, 4].

Nowadays, efforts are made to efficiently compute the synchrosqueezed version of several linear transforms such as STFT, CWT or S-transform [15, 16] and to improve the localization of strongly modulated signals using enhanced instantaneous frequency estimators [17, 18, 19]. To deal with impulses and strongly modulated signals, a new variant of the synchrosqueezing was introduced and called time-reassigned synchrosqueezing method [20]. However, this method cannot efficiently deal with mixed-content signals containing both impulsive and periodic components.

In the present paper, we propose to introduce a novel transform called the second-order horizontal synchrosqueezing aiming to improve the energy localization and the readability of the time-reassigned synchrosqueezing while remaining reversible. To this end, we use an enhanced group-delay estimator which can be mathematically related to our previous results [18].

This paper is organized as follows. In Section 2, the proper definitions of the considered transforms with their properties are presented. In Section 3, we introduce our new second-order time-reassigned synchrosqueezing transform. Section 4 presents numerical experiments involving both synthetic and real-world signals. Finally, future work directions are given in Section 5.

2 Time-reassigned synchrosqueezing in a nutshell

2.1 Definitions and properties

We define the STFT of a signal $x$ as a function of time $t$ and frequency $\omega$ computed using a differentiable analysis window $h$ as:

[TABLE]

where $j^{2}\!\!=\!\!-1$ is the imaginary unit and $z^{*}$ is the complex conjugate of $z$. A TFR also called spectrogram is defined as $|F_{x}^{h}(t,\omega)|^{2}$. Thus, the marginalization over time of $F_{x}^{h}(t,\omega)$ leads to:

[TABLE]

with $F_{x}(\omega)\!=\!\int_{\mathbb{R}}x(t)\,\mathbf{e}^{-j\omega t}\,dt$ the Fourier transform of signal $x$. Now, from Eq. (5) one can compute the Fourier Transform of $x$ as:

[TABLE]

and the following signal reconstruction formula can be deduced after applying the Fourier inversion formula:

[TABLE]

2.2 Reassignment

To improve the readability of a TFR, reassignment moves the signal energy according to: $(t,\omega)\!\mapsto\!(\hat{t}_{x}{\scriptscriptstyle(t,\omega)},\hat{\omega}_{x}{\scriptscriptstyle(t,\omega)})$, where $\hat{t}_{x}(t,\omega)$ is a group-delay estimator and $\hat{\omega}_{x}(t,\omega)$ is an instantaneous frequency estimator [9]. Both time-frequency reassignment operators $\hat{t}$ and $\hat{\omega}$ can be computed as follows in the STFT case [21, 16]:

[TABLE]

where $\mathcal{T}h(t)\!=\!th(t)$ and $\mathcal{D}h(t)\!=\!\frac{dh}{dt}(t)$ are modified versions of the analysis window $h$.

Finally, a reassigned spectrogram can be computed as $\text{RF}^{h}_{x}(t,\omega)=$

[TABLE]

The resulting reassigned spectrogram $\text{RF}_{x}(t,\omega)$ is a sharpened but non-reversible TFR due to the loss of the phase information.

2.3 Time-reassigned synchrosqueezed STFT

To overcome the problem of non reversibility, synchrosqueezing proposes to move the signal transform instead of its energy, to preserve the phase information of the original transform.

Hence, time-reassigned synchrosqueezed STFT can be defined as [20]:

[TABLE]

where $\hat{t}_{x}(t,\omega)$ corresponds to the time reassignment operator which is classically computed using Eq. (8).

The marginalization over time of the resulting transform leads to:

[TABLE]

Hence, an exact signal reconstruction from its synchrosqueezed STFT can be deduced from Eq. (13) as:

[TABLE]

3 Second-order horizontal synchrosqueezing

3.1 Enhanced group-delay estimation

Let’s consider a linear chirp signal model expressed as [18]:

[TABLE]

where $\lambda_{x}(t)$ and $\phi_{x}(t)$ respectively stand for the log-amplitude and phase and with $q_{x}=\nu_{x}+j\alpha_{x}$ and $p_{x}=\mu_{x}+j\omega_{x}$. For such a signal, it can be shown [18] that $p_{x}=\tilde{\omega}_{x}{\scriptstyle(t,\omega)}-q_{x}\,\tilde{t}_{x}{\scriptstyle(t,\omega)}$, and therefore:

[TABLE]

The proposed second-order horizontal synchrosqueezing consists in moving $F_{x}^{h}{\scriptstyle(t,\omega)}$ from the point ${\scriptstyle(t,\omega)}$ to the point ${\scriptstyle(t_{x}^{(2)},\omega)}$ located on the instantaneous frequency curve, i.e. such that $\dot{\phi}(t_{x}^{(2)})\!=\!\frac{d\phi_{x}}{dt}(t_{x}^{(2)})=\omega_{x}+\alpha_{x}t_{x}^{(2)}=\omega$. This leads to:

[TABLE]

which can be estimated by:

[TABLE]

where $\hat{q}_{x}{\scriptstyle(t,\omega)}=\hat{\nu}_{x}{\scriptstyle(t,\omega)}+j\hat{\alpha}_{x}{\scriptstyle(t,\omega)}$ is an unbiased estimator of $q_{x}$. This expression can be compared to the second-order group delay estimator introduced by Oberlin et al. [17]

[TABLE]

It can be shown using Eq.(19) that estimator $\hat{t}^{(2b)}_{x}$ is biased when $\nu_{x}\!=\!\frac{d^{2}\lambda_{x}}{dt^{2}}(t)\neq 0$.

Finally, a new second-order horizontal synchrosqueezing transform can thus be obtained using Eq.(11) by replacing the group-delay estimator $\hat{t}(t,\omega)$ by our enhanced estimator given by Eq.(20).

3.2 Theoretical considerations and computation issue

In [18, 22] we introduced two families of unbiased estimators called $(tn)$ and $(\omega n)$ involving $n$-order derivatives ($n\!\geq\!2$) with respect to time (resp. to frequency) which enable to compute Eqs. (20) and (21):

[TABLE]

with $\mathcal{D}^{n}h(t)=\frac{d^{n}h}{dt^{n}}(t)$ and $\mathcal{T}^{n}h(t)=t^{n}h(t)$. Our preliminary investigations [22] showed a slight improvement using the $(\omega 2)$ estimator in terms of accuracy in comparison to higher-order and $(tn)$ estimators.

Our implementations use the discrete-time reformulations of our previously described expressions combined with the rectangle approximation method. Thus $F_{x}^{h}[k,m]\!\approx\!F_{x}^{h}(\frac{k}{F_{s}},2\pi\frac{mF_{s}}{M})$, where $F_{s}$ denotes the sampling frequency, $k\in\mathbb{Z}$ is the time sample index and $m\in\mathcal{M}$ is the discrete frequency bin. The number of frequency bins $M$ is chosen as an even number such as $\mathcal{M}=\left[-M/2+1;M/2\right]$. It results that our method has the same computational complexity of the previously introduced second-order vertical synchrosqueezing.

The proposed method is valid for any differentiable analysis window. In our implementation111matlab code freely available at: http://www.fourer.fr/hsst, the STFT uses a Gaussian window and is also called Gabor transform. The window function is expressed as $h(t)\!=\!\frac{1}{\sqrt{2\pi}T}\,\mathbf{e}^{-\frac{t^{2}}{2T^{2}}}$ where $T$ is the time-spread of the window which can be related to $L=TF_{s}$.

4 Numerical results

4.1 Analysis of a synthetic signal

In this experiment, we consider a synthetic $500$-sample-long multicomponent real-valued signal made of two impulses, one sinusoid, one chirp and one sinusoidally modulated sinusoid. Fig. 1 compares the following TFRs: spectrogram, reassigned spectrogram, classical (frequency-reassigned) synchrosqueezing, second-order vertical synchrosqueezing, time-reassigned synchrosqueezing and second-order time-reassigned horizontal synchrosqueezing. Our computations use $M\!=\!600$, $L\!=\!8$ and a Signal-to-Noise Ratio (SNR) equal to 25 dB obtained by the addition of a Gaussian white noise. The local modulation estimator $\hat{q}_{x}^{(\omega 2)}$ is used for computing both second-order synchrosqueezing methods. The TFRs provided by the previously proposed methods are computed using the matlab implementions provided by the ASTRES toolbox [16].

The results clearly illustrate the improvement of the new second-order time-reassigned synchrosqueezing over the time-reassigned synchrosqueezing for representing the whole signal. When compared with frequency-reassigned synchrosqueezing methods, our new method has the advantage to perfectly localize the two impulses while providing a sharpened representation of the chirp and of the sinusoidally modulated sinusoid. Unfortunately, as for the time-reassigned synchrosqueezing, our method cannot localize the non-modulated sinusoid.

To assess the signal reconstruction capability, we compare in Table 1 the Reconstruction Quality Factor (RQF) of each TFR computed using [15]: $\text{RQF}=10\log_{10}\left(\frac{\sum_{n}|x[n]|^{2}}{\sum_{n}|x[n]-\hat{x}[n]|^{2}}\right)$. Thus, our results show again the advantage of the time-reassigned synchrosqueezing methods which obtain significantly higher RQF (if $M$ is chosen at least equal to the signal length) due to its theoretically exact reconstruction formula.

4.2 Draupner wave signal analysis

Now we consider a record of a possible freak wave event measured on the Draupner Platform in 1995 [23]. The signal displayed in Fig. 2(a) corresponds to the sea surface elevation deduced from the measures provided by a wave sensors consisting of a down-looking laser. The sampling frequency of this signal is $F_{s}=2.13$ Hz and its duration is 20 minutes.

4.2.1 Time-frequency representation

Fig. 2 compares the resulting TFRs provided by the STFT and its different synchrosqueezed versions (i.e. all the combination of the first- and second-order of the frequency-reassigned and time-reassigned). For our numerical results, we empirically choose $M=2660$ and $L=25$ which provide sufficiently readable results. In order to focus to the impulsive part of the signal, we have limited the analysis between 0.2 Hz and 1 Hz. As expected, the second-order time-reassigned synchrosqueezing provides the best representation to localize the 4 impulses visible in the signal. Interestingly, our results reveal the main impulse located at $t_{1}\approx 4.39$ min (also visible in Fig. 2(a)) but also 3 supplementary impulses respectively located at $t_{2}\approx 7.72$ min, $t_{3}\approx 13.36$ min and $t_{4}\approx 19.47$ min. These impulses were almost invisible in the waveform representation of the signal but have been revealed by our proposed time-frequency analysis methods.

4.2.2 Impulses detection and disentangling

Now we propose to use the synchrosqueezing signal reconstruction capability for recovering the 4 impulse signals. To this end, we compute a saliency function defined as the root mean square of the marginal over frequency band $\Omega=[0.4;1]$ Hz of the signal energy contained in its synchrosqueezing transform:

[TABLE]

A binary masked version of the transform $S_{x}^{h}(t,\omega)$ can thus be computed using $G(t)$ as:

[TABLE]

where $\Gamma$ is a defined threshold. Finally the components are extracted by applying the reconstruction formula given by Eq.(14) on $\hat{S}(t,\omega)$. Our numerical computation presented in Fig. 3 uses $\Gamma=3.37$ which corresponds to 5 times the mean value of $G(t)$. It allows us to recover the impulses locations through a peak picking and to reconstruct the corresponding waveform signal initially merged in the whole signal.

5 Conclusion and future work

A new extension of the time-reassigned synchrosqueezing called second-order horizontal synchrosqueezing was introduced for the STFT. Our experiments show a significant improvement to compute invertible and sharpened time-frequency representations of impulsive signals which cannot be addressed by vertical synchrosqueezing. Moreover, we have shown the efficiency of this technique when applied on both synthetic and real-word signals. In fact, our method helped to discover new signal components in the Draupner wave signal which could probably help to better understand the phenomenon of freak waves. Future work consist in theoretically strengthening this method, and developing new applications.

Acknowledgment

The authors would like to thank Dr. Sverre Haver who kindly agreed to answer our questions and share with us the Draupner wave signal record.