Optimization of multicomponent signals entered to the system using estimation of instantaneous frequency

Abstract

The frequency signal displays are not efficient for analyzing nonstationary signals because of their inability to represent frequency changes over time. In fact, because most of the signals are real, nonstationary, and time varying, analyzing the signals in the time–frequency domain to estimate the instantaneous frequency of a signal is inevitable. The methods of estimating the instantaneous frequency of the multicomponent signals are divided into three groups, which include the methods using signal phase derivatives that are sensitive to noise, methods that calculate the number of zero points of the signal and consider the signal frequency equal to half the frequency of the zero points and are suitable for signals that can be imagined as stationary, and methods based on time–frequency distributions and distributions such as Wigner for instantaneous frequency calculations and more for instantaneous frequency calculations on nonstationary noise signals that exhibit varied time–frequency distributions. In this article, a new hybrid algorithm is used to evaluate different distribution criteria and comparing their performance in investigating one or more features of time–frequency distributions, such as resolution and energy concentration.

Keywords

Multicomponent signal empirical mode decomposition instantaneous frequency estimation generalized likelihood ratio test Wigner–Ville distribution

1. Introduction

Frequency signal displays are not efficient because of their inability to represent frequency changes over time; they are not efficient to analyze nonstationary signals. In the time and frequency domain, there are two classical methods for signal representation that are considered separately in both methods of time and frequency variables. In the time domain, the median is calculated from the representation values in the ones. For this reason, signal frequency changes over time are not evident in the frequency domain. In time–frequency methods, it is possible to represent and analyze the signal analytically simultaneously (Biagetti et al., 2015). Transient signals (i.e., signals that start and end at specific times) can also be represented in the frequency domain using the Fourier transform (FT). The FT representation of a transient signal, s(t), is given by

\begin{array}{l} X (ω) = \int_{- \infty}^{+ \infty} x (t) e^{- j ω t} d t \\ S (t) = \frac{1}{2 π} \int_{- \infty}^{+ \infty} X (ω) e^{j ω t} d ω \end{array}

(1)

Most signals are real, non-static, and time variant. It is not possible to study such signals in the time domain or in the frequency domain separately. A time–frequency function allows the optimal analysis of signals by expressing the frequency content according to time. In the following section, we describe the time–frequency distributions (TFDs).

1.1. TFD

Time–frequency methods for analyzing non-static signals are introduced. Existing methods for analyzing non-static signals can be divided into three general categories which include instantaneous frequency (IF), linear distributions, and quadratic distributions. In the following, we will describe these distributions.

1.1.1. Instantaneous frequency

If the signal is mono-component (having only one frequency component) and the signal energy is greater than the noise contained therein, we will face

IF (t) = \frac{1}{2 π} \frac{d}{d t} ϕ (t)

(2)

where ϕ(t) is the signal phase in time domain. Instantaneous frequency signal analysis shows signal behavior in many cases but does not display all the signal details. Linear and quadratic TFDs are ways to present more detail (Biagetti et al., 2015).

1.1.2. Linear distributions

Linear TFDs map the signal from the time domain to the two-dimensional time–frequency plane. Linearity is a desirable feature in systems that work with multicomponent signals. The TFD of TF_s(t, ω) is linear if and only if the TFD of the signal S(t) = S₁(t) + S₂(t) is equal to the set of TFDs of its constituents (Franco et al., 2012)

T F_{s} (t, ω) = T F_{s_{1}} (t, ω) + T F_{s_{2}} (t, ω)

(3)

The two most commonly used linear TFDs are short-time Fourier transform (STFT) and Gabor transform. Consider the signal S(τ) with the real window and the pair ω(τ). S_ω(t, τ) is the local spectrum of the signal S(τ) at t = τ and is calculated as S_ω(t, τ) = S(τ)ω(τ − t). The STFT of the signal is the FT of the local spectrum of the signal S(τ) as $F_{s}^{ω} (t, ω) = F_{τ \to ω} {S (τ) ω (τ - t)}$ . The Gabor transform of a signal is the STFT with a Gaussian window

G T_{s} (t, ω) = \int h (τ - t) S (τ) e^{- j ω τ} d τ

(4)

where h(τ) is Gaussian window and S(t) is signal Gabor transform. If the signal S(t) is composed of the N component of S_i, then the Gabor transform of the total signals is equal to the sum of the Gabor transform of the signal components (Daubechies et al., 2011).

G T_{s} (t, ω) = \sum_{i = 1}^{N} G T_{s_{i}} (t, ω)

(5)

Time and frequency decomposition in the STFT and Gabor transform are inversely related. In other words, if the frequency resolution of the distribution increases, its time resolution decreases, and vice versa. Based on the uncertainty principle, it is impossible to achieve full time and frequency decomposition in the STFT. In general, linear TFDs cannot have high frequency and time resolution at the same time.

1.1.3. Quadratic distributions

Distributions that map a signal from the time domain to the energy density plane. Although these methods are not linear, they have high frequency and time decomposition. Because these distributions produce the signal energy density plate, there are three primary conditions:

The sum of the full-time frequency signal plane is equal to the signal energy.

The integral along the frequency of the time–frequency signal plane is equal to the signal power.

The integral part of the time–frequency signal plane is equal to the signal spectrum energy density.

It should be noted that above features (2) and (3) do not contain information about two-point (ω, t) signal energy (Oberlin et al., 2014). In the Wigner–Ville disorders (WVD), TFD forms the basis of most new distributions. In the WVD, the S(t) signal is calculated as follows

W D_{s} (t, ω) = F_{τ \to ω} {s (t + \frac{τ}{2}) s^{*} (t - \frac{τ}{2})}

(6)

where

s (t + τ / 2) s^{*} (t - τ / 2)

is a function of the signal’s instantaneous autocorrelation because of its symmetry in the Hermitian space. In other words, the WVD is a real signal. Although this distribution has high time–frequency separability, it needs to be improved because of cross-terms production that causes distortion in the time–frequency plane. Assume the signal S(t) is composed of N components, we have

S (t) = \sum_{i = 1}^{N} S_{i} (t)

(7)

If S_i,p and S_i,n are the positive and negative frequencies of the S_i component, then the WVD of the S(t) signal will be as follows

\begin{array}{l} W D_{s} (t, ω) = \sum_{i = 1}^{N} \sum_{α \in {p, n}} W D_{S_{i \cdot α}} (t, ω) \\ + \sum_{i = 1}^{N} \sum_{j = 1}^{N} \sum_{d \in {p, n}} \sum_{d \in {p, n}} W D_{S_{i \cdot α} S_{j \cdot β}}^{I} (t, ω) \end{array}

(8)

where

W D_{S_{i \cdot α}} (t, ω)

is the time–frequency transform of S_i·d (the optimal part of the TFD), and

W D_{S_{i \cdot α} S_{j \cdot β}}^{I} (t, ω)

is the external cross-term between S_i·α S_j·β equation (9) and is defined as follows

\begin{array}{l} W D_{S_{i \cdot α} S_{j \cdot β}}^{I} (t, ω) ≜ WV D_{S_{i \cdot α} S_{j \cdot β}} (t, ω) + \\ WV D_{S_{j \cdot β} S_{i \cdot α}} (t, ω) \end{array}

(9)

Adjacent and distant components can also create double cross-terms that will have undesirable effects on signal analysis. In the WVD, there is a cross-term between positive and negative frequencies. Thus, if, instead of computing the original signal distribution, it calculates the signal distribution with indefinite frequencies corresponding to the original signals, a portion of the external cross-term is omitted. Such signals are called analytical signals. The Wagner–Ville distribution is the same as the WVD, but instead of computing the original signal distribution, the corresponding analytical signal distribution is calculated (corresponding nonfunctional frequencies). The S(t) is real if and only if S(ω) = S*(−ω) or S(−ω) = S*(ω), where S(ω) is the FT of the signal. Because the real signal’s negative frequencies can be calculated from its positive frequencies, the negative frequencies can be eliminated from the actual signal repetition. In real low-pass signals, the elimination of negative frequencies can have two advantages: first, half the total bandwidth and can be sampled without any problems at half the sampling rate of the first state. In other words, the Nyquist rate is halved. Second, prevent partial cross-term (caused by the interaction of positive and negative signal frequencies in the time–frequency quadratic distribution). A sampling of the real signal S with Nyquist frequency, incidence of the reflection phenomenon is half the rate of the Nyquist. Also, a sampling of the analytical signal S with Nyquist frequency: no frequency reflection phenomenon; one way to construct an analytical signal from a real signal is to use the Hilbert transform. The signal S(t) Hilbert transform is defined as

H {s (t)} = (- j \times sgn (ω)) S (ω)

(10)

where S(ω) is FT s(t). The analytical signal with the real DC component corresponding to the signal S(t) is computed as (Jamal, et al., 2020)

Z (t) = s (t) + H {s (t)}

(11)

The WVD of the signal S(t) with the analytic signal z(t) is calculated from the relation as $W D_{s} (t, ω) = F_{τ \to ω} {z (t + τ / 2) z^{*} (t - τ / 2)}$ . The advantage of this distribution is to prevent the occurrence of frequency reflection phenomena in the sampling of real signals with Nyquist frequency. The WVD is incapable of eliminating external cross-terms between the positive signal frequency components and internal cross-terms. Kernel-based TFDs are used to reduce the cross-terms along with preserving frequency and time separability. In fact, we integrated quadratic features (high time and frequency separability) and linear (nonexistent cross-terms) distributions (Tary et al., 2014).

