Abstract

Real-world physical signals are commonly nonstationary, and their frequency details change with time and do not remain constant. Fourier transform that uses infinite sine/cosine waves as basis functions represents frequency constituents of signals but does not show the variations of the signal frequency contents over time. Multiresolution demonstration of the time-frequency domain may be achieved by the techniques that can support adjustable resolution in time and frequency. Earthquake strong motion signals are nonstationary and indicate time-varying frequency content due to the scattering from the source to the site. In this paper, we applied short-time Fourier transform, S-transform, continuous wavelet transform, fast discrete wavelet transform, synchrosqueezing transform, synchroextracting transform, continuous wavelet synchrosqueezing, filter bank synchrosqueezing, empirical mode decomposition, and Fourier decomposition methods on the near-source strong motion signals from the 7 May 2020 Mosha-Iran earthquake to study and compare the frequency content of this event estimated by these methods. According to the results that are examined by Renyi entropy and relative error, synchroextracting performed better in terms of energy concentration, and the Fourier decomposition method revealed the lowest difference between the original and reconstructed records.

1. Introduction

Time-frequency (TF) analysis (TFA) techniques have been effectively applied to characterize the time-varying properties of nonstationary signals. Fourier transform (FT), short-time Fourier transform (STFT), continuous wavelet transforms (CWT), discrete wavelet transforms (DWT), synchrosqueezing transforms (SST), synchroextracting transforms (SET), CWT synchrosqueezing (SYNCWT), filter bank synchrosqueezing (FBANK), empirical mode decomposition (EMD), and Fourier decomposition method (FDM) are among the well-known employed approaches to this problem. FT basis functions are sine and/or cosine waves of different frequency, amplitude, and phase that decompose the original signal into its harmonic components. The method may not adequately explain spectral features and local transient characteristics of nonstationary signals (e.g., strong motion earthquake signals) caused by averaging over the duration of the signal. FT analyzes the frequency components without time information. STFT can better describe local time-dependent spectral components. This method applies overlapping and sliding narrow windows centered at a definite time. The moving window extracts the local frequency contents at each time. The Fourier transforms of the divided segments are gathered to form the final time-frequency representation. Selection of the window length effects on the created spectrum. Due to the inverse relation between the window length and the related bandwidth, high resolution is not obtainable in the time and frequency domains (TFD) synchronously that limits the effective resolution of time-frequency methods. Wavelet transforms (WT) decompose the signals of interest using finite-length basis functions (mother wavelets) and their shifted, contracted, and dilated samples of the mother wavelets (so daughter wavelets) that inherit different spatial and temporal scales [1]. Many mother wavelets have been introduced for different applications (e.g., for reconstructing a signal from its wavelet transform). The resynthesizing property can be obtained if the total area under the mother wavelet becomes zero (admissibility condition). On the other hand, the finite duration of wavelets prevents propagation of the localized transient features through time by its wavelet transform [2]. Fitting wavelet basis functions to the signals of interest generate a set of wavelet coefficients (filter banks) that are under the influence of the dilated/contracted and shifted mother/daughter wavelets. The dilations/contractions build the concept of scale so that smaller scales are related to larger frequencies and larger scales to lower ones and the shifts are related to time. WT extracts time-varying characteristics of signals, due to the short length of the nonzero part of the wavelets; their coefficients at each instant of time indicate the signal features around the same time. WT applies long-length windows at low frequencies and short-length windows at high frequencies, while STFT employs a fixed window. So WTs are better correlated to time-varying features of signals. Nonstationary signals instantaneous frequency (IF) also gives information about the time-varying spectral variations. SST and SET squeeze the TF coefficients by TFA methods into the instantaneous frequency trajectory in the frequency direction. They are able to reconstruct the signals and provide better resolution [3]. Earthquake signals are nonlinear and nonstationary. It is well known that the Fourier transform is proper for stationary signal processing; the wavelet transform is appropriate for nonlinear signal analysis; and the synchrosqueezing transform better fits the nonstationary and nonlinear signals.

