# A two-step synchrosqueezing transform approach mitigating reallocation errors in time-frequency analysis

F Karkafi, M Leiber, Q Leclère, J Antoni, M El Badaoui

► To cite this version:

F Karkafi, M Leiber, Q Leclère, J Antoni, M El Badaoui. A two-step synchrosqueezing transform approach mitigating reallocation errors in time-frequency analysis. ISMA2024 Conference on Noise and Vibration Engineering, Sep 2024, Leuven, Belgium. hal-05303620

HAL Id: hal-05303620

https://hal.science/hal-05303620v1

Submitted on 8 Oct 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# A two-step synchrosqueezing transform approach mitigating reallocation errors in time-frequency analysis

F. Karkafi¹,², M. Leiber¹, Q. Leclère², J. Antoni², M. El Badaoui¹,³

¹ Safran Tech,
Rue des Jeunes Bois – Châteaufort 78772 Magny-les-Hameaux, France

² INSA Lyon, LVA, UR677,
69621 Villeurbanne, France

³ UJM-St-Etienne, LASPI, UR3059, F-42023,
Saint-Etienne, France

## Abstract

Reallocation errors are crucial in time-frequency analysis, affecting the interpretation of modulated components. The first order synchrosqueezing transform (SST-1) reallocates signal energy around instantaneous frequencies but introduces reallocation errors. The second order SST (SST-2) rectifies these errors but amplifies noise in practical applications. This paper proposes a two-step method to address this issue. First, a differentiable short-time Fourier transform optimizes window lengths via gradient descent, dynamically determining optimal lengths based on predefined criteria. This pre-processing step corrects the local group delay. Then, SST-1 corrects frequency reallocation errors. This combined approach achieves the same correction as SST-2 without noise amplification. The methodology is validated by estimating the instantaneous angular speed of a rotating machine, tested on both simulated and experimental data. The results demonstrate the approach's effectiveness in mitigating reallocation errors and overcoming phase noise challenges.

## 1 Introduction

In contemporary industrial applications, vibration analysis is vital for the condition monitoring (CM) of rotating machinery [1, 2]. However, a pressing challenge has commanded considerable attention over the past two decades: the presence of speed fluctuations within recorded vibration signals [3, 4, 5]. Analyzing or processing these signals efficiently requires the use of local transformations, called time-frequency representations (TFRs). One of the most popular TFRs is the short-time Fourier transform (STFT), which essentially applies a window to the signal around the time of interest before performing the discrete Fourier transform (DFT). While the STFT is invertible and suitable for processing non-stationary signals, it has a major drawback: the time-frequency (TF) resolution is limited by the choice of the window size. This trade-off is dictated by the Heisenberg-Gabor uncertainty principle [6]. A small temporal window ideal for precise time localization, leads to poor frequency resolution, and vice versa. Additionally, the fixed window size is unsuitable for signals with varying regimes, where varying characteristics within the same signal require different resolutions. As a result, the STFT struggles to adapt and provide clear and readable TFRs for such complex signals. Many works in the past decades have tackled this limitation by using quadratic TFRs, such as the Wigner-Ville distribution [7], which offer high resolution in both time and frequency domains but exhibit strong interference, hampering the representation.

The synchrosqueezing transform (SST) is yet another approach, whose initial goal is to improve the readabil- ity of the TFR using a reassignment procedure. This improvement allows for the retrieval of the modes of a multicomponent signal (MCS) from the reassigned transform [8], which is of great practical interest. SST has undergone numerous extensions: it has been adapted to the STFT, resulting in the STFT-based SST [9], with
mathematical analysis provided in [10, 11, 12]. Additionally, a multi-taper approach has been proposed to enhance both the resolution and robustness of SST, as seen in the ConceFT method [13], and the analysis of non-harmonic waves has been explored in [14]. Despite these advancements, one major limitation associated with SST is its inability to handle MCSs with modes exhibiting strong frequency modulation, which are very common in many fields of practical interest. To address the limitations in handling signals with significant frequency modulation, an adaptation of SST, termed the second-order SST (SST-2), was introduced in [12], with its theoretical foundations established in subsequent work [15]. This modified transform has been utilized in demodulation algorithms [16] and extended to handle modes with rapidly oscillating phases [17]. The SST-2 method calculates the ratio of gradients of both the channelized instantaneous frequencies and the local group delays. However, SST-2 assumes a clean phase component, which is not typical in complex industrial systems such as rotating machinery and aircraft engines. Consequently, the computation of gradients ratio tends to amplify unaccounted phase noise, often even at a stronger level than the reallocations errors affecting the first order SST (SST-1). To improve STFT-based representations, there has been interest in finding optimal parameters, such as the window length [18, 19]. One approach is the differentiable short-time Fourier transform (DSTFT) [20, 21, 22]. It enables optimizing STFT parameters using gradient descent based on user-defined criteria. This paper utilizes the DSTFT by employing a time-varying sliding window length to correct local group delay (time reallocation) concerning underlying amplitude and frequency modulation. The optimization aims to improve the signal-to-noise ratio (SNR) of target frequencies while attenuating unwanted ones. Following the optimization, SST-1 is applied to rectify frequency reallocation errors. The combined output of these two steps achieves a similar correction to SST-2 but avoids noise amplification. Finally, the assessment is done by comparing instantaneous frequency (IF) estimation accuracy. Section 2 delves into the theoretical foundations underpinning the frequency reallocation framework, followed by a detailed exposition of the differentiable optimization of the STFT in Section 3. The evaluation in Section 4 encompasses both simulated and real-world experimental data, assessing the performance of the proposed method in the IF estimation case.

