**1. Introduction**

High-rate dynamic systems are defined as engineering systems experiencing high-amplitude disruptions (acceleration >100 g*n*) within a very short duration (<100 ms) [1]. Examples of high-rate systems include blast mitigation mechanisms, hypersonic aircraft, and advanced weaponry. Generally, these systems experience rapid changes in their dynamics that can cause malfunctions and sudden failures. The capability to conduct real-time identification of such changes combined with real-time adaptation is critical in ensuring their continuous operation and safety [2]. In terms of system identification, the high-rate problem consists of being able to identify and quantify changes in dynamics over a very short period of time for dynamics containing (1) large uncertainties in the external loads; (2) high levels of non-stationarities and heavy disturbance, and (3) unmodeled dynamics from changes in system configuration.

Work directly addressing the problem of high-rate state estimation, or system identification is limited. In previous work, the authors have proposed an adaptive sequential neural network with a self-adapting input space enabling fast learning of nonstationary signals from high-rate systems [3]. Although the data-based technique showed great promise at high rate state estimation, it did not provide insight into the system's physical characteristics, as it is generally the case with data-based techniques. Physics-driven methods, such as those borrowing on model reference adaptive system (MRAS) theory, showed great promise in handling nonlinearities, uncertainties, and perturbations [4,5]. MRAS was applied to the problem of high-rate state estimation in [6], where the position of a moving cart was accurately identified under 172 ms through a time-based adaptive algorithm used in reaching the reference model with an average computing time of 93 *μs* per step, obtained through numerical simulations conducted in MATLAB. A frequency-based approach was proposed in [7] to identify the position of that same cart. Their algorithm consisted of extracting the dominating frequency through a Fourier transform over a finite set of data and matching that frequency to a set of pre-analyzed finite element models. The authors applied their algorithm experimentally using a field-programmable gate array (FPGA), and were able to accurately identify the position of the cart within 202 ms with a 4.04 ms processing time per step.

A net advantage of frequency-based methods over time-based methods is that they do not typically rely on the tuning of parameters such as adaptive gains. However, they are harder to apply in real time because they are inherently batch processing techniques. It follows that, to enable applications to high-rate systems, one must integrate a temporal approach to the frequency technique in order to extract the required real-time information, a method known as time-frequency representation (TFR). For example, this was done in [7] through the use of a non-overlapping sliding window of 198 ms length. The objective of this paper is to explore the applicability of TFRs to the high-rate problem.

TFRs are widely used for the detection and quantification of faults through vibration-based data [8]. Frequency domain characteristics, such as frequencies, damping ratios, energy in different frequency ranges, and time-frequency domain characteristics, such as time-frequency spread [9], can be used as key features to conduct structural health monitoring [9]. A number of TFRs for instant frequency recognition have been proposed. Popular approaches include linear non-parametric methods, such as short-time Fourier transform (STFT), wavelet transform (WT), and Wigner–Ville distribution (WVD) [9,10]. The application of these methods results in a trade-off between time and frequency resolutions [11]. An adaptive non-parametric method have also been proposed, including the Hilbert–Huang transform (HHT) [12–14], the Cauchy continuous wavelet transform (CCWT) [15], the instantaneous frequencies (IF) re-assignment methods, synchrosqueezing transform (SST) [16], and the multi-SST (MSST) [17].

In this paper, five of these methods are selected, namely the STFT, WT, WVD, SST, and MSST, and their real-time applicability are compared with a specific focus on weakly time-varying systems, here defined as systems with continuously time-varying frequencies, in opposition to abrupt changes (steps, jumps, shifts) in frequencies [18,19]. The objective is to provide a path in designing the next generation of real-time high-rate algorithms. The quantification of performance is conducted on experimental data from the DROPBEAR (Dynamic Reproduction of Projectiles in Ballistic Environments for Advanced Research) testbed [20], which includes the moving cart dynamics used in [6,7]. The remainder of the paper is organized as follows. First, the background on the five TFR methods is presented—after DROPBEAR is introduced and the performance of the methods is analyzed numerically on two different sets of experimental data. Lastly, the performance of each method is compared, and key conclusions are drawn.

#### **2. Time-Frequency Response Methods**