On 7 May 2020, a shallow-depth (11 km) medium-sized Mn5.0 earthquake occurred around the suburbs of Tehran near Damavand city and Mosha village. The event was felt in a great region including neighboring provinces: Tehran, Qom, Alborz, and Semnan. The event left two deaths and some injured people around the Mosha village and made a lot of people stay out of their residences on that night. The epicentral area is across the Central Alborz (Figure 1) that is regionally under the influence of the convergence between the Persian Gulf and Eurasian plate that is reflected in high mountain ranges along the Alborz Mountains [4]. North Tehran reverse fault and Eastern Mosha fault are near the epicenter. The Eastern Mosha fault that is closer to the epicentral area is a reverse strike-slip fault, and according to IRSC, this quake was estimated to be Mw4.9 that showed oblique-slip motion. It was followed by several small-sized aftershocks. The general construction quality in Iran involving this region is considered to be low, and the residences are usually constructed without regarding and/or effective application of seismic design codes. Even medium-sized earthquakes (e.g., 2003 Bam earthquake) leave considerable destructions and causalities in Iran [5]. Studying the earthquakes such as Mosha can help improve not only our knowledge of local seismicity but also seismic design codes and consequently construction qualities.

2. Methods

FT in the time domain provides the frequency-amplitude representation of signals. FT coefficients are calculated as inner products of basis functions (sine/cosine waves) of infinite duration. Any sudden change in the time domain will propagate over the entire frequency domain. The forward and inverse FT is formulated as follows:

The discrete version of FT for a signal of length N and its inverse is suggested as follows:

Although FT presents the signal frequency content, the change of frequencies in time is not provided by this spectrum.

2.1. Short-Time Fourier Transform

The Fourier spectrum of each signal provides a complete frequency-magnitude relation, while one cannot understand the temporal arrangement of the same frequencies. To overcome this problem, we can divide the signal into a series of short-time overlapping segments, apply discrete Fourier transform (DFT) on each segment, move to the next segment, and reapply DFT to the overlapped segment. The segment size should be so small so that the signal frequency content does not change over time for that segment (signal keeps to be stationary over that segment). Let x [n] be a real-time limited signal; then the windowed version of the signal is as follows:and the STFT and inverse STFT (ISTFT) of the signal will bewhere m∈{1, 2, …, M} is the local time index, 1 ∈ {1, 2, …, L − 1} is the segment index, k ∈ {1, 2, …, K} is the frequency bin index, and h is the hop size (the distance that each window moves in the unit of samples). X (k, l) is the localized STFT spectrum of x [n] that indicates the time-frequency representation (TFR) of the signal about the frequency bin k and the time index lh. Finally, to recover the original signal, all the signal segments in the time domain should undergo an overlap-and-add (OLA) process [6]. For STFT and ISTFT estimations, we followed suggestions of [6] and applied the same considerations (Figure 2). The window length applied for STFT is 1,024 samples. One of the vulnerabilities of STFT is its limited frequency-time resolution originated from the predetermined window size.

2.2. S-Transform

The standard S-transform (ST) and its inverse [7] for signal x (t) are defined as follows:for which a wavelet function can be proposed as follows:

Therefore, it may be considered a special case of CWT. Then the transform bypasses the window size problem of STFT while holding some WT properties. The permanent trend of the standard transform wavelet function restricts its applications, so several authors have suggested different forms for the S-transform [8]. On the other hand, S-transform may be viewed as a type of windowed Fourier transform, where the length of the centered window at τ is inversely proportional to f. Discrete ST transforms convert an n point signal into an n2 point one, so they are highly redundant. Discrete orthonormal Stockwell transform (DOST) has been developed to overcome this problem. DOST offers multiple scales, an orthonormal basis, and an absolutely referenced phase. Discrete cosine-based DOST (DCST) provides a real-valued transform [9]. We used the DCST method (Figure 3) as suggested by [10] to transform the records of the Mosha earthquake.

2.3. Continuous Wavelet Transform

Multiresolution analysis decomposes each signal at different frequencies with different resolutions. At low frequencies, the method produces lower time resolution and higher frequency resolution. At high frequencies, it provides higher time resolution and lower frequency resolution. Like STFT and ST, continuous wavelet transform multiplies the signal with a mother wavelet ψ (t) as follows:where s is a scale factor, ∗ shows complex conjugate, τ is a time shift, 1/√s guaranties that energy will be the same for all s, and c is a constant that is related to the mother wavelet. CWT (τ, s)x is a matrix of coefficients including time-scale details of x (t) that is the sum over all time points of the signal multiplied by the shifted and scaled versions of the function ψ [11]. Mother wavelets are usually vibrating signals with zero mean and unit energy. Many different mother wavelets suitable for different applications have been defined in the literature. In this section, we employed Morlet wavelet [12] to apply on the strong motions of this event (Figure 4), the spacing between discrete scales: 0.1 and the mother wavelet parameter: 6. Roughly, one can classify wavelets to real and analytic: real wavelets are usually employed to find sharp signal transitions, while analytic wavelets detach amplitude and phase information [13]. Based on (7), when s increases, mother wavelet contracts and vice versa. So small scale refers to dilated, and large scale refers to contracted versions of signals. The resolution of continuous-time signals does not change with scale change because the change in scale can be reversed. In discrete-time signals, decreasing the scale leads to oversampling that may be undone, so it does not change the resolution. But increasing the scale leads to undersampling and reduces the resolution [14].