1.2. Kernel-based TFDs

Internal terms can be reduced by the convolution of WVD with a smoothing core. This type of distribution is considered to be a quadratic distribution class, Li and Liang (2012), and is defined as follows

ρ_{s} (t, ω) = γ (t, ω) \times WV D_{s} (t, ω)

(12)

where ρ_s is the core-based TFD of the S(t) signal with a two-dimensional kernel γ(t, ω). The convolution action in the above relation reduces the frequency separability. In other words, this distribution reduces the cross-terms to reduce the time–frequency separability. Generally, the design of the TFD cores is conducted in four areas, including time–frequency, time lag, Doppler frequency, and Doppler lag. The time-lag kernel is the inverse FT of the kernel and time–frequency domain. As seen in equation (12), a Doppler-independent (DI) kernel is indeed independent of Doppler in all four domains. In the Doppler-lag domain, it is a function of lag alone. Equation (12) shows that a DI kernel causes smoothing of the WVD in the frequency direction only. The equation (12) can be defined as

G (t, τ) = F_{τ \to ω}^{- 1} {Y (t, ω)}

, which shows that a quadratic TFD with a DI kernel is a windowed WVD; the windowing is applied in the lag direction before Fourier transformation from lag to frequency. The kernel of time–frequency domain is defined as

ρ_{s} (t, ω) = \int_{R} G (t, τ) \times k_{z} (t, τ) e^{- j ω τ} d τ

(13)

where k_z(t, τ) is instantaneous autocorrelation of analytical signal z(t) as

k_{z} (t, τ) = z (t + τ / 2) z^{*} (t - τ / 2)

, and similarly

k_{x} (t, τ) = x (t + τ / 2) x^{*} (t - τ / 2)

. Doppler frequency domain to simplify kernel is calculated as g(v, ω) = F_t→v{γ(t, ω)}. Therefore, equation (12) is rewritten as follows

ρ_{s} (t, ω) = \int_{R} g (v, ω) \times k_{z} (v, ω) e^{j ω t} d t

(14)

where k_z(v, ω) is spectral autocorrelation function from analytical signal z(t) of signal s(t) and is defined as

K_{z} (v, ω) = x (ω + v / 2) x^{*} (ω - v / 2)

. Because the core in the Doppler-delay domain can be calculated from one of the two relations, g(v, τ) = F_t→v{G(t, τ)} and

g (v, τ) = F_{t \to w}^{- 1} {G (v, w)}

, the equation (12) can be written as follows

ρ_{s} (t, ω) = \int_{R} \int_{R} g (v, w) \times A_{z} (v, τ) e^{- j ω τ + 2 π j t v} d τ d v

(15)

where A_z(v, τ) is function ambiguity of analytical signal z(t).

1.3. The evaluation criteria

In this section, some criteria for evaluating signal TFDs are given. The need for these criteria in time–frequency methods is essential for two reasons, comparison of different distributions and automatic and adaptive selection of the parameters of a distribution in a particular application. These criteria examine one or more features of the TFD, such as separability and energy concentration. Because the time–frequency plane of the quadratic distributions are all signal energy planes, the energy concentration can be defined as the number of nonzero points of distribution. The smaller the number of points, the greater the concentration of signal energy (Coifman and Wickerhauser, 1992).

1.3.1. Shannon entropy

This criterion measures the energy concentration and the separability of the TFD. To represent the time–frequency TF(n, ω), the following relation is calculated

E_{Shannon}^{TF (n, ω)} = \sum_{n} \sum_{ω} TF (n, ω) lo g_{2} (TF (n, ω))

(16)

where the TF{·} can be selected by any TFDs. Equation (16) is not suitable for the representation of time–frequency with a negative value such as WVD and uses the following equation

E_{Shannon}^{TF (n, ω)} = \sum_{n} \sum_{ω} | TF (n, ω) | lo g_{2} (| TF (n, ω) |)

(17)

where the higher the Shannon entropy, the lower the distribution energy’s separability and concentration, and vice versa.

1.3.2. Renyi entropy

If a signal is made of many basic components, its time–frequency plane will be highly complex (containing a lot of information). This criterion determines the complexity of the TFD and is derived from the following equation as

Reny i_{α}^{TF (n, ω)} = \frac{lo g_{2} (\sum_{n} \sum_{ω} T F^{α} (n, ω))}{1 - α}

(18)

where α is entropy order, and usually greater than one; energy-unfocused entropy is generally calculated as

Reny i_{α}^{TF (n, ω)} = \frac{1}{1 - α} lo g_{2} (\frac{\sum_{n} \sum_{ω} T F^{α} (n, ω)}{\sum_{n} \sum_{ω} TF (n, ω)})

(19)

where the lower the amount of Renyi entropy, the higher the quality and the lower the complexity of the TFD, and vice versa. The less complexity of a TFD can be seen as equivalent to high separability and greater energy concentration. Given that the signal’s IF at any time calculates the largest frequency component of the signal in terms of amplitude, it cannot show the details of the components of the signals. Hence, the need for more representative analyzes. Linear TFDs display more details of the signal. Although this method has a good frequency and time separability, it does not have good frequency and time ones. Hence, time–frequency signal distributions are considered as a significant concern. Quadratic time–frequency methods are for high frequency and time separability but also with cross-terms. Therefore, kernel-based methods are designed and used to reduce cross-terms by preserving time–frequency separability. The simulation of the evaluation criteria for comparing different distributions and automatic and adaptive selection of distribution parameters is shown in Figure 1.

Figure 1.

Evaluation of Renyi entropy criteria (solid line) and Shannon entropy criteria (dashed line) for orientation: comparison of different distributions, and automatic and adaptive selection of parameters of distribution in the specific application.

2. Signal analysis using FT

Suppose that the number of sources P (P is known) is dotted in space and that the RF waves are propagating flat, we want to estimate the signals to enter and extract the signals entering these known locations. For this, we assume that the signal is received at the M location of the system. Signals are collected after propagation with an additive white Gaussian noise (AWGN) and are received by each of the locations in the system, Pham and Meignen (2017). The signal received by the mth receiver of the array at time t will be as follows