This section gives an introductory background on the TFR methods under comparison. These include the following traditional TFR methods: STFT, WT, WVD, and TF re-assignment methods: SST and MSST.

#### *2.1. Traditional TFR Methods*

#### 2.1.1. Short-Time Fourier Transform

The Fourier transform of a sequence of an *N* discretely sampled signal *x*[*n*] to the discrete frequency domain *X*[*k*] is taken as:

$$X[k] = \sum\_{n=0}^{N-1} x[n]e^{-j2\pi nk/N} \quad n = 0, 1, 2, \ldots, N-1 \tag{1}$$

where *j* is the imaginary unit, and *k* is the corresponding frequency. To integrate a temporal notion, the Fourier transform can be applied over short time segments through a moving window, where the signal can be assumed stationary between two consecutive segments, a method known as STFT. The local Fourier spectrum of each segment can be generated around the position of the window, and the frequency's temporal variation can be observed locally. This is done by modifying Equation (1) as follows [21,22]:

$$\text{STFT}\_{\text{x}}[m,k] = \sum\_{n=0}^{L-1} \text{x}[n]g[n-m]e^{-j2\pi nk/N} \tag{2}$$

where *g*[*n*] is the windowing function of length *L*, and *m* denotes a time shift. With the STFT, a narrower window will improve the time domain resolution but will result in a lower frequency domain resolution. Conversely, a wider window will improve the frequency domain resolution but will result in a lower time domain resolution.

### 2.1.2. Wavelet Transform

The WT method is known for its superior spectral resolution by overcoming the STFT's requirement of predefining a window length. It provides a linear time-frequency representation based on a preselected mother wavelet *ψ* using simultaneous dilation and translation operations. The discretized version of the continuous WT of a signal *x*[*n*] is written [23,24]:

$$\text{WT}\_{\Psi}[m,k] = \frac{1}{\sqrt{c}} \sum\_{n=0}^{L-1} \mathbf{x}[n] \Psi \left[ \left( \frac{n}{c} - m \right) T \right] \tag{3}$$

where *c* is scaling factor. Generally, the WT method has better temporal resolution and lower frequency resolution for higher frequency contents, and better frequency resolution and lower temporal resolution for lower frequency contents. The resulting time-frequency signal transform may be blurred and cannot achieve high resolution simultaneously in time and frequency.

#### 2.1.3. Wigner–Ville Distribution

The WVD method is an approach based on the quadratic energy density obtained through an instantaneous autocorrelation function [25]:

$$\text{WVD}\_{\text{x}}[m\_{\text{}}k] = \sum\_{\text{n}} \text{x} \left[n + m/2\right] \text{x}^{\ast} \left[n - m/2\right] e^{-j2\pi mk/N} \tag{4}$$

where the asterisk denotes the complex conjugate. Challenges in applying the WVD include interferences and negative values [26]. There are times when cross-terms produce oscillatory interference with multiple frequency components, and the magnitude of the interference may range from extremely negative to extremely positive values.

#### *2.2. Time-Frequency Reassignment*

#### 2.2.1. Synchrosqueezing Transform

The SST method is time-frequency reallocation method that yields finer time-frequency representations for a non-stationary multi-component signal. With SST, the energies of the time-frequency coefficients are reassigned to achieve higher energy around the trajectories of IF, resulting in more accurate tracking of these IFs [27]. The SST representation is conducted using:

$$\text{SST}\_{\Phi}[m,k] = \frac{1}{\Delta k} \sum\_{c\_{m}} \text{WT}\_{\Psi}[m,k] \varepsilon^{-3/2} \Delta c\_{m} \tag{5}$$

where *cm* is the discrete scale for which the wavelet decomposition WT*ψ*[*m*, *k*] is computed, and Δ*cm* = *cm*−<sup>1</sup> − *cm* is a scaling step. The SST results in a more concentrated profile and unique IF curves but is inherently more computationally intensive comparing to the WT method. A limitation of the SST method is that it assumes a weakly time-varying system.

#### 2.2.2. Multi-Synchrosqueezing Transform

The MSST method is an improvement of the SST that results in a sharper energy concentration of the TFR [17]. It does not require any other redundant parameter or a priori information to demodulate the frequency-modulated signals, and thus it can be applied beyond weakly time-varying systems. MSST is formulated to post-processes STFT:

$$\text{MSST}^1[m,k] = \text{STFT} - \text{SST}\_x[m,k] = \sum\_{n=0}^{N-1} \text{STFT}\_x[m,k] \delta[k - \hat{k}[m,k]] \tag{6}$$

Multiple iterations of the process can be conducted with

$$\text{MSST}^{N}[m,k] = \sum\_{n=0}^{N-1} \text{MSST}^{[N-1]}[m,k]\delta[k-\hat{k}[m,k]] \tag{7}$$

where MSST*N*[*m*, *k*] is the SST at the *Nth* iteration for *N* > 2. Increasing *N* will yield better IF estimates, but at the cost of higher computational time.

#### **3. Methodology**

This section introduces the research methodology. It includes a description of the experimental setup, data collection process, and TFR performance investigation methodology.

#### *3.1. Experiment Setup*

DROPBEAR is an experimental testbed designed to conduct reproducible high-rate dynamic responses [20]. Shown in Figure 1, the testbed features a cantilever steel beam with a rolling cart moving along the beam. The movable cart functions as a variable pin, used to mimic sudden or gradual changes in stiffness due to a change in boundary conditions. The beam was equipped with one PCB 353B17 accelerometer connected to the beam located 300 mm away from the clamp to measure the beam's response to the moving cart. A modal hammer PCB 086C01 was used to excite the beam at the tip. Figure 1 also shows an electromagnet that can be used to simulate a sudden change in mass through a controlled drop. That feature was not utilized in generating test data used in this paper.

**Figure 1.** DROPBEAR testbed: (**a**) picture and (**b**) schematic of the setup.

Two types of experiments were conducted. First, the beam was excited with the cart fixed at various locations (i.e., "fixed cart configuration"): positions A, B, C, D (respectively 50 mm, 100 mm, 150 mm, and 200 mm away from the clamp). For each location, the cart was maintained in place for 2 s and the beam impacted at 0.5, 2.5, 4.5, and 6.5 s (Figure 2a). Data from the accelerometer were sampled at 1 kHz and used to compute frequency response functions (FRFs) using the *H*<sup>1</sup> estimation method [28] plotted in Figure 2c and extract the fundamental natural frequencies at 26.5, 31, 38, and 47.5 Hz (Positions A, B, C, and D, respectively). The fixed cart configuration was used to examine the performance of each TFR method over the entire dataset, without the use of sliding windows (except for the STFT). Second, the cart was moved (i.e., "moving cart configuration") starting at 50 mm from the clamp at 0.5 s to 200 mm from the clamp over 1 s, stayed for 1.39 s, and then returned to the initial position at 4.26 s. The acceleration was sampled at 25 kHz. The time series response of the beam is plotted in Figure 2b. The moving cart configuration was used to examine the real-time applicability of methods through the use of sliding or non-overlapping sliding windows. The type and size of windows, along with the type of wavelets, were selected heuristically to obtain the best overall performance.

**Figure 2.** (**a**) signal from the fixed cart configuration at 50, 100, 150, and 200 mm; (**b**) signal from moving cart configuration between 50 mm to 200 mm; and (**c**) frequency response functions (FRFs) for DROPBEAR with fixed cart positions at 50, 100, 150, and 200 mm.

## *3.2. Performance Analysis*

Analyses were performed in MATLAB 2019b, with an Intel(R) Core(TM) i7-7700 CPU 3.6 GHz. The performance of the TFR methods in the fixed cart configuration is assessed as a function of four performance metrics (*J*<sup>1</sup> to *J*4). Metric *J*<sup>1</sup> is the mean absolute error between the estimated *ω*ˆ and real *ω*

frequency over the length *n* of the signal. Metric *J*<sup>2</sup> is the energy concentration of the TFR using Renyi entropy [29]. Metric *J*<sup>3</sup> is the computation time per window. Metric *J*<sup>4</sup> is the mean convergence time when the estimation error falls and remains within an error threshold, here taken as 5%. Metrics *J*<sup>1</sup> and *J*<sup>2</sup> can be expressed mathematically as follows:

$$\begin{aligned} J\_1 &= \sum\_{i=1}^n \frac{|\hat{\omega}\_i - \omega\_i|}{n} \\ J\_2 &= \iint \log|\text{TFR}(t, \omega)^3| dt d\omega \end{aligned} \tag{8}$$