2.4. Discrete Wavelet Transform

A sequence x (n) of N points is disintegrable into N + 1 levels (j = −1, 0, 1, …, N − 1), where each level j contains i = 2j evenly distanced (at 2N/j) overlapped wavelets [1]. So for every wavelet i at level j, the discrete wavelet transform (DWT) of x (n) is as follows:

Reference [13] expanded a DWT algorithm that exerts a sequence of high- and low-pass filters to estimate DWT (i, j). The high-passed filter part explains the high details and produces 2N − 1 coefficients. The low-passed part will be forwarded to the subsequent repetition to generate high-passed (to make 2N − 2 coefficients) and low-passed sections. These iterations continue to calculate all the 2N coefficients. The strong motion signals in this study have been analyzed with Daubechies wavelet (Figure 5).

2.5. Synchrosqueezing Transform

Synchrosqueezing transform (SST) is developed by merging TFR techniques and reassignment. It is supposed that SST will ameliorate the energy concentration of TFR by exerting a postprocessing reallocation to the original TFR. SST calculates instantaneous frequency (IF) from the derived time-frequency coefficients and continues with squeezing and reassignment phases. Multiple TFA methods are applicable in the first phase: STFT, S-transform, CWT, and so on [8]. The reassignment method (RS) reappoints the TF spectrogram into the IF path along the 2D time-frequency direction, while SST compresses the TF coefficients into the instantaneous frequency path just in the frequency direction. RS attains high resolution, but the method is not able to do signal reconstruction. Although SST gives lower resolution, the procedure is capable of reconstructing (Figure 6) the analyzed signal [3].

2.6. Synchroextracting Transform

Almost all physical signals involve noise that may propagate into the TFR plane. This unanticipated noise can also spread into the SST results and increase noise powerfulness.

To overcome this feature, the synchroextracting transform (SET) has been suggested. Theoretically, the synchroextracting operator extracts the TF coefficient just in the IF path of interest and rectifying the rest of the TF coefficients. This method employs the time-frequency coefficients containing the maximum value so as to reduce the effect of noises on the time-frequency results (Figure 7). While SST hires squeezing, this method uses extracting (i.e., synchroextracting transform) with better energy concentration properties [3].

2.7. CWT Synchrosqueezing

Synchrosqueezing is able to draw out and determine the time-varying frequency content of the signal components, while the method can reconstruct the same signal after its decomposition. Reference [15] showed that continuous wavelet synchrosqueezing (SYNCWT) offers better precision in time and frequency domains. They presented that synchrosqueezing is invariable to the selection of the basis wavelet. However, the practical differences are caused by the wavelet’s relative concentrations in frequency time. The continuous wavelet synchrosqueezing decomposition and reconstruction of signals of this earthquake has been performed by the Morlet wavelet (Figure 8).

2.8. Filter Bank Synchrosqueezing

Using filter banks is one of the solutions proposed to help solve time and frequency trade-offs. Reference [16] developed SST in order to employ filter banks to provide a sharp time-frequency representation combined with reconstruction property. They expanded reconstructable synchrosqueezing filter banks (FBANK) and applied them by utilizing the nonuniform ERBlet filter banks [17]. In this contribution, we exercised the same method on the Mosha earthquake accelerograms (Figure 9). For a detailed explanation and discussion of the method, the readers would refer to [16] and the references therein.

2.9. Empirical Mode Decomposition

Empirical mode decomposition (EMD) is a method for analyzing nonlinear and nonstationary signals that has been applied in several research areas such as medical, meteorological, geophysical, and so on. Unlike many other decomposition methods that implement predefined fixed basis for their modeling, this technique decomposes the given signal into a group of the finite number of narrow-banded intrinsic mode functions (IMFs) and residual through a sifting procedure. We performed EMD as described in [18, 19] such that the stopping criterion is similar to the one proposed in [20]. For better presenting the results, we have just displayed the first three IMFs in Figure 10.

2.10. Fourier Decomposition Method