## 2 Frequency Reallocation

Under the assumption that a rotating machinery vibration signal can be represented as the sum of AM/FM components, $x(t)$ can be expressed as:

$$x(t) = \sum_{k} A_{k}(t)e^{j\alpha_{k}\phi(t)} + \nu(t), \qquad (1)$$

where $\alpha_k$ represents the $k^{th}$ harmonic order, $A_k(t)$ the instantaneous amplitude of harmonic $k$, $\nu(t)$ is an added noise of variance $\sigma_\nu^2$ such that the SNR of the $k^{th}$ harmonic would be $A_k^2(t)/2\sigma_\nu^2$, $\phi(t)$ the instantaneous angle of rotation expressed in [rad], and its derivative with respect to the time $\omega(t) = \frac{1}{2\pi} \frac{d\phi(t)}{dt}$ is the instantaneous angular speed in [Hz]. To address the requirement for time-localized frequency descriptors, the STFT of $x(t)$ is employed and computed by convolving the signal with a sliding window function $g$, defined as:

$$S_x^g(\tau, f) = \int_{\mathbb{R}} x(t)g(t - \tau)e^{-2\pi j f (t - \tau)} d\tau, \qquad (2)$$

where $\tau$ represents the tims shift and $f$ is the frequency variable. The STFT squared magnitude, $|S_x^g(\tau, f)|^2$, when represented in the TF plane, is referred to as the spectrogram of the signal $x(t)$. The frequency resolution of the STFT is determined by the window duration: longer windows provide better frequency resolution $\Delta f$ but poorer time resolution $\Delta \tau$, and vice versa. This trade-off is governed by the Heisenberg-Gabor uncertainty principle, which states that $\Delta \tau \Delta f \ge \frac{1}{4\pi}$. This principle highlights the fundamental limitation in simultaneously achieving high resolution in both time and frequency domains.
## 2.1 First Order SST

To refine the precision of IF representations within signal components, the concept of frequency reallocation is introduced. This involves redistributing frequency components within the signal's TFR. The STFT filters the signal across frequency bands, decomposing it into multiple frequency component signals. The IFs of the signal components can be calculated using Rihaczek's equations [23]. For a more in-depth exploration of this concept, Nelson's work [24] is a useful reference, where the term channelized instantaneous frequency (CIF) is used to describe the simultaneous IFs computed within each time instant of the STFT. The STFT, $S_x^g(\tau, f)$, is a complex quantity characterized by its slowly varying magnitude $X(\tau, f)$ and phase $\Phi(\tau, f)$, defined as:

$$S_x^g(\tau, f) = X(\tau, f)e^{j\Phi(\tau, f)}. \qquad (3)$$

Ultimately, the channelized IF $\hat{\omega}_x(\tau, f)$ can be deduced from Eq (3) as the time derivative of the phase $\frac{1}{2\pi}\partial_\tau \Phi(\tau, f)$ where $\partial_a Z(a, b) = \frac{\partial Z(a, b)}{\partial a}$. Therefore, $\hat{\omega}_x(\tau, f)$ can be computed as:

$$\hat{\omega}_x(\tau, f) = \operatorname{Im} \left\{ \frac{\partial_\tau S_x^g(\tau, f)}{2\pi S_x^g(\tau, f)} \right\}. \qquad (4)$$

Now that the STFT frequency operator is defined, one can conduct the signal phase $\phi(t)$ study to compute the CIF corresponding to $\phi(t)$ to reallocate the signal frequency components. Referring to Eqs (1) and (2), and assuming that the signal envelope $A_k(t)$ is mostly constant over the time interval spanned by $g(t - \tau)$, $S_x^g(\tau, f)$ can be expressed as