x_{m} (t) = \sum_{i = 1}^{p} s_{i} ((t - τ_{m} (θ_{i}) + e_{m} (t)

(20)

for m = {1, 2, …, M} and i = {1, 2, …, p}, and s_i(t) is signal emitted by ith source, τ_m(θ_i) is delay of receiving the ith source signal through the point of mth in the system, n_m(t) is noise received by mth point in the system, and θ_i is the angle defining the ith source location. If we sampled the signals received in the system, the equation (20) is expressed in the discrete-time domain by converting t to n, and converting τ_m(θ_i) to

d_{m_{i}}

x_{m} [n] = \sum_{i = 1}^{p} s_{i} (n - d_{m i}) + n_{m} [n]

(21)

where

d_{m_{i}}

is ith signal delay on the mth location in the system. Also,

d_{m_{i}}

is a function of the known location of the received signal and the location of the signal sources, as

d_{m_{i}} = f (θ_{i}, ϕ_{m})

, where θ_i is ith source location, and ϕ_m is mth point location in the system. Purpose of estimating vectors

Θ = {[Θ_{1}, Θ_{2}, \dots, Θ_{p}]}^{T}

, and

{[s_{1} (n), s_{2} (n), \dots, s_{P} (n)]}^{T}

by receiving the system received signals

{[x_{1} (n), x_{2} (n), \dots, x_{m} (n)]}^{T}

. For this purpose, it is assumed as follows:

The noise at each of the receiving points in the white and Gaussian system is independent of each other at different points and signals.

The number of sources of P is known.

The signals have a wide bandwidth and are centered around a known frequency ω₀.

There is enough distance between the two signal sources so that they can be separated.

The zero-frequency component of the signals is not discussed.

In hypotheses, signals are broad bandwidth and can be random or deterministic processes. Here, we study this topic based on STFT. We know it is necessary for the analysis based on the FT. If the signal x(n) is present over time, then the FT is noncausal. In addition, if the signal statistics change over time, meaning nonstationary sources, the Fourier analysis method is not feasible, but with the help of STFT, both of the above axes will be eliminated (Oberlin et al., 2015). With this method, we can analyze the signals whose spectrum changes over time in the frequency domain and estimate the parameters of the signals measured for a limited time. In STFT method, $X (n, ω) = \sum_{r = - \infty}^{\infty} x (r) w (n - r) e^{- j ω r}$ and ω(n) is window. The discrete form of X(n, ω) obtained by sampling the above equation at N is called discrete STFT (DSTFT) and is as follows

x (n, k) = \sum_{r = - \infty}^{\infty} x (r) w (n - r) e^{- j 2 π k r / N}

(22)

By holding the FT at a finite time, the signal Y_n(m), that means X(n, ω), can be written as $x (r) w (n - r) = (1 / 2 π) \int_{- π}^{π} x (n, ω) e^{j ω r} d w$ . For n = r, we have $x (n) = (1 / 2 π ω_{0}) \int_{- π}^{π} x (n, ω) e^{j ω r} d w$ . Unlike x(n, ω), which is reversible, X(n, k) or DSTFT is only recursive under certain conditions as

x_{m} (n, k) = \sum_{i = 1}^{p} s_{i} (n - d_{m_{i}}, k) e^{- j 2 π d_{m_{i}} k / N} + E_{m} (n, k)

(23)

Assume ω = e^−j2πk/N, and $[d_{m_{1}}, d_{m_{2}}, \dots d_{m_{P}}]$ is signal source delay vector at mth point of the system. Because x_m(n, k) = a_m(k).S_m(n,k) + E_m(n,k), where $a_{m} (k) = {[W^{d_{m_{1}}} W^{d_{m_{2}}} \dots W^{d_{m_{P}}}]}^{T}$ and $X (n, k) = {[X_{1} (n, k), X_{2} (n, k), \dots, X_{m} (n, k)]}^{T}$ , it can be written as

\begin{array}{l} x (n, k) = diag (a_{1}^{T} (k), a_{2}^{T} (k), \dots, a_{M}^{T} (k)) \\ diag (S_{1} (n, k), S_{2} (n, k), \dots, S_{M} (n, k)) + E (n, k) \end{array}

(24)

where diag is the diagonal matrix whose elements are diagonal vectors in parentheses with length P, and E(n, k) is equal to

E (n, k) = {[E_{1} (n, k), E_{2} (n, k), \dots, E_{m} (n, k)]}^{T}

(25)

Considering the large dimension of equation (24), we consider the conditions that S_m(n, k) of m from one to M equals a good approximation so that we can write to the signal source ith as

S_{i} (n, k) \approx S_{i} (n - d_{1, i}, k) \approx \dots \approx S_{i} (n - d_{M, i}, k)

(26)

The estimation of the signal s(n, k) is obtained by assuming the availability of resources (with matrix D) based on maximum likelihood estimation (MLE) as follows

\hat{S} (n, k) = {[A {(D, K)}^{H} A (D, K)]}^{- 1} A {(D, K)}^{H} X (n, k)

(27)

To determine the estimates that are a function of vectors θ and ϕ, that is, the locations of the sources and the receiver, we use the statistical gradient method based on the logarithm of the conditional probability density function. We define the θ_is such that the function L defined in $L (n, θ) = - \sum_{k = 1}^{N - 1} 1 / σ N^{2} (k) | X (n, k) - A (D (θ, k) S {(n, k)}^{2} |$ , which is related to the logarithm of the probability density function, is maximized. We need to find θ_i signals that equal $\partial L (n) / \partial θ_{i} = 0$ . Then, we can write with the static gradient method and given that D is a function of the (ϕ, θ) vectors

\begin{matrix} {\hat{θ}}_{i} (n + 1) = {\hat{θ}}_{i} (n) - μ \nabla θ_{i} L (n) = {\hat{θ}}_{i} (n) \\ + μ \sum_{k = 1}^{N - 1} \frac{S_{i} {(n, k)}^{H}}{σ N^{2} (k)} \nabla_{θ_{i}} \cdot A {(D (θ, k), k)}^{H} \\ \cdot [X (n, k) - A (D (θ, k) \cdot k) \cdot {\hat{S}}_{i} (n, k)] \end{matrix}

(28)

Therefore, the adaptive algorithm, after the initial selection of the parameter μ and determination of the initial values of the input angle θ_i, finalizes the input angles of the signals using the following methods: (1) Using equation (27), s(n, k) and given the presence of X(n, k) is the STFT of the data and the convergence of A(D, k) is obtained by estimating $\hat{S} (n, k)$ . (2) Find out the initial estimate of the angle of entry as well as calculate $\hat{S} (n, k)$ from the previous step, using equation (28) to obtain the new value of $\hat{θ}$ , and then repeat one step until the algorithm converges.

3. Signal processing theories and methods

3.1. Empirical mode decomposition

The empirical mode decomposition (EMD) is an adaptive tool for analyzing nonlinear and nonstationary signals that separate the component of the signal based on the local behavior of the signal. The signal will be decomposed into a set of mono-component functions called empirical state functions using the EMD method. An intrinsic state function is similar to a harmonic function. Like a harmonic function, it does not have a fixed amplitude and frequency and has different frequencies with different amplitudes. At each stage of the signal decomposition into its frequency components, the high-frequency components are separated first, and this process continues until the components with the lowest frequency remain. In other to consider a waveform as an empirical mode function, the following two conditions must be met simultaneously: the number of extremes and the number of zero-cross is equal to or equal to at most one number. The mean value of the upper and lower curves of the curve at each point is equal to zero (Wang et al., 2014). The empirical mode functions will be obtained in the following steps from a time series:

Determine the local maximum and minimum point of the input signal;

Obtain the top curve envelope by connecting the local maximum points of the time series and repeat this operation with the local minimum points to create the bottom curve envelope by cubic spline.

Calculate the mean of the top and bottom envelope of the curve (m₁(t)).

Subtract the upper and lower bandwidth of the input signal and generate the first signal component as

h_{1} (t) = x (t) - m_{1} (t)

(29)

where m(t) is the mean of the top and bottom envelope of the curve, h(t) is the first signal component, and x(t) is the input signal. If h₁(t) satisfies both conditions of the empirical case function simultaneously, it is known as the first mode function (mif₁(t)) of the signal. Otherwise, h₁(t) is assumed to be the main function and steps (1)–(4) repeated until h₁(t) can be defined as h₁(t) = X(t) − m₁(t).

5. In the fifth step, k orders are repeated until h_1k(t) has both conditions related to the empirical mode function. In this case, c₁ = h_(1,k) is taken as the first empirical mode function.

6. Separate the first large mode function from the signal x(t) and obtain the first remain u₁ as u₁(t) = x(t) − c₁(t).

7. Consider the signal u₁(t) as the main signal and repeat the first to seventh steps to obtain the second empirical mode function. The steps above are repeated n times to obtain the n empirical mode function. This algorithm stops when u_n(t) of each signal is single frequency and has no other extractable frequency component.

3.2. Cluster empirical mode analysis

Mixing modes is perhaps the most common problem we encounter when working with the algorithm $x (t) = \sum_{i = 1}^{n} c_{i} + u_{n}$ . In decomposition of empirical mode, a particular signal may not be separated into identical inherent mode functions. The EMD method is capable of separating modes with a specific frequency interval and amplitude. The cluster EMD method, introduced with the idea of adding AWGN at all signal degradation steps to solve this problem, will effectively eliminate state mixing. All data are almost noisy. To generalize this idea, noise is added to the input signal repeatedly. Although adding noise may result in a lower signal-to-noise ratio (SNR), the added AWGN can effectively facilitate the process of EMD decomposition by the uniform distribution. Low SNR values do not affect the decomposition method (Wu, 2011). The EMD method can be described in the following steps:

Add AWGN to the data as x_j(t) = x(t) + D(t), where j = 1, 2, …, M, M is number of preset attempts, and D is the amplitude of AWGN.

Data analysis with noise added by the EMD method as

x_{j} (t) = \sum_{i = 1}^{N_{j}} c_{i j} + u_{N j}

(30)

where c_ij represents the ith empirical mode function from ith attempt, u_Nj represents remaining from j − th attempt, and NJ is the number of empirical mode function for j − th attempt.

3. The first and second steps will be repeated several times with different white noise.

4. The average empirical mode functions are calculated from the third step and are obtained as final empirical mode functions as

c_{i} (t) = \frac{1}{M} (\sum_{i = 1}^{M}), i = 1,2, \dots, I

(31)

where I is the minimum number of empirical mode functions among all attempts. The number of attempts and the amplitude of the added noise is two determinants of the cluster EMD method’s performance. For this to work well, the added noise amplitude must not be too small, and empirically, it is usually considered to be a standard deviation of 0.2, and the number of attempts at the optimal value is usually 100.

3.3. Modified EMD method for noisy signal analysis in time–frequency domain

The Hilbert–Huang transform (HHT) method has been used in recent years to process received signals, which results in noise-free signals (Coifman and Wickerhauser, 1992). Therefore, it can be applied to the noisy signals by modifying this method and consequently modifying the EMD method to decompose the signals into single components. The HHT method is based on nonlinear and nonstationary signal analysis based on EMD and Hilbert transform. Standard deviation in the EMD method using the HHT method is as follows

{SD}_{k} = \sum_{t = 0}^{T} (\frac{{| h_{1 (k - 1)} (t) - h_{1 k} (t) |}^{2}}{h_{1 (k - 1)}^{2} (t)})

(32)

such that the stop condition is 0.2–0.3. By modifying the EMD method, instead of the stop criterion proposed by Huang, we use the following proposed stop criterion

{ST}_{k} = \frac{\sum_{t = 0}^{T} {| h_{1 (k - 1)} (t) - h_{1 k} (t) |}^{2}}{\sum_{t = 0}^{T} h_{1 (k - 1)}^{2} (t)}

(33)

In the stop criterion used by Huang, the expression inside the bracket is first calculated, and if the face and denominator of the fraction in a time step are very small, the error in the program is $0 / 0$ and is applied ∑ to the desired interval (t = 0: T), the correct result is not achieved. Because of the complexity of the experimental signals due to the presence of noise, such a problem causes a computational error and a stopping of the program process. In the proposed stopping criterion in equation (12), the operator ∑ is applied separately to the face and the denominator at any given time interval. Then, the two values are divided and the problem does not arise, yielding better results. In the second step, we also use the signal decomposition process, which is designed to create upper and lower coverage curves instead of using cubic spline smoothly (Franco et al., 2012). The cubic spline for n points is defined as $s_{i} (x) = y_{i} + b_{i} (x - x_{i}) + c_{i} {(x - x_{i})}^{2} + d_{i} {(x - x_{i})}^{3}$ ; x ∈ (x_i, x_i+1). To calculate the unknown coefficient, two conditions of continuous and smoothness are used.

3.3.1. Continuous conditions

The function must be continuous, so for the inner points, that is, the curve fitting points as S_i(x_i+1) = S_i + 1(x_i+1) = y_i+1, i = 1, 2, …, n − 2, and for the endpoints x_n, we have S_n−1(x_n) = y_n.

3.3.2. Smoothness conditions

There should be one in the interior of the slope of the curves. That is, the first derivative of these curves is equal to x at these points. So for inner points, we have S^′_i(x_i + 1) = S^′_i + 1(x_i + 1); i = 1, 2, …, n − 2. The curves inside the concave points should be the same. That is, the second derivative of these curves is equal to x at these points. So for inner points, we have $S_{i}^{''} (x_{i + 1}) = S_{i + 1}^{''} (x_{i + 1})$ ; i = 1, 2, …, n − 2. As can be seen, the cubic spline is a flexible method, and without any curve control created by this method, it will pass through all the data, whereas the measured signals are subject to instantaneous changes and unwanted variations. Under these conditions, the cubic spline method will follow these fluctuations without the slightest flexibility and will not affect noise reduction effects. To fix this, we use the smooth spline for n points (x₁, y₁), (x₂, y₂), …, (x_n, y_n) as follows

S (x) = β_{0} + β_{1} x + β_{2} x^{2} + β_{3} x^{3} + \sum_{i = 1}^{n} θ_{j} {(x - x_{i})}^{3}

(34)

where S is a third-degree polynomial in each interval (x_i, x_(i+1)) and has the first and second derivative continuous. To define this function, we need to calculate the number of n + 6 computation parameters. To do this, the functions must minimize as

p \sum_{i = 1}^{n} {(y_{i} - S (x_{i}^{2}))}^{2} + (1 - p) \int_{x_{1}}^{x_{n}} {(d^{2} S / d x^{2})}^{2} d x

. The balance between the two sentences is determined by the smooth p parameter and determines how smooth or close the data are to the spline. The smoothing parameter is defined as 0 and 1, where p = 0 is the straight line with the least squares method, and p = 1 corresponds to the cubic spline. Experimental and simulation results at p = 1/(1 + h³/6) give the best results where h is the mean distance of points. By minimizing the above expression, the function S is obtained. Smooth spline produces smoother results by providing control over how the curve passes through the data, and it is more resistant to unwanted variations and instantaneous changes because of the presence of noise and reduces the effects of noise. For more detailed analysis, modified EMD and EMD can be compared in simulations.

3.4. Multicomponent signal detection using the Neyman–Pearson criterion

It assumes that the noisy signal is entered into the system. The waveform or signal received (y) consists of the signal and noise or noise alone (Daubechies et al., 2011) as

x [n] = y [n] = \sum_{i = 1}^{p} a_{L} s_{i L} + w [n]

(35)

where the received signal is a set of randomly scaled component p, signal space {s₁, s₂, …, s_N} is orthogonal, and N is subspace dimension. The number p denotes the number of components of the signal considered as signal components: {s_i1, s_i², …, s_ip}. Assume that the maximum p is smaller than N and denote them by B, where B ≤ N. The following three conditions are considered to continue:

If p > 0, the signal is detected.

If p > 0, the signal power estimation is used to determine the actual number of signal components as P ∈ 1, 2, …, B.

If p = p₀ > 0, the signal components’ estimation is conducted by identifying the signal components present and the availability of p₀. The equation (35) can be represented as an optimal structure by $M = \sum_{p = 1}^{B} (\begin{matrix} N \\ P \end{matrix})$ equation. Considering hypothesis assumptions of likelihood to estimate and detect multicomponent signals with the components of {s_i1, s_i2, …, s_ip} under i_p ∈{1, 2, …, N}, p = {1, 2, …, B}, and M is likelihood ratio.

We use the weighted maximum likelihood estimator to implement a weighted generalized likelihood ratio test (WGLRT) with a randomized threshold. The optimal rank selection structure uses the weighted averages of B with simulated ratios of $(\begin{matrix} N \\ P \end{matrix})$ . Each mean is associated with a constant p of the signal components. The optimal structure detector compares the weighted average of all ratios similar to the final M with a threshold. In each of the three conditions stated, weights and thresholds of detection shall be determined by solutions to nonlinear optimization problems and false alarm rate α. According to the Neyman–Pearson (NP) criterion, the optimal detector is a detector that maximizes the detection probability for the likelihood of constant false alarm as S = ∑_{i = 1}Na_ie ^jϕi, S = Ae ^jϕS_k = a_ke ^jpk, and k = 0, 1, …, N − 1. The detection problem is equivalent to testing the hypothesis model, where H₀ = y = n and H₁ = y = s + n are no signal hypothesis, and signal presence hypothesis, respectively. Detector test according to NP criterion as

\nabla (y) = \frac{f_{y} (y | H_{1}) \geq T : H_{1}}{f_{y} (y | H_{0}) \leq T : H_{0}}

(36)

where f_y(y|H₁) is the probability density function of the received vector providing the signal exists, and f_y(y|H₀) is the probability density function of the received vector providing no signal. In the equation (36), T is the threshold level, which is determined by the optimal value of P_Fa and ∇(y) is the final correct function. If the probability density functions of the desired signal vectors and the noise are f_s(s) and f_n(n), then the equation (34) becomes

\begin{array}{l} \nabla (y) = \frac{\int f_{n} (y - s) f_{s} (s) d s \geq T : H_{1}}{f_{n} (y) \leq T : H_{0}} \\ \Rightarrow GLRT = \nabla (y | s_{m L}) = \frac{F_{n} (y - s_{M L}) \geq T : H_{1}}{F_{n} (y) \leq T : H_{0}} \end{array}

(37)

where instead of integrating the random parameter into the generalized likelihood ratio test (GLRT) estimation, it replaces the optimal MLE vectors. When no statistical information is available from the signal, the generalized likelihood ratio (GLR) overlaps with the average likelihood ratio (ALR). The accuracy and performance of the GLR depend directly on the accuracy and performance of the ALR estimation. Therefore, given MLE’s asymptotic properties, it can be said that the GLR test is asymptotically consistent with ALR (Feldman, 2006). Mono-component signals’ behavior differs from that of the multicomponent, in contrast to the time–frequency domain characterized by a single edge in a concentrated energy region (Liutkus et al., 2017). For the real signal s(t), the analytical equivalent of z(t) is defined as follows

z (t) = s (t) + j H {s (t)} = α (t) e^{j ϕ (t)}

(38)

where H{s(t)} is the Hilbert transform, s(t)α(t) is the instantaneous amplitude of the signal, and ϕ(t) the signal phases. The IF involves the frequency variations of the signal over time. In the frequency-modulated (FM) signal, IF often involves low values (Feldman, 2006). The IF of the mono-component signal z(t) is equal to the first derivative of its instantaneous phase and is expressed as ω(t) = ϕ^′(t). In addition, the maximum edge value for the IF estimation of the z(t) signal is defined as follows (Liutkus et al., 2017)

w (t) = \underset{f}{argmax} | TF D_{z} (t, f) |

(39)

where TFD_z(t, f) is the time distribution of the signal frequency z(t) (Liutkus et al., 2017). In other words, the analytic multicomponent signal x(t) can be modeled as a sum of two or more than two mono-component signals

x (t) = \sum_{m = 1}^{M} z_{m} (t) = \sum_{m = 1}^{M} α_{m} (t) e^{j ϕ_{m} (t)}

(40)

where z_m(t) is the IF of any mono-component signal. In equation (40), M is the number of signal components; α_m(t), the instantaneous amplitude of the component Mth; and ϕ_m(t), the instantaneous phase. When calculating the FT of the signal s(t) in equation (38), α(t) must be a low-frequency function with a spectrum that has no overlap with the spectrum e^jϕ(t) (Biagetti et al., 2015). To determine the IF of the multicomponent signals, the component separation instruction must be performed faster than the IF estimation of the extracted signal components. However, when dealing with multicomponent signals, their TFD often includes a cross-term that significantly disrupts the discussion of time–frequency signals (Liutkus et al., 2017). Therefore, the component separation procedure is relatively more difficult. Therefore, selecting the appropriate TFD plays a key role in the efficiency and efficiency of the extraction of signal components. Various interference attenuation reduced interference distributions (RID) have been proposed to have high-resolution frequency time signals, such as modified B distribution (MBD) and the Bessel kernel-based RID Pham and Meignen (2017). Measurement and for time–frequency resolution and component decomposition has been proposed in Herrera et al. (2015). Methods for extracting signal components for two or more independent signals are often statistically expressed as blind source separation (BSS). The blind expression indicates that neither the component structure nor the source signals are present. They are not known in advance (Liutkus et al., 2017). So, the main form of BSS is to determine the main waveform of resources. When only their combination is available (Tary et al., 2014) because of a wide range of potential applications, BSS is looking to itself. The results are based on multiple BSS techniques that can be classified into time domain methods Iatsenko et al. (2015). The IF estimation methods for noisy signals can be divided into two categories: For the signal in the presence of multiplying noise or a time-varying amplitude signal, the Wigner–Ville spectrum, or the Wigner–Ville polynomial distribution can also be used as proposed in Wu (2011). For the FM polynomial signals in the presence of high noise and high SNR, the Wigner–Ville polynomial distribution is proposed based on the IF estimation method (Oberlin et al., 2015). A simple flowchart of the new multicomponent IF estimation method is shown in Figure 2. As can be seen, the multicomponent IF estimation using the proposed method of extracting the modified components is quite evident. In the component extraction procedure, the first step is to select TFD. When faced with multicomponent signals, TFD selection plays an essential role because of the presence of unwanted cross-terms. It should be noted that these cases always disrupt the time–frequency signal display. The best known TFD is a linear mono-component with WVD of analytical FM signal that can be defined by equation (41) (Amin et al., 2019)

WV D_{z} (t, f) = \int_{- \infty}^{+ \infty} z (t + \frac{τ}{2}) \cdot z^{*} (t - \frac{τ}{2}) \cdot e^{- j 2 π f τ} d τ

(41)

Figure 2.

Simple flowchart of the new instantaneous frequency estimation algorithm.

One of the WVD major disadvantages of multicomponent or mono-component by nonlinear IF is the presence of interference and frequency resolution losses (Liutkus et al., 2017). To reduce the cross-term in WVD, the instantaneous signal autocorrelation function can be used as $z (t + τ / 2) \cdot z^{*} (t - τ / 2)$ , which can be windowed in the direction of lag τ. Of course, before getting the FT as

PWV D_{z} (t, f) = \int_{- \infty}^{+ \infty} h (τ) z (t + \frac{τ}{2}) \cdot z^{*} (t - \frac{τ}{2}) \cdot e^{- j 2 π f τ} d τ

(42)

where h(τ) is a symmetric time window. In fact, this windowing is equivalent to smoothing the distribution along the frequency axis (

PWV D_{z} (t, f) = \int_{- \infty}^{+ \infty} H (v - ξ) WV D_{z} (t, ξ) d ξ

), where H(v) is the FT of h(t). Pseudo WVD is defined as follows for the desired smoothing of time and frequency independently of WVD (Sharma et al., 2017)

SPWV D_{z} (t, f) = \int_{- \infty}^{+ \infty} h (τ) \int_{- \infty}^{+ \infty} g (s - t) z^{*} (s - \frac{τ}{2}) e^{- j 2 π f τ} d s \cdot d τ

(43)

where h(τ) is the function of frequency smoothing window in the time domain, and g(t) is time smoothing function. The efficiency of the IF method presented in the current study depends on the proper selection of TFD. Therefore, high-resolution TFD should be used to reduce interference. There are several TFDs that have such features, some of which are defined in Wang et al. (2014). The MBD method can be used to reduce the cross-term and enhance the resolution (Daubechies et al., 2011)

\begin{array}{l} MB D_{z} (t, f) = \int \int_{- \infty}^{+ \infty} \frac{\cos h^{- 2 β} (t - τ)}{\int_{- \infty}^{+ \infty} \cos h^{- 2 β}} \\ \cdot z (u + \frac{τ}{2}) z^{*} (u - \frac{τ}{2}) e^{- j 2 π f τ} d τ \end{array}

(44)

where β parameter (0 ≤ β ≤ 1) is used to control the resolution of the cross-term distribution and eliminate (Wang et al., 2014). In general, there is a compromise between the TFD and MBD characteristics, which makes it less usable (Li and Liang, 2012). Also, MBD is a way to find the appropriate TFD to continue IF estimation (Li and Liang, 2012). With these features in mind, we try to extend the results obtained using MBD and compare it with other methods. One method is RID with kernel based on Bessel function RIBD (Coifman and Wickerhauser, 1992), which is used to eliminate cross-term with good performance in time–frequency resolution and independent window in τ domain v (Wu, 2011)

RID B_{z} (t, f) = \int_{- \infty}^{+ \infty} h (τ) R_{z} (t, τ) e^{- j 2 π f τ} d τ

(45)

where

R_{z} (t, τ) = \int_{- | τ | / 2}^{| τ | / 2} g (v) / | τ | (1 + \cos 2 π v / τ) x (t + v + | τ | / 2) x^{*} (t + v - | τ | / 2) d v

, t is time, f is frequency, z* is complex conjugation, and h(τ) and g(v) are time and frequency smoothing window, respectively.

3.5. Improved IF estimation method using sliding window

First, consider a nonstationary discrete multicomponent signal in the presence of additive noise as y(n) = x(n) + ϵ(n), where

x (n) \sum_{m = 1}^{M} z_{m} (n) = \sum_{m = 1}^{M} α_{m} (n) e^{j ϕ_{m} (n)}

(46)

where, M is the number of signal components, mixed Gaussian white noise with independent real and imaginary parts with mean zero and variance of

d_{ϵ}^{2} / 2

, and instantaneous amplitude of component mth a_m(n) and instantaneous phase ϕ_m(n). The component TF can be estimated from the TFD signal as

{\hat{ω}}_{m} (n, h) = \underset{k}{argmax} | TF D_{m} (n, k, h) |

(47)

where TFD_m(n, k, h); TFD comprises the component extracted from the multicomponent signal calculated using the h-length window. Using the experiments in Wang et al. (2014), the IF estimation error is calculated as follows

Δ {\hat{ω}}_{m} (n, h) = ω_{m} (n) - {\hat{ω}}_{m} (n, h)

(48)

| Δ {\hat{ω}}_{m} (n, h) | \leq | bia s_{m} (n, k) | + k f_{m} (h)

(49)

The probability P(k), where k is the quantile standard Gaussian distribution and d_m(h) is the standard deviation of the component estimation error as follows

d_{m} (h) = \sqrt{\frac{d_{ϵ}^{2}}{2 {| A_{m} |}^{2}} (1 + \frac{d_{ϵ}^{2}}{2 {| A_{m} |}^{2}}) \frac{T}{h^{3}} \frac{E}{F^{2}}}

(50)

where,

E = \sum_{i = 1}^{h} w (i) {(i / h)}^{2}

F = \sum_{i = 1}^{h} {[w (i)]}^{2} {(i / h)}^{2}

, and w(n) is a symmetric window with length h. Values of F and E depend on the type of window used to calculate TFD. For example, in the rectangular window E = F = 1/12 as seen in Wang et al., 2014, |bias_m(n, k)| ≤ kd_m(h). Therefore, the result of equation (49) is as follows

| Δ {\hat{ω}}_{m} (n, h) | \leq 2 k d_{m} (h)

(51)

Equations (48) and (51) state that w_m(n) belongs to the safe interval D_m(n·L) = (L_m(n·L)·U_m(n·L)). With a probability of greater P(k), always P(k) results in near one. The upper range U_m(n·L) and the lower L_m(n·L) are defined as

\begin{matrix} U_{m} (n, L) = [{\hat{ω}}_{m} (n, h) + 2 k α_{m} (h)] \\ L_{m} (n, L) = [{\hat{ω}}_{m} (n, h) - 2 k α_{m} (h)] \end{matrix}

(52)

where L is the number of sequences h in the set of window width increments H = {h₁|h₁ < h₂ < … < h_j}. The IF estimation method proposed in Feldman (2006) is a sequence of TFD that calculates H for each window width. In this study, using the proposed method H = {h₁|h₁ = k_L−1 + 2}, similar to what Sharma et al. (2017) has stated, we consider

h_{1} = N / 8

and

h_{j} = N / 2 + 1

. Then, the component separation and extraction procedure are performed. The IF estimation is then computed using equation (47), for each signal component, which follows the D_m(n, L) time interval calculations for each time constant nT (T is time or sampling period), and window width h. This adaptive approach always follows the conflicts between D_m(n·L) and D_m(n, L − 1). Two-time intervals and opposite zero denote the best window width for any nT time constant as

D_{m} (n, L) \cap D_{m} (n, L - 1) \neq 0

(53)

Based on the equation (51) inequality, conditions of these two methods are not fully available, causing the bias estimation to be too large compared with the variance and not producing the 50 unequal conditions. According to the material, k plays an important role in current calculation tasks, including proper window size calculations as well as estimation accuracy (Gianfelici et al., 2007). Various computational methods have been used in this area. As noted in Jones and Baraniuk (1994), smaller k values result in windows of much shorter width, when k values increase (P(k) → 1). The window width is higher than usual and affects the accuracy of estimation. One way to improve the application of the window’s appropriate width is to follow the interval overlap. For this purpose, we introduce C_m(n, L) as the overlap value between two intervals

C_{m} (n, L) = | D_{m} (n, k) \cap D (n, L - 1) |

(54)

To have an overlap size of two finite intervals, C_m(n, L) must be normalized by the current interval size

O_{m} (n, L) = \frac{c_{m} (n, L)}{| D_{m} (n, L) |}

(55)

where the value of O_m(n, L), unlike C_m(n, L), always belongs to the interval [0, 1], so it is easy to apply the threshold value of O_c as an additional criterion to select the most appropriate window width as O_m(n, l) > O_c, where

O_{m} (n, l) = {\begin{array}{l} 0, & c_{m} (n, L) = 0 \\ 1, & C_{m} (n, L) = | D_{m} (n, L) | \\ 〈 0,1 〉 & otherwise \end{array}

(56)

4. Results and discussion

Time–frequency analysis (TFA) is potentially applied to nonstationary signals. Various TFA methods, including STFT, quadratic EMD distributions based on HHT, and continuous wavelet conversions, have been proposed in recent years to deal with nonstationary signals. Figure 3 illustrates a flowchart of the IF estimation process using the demodulated signal extraction method (Daubechies et al., 2011). A mono-component signal is described in the (t, f) domain by one single ridge, corresponding to an elongated region of energy concentration. Furthermore, interpreting the crest of the ridge as a graph of IF versus time, we require the IF of a mono-component signal to be a single-valued function of time. Such a mono-component signal s(t) has an analytic associate of the form

z (t) = a (t) e^{j ϕ (t)}

(57)

where a(t), known as the instantaneous amplitude, is real and positive; ϕ(t), known as the instantaneous phase, is differentiable; and a(t) and e^jϕ(t) are spectrally disjoint. In addition, s(t) itself is real and asymptotic with amplitude modulation a(t), it can also be expressed approximately as x(t) = a(t)cosϕ(t). A multicomponent signal may be described as the sum of two or more mono-component signals such that (Khan and Ali, 2020)

z (t) = \sum_{i = 1}^{M} a_{i} (t) e^{j ϕ (t)}

(58)

Figure 3.

Flowchart of the instantaneous frequency estimation process using the demodulated signal extraction method.

The decomposition into components is not necessarily unique. This model allows the extraction and separation of components from a given multicomponent signal using (t, f) filtering methods (Khan and Ali, 2020). Given a real signal s(t), we can construct the complex signal

z (t) = s (t) + j H {s (t)} = α (t) e^{j ϕ (t)}

(59)

where z(t) is analytic, and ϕ(t) = 2πf_ct + ψ. Gabor complex signal and complex signal spectral are defined in equations (60) and (61), respectively, as

z (f) = \int_{- \infty}^{+ \infty} z (t) e^{- j 2 π f t} d t

(60)

z (f) = \int_{- \infty}^{+ \infty} α (t) e^{[ϕ (t) - 2 π f t]} d t

(61)

when the above integral is maximized when the derivative ψ(t) = ϕ(t) − 2πf_ct is zero (Stankovic et al., 2014)

\begin{array}{l} \frac{d ψ}{d t} = 0 \to \frac{d}{d t} [ϕ (t) - 2 π f_{c} t] = 0 \\ \to f_{c} = \frac{1}{2 π} \frac{d ϕ (t)}{d t} \Rightarrow f_{i} (t) = \frac{1}{2 π} \frac{d}{d t} [\arg z (t)] \end{array}

(62)

Then, E{F} = E{F_i} means mean frequency is equal to the mean time of IF (Stankovic et al., 2014)

E {F^{n}} = \frac{\int_{- \infty}^{+ \infty} f {| z (F) |}^{2} d f}{\int_{- \infty}^{+ \infty} f {| z (F) |}^{2} d f}

(63)

E {F_{i}} = \frac{\int_{- \infty}^{+ \infty} f_{i} (t) {| z (t) |}^{2} d t}{\int_{- \infty}^{+ \infty} {| z (t) |}^{2} d t}

(64)

where E{F_i} is time averaging and E{Fⁿ} is frequency averaging. Using WVD time and frequency signal distribution, the IF can be calculated as follows

ω (t, f) = \int_{- \infty}^{+ \infty} z (t + \frac{τ}{2}) z^{*} (t - \frac{τ}{2}) e^{- j 2 π f τ} d τ

(65)

where ω(t, f) is WV given function. And

f_{i} (t) = \frac{\int_{- \infty}^{+ \infty} f ω (t, f) d f}{\int_{- \infty}^{+ \infty} ω (t, f) d f}

(66)

This is a new method using WVD. From equations (64) and (65), we have

\begin{array}{l} E {f_{i}} = \frac{\int_{- \infty}^{+ \infty} [\int_{- \infty}^{+ \infty} f ω (t, f) d f / \int_{- \infty}^{+ \infty} ω (t, f) d f]}{\int_{- \infty}^{+ \infty} {| z (t) |}^{2}} \\ = \frac{\int_{- \infty}^{+ \infty} \int_{- \infty}^{+ \infty} f ω (t, f) {| z (t) |}^{2} d f d t}{\int_{- \infty}^{+ \infty} \int_{- \infty}^{+ \infty} ω (t, f) {| z (t) |}^{2} d f d t} \\ = \frac{1}{2 π L} \int_{- \infty}^{+ \infty} \int_{- \infty}^{+ \infty} ω (t, f) {| z (t) |}^{2} d f d ϕ (t) \\ = \frac{1}{2 π L} \int_{- \infty}^{+ \infty} \int_{- \infty}^{+ \infty} ω (t, f) {| z (t) |}^{2} d f [2 π f_{i} (t) d t] \\ = \frac{1}{L} \int_{- \infty}^{+ \infty} \int_{- \infty}^{+ \infty} f_{i} (t) {| z (t) |}^{2} d t ω (t, f) d f \end{array}

(67)

where

L = \int_{- \infty}^{+ \infty} \int_{- \infty}^{+ \infty} ω (t, f) {| z (t) |}^{2} d f d t

. Also, equation (65) can be rewritten as follows

f_{i} (t) = \frac{1}{2 π} \cdot \frac{d ϕ (t)}{d t} = \frac{1}{2 π} \sum_{i = 2}^{n} [\frac{ϕ (t_{i}) - ϕ (t_{i - 1})}{t_{i} - t_{i - 1}}]

(68)

and

S (t_{i}) = \sum_{i = 2}^{n} x_{1} (t_{i} - 1) x_{1} (t_{i})

(69)

The average of f_i(t) is obtained from the following equation as

f_{i} ave (t) = \frac{1}{n} \sum_{i = 1}^{n} f (t_{i} - i)

(70)

To achieve the signal phase, we need to define the imaginary part of the signal, which is achieved by using the following equation as

\begin{matrix} z_{1} (t) = x_{1} (t) + j H {x_{1} (t)} = \\ A_{1} (t) e^{- j ϕ_{1} (t)} = Arctan [\frac{H {x_{1} (t)}}{x_{1} (t)}] \end{matrix}

(71)

where x₁(t) = A₁(t)cos ϕ₁(t). Frequency-modulated demodulated final part is defined as D(t) = d_m(t), where d_m(t) = d_m−1(t)/e_m(t). In this way, the instantaneous amplitude of the signal will be determined by the following equation as

A_{1} (t) = \frac{x_{1} (t)}{D_{1} (t)}

(72)

where according to equations (71) and (72), we will have

z_{1} (t) = \frac{x_{1} (t)}{D_{1} (t)} e^{- j ϕ_{1} (t)}

(73)

We define $x_{i}^{(d)} (t)$ of the difference between s_i(t) and s_i−1(t). Therefore, signal x(t) can be rewritten as follows

x (t) = \sum_{i = 1}^{n} x_{i}^{(d)} (t) = \sum_{i = 2}^{n} s_{i} (t) - s_{i - 1} (t)

(74)

where

\begin{matrix} s_{i} (t) = \sin (ϕ_{l i} (t)) H {x (t) \cos (ϕ_{l i} (t))} \\ - \cos (ϕ_{l i} (t)) H {x (t) \sin (ϕ_{l i} (t))}, i = 1,2, \dots, n \end{matrix}

(75)

Equation (75) can also be expressed as follows

\begin{array}{l} s_{i} (t) = x (t) \sin (ϕ_{l i} (t)) H {\cos (ϕ_{l i} (t))} \\ - x (t) \cos (ϕ_{l i} (t)) H {\sin (ϕ_{l i} (t))} \\ = x (t) \times (\sin (ϕ_{l i} (t)) H {\cos (ϕ_{l i} (t))} \\ - \cos (ϕ_{l i} (t)) H {\sin (ϕ_{l i} (t))}) \end{array}

(76)

where H{cos(ϕ_li(t))} = sin(ϕ_li(t)) and H{sin(ϕ_li(t))} = cos(ϕ_li(t)). The signal s(t) is analytic with a real DC component, if

H {X (f)} = (- j sgn f) X (f)

(77)

where X(f ) and H{X(f )} are the FTs of x(t) and H{x(t)}, respectively, and where

sgn (s_{i} (t)) = {\begin{array}{l} 1 & s_{i} (t) > 0 \\ 0 & s_{i} (t) = 0 \\ - 1 & s_{i} (t) < 0 \end{array}

(78)

Assume that ϕ_li(t) = ω_i(t), therefore

\begin{matrix} ω_{i} (t) = Arctan H {\tan ω_{i} (t)} \\ = H {Arctan (\tan ω_{i} (t))} \Rightarrow ω_{i} (t) = H {ω_{i} (t)} \end{matrix}

(79)

where as a result f_i(t) = H{f_i(t)} =

\int_{- \infty}^{+ \infty} f_{i} (τ) / π (t - τ) d τ

, and f_i(t) for i = 1, 2 is f₁(t) + f₂(t). According to equation (77), we have

\begin{matrix} f_{i} (t) = H {f_{0} (1 + m g (t))} \\ = f_{0} + f_{0} m H {g (t)} = f_{0} [1 + m H {g (t)}] \end{matrix}

(80)

and the signal phase is defined as follows

\begin{matrix} ϕ (t) = arg (z (t)) = \arctan \frac{H {x (t)}}{x (t)} + arg a \\ = \arctan \frac{H {A (t) D (t)}}{A (t) D (t)} = \arctan \frac{H {A (t)}}{| z (t) |} \end{matrix}

(81)

where f_i(t) is written consequently as

f_{i} (t) = \frac{d}{d t} [\arctan \frac{H {A (t)}}{| z (t) |}]

(82)

Similarly, z(t) can be written as $A_{1} (t) e^{- j ϕ_{1} (t)} + A_{2} (t) e^{- j ϕ_{2} (t)}$ , where ϕ₁(t) = (ω₀ − θ) (t), ϕ₂(t) = (ω₀ + θ) (t), and ω₀ = ϕ(t) ± θ(t). On the other hand, given that z(t) = x(t) + jy(t), and equation (62), f_i(t) is rewritten as follows

\begin{matrix} f_{i} (t) = \frac{x' (t) y (t) - x (t) y' (t)}{x^{2} (t) + y^{2} (t)} \\ = \frac{1}{2 π} \frac{d ϕ_{i} (t)}{d t} ≃ \frac{ϕ_{i} (t + T) - ϕ_{i} (t)}{2 π t} \end{matrix}

(83)

where ϕ_i(t) = arg(z(t)), so

\begin{array}{l} f_{i} (t) ≃ \frac{\arg z (t + T) - \arg (z (t))}{2 π t} \\ = \frac{1}{2 π t} \arg [\frac{z (t + T)}{z (t)}] \end{array}

(84)

f_{i} (t) ≃ \frac{1}{2 π T} Arctan [\frac{x (t) y (t + T) - x (t + T) y (t)}{x (t) y (t + T) + x (t + T) y (t)}]

(85)

If A = z(t + T) = x(t + T) + jy(t + T) and B = z(t) = x(t) + jy(t), define A and B as

{\begin{array}{l} A = a_{1} (t + T) + j b_{1} (t + T) \\ B = a_{2} (t) + j B_{2} (t) \end{array}

(86)

given that

\arg [\frac{A}{B}] = [\frac{a_{2} b_{1} - a_{1} b_{2}}{{(a_{1} + a_{2})}^{2} + {(b_{1} + b_{2})}^{2}}]

(87)

Given a quadrilateral circle, the equation (85) can be rewritten with approximation as follows

f_{i} (t) ≃ \frac{1}{2 π T} \times \frac{s (t) H {s (t + T)} - s (t + T) H {s (t)}}{{(s (t) + s (t + T))}^{2} + {(H {s (t)} + H {s (t + T)})}^{2}}

(88)

Because ϕ_Li = H{ϕ_Li}, we will have

s (t) = {s (t)} \to s (t + T) ≃ H {s (t + T)}

(89)

According to equation (89), it can be written as

\begin{matrix} f_{i} (t) ≃ \frac{1}{2 π T} \frac{s (t) [H {s (t + T) - s (t + T)]}{{(s (t) + s (t + T))}^{2} + {(s (t) + H {s (t + T})}^{2}} \\ = \frac{1}{2 π T} \frac{s (t) [(1 / T - 1) s (t + T)]}{{(s (t) + s (t + T))}^{2} + {(H {s (t)} + H {s (t + T)})}^{2}} \\ \frac{1}{2 π T^{2}} \frac{(1 - T) S (T) S (t + T)}{{(s (t) + s (t + T))}^{2} + {(H {s (t)} + H {s (t + T)})}^{2}} \end{matrix}

(90)

In Figure 4, the three main signal frequency components are extracted to estimate the IF, and by obtaining the computational error shown in Figure 5, the accuracy of the distributive separation is increased. The theoretical results are always closer to the practical. In Figure 6, signal representation using evaluation criteria with different SNR is shown that it is applied for three modes. In IF calculations for static signals, the instantaneous phase derivative can be used, but in nonstationary signals, this operation causes a computational error and can make it difficult to estimate the measurement parameter accurately. In the current approach, the goal is not to optimize the input signal, but to use mathematical functions to optimize the IF of nonstationary signals to reduce computational time. Analysis, evaluation and compensation of component-errors, was simulated using SNR, the optimal variance for the least complexity of subsystem from evaluation criterion. Because of the proximity of SNR in mode-1 and mode-2, we see fewer fluctuations than mode-3. As a result, separating the components in two modes with close SNR will reduce the component error in the cross-term. Figure 7 shows the EMD method’s algorithm separating the signal components from which the proposed method is adapted. The proposed method is compared with the EMD method, and according to the obtained results, the relative stability of the new method is clearly evident. In the block diagram of Figure 8, the proposed method algorithm is displayed. The proposed method in the structure is very similar to the experimental mode analysis method, with the difference that to increase the accuracy in the analysis of signal components, a step has been added to the mean calculations that despite the simulation time, we will always achieve the expected results. In the proposed method, two stop criteria SD_k and ST_k were used using equations (32) and (33). The stop condition was considered less than 0.1 for both parameters SD_k and ST_k. The Figure 9 shows the process of decomposing the input signal components into the system using the proposed algorithm. The relationships required to implement and simulate the proposed method are presented below. Consider the input signal x(t) varying with time t to the subsystem is defined by equation (91) as

x (t) = c (t) + r (t)

(91)

where residual term r(t) is the trend or the baseline signal. Assume that c(t) is an oscillatory function with mean zero, its integral in the interval of two local maximum (or minimum) points t_k and t_k+2 should be equal to zero

E [c (t)] = \int_{t_{k}}^{t_{k + 2}} c (t) d t = 0

(92)

Figure 4.

Extraction of the three main signal components to estimate the instantaneous frequency.

Figure 5.

Simulation of computational error resulting from the component separation.

Figure 6.

Signal representation capability using evaluation criteria.

Figure 7.

Flowchart representing the algorithm of empirical mode decomposition technique.

Figure 8.

Flowchart representing of the proposed algorithm.

Figure 9.

Part of the signal x(t) analyzed.

Therefore, from equations (94) and (95), we can conclude

\begin{array}{l} \int_{t_{k}}^{t_{k + 2}} x (t) d t = \int_{t_{k}}^{t_{k + 2}} c (t) d t + \int_{t_{k}}^{t_{k + 2}} r (t) d t \\ = \int_{t_{k}}^{t_{k + 2}} x_{mean} (t) d t = r_{0} (t_{k + 2}) - r_{0} (t_{k}) \end{array}

(93)

where r₀(t) is the original function of r(t) and a real-valued function defined on interval (t_k, t_k+2); according to the Lagrange differential theorem of mean, it is known that there is at least one point t₁ between t_k and t_k+2 to make

R (t_{k + 2}) - R (t_{k}) = (t_{k + 2} - t_{k}) R' (t_{1}) = (t_{k + 2} - t_{k}) r (t_{1})

(94)

Using equations (93) and (94), it can be rewritten as

\int_{t_{k}}^{t_{k + 2}} x (t) d t = (t_{k + 2} - t_{k}) r (t_{1})

(95)

Similarly, it can be derived that there is a point t_m between t_k+1 and t_k+3 to make

\int_{t_{k + 1}}^{t_{k + 3}} x (t) d t = (t_{k + 3} - t_{k + 1}) r (t_{m})

(96)

Choosing a point C located at the central point from t_k to t_k + 1, as shown in Figure 10, can derive new relations from equations (97) and (98), respectively, given by

\begin{array}{l} \int_{t_{k + 1}}^{t_{k + 2}} x (t) d t = \int_{t_{k}}^{t_{k} + Δ t} x (t) d t + \int_{t_{k} + Δ t}^{t_{k + 2} + Δ t} x (t) d t \\ - \int_{t_{k + 2}}^{t_{k + 2} + Δ t} x (t) d t = (t_{k + 2} - t_{k}) r (t_{1}) \end{array}

(97)

\begin{array}{l} \int_{t_{k}}^{t_{k + 3}} x (t) d t = \int_{t_{k + 1} - Δ t}^{t_{k + 1}} x (t) d t + \int_{t_{k + 1} - Δ t}^{t_{k + 3} - Δ t} x (t) d t \\ + \int_{t_{k + 3} - Δ t}^{t_{k + 3}} x (t) d t = (t_{k + 3} - t_{k + 1}) r (t_{m}) \end{array}

(98)

where Δt = (t_k+1− t_k)/2. To introduce the values of extrema at four points t_k, t_k+1, t_k+2, and t_k+3, one can use the following approximate relations

\int_{t_{k}}^{t_{k} + Δ t} x (t) d t - \int_{t_{k + 2}}^{t_{k + 2} + Δ t} x (t) d t \approx - Δ t (x_{k + 2} - x_{k})

(99)

- \int_{t_{k + 1} - Δ t}^{t_{k + 1}} x (t) d t + \int_{t_{k + 3} - Δ t}^{t_{k + 3}} x (t) d t \approx Δ t (x_{k + 3} - x_{k + 1})

(100)

\frac{\int_{t_{k} + Δ t}^{t_{k + 2} + Δ t} x (t) d t}{t_{k + 2} - t_{k}} \approx \frac{\int_{t_{k + 1} - Δ t}^{t_{k + 3} - Δ t} x (t) d t}{t_{k + 3} - t_{k + 1}}

(101)

t_{m} - t_{1} \approx 2 Δ t

(102)

Figure 10.

Approximate to integrals over different intervals.

According to approximation r(t) ≅ x_mean(t) and equation (102), we can write as follows

r (t) \approx \frac{x_{\max} (t) + x_{\min} (t)}{2}

(103)

Using equation (103), r₁(t) and r₂(t) are obtained according to equations (104) and (105), respectively, as

r (t_{1}) = \frac{r (t) + x_{\max} (t)}{2}

(104)

r (t_{2}) = \frac{r (t) + x_{\min} (t)}{2}

(105)

Using equation (95) and based on the Lagrange differential theorem of mean, at least one point between τ_k and τ_k+2 can be considered, and equation (103) can be obtained. With due regard to the behavior of r(t), x_max(t), and x_min(t) uniformly varying with time over the interval (t_k, t_k+3), one has

r (t) = \frac{r (t_{1}) + r (t_{2})}{4} \approx x_{mean} (t)

(106)

Using equations (93) and (96), it can be rewritten as follows

\int_{τ}^{τ_{k + 2}} x (t) d (t) = \frac{1}{4} \int_{τ}^{τ_{k + 2}} r (t_{1}) + r (t_{2}) d t

(107)

where r(t₁) and r(t₂) are obtained from the following equations as

r (t_{1}) = \frac{x_{mean} (t) + x_{\max} (t)}{2} = \frac{3}{4} x_{\max} (t) + \frac{1}{4} x_{\min} (t)

(108)

r (t_{2}) = \frac{x_{mean} (t) + x_{\min} (t)}{2} = \frac{3}{4} x_{\min} (t) + \frac{1}{4} x_{\max} (t)

(109)

As can be seen from the above equations, instead of using x_mean(t) in signal component decomposition relations, as discussed in the EMD method, three-quarters of the signal cap can be used to decompose the components, which with these changes is always faced with increasing computational iterations and consequently increasing computational time, but has far less phase error than the empirical mode method, which is one of the advantages of the proposed method over the experimental mode method. Figure 11 shows the input signal from the linear combination of three sine functions in the time domain to extract frequency components and simulate IF. Figure 12 compares the simulation of the new method and the EMD method.

Figure 11.

Linear combination of three input signals to the subsystem considering frequency components f_h = 200, f_i = 45, f_L = 150, f_iL = 40, f_LL = 100, and ϵ = 1.

Figure 12.

The left column shows the results of applying the proposed method, while the right column shows the results of applying the experimental mode analysis method.

It should be noted that this simulation does not necessarily mean that the proposed method can always be better than the EMD method or other existing methods, but in nonlinear conditions and considering the design parameters in real conditions, it always has less phase error that is shown in Figure 13. Figure 14 compares the two proposed methods and the EMD method by simulating the IF of each component and the IF of the total input signals to the subsystem. It is shown that both methods have similar results to each other.

Figure 13.

Comparison of phase error of the proposed method and the empirical mode decomposition method.

Figure 14.

Simulation of instantaneous frequency of each component of input signals in the proposed method.

5. Conclusion

According to the mentioned theoretical and practical issues in this research, a new method of signal component analysis was presented to analyze the input signals to the subsystem more accurately based on the analysis of signal component, adding a step to mean and application calculations. Three-quarters of signal compression and its comparison with the EMD method is one of the most prominent empirical and efficient methods available. The signal can be subdivided into its constituents according to different component extraction methods to calculate the IF of a multicomponent signal. The existing method is then used to calculate the mono-component IF, and this method is applied to multicomponent signals. The component that overlaps the components is inefficient. Therefore, a hybrid algorithm based on evaluating different distribution criteria and comparing their performance in solving one or more features is used to investigate the frequency distribution such as separability and energy concentration. In general, the methods of estimating the IF of multicomponent signals are divided into three groups which include methods using signal phase derivative, which is sensitive to noise; methods that calculate the number of zero points of the signal and consider the signal frequency equal to half the frequency of the zero points and are suitable for signals that can be thought of as stationary; methods based on TFD and distributions such as WVD for IF calculation; and more for nonstationary noise signals for IF calculation. In this article, all three methods have been investigated, and a new method has been proposed concerning the obtained estimates, which far outweigh the drawbacks of other methods. According to the results, the proposed method always has a lower phase error than the EMD method.

Footnotes

Declaration of Conflicting Interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) received no financial support for the research, authorship, and/or publication of this article.

ORCID iDs

Milad Daneshvar

Pouria Salehi

References

Amin

Zhang

Himed

(2019) Improved instantaneous frequency estimation of multi-component FM signals. In: 2019 IEEE radar conference, Boston, USA, 22–26 April 2019. Boston, MA: IEEE.

Biagetti

Crippa

Curzi

, et al. (2015) Analysis of the EMG signal during cyclic movements using multicomponent AM-FM decomposition. IEEE Journal of Biomedical and Health Informatics 19(5): 1672–1681.

Coifman

Wickerhauser

(1992) Entropy-based algorithms for best basis selection. IEEE Transactions on Information Theory 38(2): 713–718.

Daubechies

(2011) Synchrosqueezed wavelet transforms: an empirical mode decomposition-like tool. Applied and Computational Harmonic Analysis 30(2): 243–261.

Feldman

(2006) Time-varying vibration decomposition and analysis based on the Hilbert transform. Journal of Sound and Vibration 295(3): 518–530.

Franco

Gumery

Vuillerme

, et al. (2012) Synchrosqueezing to investigate cardio-respiratory interactions within simulated volumetric signals. In: Proceedings of the 20th European signal processing conference (EUSIPCO), Bucharest, Romania, 27–31 August 2012, pp. 939–943. Bucharest, Romania: IEEE.

Gianfelici

Biagetti

Crippa

, et al. (2007) Multicomponent AM-FM representations: an asymptotically exact approach. IEEE Transactions on Audio, Speech and Language Processing 15(3): 823–837.

Herrera

Tary

van der Baan

, et al. (2015) Body wave separation in the time-frequency domain. IEEE Geoscience and Remote Sensing Letters 12(2): 364–368.

Iatsenko

McClintock

PVE

Stefanovska

(2015) Linear and synchrosqueezed time-frequency representations revisited: overview, standards of use, resolution, reconstruction, concentration, and algorithms. Digital Signal Processing 42: 1–26.

10.

Jamal

Nabeel

Ali

, et al. (2020) Multi-component instantaneous frequency estimation using signal decomposition and time-frequency filtering. Signal, Image and Video Processing 5: 1–8.

11.

Jones

Baraniuk

(1994) A simple scheme for adapting time-frequency representations. IEEE Transactions on Signal Processing 42(12): 3530–3535.

12.

Khan

Ali

(2020) A robust and efficient instantaneous frequency estimator of multi-component signals with intersecting time-frequency signatures. Signal Processing 177: 1–6.

13.

Liang

(2012) A generalized synchrosqueezing transform for enhancing signal time-frequency representation. Signal Processing 92(9): 2264–2274.

14.

Liutkus

Stoter

Kitamura

, et al. (2017) The 2016 signal separation evaluation campaign, latent variable analysis and signal separation. In: 13th international conference, LVA/ICA, Grenoble, France, 22–26 April 2019, pp. 323–332. IEEE: Boston, MA.

15.

Oberlin

Meignen

Perrier

(2014) The fourier-based synchrosqueezing transform. In: IEEE international conference on acoustics, speech and signal processing (ICASSP), Florence, Italy, 4–9 May 2014. IEEE.

16.

Oberlin

Meignen

Perrier

(2015) Second-order synchrosqueezing transform or invertible reassignment? towards ideal time-frequency representations. IEEE Transactions on Signal Processing 63(5): 1335–1344.

17.

Pham

Meignen

(2017) High-order synchrosqueezing transform for multicomponent signals analysis-with an application to gravitational-wave signal. IEEE Transactions on Signal Processing 65(12): 3168–3178.

18.

Sharma

Vignolo

Schlotthauer

, et al. (2017) Empirical mode decomposition for adaptive AM-FM analysis of speech: a review. Speech Communication 88: 39–64.

19.

Stankovic

Djurovic

Stankovic

, et al. (2014) Instantaneous frequency in time frequency analysis: enhanced concepts and performance of estimation algorithms. Digital Signal Processing 35: 1–13.

20.

Tary

Herrera

Han

, et al. (2014) Spectral estimation-what is new? what is next?. Reviews of Geophysics 52(4): 723–749.

21.

Wang

Chen

Cai

, et al. (2014) Matching demodulation transform and synchrosqueezing in time-frequency analysis. IEEE Transactions on Signal Processing 62(1): 69–84.

22.

(2011) Adaptive Analysis of Complex Data Sets. PhD Thesis, Princeton University, New Jersey, USA.