Although Fourier-based methods are well known for spectrum analysis, the Fourier decomposition method (FDM) produces a number of band-limited Fourier intrinsic band functions (FIBFs) suitable for analysis of nonlinear and nonstationary signals [21, 22].

FDM is a generalized Fourier expansion with variable amplitudes and frequencies of signals that apply Fourier-based adaptive signal decomposition technique. The FIBFs of FDM build the basis of the algorithm, which is adaptive, local, orthogonal, and complete. The instantaneous frequencies of the FIBFs yield time, frequency, and energy distribution of time series. The FIBFs are zero-mean orthogonal functions that admit analytic representation. The analytic FIBFs are monocomponent signals comprising the sum of zero-mean sinusoidal functions of consecutive frequency bands. Figure 11 shows DCT-based FDM representation of the employed strong motion record produced by decomposing data into 10 equal energy FIBFs: FIBF0 (0–4.4189) Hz, FIBF1 (4.4189–6.3354) Hz, FIBF2 (6.3354–9.2163) Hz, FIBF3 (9.2163–9.5947) Hz, FIBF4 (9.5947–10.4736) Hz, FIBF5 (10.4736–10.9131) Hz, FIBF6 (10.9131–11.6455) Hz, FIBF7 (11.6455–12.2681) Hz, FIBF8 (12.2681–13.2446) Hz, and FIBF9 (13.2446–100) Hz. In Figure 11, just three samples of the calculated FIBFs are shown.

3. Results

Common TFA methods extend any 1D signal into a 2D TF plane in order to look at the signal time-varying attributes. Meanwhile, TFA methods suffer time-frequency uncertainty that leads to blurry images and imprecise TFA presentations. High-resolution methods improve TFA representation of the time-varying specifications. Earthquake time series is nonlinear and nonstationary. We have applied several TFA methods to decompose and recover earthquake accelerograms in order to compare the methods. For comparing the efficiency of the results, we utilized two measures.

The first indicator is Renyi entropy (9) that is an index to assess the energy distribution of the TFRs (Table 1). This indicator definition is as follows:where α is the order of the indicator (default α = 3). Lower values signify that the utilized TFR method produces more concentrated TFR and supply detailed specifications of the analyzed signal features [23].

The second measure is a relative error defined as follows:where norm is the second norm of the signal. The relative error may be considered an illustration of the difference between the original and reconstructed records.

We have applied short-time Fourier transform (Figure 2), discrete cosine-based S-transform (Figure 3), continuous wavelet transform (Figure 4), fast discrete wavelet transform (Figure 5), synchrosqueezing transform (Figure 6), synchroextracting transform (Figure 7), continuous wavelet transform synchrosqueezing (Figure 8), ERBlet filter bank synchrosqueezing (Figure 9), empirical mode decomposition (Figure 10), and Fourier decomposition method (Figure 11) to the vertical, lateral, and transverse components of the accelerograms of the 7 May 2020 Mosha earthquake.