In terms of performance assessment, small values for *J*1, *J*3, and *J*<sup>4</sup> are desired, while a high value for *J*<sup>2</sup> is desired. The performance of the TFR methods in the moving cart configuration is only assessed using performance metrics *J*<sup>1</sup> and *J*3, given that the problem of interest is frequency tracking.

#### **4. Results and Discussion**

This section presents results from the fixed and moving cart configurations and discusses the performance of the different TFR methods for applicability to high-rate system identification.

#### *4.1. Fixed Cart Configuration*

A parametric study was first conducted to study the influence of TFR parameters and select the optimal parameters in comparing performance across TFRs. The study starts with the STFT, where the window length and type are investigated. In the investigation, the window size ranged from 128 to 768 at an interval size of 64, window length overlaps were a half and a quarter of the window length, and windowing functions were Hanning, Gaussian, and Blackman, as shown in Figure 3. Figure 4a plots the results. To facilitate the comparison, the four metrics were normalized to the highest value at 1. It can be observed that both *J*<sup>1</sup> and *J*<sup>2</sup> converge after a window length of 512. In addition, under the *J*<sup>2</sup> and *J*<sup>3</sup> metrics, the half overlapped window performs better than the quarter overlapped windows. Under *J*4, the half overlapped Gaussian window performs similarly to the quarter overlapped windows. However, the other half overlapped window functions generally perform worse, and the window length of 512 appears to yield optimal results. From these results, a window length of 512 with a half-overlapping Gaussian windowing function was selected as the optimal set of parameters.

**Figure 3.** Window functions used in STFT.

The parametric study for the WT consisted of evaluating the performance under different wavelets, including the Morse, Morlet, and Bump wavelets [30,31]. Figure 4a plots the results. Results show a similar performance across all wavelet types, with slightly better performance for the Morlet wavelet observable under *J*3. The Morlet wavelet was selected as the optimal parameter For the MSST, the number of 2 to 5 iterations were studied. Results were also similar across all iteration numbers, and thus a total of two iterations were selected as the optimal parameter.

**Figure 4.** Parametric investigation for (**a**) STFT; and (**b**) WT.

Results from the fixed cart configuration experiment using the selected optimal parameters are shown in Figure 5. Figure 5a is the time series response over 2 s under each location. Figure 5b–f show the time-frequency content extracted by the STFT, WT, WVD, SST, and MSST methods, respectively. Ridge detection is used to identify the first natural frequency (shown as a solid red line) by extracting the maximum-energy time-frequency ridge of the spectrograms [32]. The measured frequencies from the FRF (Figure 2c) are also shown in blue dashed lines. Table 1 lists performance under metrics *J*1–*J*4.

Visual observation of the time-frequency plane (Figure 5b–f) shows that the WVD provides the best frequency identification at position A, followed by the SST, while the WT yielded high variance in results. This may be attributed to the weaker acceleration response from the beam compared to other positions (Figure 5a). Over other positions, all methods appear to have identified a stable frequency, with the SST and MSST converging faster than the other methods, followed by the STFT. An examination of the performance metrics listed in Table 1 confirms these observations. The SST provided the most precise estimate of frequency, while the WT was the least precise (*J*1). The SST and MSST methods showed the highest energy concentrations (*J*2), as expected from time-frequency reassignment methods. However, their computation time (*J*3) was significantly higher compared to the other methods, with the STFT and WT being the fastest. In terms of convergence (*J*4), the SST and MSST were the fastest, while the STFT was the slowest. From the static cart experiment results, it appears that the WT is one of the most applicable methods to the high-rate problem given its fast convergence speed and lower computation time, although at the cost of precision on the estimation of frequency over weak signals.

**Figure 5.** Time-frequency planes from the fixed cart configuration obtained using the (**a**) signal from the fixed cart configuration at 50, 100, 150, and 200 mm; (**b**) STFT; (**c**) WT; (**d**) WVD; (**e**) SST; and (**f**) MSST method with extracted frequencies (red solid lines) and estimated true frequencies (dashed blue lines).


**Table 1.** Fixed cart configuration time-frequency analysis comparison.