$$S_x^g(\tau, f) = A_k(\tau) \int e^{j\phi(t)} g(t - \tau) e^{-j2\pi f(t-\tau)} dt. \qquad (5)$$

The phase $\phi(t)$ can be expanded in a Taylor series around $\tau$. Considering the case of SST-1, the transform assumes that the IF is slowly changing, such that the second-order derivative $\phi''(\tau) \approx 0$. Thus, the first-order phase is defined as

$$\phi(t) = \phi(\tau) + (t - \tau)\phi'(\tau). \qquad (6)$$

Regardless of the considered order, the goal is to estimate the CIF expressed in Eq (4). Therefore, after integrating the first-order phase from Eq (6) in Eq (5), and applying the variable change $v = t - \tau$, one can conclude that $\operatorname{Im}\left\{\frac{\partial_\tau S_x^g(\tau, f)}{2\pi S_x^g(\tau, f)}\right\} = \frac{\phi'(\tau)}{2\pi}$, representing the first-order CIF (CIF-1) $\hat{\omega}_x^{[1]}(\tau, f)$. To effectively estimate $\hat{\omega}_x^{[1]}(\tau, f)$, referring to Eq (2), $\partial_\tau S_x^g$ can be expressed as:

$$\partial_\tau S_x^g = 2\pi j f S_x^g(\tau, f) - S_x^{g'}(\tau, f), \qquad (7)$$

where $g'$ is to the time-derivative of the sliding window. Consequently, by replacing $\partial_\tau S_x^g(\tau, f)$ in Eq (4), the speed operator takes the following expression:

$$\hat{\omega}_x^{[1]}(\tau, f) = f - \frac{1}{2\pi} \operatorname{Im} \left\{ \frac{S_x^{g'}(\tau, f)}{S_x^g(\tau, f)} \right\}. \qquad (8)$$

After computing the signal CIF-1 $\hat{\omega}_x^{[1]}(\tau, f)$, the synchrosqueezing process commences. The STFT-based SST of the signal $x$ is essentially defined by the threshold $\gamma$. This parameter enables the adaptation of the SST to the specific characteristics of the signal and the analysis requirements. The synchrosqueezing process is mathematically defined as follows:
$$ \tilde{T}_x^\gamma(\tau, \omega) = \int_{|S_x^g(\tau, f)|>\gamma} S_x^g(\tau, f) \delta\left(\omega - \hat{\omega}_x^{[1]}(\tau, f)\right) df, \quad (9) $$

where $\delta$ denotes the Dirac distribution. The SST formulation, theoretical foundation and further definitions were settled in [25].

## 2.2 Second Order SST

The applicability of SST-1 is restricted to a class of MCSs composed of slightly perturbed purely harmonic modes, also known as small speed fluctuations. To overcome this limitation, an extension was introduced based on a more accurate IF estimate, called SST-2 [12, 25]. To achieve this, one must consider the speed fluctuations within the computed CIF-1 and refine its estimation. This can be done by adding the previously ignored speed fluctuations $\phi''(\tau)$ such that

$$ \phi(t) = \phi(\tau) + (t - \tau)\phi'(\tau) + \frac{(t - \tau)^2}{2}\phi''(\tau). \quad (10) $$

Again, this means that the third order term $\phi'''(\tau)$ is neglected assuming that $\phi''(\tau)$ is slowly varying. After integrating the second-order phase from Eq (10) in Eq (5), $S_x^g(\tau, f)$ becomes

$$ S_x^g(\tau, f) = A_k(\tau) \int e^{j\left(\phi(\tau) + (t-\tau)\phi'(\tau) + \frac{(t-\tau)^2}{2}\phi''(\tau)\right)} g(t-\tau)e^{-j2\pi f(t-\tau)}dt. \quad (11) $$

On the one hand, computing the time derivative $\partial_\tau S_x^g$ can be expressed as:

$$ \partial_\tau S_x^g(\tau, f) = j\phi'(\tau)S_x^g(\tau, f) + A_k(\tau)e^{j\phi(\tau)}j\phi''(\tau)I(\tau, f), \quad (12) $$

where $I(\tau, f) = \int vg(v)e^{jv(\phi'(\tau)+\frac{v}{2}\phi''(\tau)-2\pi f)}dv$ with $v = t - \tau$. On the other hand, computing the frequency derivative $\partial_f S_x^g$ is expressed as

$$ \partial_f S_x^g(\tau, f) = -2\pi j A_k(\tau) e^{j\phi(\tau)} I(\tau, f). \quad (13) $$

Therefore, defining $I(\tau, f)$ in function of $\partial_f S_x^g(\tau, f)$, substituting it in Eq (12) and taking the imaginary part results in

$$ \mathrm{Im} \left\{ \frac{\partial_\tau S_x^g(\tau, f)}{S_x^g(\tau, f)} \right\} = \phi'(\tau) - \phi''(\tau) \mathrm{Im} \left\{ \frac{\partial_f S_x^g(\tau, f)}{2\pi S_x^g(\tau, f)} \right\}. \quad (14) $$

An analogous relation applies to the STFT: each frequency component at a specific time can be treated as a frequency “signal” whose group delay can be computed using Rihaczek’s equations [23], thus yielding new corrected event times $\hat{\tau}_x(\tau, f)$ for each time-frequency cell. Following insights from Nelson [24] and based on Eqs (2) and (3), many studies define the local group delay (LGD) of a signal as $\hat{\tau}_x(\tau, f) = \tau + \frac{1}{2\pi}\partial_f\Phi(\tau, f)$ [26, 27]. Therefore, referring to the same equations, $\hat{\tau}_x(\tau, f)$ can be computed as [26, 27, 28, 29]:

$$ \hat{\tau}_x(\tau, f) = \tau + \mathrm{Im} \left\{ \frac{\partial_f S_x^g(\tau, f)}{2\pi S_x^g(\tau, f)} \right\}. \quad (15) $$

To effectively estimate $\hat{\tau}_x(\tau, f)$, referring to Eq (2), $\partial_f S_x^g(\tau, f)$ can be expressed as:

$$ \partial_f S_x^g = -2\pi j \int x(t)(t - \tau)g(t - \tau)e^{-2\pi j f(t - \tau)}dt = -2\pi j S_x^{\tau g}(\tau, f). \quad (16) $$
Therefore, the LGD can be concluded as

$$ \hat{\tau}_x(\tau, f) = \tau - \text{Re} \left\{ \frac{S_x^{\tau g}(\tau, f)}{S_x^g(\tau, f)} \right\}, \qquad (17) $$

and replacing the LGD expression from Eq (15) in Eq (14) results in

$$ \hat{\omega}_x^{[2]}(\tau, f) = \phi'(\tau) + \phi''(\tau)(\tau - \hat{\tau}_x(\tau, f)). \qquad (18) $$

As seen in the first order case, $\phi'(\tau)$ is none other than CIF-1 $\hat{\omega}_x^{[1]}(\tau, f)$. In this context, $\phi''(\tau)$ is the signal speed fluctuations that are represented by the time-difference of CIF-1 within their LGD. Thus, the frequency modulation operator $\hat{q}(\tau, f)$ corresponding to $\phi''(\tau)$ is defined as follows [12, 29]:

$$ \hat{q}(\tau, f) = \frac{\partial_{\tau} \hat{\omega}_{x}^{[1]}(\tau, f)}{\partial_{\tau} \hat{\tau}_{x}(\tau, f)}. \qquad (19) $$

Consequently, whenever the LGD is varying along time, CIF-2 refines CIF-1 as follows

$$ \hat{\omega}_x^{[2]}(\tau, f) = \begin{cases} \hat{\omega}_x^{[1]}(\tau, f) + \hat{q}(\tau, f) \times (\tau - \hat{\tau}_x(\tau, f)) & \text{if } \partial_{\tau} \hat{\tau}_x(\tau, f) \neq 0 \\ \hat{\omega}_x^{[1]}(\tau, f) & \text{otherwise.} \end{cases} \qquad (20) $$

The SST-2 formulation, theoretical foundation, and further definitions were also settled in [25]. This formulation was extended in other papers to higher orders by expanding the Taylor series even more. However, it is valid only if the phase $\Phi(\tau, f)$ is accurately extracted. Referring to Eq (3), if $\Phi(\tau, f)$ contains significant random noise, the computation of $\hat{q}(\tau, f)$ via the ratio of both frequency and delay operators differences will amplify the noise, often at a stronger level than the reallocation limitations of SST-1. Therefore, this paper considers correcting the reallocation errors by adapting the window length through differentiable STFT optimization.

## 3 Differentiable STFT

The differentiable STFT (DSTFT) has been introduced in recent studies [20, 21, 22]. It extends the STFT by enabling gradient-based optimization of its parameters, such as the window and hop lengths. This paper focuses exclusively on optimizing the window length. The core concept of DSTFT lies in decomposing the window length $M$, into an integer window support $N$ (representing the upper bound $M \le N$) and a continuous time resolution $\Theta$, which is differentiable and learnable [20]. Consequently, a digital transform is necessary to match this analysis. Adaptive DSTFT can be defined by employing a continuous window length for each frame $i$ of the STFT [21]:

$$ S_x^{g_N}[i, k] = \frac{1}{F_s} \sum_{n=0}^{N-1} x[n] g_N \left( \frac{n}{F_s} - \zeta_i, \Theta_i \right) e^{-2\pi j \eta_k \frac{n}{F_s}}, \qquad (21) $$

with $n \in \mathbb{Z}$ the sample index, $\zeta$ is the discrete time-shift parameter ranging from $\frac{T}{2}$ to $(\frac{N}{F_s} - \frac{T}{2})$ with a time step of $\Delta\zeta = T(1 - \gamma_{overlap})$ such that $T = \frac{M}{F_s}$ is the window duration, $M$ is the window length, $F_s$ is the sampling frequency and $\gamma_{overlap}$ is the overlap factor, $i \in \mathbb{Z}$ corresponds to the $i^{th}$ sliding window starting at $\zeta_i$, $\eta_k \in \mathcal{F}$ the discrete frequency of the DFT where $\eta_k = k\Delta\eta$, $k \in \mathbb{Z}$ being an integer increment ranging from 0 to $\lfloor \frac{F_s}{2\Delta\eta} \rfloor$, the frequency step $\Delta\eta = \frac{1}{Tz}$ where $z$ is the zero-padding coefficient and $\mathcal{F}$ represents the grid of frequency bins, $g_N(\zeta, \Theta_i)$ represents a tapering function that is continuous and differentiable with respect to both of its arguments $(\zeta, \Theta)$ and possesses compact support. This means that $g_N(\zeta, \Theta) = 0$ for all $\zeta \notin [0, \Theta]$. In this work, the analysis window is defined by the Hann function expressed as:
$$g(\varsigma, \Theta) = \begin{cases} \frac{1}{2} \left(1 - \cos\left(\frac{2\pi\varsigma}{\Theta}\right)\right) & \text{for } 0 \le \varsigma \le \Theta \\ 0 & \text{else.} \end{cases} \qquad (22)$$

As a result, the SST-1 formulation in (9) becomes:

$$\tilde{T}_x^\gamma[i, l] = \sum_{k: \{|\hat{\omega}_x^{[1]}[i,k] - \eta_l| \le \frac{\Delta\eta}{2} \cap |S_x^{gN}[i,k]| > \gamma\}} S_x^{gN}[i, k], \qquad (23)$$

where $l \in \mathbb{Z}$ being an integer increment ranging from 0 to $\lfloor \frac{F_s}{2\Delta\eta} \rfloor$. The performance of this synchrosqueezed representation (SSR) will be assessed by extracting the IF. In this paper, the methodology involves optimizing the STFT window lengths via gradient descent based on a given criterion. Specifically, the goal is to adapt the window lengths so they are well-distributed over the entire signal, ensuring the capture of all varying regimes. This adaptation aims to minimize energy leakage while achieving better concentration of energy in the temporal plane. To achieve this in each frame, the criterion is set to maximize the empirical kurtosis of each frame $i$:

$$\mathcal{K}(\Theta_i) = \frac{\sum_k |S_x^{gN, \Theta_i}[i, k]|^4}{\left(\sum_k |S_x^{gN, \Theta_i}[i, k]|^2\right)^2}. \qquad (24)$$

# 4 Results

This section presents a comprehensive assessment of the proposed methodology through the analysis of both simulated and experimental data. Considering diverse contexts and challenging scenarios offers a comprehensive perspective on the adaptability of the proposed methodology. The performance of the approach is evaluated by highlighting its IF estimation accuracy and comparing it against two other methods. The evaluation scheme includes:

1. Fixed window length STFT + SST-1

2. Fixed window length STFT + SST-2

3. Proposed methodology: DSTFT + SST-1

## 4.1 Evaluation on Simulated Data

To evaluate the performance of the proposed methodology, a simulation example is first presented. The nonstationary signal was generated with a non-constant speed profile varying between 10 and 60 Hz. The signal was simulated with a sampling frequency of 8 kHz, and a total length of $N = 160,000$ samples, resulting in a virtual acquisition time of 20 seconds. The signal $x[n]$ consists of multiple harmonics, including four key components $\{\alpha_1, \alpha_2, \alpha_3, \alpha_4\}$: one with an order of 1, a second one that appears at 7.5, a third one at 15, and a fourth one at 25. These components were intentionally selected to avoid interference with the first harmonic. To further emulate complex conditions, white Gaussian noise was introduced for each harmonic phase with corresponding variances $\nu_{\phi_1}^2 = \nu_{\phi_2}^2 = 1$, $\nu_{\phi_3}^2 = 5$ and $\nu_{\phi_4}^2 = 10$. Moreover, a background white Gaussian noise $\nu$ was added to pollute all the harmonics within the signal with a SNR of 5 dB. The amplitudes of all harmonic components were unitary, contributing to the overall characteristic of the simulated signal. The signal was finally convolved with a second-order impulse response $h[n]$ to encompass the effect of a resonant transmission path. Thus the final form:
$$x[n] = \left( \sum_{k=1}^{4} A_k[n] \cos(\alpha_k \times \theta[n] + \phi_k[n]) + \nu[n] \right) \otimes h[n], \quad (25)$$

where $n = 0, 1, ..., N - 1$ is the sample index, $k$ is the harmonic number, $A_k[n]$ is the amplitude of harmonic $k$, $\theta[n]$ is the varying phase, $\alpha_k$ is the harmonic order and the symbol $\otimes$ is used to describe the convolution. Firstly, DSTFT is initiated as a standard STFT using a sliding Hann window of length 800, an overlap factor of 80% and $N_{FFT} = 2^{13}$. After 100 iterations, DSTFT rectified the window length distribution by maximizing the kurtosis of each frame's spectrum accordingly. A zoom of the resulting spectrograms and window lengths distribution are visualized in figure 1.

<figure><img src="image_2.png" /><figcaption>Figure 1: Comparison of spectrograms and window length distributions obtained using different methods. Top left: Zoomed-in spectrogram computed using the mean of window length distribution (bottom left). Top right: Zoomed-in spectrogram resulting from DSTFT optimization using window length distribution (bottom right) found with DSTFT (dashed black lines) overlapped with a scaled version of the actual $1/\tilde{\theta}^2$ (green lines) and smoothed version of it (purple lines).</figcaption></figure>

The standard spectrogram is computed using the mean of the resulting window lengths distribution. In the same figure, it is shown that the window length distribution resulting from the criterion optimization closely follows a scaled version of $1/\tilde{\theta}^2$, where $\tilde{\theta}$ represents the actual angular acceleration. Due to the challenge of capturing the smallest acceleration fluctuations, which requires much smaller windows (ultimately for each instant of time), a smoother version of $1/\tilde{\theta}^2$ is presented in the figure. This smoother version shows that the window distribution aligns with a slower variation of the true $1/\tilde{\theta}^2$. On the one hand, this representation reflects the maximum feasible window lengths change attempted during the simulation. On the other hand, it demonstrates the effectiveness of the optimization process in adapting to the signal's overall varying regimes. The SSRs of the three methods are visualized in figure 2. It can be seen that, despite of its attempt to correct reallocation errors, SST-2 seems to amplify the phase noise. Additionally, the zoomed-in version of each SSR illustrates the capability of the proposed methodology to accurately track both low and high-speed fluctuations. Figure 3 depicts the estimated speed obtained by computing the expected value of each SSR within the frequency band [10, 65]Hz. Zoomed-in versions showcase both stationary and fast-varying segments of the speed, highlighting the adaptive nature of the proposed methodology. For a comprehensive evaluation of each method's performance, mean and median absolute errors are presented in Figure 4. The histogram demonstrates that the proposed methodology surpasses the performance of other conventional
methods.

<figure><img src="image_3.png" /><figcaption>Figure 2: Illustrations of the SSRs obtained from three different methods. Top row: SSRs corresponding to STFT+SST-1 (left), STFT+SST-2 (middle), and DSTFT+SST-1 (right). Zoomed-in versions focusing on the time interval [0, 3]s are provided below each SSR for detailed analysis.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Figure 3: Comparison of estimated IF using different methods with a zoomed-in version focusing on the time intervals [0.6, 1.6]s and [8.9, 9.9]s.</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>Figure 4: Mean and median absolute errors for each method of the simulated data with respect to the actual speed.</figcaption></figure>

## 4.2 Evaluation on Experimental Data

The vibration analysis case originates from the Safran contest at the Surveillance 8 conference as documented in Ref [30]. This dataset comprises measurements of vibrations and tachometer readings obtained during a ground test campaign on a civil aircraft engine with two bearings exhibiting damage. In Figure 5, one can identify the aircraft engine under observation, along with the identified locations of the affected bearings and the positions of the sensors installed for data collection and the gearbox configuration. It encompasses two primary shafts and an accessory gearbox containing components such as filters, alternators, pumps, and a starter. The high-pressure shaft (HP) is connected to the accessory gearbox via a radial drive shaft RDS and a horizontal drive shaft HDS.

<figure><img src="image_6.png" /><figcaption>Figure 5: General overview of the engine and the accessory gearbox kinematics.</figcaption></figure>

Similar to the simulation case, DSTFT is initiated as a standard STFT with a sliding Hann window of length 4000, an overlap factor of 80%, and $N_{FFT} = 2^{16}$. Following 100 iterations, DSTFT rectifies the window length distribution by maximizing the kurtosis of each frame's spectrum. A detailed view of the resulting spectrograms and window length distribution is presented in Figure 6.

Notably, the standard spectrogram is computed using the mean of the resulting window length distribution. An observation reveals that window lengths increase during quasi-stationary regimes and decrease during fast-varying ones. The results of the three comparative methods are illustrated in Figure 7. This figure
<figure><img src="image_7.png" /><figcaption>Figure 6: Comparison of spectrograms and window length distributions obtained using different methods. Top left: Zoomed-in spectrogram computed using the mean of window length distribution (bottom left). Top right: Zoomed-in spectrogram resulting from DSTFT optimization using window length distribution found with DSTFT (bottom right).</figcaption></figure>

displays the estimated speed derived from computing the expected value of each SSR within the frequency band [175, 250]Hz.

<figure><img src="image_8.png" /><figcaption>Figure 7: Comparison of estimated IF using different methods with a zoomed-in version focusing on the time intervals [10, 20]s and [45.5, 46.5]s.</figcaption></figure>
Zoomed-in versions provide a detailed view of both stationary and fast-varying segments of the speed, thereby emphasizing the adaptive nature of the proposed methodology. To comprehensively evaluate the performance of each method, mean and median absolute errors are showcased in Figure 8. The histogram underscores that the proposed methodology outperforms the tested ones.

<figure><img src="image_9.png" /><figcaption>Figure 8: Mean and median absolute errors for each method of the experimental data with respect to the actual speed.</figcaption></figure>

# 5 Conclusion

This paper introduces an innovative approach to mitigate reallocation errors in time-frequency analysis by combining the SST-1 with the DSTFT. The primary goal is to address the challenges posed by varying regimes in signals, particularly in the presence of phase noise, without the amplification issues commonly observed in SST-2. The proposed approach using gradient-based optimization for optimal parameters demonstrated effectiveness through both simulated and experimental data. In the simulation case, the methodology showcased its ability to accurately track low and fast-varying speed fluctuations and provided a more concentrated TFR. The experimental case, involving vibration analysis of a civil aircraft engine, further validated the methodology's robustness and adaptability. The results consistently showed that the proposed method outperformed traditional approaches, both in terms of mean and median absolute errors in speed estimation. Future work will focus on exploring the use of different window shapes and alternative criterion tailored to various applications.

# Acknowledgements

F. Karkafi gratefully acknowledges the European Commission for its support of the Marie Skłodowska Curie program through the ETN MOIRA project (GA 955681).

# References

[1] R. B. Randall and J. Antoni, "Rolling element bearing diagnostics—a tutorial," *Mechanical Systems and Signal Processing*, vol. 25, no. 2, pp. 485–520, 2011.
[2] S. Braun, *Mechanical signature analysis: theory and applications*, 1986.

[3] D. Remond, J. Antoni, and R. B. Randall, "Editorial for the special issue on instantaneous angular speed (ias) processing and angular applications," *Mechanical Systems and Signal Processing*, vol. 44, no. 1-2, 2014.

[4] C. Peeters et al., "Review and comparison of tacholess instantaneous speed estimation methods on experimental vibration data," *Mechanical Systems and Signal Processing*, vol. 129, pp. 407-436, 2019.

[5] P. Cédric, "Advanced signal processing for the identification and diagnosis of the condition of rotating machinery," Thesis, Université de Lyon ; Vrije universiteit Brussel, Nov. 2019.

[6] K. Grochenig, *Foundations of Time-Frequency Analysis*, ser. Applied and Numerical Harmonic Analysis. Birkhäuser Boston, 2001.

[7] P. Flandrin, *Time-Frequency/Time-Scale Analysis*. Academic Press, 1998, vol. 10.

[8] I. Daubechies, J. Lu, and H.-t. Wu, "Synchrosqueezed wavelet transforms: an empirical mode decomposition-like tool," *Applied and Computational Harmonic Analysis*, vol. 30, no. 2, pp. 243-261, 2011.

[9] G. Thakur and H.-t. Wu, "Synchrosqueezing-based recovery of instantaneous frequency from nonuniform samples," *SIAM Journal on Mathematical Analysis*, vol. 43, no. 5, pp. 2078-2095, 2011.

[10] H.-T. Wu, "Adaptive analysis of complex data sets," Ph.D. dissertation, Princeton University, Princeton, NJ, USA, 2011.

[11] F. Auger, P. Flandrin, Y.-T. Lin, S. McLaughlin, S. Meignen, T. Oberlin, and H.-T. Wu, "Time-frequency reassignment and synchrosqueezing: an overview," *IEEE Signal Processing Magazine*, vol. 30, no. 6, pp. 32-41, 2013.

[12] T. Oberlin, S. Meignen, and V. Perrier, "The fourier-based synchrosqueezing transform," in *2014 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)*. IEEE, 2014, pp. 315-319.

[13] I. Daubechies, Y. G. Wang, and H.-T. Wu, "Concept: concentration of frequency and time via a multitermed synchrosqueezed transform," *Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences*, vol. 374, 2016.

[14] H.-T. Wu, "Instantaneous frequency and wave shape functions (i)," *Applied and Computational Harmonic Analysis*, vol. 35, pp. 181-199, 2013.

[15] R. Behera, S. Meignen, and T. Oberlin, "Theoretical analysis of the second-order synchrosqueezing transform," *Appl. Comput. Harmon. Anal.*, vol. 45, pp. 379-404, 2018.

[16] S. Meignen, D.-H. Pham, and S. McLaughlin, "On demodulation, ridge detection and synchrosqueezing for multicomponent signals," *IEEE Trans. Signal Process.*, vol. 65, no. 8, pp. 2093-2103, 2017.

[17] D.-H. Pham and S. Meignen, "High-order synchrosqueezing transform for multicomponent signals analysis – with an application to gravitational-wave signal," *IEEE Trans. Signal Process.*, vol. 65, pp. 3168-3178, 2017.

[18] C. Peeters, A. Jakobsson, and J. Helsen, "Improved time-frequency representation for non-stationary vibrations of slow rotating machinery," *PHM Society European Conference*, vol. 7, pp. 401-409, 2022.

[19] G. Protopapadakis, C. Peeters, Q. Leclère, and J. Helsen, "Enhancing the performance of the multi-order probabilistic approach in angular speed estimation through adaptive window selection," in *Surveillance, Vibrations, Shock and Noise*. Toulouse, France: Institut Supérieur de l'Aéronautique et de l'Espace [ISAE-SUPAERO], July 2023, https://hal.archives-ouvertes.fr/hal-04165643.
[20] M. Leiber, A. Barrau, Y. Marnissi, and D. Abboud, "A differentiable short-time fourier transform with respect to the window length," in *2022 30th European Signal Processing Conference (EUSIPCO). IEEE, 2022*.

[21] M. Leiber, Y. Marnissi, A. Barrau, and M. El Badaoui, "Differentiable adaptive short-time fourier transform with respect to the window length," in *ICASSP 2023-2023 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP). IEEE, 2023*.

[22] M. Leiber, Y. Marnissi, A. Barrau, and M. El Badaoui, "Differentiable short-time fourier transform with respect to the hop length," in *2023 IEEE Statistical Signal Processing Workshop (SSP). IEEE, 2023*.

[23] A. Rihaczek, "Signal energy distribution in time and frequency," *IEEE Transactions on Information Theory*, vol. 14, no. 3, pp. 369-374, 1968.

[24] D. J. Nelson, "Cross-spectral methods for processing speech," *J. Acoust. Soc. Am.*, vol. 110, pp. 2575-2592, 2001.

[25] R. Behera, S. Meignen, and T. Oberlin, "Theoretical analysis of the second-order synchrosqueezing transform," *Appl. Comput. Harmon. Anal.*, vol. 45, pp. 379-404, 2018.

[26] K. Kodera, C. De Villedary, and R. Gendrin, "A new method for the numerical analysis of non-stationary signals," *Phys. Earth Planet. Interiors*, vol. 12, no. 2-3, pp. 142-150, Aug. 1976.

[27] K. Kodera, R. Gendrin, and C. Villedary, "Analysis of time-varying signals with small bt values," *IEEE Transactions on Acoustics, Speech, and Signal Processing*, vol. 26, no. 1, pp. 64-76, 1978.

[28] F. Auger, E. Chassande-Mottin, and P. Flandrin, "Making reassignment adjustable: The levenberg-marquardt approach," in *2012 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)*, 2012, pp. 3889-3892.

[29] S. Meignen, T. Oberlin, and D.-H. Pham, "Synchrosqueezing transforms: From low- to high-frequency modulations and perspectives," *Comptes Rendus. Physique*, vol. 20, no. 5, pp. 449-460, 2019. [Online]. Available: https://hal.archives-ouvertes.fr/hal-02573644

[30] J. Antoni, J. Griffaton, H. André, L. A. Valencia, F. Bonnardot, O. Cardona-Morales, G. Castellanos-Dominguez, A. Daga, Q. Leclère, and C. Vicuña, "Feedback on the surveillance 8 challenge: vibration-based diagnosis of a safran aircraft engine," *Mech. Syst. Signal Process.*, vol. 97, pp. 112-144, 2017.