We have selected the first three close accelerographic stations to the mainshock (Mosha3 station 13 km, Rodhen1 station 22 km, and Arjmand1 station 37 km from the epicenter). The sampling frequency is 200 Hz. The results for Mosha3 (the closest station) station are presented in this text, and the results for Arjmand1 and Rodhen1 stations are presented in a separate file attached as supplementary material (Figures s1 to s20 and Tables s1 to s4). Based on BHRC (http://www.bhrc.ac.ir; last accessed 1 August 2020), the peak ground acceleration (PGA) of nearly equal to120 cm/sec2 was recorded on the transverse component of Rodhen1 station. The maximum PGA of nearly equal to 112 cm/sec2 was recorded on the transverse component of Mosha3 station, and PGA of nearly equal to 13 cm/sec2 was recorded on the lateral and transverse components of Arjmand1 station. West-ward stations generally recorded greater accelerations than east-ward stations. The reconstructed records results are compared by employing Renyi entropy (Table 1) and relative error (Table 2). According to the tables, for all 3 station records, synchroextracting transform presents the lowest Renyi index (so a better energy concentration characteristics), and the Fourier decomposition method showed the lowest relative error (so the lowest disagreement between the original and reconstructed records). In terms of the Renyi index, after synchroextracting transform, the short-time Fourier transform, ERBlet filter bank synchrosqueezing, CWT synchrosqueezing, continuous wavelet transform, synchrosqueezing transform, and discrete cosine S-transform showed better performance. The Renyi entropy calculation was not applicable to the fast discrete wavelet transform. Regarding relative errors, the performance of empirical mode decomposition, discrete cosine S-transform, ERBlet filter bank synchrosqueezing, fast discrete wavelet transform, synchrosqueezing transform, CWT synchrosqueezing, short-time Fourier transform, continuous wavelet transform, and synchroextracting transform was better.

Considering that FDM delivers the best results in terms of relative error, we have depicted marginal spectrum and instantaneous energy density as explained in [21]. Figure 12 represents marginal spectrum by FDM for vertical, lateral, and transverse components of the Mosha3 station that is the closest station to the epicenter. The marginal spectrum provides a tool to evaluate the total energy or amplitude contribution from each frequency. Based on Figure 12, almost all of the energy of the vertical record lies within 40 Hz and for the horizontal data is located within 20 Hz.

Figure 13 demonstrates the instantaneous energy density calculated for this quake. All three components show a bimodal character in delivering energy density. The higher density on the vertical component may be taken as an indication of a stronger dip-slip component around this station. Most of the vertical energy density has been delivered within 10 seconds, while for horizontal components, it took nearly 20 seconds.

4. Discussion and Conclusion

Practically, signals may be classified as deterministic and random signals. Deterministic signals are analytically known, and their values are computable with enough certainty. Random signals are divided into stationary and nonstationary signals. Statistical characteristics of stationary signals do not change with time, but statistical features of nonstationary signals vary with time. Time-domain presentation of signals indicates energy and duration of the components, while their frequency representation indicates repetitive characteristics of the components [24]. The traditional tool to study signals in the frequency domain is FT. One of the main drawbacks of FT is losing time information after transforming to frequency domain. TFA tools are powerful methods to transform one-dimensional signals from time domain into two-dimensional time-frequency domain, in order to display how the energy of the signal is distributed over time-frequency domain simultaneously. One can roughly classify TFAs into linear and quadratic transforms. The linear transforms (e.g., STFT and WT) decompose a signal into atomic parts in order to simultaneously and suitably localize each part in time and frequency. One of the vulnerabilities of STFT is limited frequency-time resolution originated from the predetermined window size. Spectrogram (square magnitude of STFT) is able to show the frequency content of the signals, but its estimated frequency content depends on the selected window function. The changing size of the mother/daughter wavelets offered by WT techniques provides more flexibility to surpass the trade-off between frequency and time resolution. But the finite length of the method causes spectral contamination [8]. Quadratic transforms (e.g., Wigner–Ville transform) explain the energy of signal over two descriptive variables: frequency and time. Linear TFAs suffer from the poor resolution that depends on the windows used to analyze the signal and the quadratic transforms suffer interferences for monocomponent and cross-terms in multicomponent signals. Fourier-based transform (i.e., FT and STFT) has fixed frequency-time resolution, and wavelet-based transforms (e.g., CWT and DWT) have a changeable but limited frequency-time resolution. The variable time-frequency resolution of the wavelet transforms qualifies them for extracting details from the broadband and very broadband signals. The SST-based transforms resemble to improve time-frequency resolution by reassigning instantaneous frequencies. Engaging with instantaneous frequencies makes SST-based methods fit better with signals including narrowband contents. Accordingly, it seems that there is no TFA method that could be properly applied to all physical signals. Based on the results that have been derived in this study, synchroextracting transform performed better in terms of energy concentration estimated by the Renyi index. The Fourier decomposition method appears to be better regarding the difference between the original and synthesized records derived by the relative errors. Based on Figures 2–9 (second rows), most of the event’s energy has been released in nearly 10 seconds after the arrival of the primary wave (P-wave) and was delivered relatively in lower frequencies. Short-time Fourier transform spectrograms represent more energy in lower periods relative to the other methods. Marginal spectrum density by FDM shows that most of the energy is located within 40 Hz for the vertical record and within 20 Hz for the horizontal records. On the other hand, vertical record displays stronger energy that is delivered in a shorter time, while horizontal components transmitted smaller energy in a longer time.

Data Availability

The strong motion data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest

The author declares that he has no conflict of interest.

Acknowledgments

The author acknowledges IRSC (irsc.ut.ac.ir) and BHRC (http://www.bhrc.ac.ir) for providing seismicity data and accelerograms. The author has adopted codes and/or scripts from some of the authors referenced within the article. GMT [25] was used for drawing Figure 1.

Supplementary Materials

Supplementary materials include Figures s1 to s20 and Tables s1 to s4 that present the results of the calculations for the accelerograms of the Arjmand1 and Rodhen1 stations. (Supplementary Materials)