RRCNN: A novel signal decomposition approach based on recurrent residue
  convolutional neural network
The decomposition of non-stationary signals is an important and challenging
task in the field of signal time-frequency analysis. In the recent two decades,
many signal decomposition methods led by the empirical mode decomposition,
which was pioneered by Huang et al. in 1998, have been proposed by different
research groups. However, they still have some limitations. For example, they
are generally prone to boundary and mode mixing effects and are not very robust
to noise. Inspired by the successful applications of deep learning in fields
like image processing and natural language processing, and given the lack in
the literature of works in which deep learning techniques are used directly to
decompose non-stationary signals into simple oscillatory components, we use the
convolutional neural network, residual structure and nonlinear activation
function to compute in an innovative way the local average of the signal, and
study a new non-stationary signal decomposition method under the framework of
deep learning. We discuss the training process of the proposed model and study
the convergence analysis of the learning algorithm. In the experiments, we
evaluate the performance of the proposed model from two points of view: the
calculation of the local average and the signal decomposition. Furthermore, we
study the mode mixing, noise interference, and orthogonality properties of the
decomposed components produced by the proposed method. All results show that
the proposed model allows for better handling boundary effect, mode mixing
effect, robustness, and the orthogonality of the decomposed components than
existing methods.
Feng Zhou, Antonio Cicone, Haomin Zhou
RRCNN: A novel signal decomposition approach based on recurrent residue convolutional neural network

###### Abstract

The decomposition of non-stationary signals is an important and challenging task in the field of signal time-frequency analysis. In the recent two decades, many signal decomposition methods led by the empirical mode decomposition, which was pioneered by Huang et al. in 1998, have been proposed by different research groups. However, they still have some limitations. For example, they are generally prone to boundary and mode mixing effects and are not very robust to noise. Inspired by the successful applications of deep learning in fields like image processing and natural language processing, and given the lack in the literature of works in which deep learning techniques are used directly to decompose non-stationary signals into simple oscillatory components, we use the convolutional neural network, residual structure and nonlinear activation function to compute in an innovative way the local average of the signal, and study a new non-stationary signal decomposition method under the framework of deep learning. We discuss the training process of the proposed model and study the convergence analysis of the learning algorithm. In the experiments, we evaluate the performance of the proposed model from two points of view: the calculation of the local average and the signal decomposition. Furthermore, we study the mode mixing, noise interference, and orthogonality properties of the decomposed components produced by the proposed method. All results show that the proposed model allows for better handling boundary effect, mode mixing effect, robustness, and the orthogonality of the decomposed components than existing methods.

keywords: Empirical mode decomposition; Adaptive signal decomposition; Signal local average; Convolutional neural network; Residual network +
Footnote †: journal: Pattern Recognition

## 1 Introduction

With the development of technology, many every day-signals that exhibit nonlinearity and non-stationarity, such as human speech, radar systems, and seismic waves, can be accurately captured. It is well known that
decomposing and exploring features of this kind of signals is quite challenging due to their nonlinear and non-stationary characteristics.

In the past two decades, many studies have emerged for processing non-stationary signals. One of the most representative works is the empirical mode decomposition (EMD) algorithm along with the Hilbert spectrum analysis proposed by Huang et al. in 1998 [1]. Because EMD is fully data-driven, and can adaptively decompose a signal into several intrinsic mode functions (IMFs), it has already shown its usefulness in a wide range of applications, including semantic recognition [2], alcoholism identification [3], and stock trend prediction [4]. Despite its remarkable success, it still lacks mathematical foundations and is sensitive to noise and sampling. This sparked many efforts to improve the EMD. The improvements share the same feature: a signal is decomposed into several simpler components, and then a time-frequency analysis method is applied to each component separately. These signal decomposition methods can be mainly achieved in two ways: by iteration or by optimization.

Methods based on iteration include many techniques, such as moving average, partial differential equation (PDE) and filter. For instance, Smith presented a new iteration method, based on the local average, to decompose the non-stationary signals into a set of functions [5]. Delechelle et al. proposed a new approach that resolves one major problem in the EMD, that is, the mean envelope detection of a signal, in virtue of a parabolic PDE [6]. Hadji et al. used the differential calculus on envelopes, which makes them prove that iterations of the sifting process are well approximated by the resolution of PDE [7]. Hong et al. introduced a novel sifting method based on the concept of the local integral mean of a signal [8]. And Cicone et al. studied the method based on iterative filtering to compute the local average, which is utilized to replace the mean of the upper and lower envelopes in the sifting procedure of the EMD [9; 10]. Tu et al. proposed the iterative nonlinear chirp mode decomposition (INCMD) [11] under the framework of the variational nonlinear chirp mode decomposition.

On the other hand, there are methods based on optimization. Peng et al. designed an adaptive local linear operator-based optimization model to decompose a signal into several local narrow band signals [12]. Oberlin et al. proposed an optimization model in computing the mean envelope to replace the original one in EMD [13]. Inspired by the compressed sensing theory, Hou et al. studied a new adaptive data analysis method, which can be seen as a nonlinear version of compressed sensing and provides a mathematical foundation of the EMD method [14]. Flandrin et al. proposed a convex optimization procedure in order to replace the sifting process in the EMD, which follows the idea of texture-geometry decomposition with further specific EMD features such as quasi-orthogonality and extrema-based constraints [15; 16]. Dragomiretskiy et al. put forward the variational mode decomposition (VMD), whose goal is to decompose a signal into a discrete number of modes, that have specific sparsity properties while reproducing the input [17]. Rehman et al. generalized the VMD method to multivariate or multichannel data [18]. And Zhou et al. presented a new mathematical framework by finding the local average based on the local variational optimization model [19].

In addition, there are some methods that cannot be classified into the above two categories. For instance,
Daubechies et al. proposed the method, called synchrosqueezed wavelet transforms, by combining the wavelet analysis and reallocation method [20]. Gille presented the approach, called empirical wavelet transform (EWT), to build adaptive wavelets [21], whose main idea is to extract the different modes by designing an appropriate wavelet filter bank. Singh et al. studied the adaptive Fourier decomposition method (FDM) based on the Fourier theory, which decomposes any data into a small number of "Fourier intrinsic band functions" [22; 23]. And Wang et. extended the adaptive FDM to the multi-channel case [24].

According to the works described above, we find that whether the method is based on iterative or optimization, calculating the local average of a given signal is very critical. For example, in EMD [1], the mean of the upper and lower envelopes are used to measure the local average of the signal; the local variational optimization model is constructed to compute the local average in [19]; and in the iterative filtering method [9; 10], the low-pass filter is employed to find the local average. Although there exist many studies on the characterization of the local average, it is basically impossible to find a method suitable for all signals from a practical point of view. Discussing the local average customized according to the type of signal, it not only provides a new research perspective, but also is likely to become the trend in the near future in signal processing for non-stationary data.

In recent years, thanks to the remarkable results obtained in fields of research like image and natural language processing, the usage and application of deep learning methods have spread widely in an ample variety of research fields, like image processing [25] and natural language processing [26]. In signal processing, deep learning models have been used, so far, to achieve various goals, such as: noise removal [27; 28], forecasting [29; 30; 31], and detection [32; 33]. However, to the best of our knowledge, not a single method has been proposed so far in the literature, which allows to decompose a given non-stationary signal into simple oscillatory components, like the IMFs, which is solely based on deep learning techniques.

For this reason, in the current work we propose an innovative signal decomposition algorithm, named recurrent residual convolutional neural network (RRCNN), which is based on deep learning models. In the RRCNN method, in fact, we first characterize the local average of a signal under the framework of deep learning, and then use it to handle signal decomposition. Specifically, the 1-Dimensional (1-D) convolutional neural network is primarily designed to adaptively compute the local average of the input signal, which is similar to the moving average method and the filter operation in the iterative filtering method, except that the weights in 1-D convolutional neural network are not fixed, but are learned adaptively during the training phase according to the input signals. Moreover, both the residual and recurrent structures are employed to amend the computed local average, which is consistent with the role of the inner loop process in many existing iterative-based signal decomposition methods. After the local average model is derived, it is cascaded in series to realize the decomposition model of the signal, whose function is equivalent to the outer loop structure of the existing iterative-based signal decomposition methods.

Although the proposed method looks similar to those iterative-based decomposition methods that contain a two-loop structure, the use of the deep learning techniques makes this method have the following peculiarities:
1. Unlike the moving average method and the filter operation in the iterative filtering method, the convolutional filter weights that appear in the proposed RRCNN model, are not fixed in advance, but are learnt adaptively in the training phase according to the inputs.
2. Since the proposed RRCNN model is constructed under the framework of deep learning, it makes RRCNN more flexible and adaptive in finding the local average and achieving the decomposition for a given signal. In particular, the nonlinear activation function can be added after the convolutional operation to increase the expression ability. The loss function also can be customized according to the requirements that the ground truths usually have in the specific application.
3. Several artificial signals are constructed to verify the performance RRCNN in terms of local average characterization, noise interference, mode mixing, and orthogonality. Furthermore, we compare the RRCNN model with the state-of-the-art methods. In addition, we also use the solution of the Duffing and Lorenz equations, and the real data of the length of day (LOD) to evaluate the approximation ability of RRCNN to the existing models.
4. Generally speaking, the RRCNN model takes a certain amount of time in the training phase, which is a commonality of deep learning-based models. However, once the model training is completed, the computational efficiency in the prediction stage is relatively fast, especially it can use the parallelization mechanism to predict multiple signals at the same time, which is not available in most existing methods.
5. RRCNN has the limitations brought from the supervised model. For example: in the training phase, each input signal needs to know its label in advance.

The rest of the paper is organized as follows. We review the iterative filtering method and provide its algorithm in Section 2.1. And the concept of \(\beta\)-smooth function and its properties are given in Section 2.2, which are used for proving the convergence of the proposed model. In Section 3, the new local average method and the derived signal decomposition method, collectively called the RRCNN, are proposed. Moreover, the training process and convergence analysis of RRCNN are given in this section. In Section 4, we study a series of examples to evaluate the performance of RRCNN compared with the existing methods. Finally, we give the conclusion in Section 5.

## 2 IF and \(\beta\)-smooth function

### IF

The iterative filtering (IF) [9] is a recurrent algorithm that decomposes a nonlinear and non-stationary signal into a number of IMFs. The main idea of IF is the subtraction of local moving averages from the signal iteratively, where the local moving averages are calculated through convolutions with low-pass filters. Alg. 1 shows the
detailed steps, where the parameter \(l_{n}\), called the filter length, is important in the IF method, and is determined by the information contained in the signal itself; \(w_{n}(\cdot)\) represents the low-pass filter function.

```
Data: Given a signal \(x(t)\) Result:\(IMF\) while the number of extrema of \(x\geq 2\)do
1\(n=0\);
2\(x_{1}(t)=x(t)\);
3while the stopping criterion is not satisfied do
4 compute the filter length \(l_{n}\) and filter weight function \(w_{n}\) for \(x_{n}\);
5\(x_{n+1}(t)=x_{n}(t)-\int_{-l_{n}}^{l_{n}}x_{n}(t+y)w_{n}(y)dy\);
6\(n=n+1\);
7
8 end while
9\(IMF=IMF\cup\{x_{n}\}\);
10\(x(t)=x(t)-x_{n}(t)\);
11
12 end while
13\(IMF=IMF\cup\{x\}\).
```

**Algorithm 1**Iterative filtering (IF)

### \(\beta\)-smooth function and some of its properties

We first introduce the concepts of \(L\)-Lipschitz continuous and \(\beta\)-smooth for a function from [34].

**Definition 2.1**.: _A function \(f\) is said to be \(L\)-Lipschitz continuous if for all \(x,y\in\mathcal{X}\), \(\|f(x)-f(y)\|\leq L\|x-y\|\), where \(\mathcal{X}\) denotes the convex domain of \(f\), and \(L\) is called the Lipschitz constant._

**Definition 2.2**.: _A continuously differentiable function \(f\) is \(\beta\)-smooth if the gradient \(\nabla f\) is \(\beta\)-Lipschitz, that is if for all \(x,y\in\mathcal{X}\), \(\|\nabla f(x)-\nabla f(y)\|\leq\beta\|x-y\|\), where \(\mathcal{X}\) is the convex domain of \(f\)._

Then, for a unconstraint optimization problem, if its objective function is \(\beta\)-smooth, we can prove that the sequence generated by the gradient descent algorithm converges to a stationary point when the learning rate is small enough. The details can be found in Theorem 2.1.

**Theorem 2.1**.: _Let \(f\) be a \(\beta\)-smooth function and \(f^{*}=\min f(x)>-\infty\). Then the gradient descent algorithm with a constant learning rate \(\lambda<\frac{2}{\beta}\), i.e., \(x^{(k+1)}=x^{(k)}-\lambda\nabla f(x^{(k)}),\) converges to a stationary point, i.e., the set \(\{x:\nabla f(x)=\mathbf{0}\}\)._

Proof.: According to the gradient descent algorithm, i,e.,

\[x^{(k+1)}=x^{(k)}-\lambda\nabla f(x^{(k)}), \tag{1}\]

as \(f\) is \(\beta\)-smooth, we have

\[f(x^{(k+1)}) \overset{(a)}{\leq}f(x^{(k)})+\nabla f(x^{(k)})(x^{(k+1)}-x^{(k) })+\frac{\beta}{2}\|x^{(k+1)}-x^{(k)}\|^{2}\] \[\overset{(b)}{=}f(x^{(k)})-\lambda\|\nabla f(x^{(k)})\|^{2}+ \frac{\beta\lambda^{2}}{2}\|\nabla f(x^{(k)})\|^{2}\] \[=f(x^{(k)})-\lambda(1-\frac{\beta\lambda}{2})\|\nabla f(x^{(k)})\| ^{2},\]
where the inequality (a) follows from Lemma 3.4 in [34], and the equality (b) is obtained from Eqn. (1). Due to \(\lambda<2/\beta\), it becomes

\[\|\nabla f(x^{(k)})\|^{2}\leq\frac{f(x^{(k)})-f(x^{(k+1)})}{\lambda(1-\frac{ \beta\lambda}{2})}.\]

Next, we have

\[\sum_{k=0}^{K}\|\nabla f(x^{(k)})\|^{2}\leq\frac{1}{\lambda(1-\frac{\beta \lambda}{2})}\sum_{k=0}^{K}(f(x^{(k)})-f(x^{(k+1)}))=\frac{f(x^{(0)})-f(x^{(K+1 )})}{\lambda(1-\frac{\beta\lambda}{2})}\leq\frac{f(x^{(0)})-f(x^{*})}{\lambda( 1-\frac{\beta\lambda}{2})},\]

where \(x^{*}\) denotes the global optimization point. Taking the limit as \(K\rightarrow+\infty\), we have \(\sum_{k=0}^{+\infty}\|\nabla f(x^{(k)})\|^{2}\leq+\infty\). Hence, \(\lim_{k\rightarrow+\infty}\nabla f(x^{(k)})=0\) is obtained.

## 3 RRCNN inner loop block and RRCNN

### RRCNN inner loop block

The main operation in IF is the computation of moving average, which is essentially realized by the convolution operation, where the filter length depends on the given signal, and the filter weights are mainly given by some empirical functions selected artificially a priori. Therefore, it is very natural to convert the convolution operation into a 1-D convolutional neural network model, where both the filter length and the filter weights can be learnt adaptively according to the input signals given in advance. Furthermore, some ingenious mechanisms in deep learning, such as the nonlinear activation function, the residue learning [35], etc., can be adopted to make it more flexible. The structure we design to mimic the inner while loop of Alg. 1, is graphically depicted in Fig. 1. Since it mainly contains the recurrence mechanism, the convolutional layer and the subtraction operation, we call it the recurrent residual convolutional neural network (RRCNN) inner loop block.

```
Data:\(X\in\mathbb{R}^{N}\) Result: the local average of \(X\)
1 Initialize \(i=0\) and \(X^{(0)}=X\);
2while\(i<S\)do
3 The input \(X^{(i)}\) goes through the first 1-D convolutional layer, i.e., \(X_{C1}:=Conv1D(X^{(i)},W_{1}^{(i)},padding=True,activation=\tanh)\);
4 Transfer \(X_{C1}\) to the second convolutional layer, i.e., \(X_{C2}:=Conv1D\) (\(X_{C1},\tilde{W}_{2}^{(i)},padding=True\)), where \(\tilde{W}_{2}^{(i)}:=softmax(W_{2}^{(i)})\);
5\(X^{(i+1)}=X^{(i)}-X_{C2}\);
6\(i=i+1\);
7
8 end while
9\(\hat{Y}=X^{(S)}\) and \(X-\hat{Y}\) are the IMF and the local average respectively.
```

**Algorithm 2**RRCNN inner loop block
As shown in Fig. 1, the inner loop mainly consists of a judgment-loop structure and a residue operation, and the judgment-loop structure is formed of two convolutional layers and a residual operation. Suppose \(X\in\mathbb{R}^{N}\) (the vectors in this article are column vectors by default unless otherwise specified.) denote the input, the output of the RRCNN inner loop block, called \(\hat{Y}\in\mathbb{R}^{N}\), is computed as the following Alg. 2. And \(X-\hat{Y}\) is the local average of the input signal obtained by the RRCNN inner loop block.

Mathematically, the output of the RRCNN inner loop block can be expressed as: \(\hat{Y}=F(X,\mathbf{W})\), where \(\mathbf{W}\) denotes the undetermined weights in the RRCNN inner loop block, and the function \(F\) represents the structure of RRCNN inner loop block, which is composed of the operators including convolution, nonlinear activation and subtraction. The detailed process of \(F\) can be formulated as: \(F(X,\mathbf{W})=f(X^{(S-1)},\mathbf{W}^{(S-1)})\), where \(S\) represents the number of recursion in the RRCNN inner loop block, \(\mathbf{W}^{(S-1)}\) is the undetermined weights in the \((S-1)\)-th recursion, and the function \(f\) and \(X^{(S-1)}\) are defined as:

\[\begin{cases}f(X^{(i)},\mathbf{W}^{(i)})=\tanh(X^{(i)}*W_{1}^{(i)})*\tilde{W} _{2}^{(i)},\\ X^{(i+1)}=X^{(i)}-f(X^{(i)},\mathbf{W}^{(i)}),\end{cases} \tag{2}\]

where \(X^{(0)}=X\), \(\mathbf{W}^{(i)}\) is the undetermined weights in the \(i\)-th recursion that it includes the weights, denoted as \(W_{1}^{(i)}\in\mathbb{R}^{K_{1}}\) and \(W_{2}^{(i)}\in\mathbb{R}^{K_{2}}\), \(\tilde{W}_{2}^{(i)}:=softmax(W_{2}^{(i)})=\{\frac{exp(W_{2}^{(i)})}{\sum_{k} exp(W_{2}^{(i)})}\}_{l=1}^{K_{2}}\), \(*\) is the 1-D convolution operation, and \(i=0,1,\ldots,S-2\).

It is worth pointing out that the roles of the two 1-D convolutional layers in each judgment-loop are different. The role of the first convolutional layer, which is configured with a non-linear activation function (we select tanh in this work), is to enhance the nonlinear expression ability of the method. Whereas, the purpose of the second convolutional layer is to make the result more reasonable to describe the local average of the signal. Therefore, the non-negativity and normalization restrictions of its weights are added; and there is no nonlinear activation function configured with it. The use of padding in the two layers is to ensure that the length of the output is consistent with the input. We will discuss the training details of the RRCNN inner loop block in the following

Figure 1: Graphic illustration of the RRCNN inner loop block.
section.

### Rrcnn

After the RRCNN inner loop block for the identification of the local average of a given signal is constructed, we can cascade a finite number of RRCNN inner loop blocks together to derive the signal decomposition, which is called RRCNN, and is shown in Fig. 2. According to it, an input signal \(X\in\mathbb{R}^{N}\) (also denoted as \(X_{0}\)) can be decomposed into \(M\) IMFs \(\hat{\mathbf{Y}}:=\{\hat{Y}_{i}\}_{m=1}^{M}\) (each \(\hat{Y}_{i}\in\mathbb{R}^{N}\)) and a residue \(X_{M}\in\mathbb{R}^{N}\). The detailed steps of RRCNN are listed in Alg. 3.

The output of RRCNN can be formulated as:

\[\begin{cases}\hat{Y}_{m}=F(X_{m-1},\mathbf{W}_{m}),\\ X_{m}=X_{m-1}-\hat{Y}_{m},\end{cases} \tag{3}\]

where \(m=1,2,\ldots,M\), \(X_{0}=X\), \(F(X_{m-1},\mathbf{W}_{m})\) is the \(m\)-th RRCNN inner loop block whose purpose is to extract an IMF from \(X_{m-1}\), and \(\mathbf{W}_{m}\) denotes the undetermined weights of the \(m\)-th RRCNN inner loop block.

``` Data:\(X\in\mathbb{R}^{N}\), and the number of IMFs \(M\) Result: the IMFs and residue of \(X\) Initialize \(m=1\), \(X_{0}=X\); while\(m\leq M\)do

Compute the \(m\)-th IMF and the local average for the input \(X_{m-1}\), denoted as \(\hat{Y}_{m}\) and \(X_{m}\) respectively, according to the RRCNN inner loop block; \(m=m+1\);

end while\(\hat{\mathbf{Y}}=\{\hat{Y}_{m}\}_{m=1}^{M}\) and \(X_{M}\) are the resulting IMFs and residue of RRCNN.

Figure 2: Graphic illustration of the RRCNN.
All the generated IMFs are concatenated as the outputs of the RRCNN model. The errors between the outputs, i.e., \(\hat{\mathbf{Y}}\), and the labels that are composed of the true first \(M\) IMFs, denoted as \(\mathbf{Y}\in\mathbb{R}^{N\times M}\), are computed by the loss function. For example, the loss function can be expressed as:

\[L(\hat{\mathbf{Y}},\mathbf{Y})=\|\hat{\mathbf{Y}}-\mathbf{Y}\|_{F}^{2}, \tag{4}\]

where the errors are measured by mean square error (MSE), and \(\|\cdot\|_{F}\) denotes the Frobenius norm. In the RRCNN model equipped with the loss function as in Eqn. (4), the computational complexity of the forward process of RRCNN is mainly attributed to the computation of the convolutional layer, which is \(O(N\cdot K)\), where \(N\) and \(K\) denote the length of the signal and the size of the convolutional filter, respectively.

The loss can be customized according to the characteristic of the decomposition task. For example, if the 3rd IMF are smooth, the quadratic total variation term, expressed as

\[QTV(\hat{\mathbf{Y}}_{\Omega_{1}}):=\sum_{m\in\Omega_{1}}\sum_{t=1}^{N-1}( \hat{Y}_{(t+1),m}-\hat{Y}_{t,m})^{2},\]

can be added to the loss function, where \(\Omega_{1}\) represents the set of subscripts of those smooth components (here \(\Omega_{1}=\{3\}\)), i.e.,

\[L(\hat{\mathbf{Y}},\mathbf{Y})=\|\hat{\mathbf{Y}}-\mathbf{Y}\|_{F}^{2}+\eta QTV (\hat{\mathbf{Y}}_{\Omega_{1}}), \tag{5}\]

where \(\eta\geq 0\) is a penalty parameter, \(\hat{\mathbf{Y}}\) is the dependent variable of the function \(F(\cdot,\cdot)\), and its independent variables are \(X\) and \(\mathbf{W}\) respectively.

Moreover, if the 2nd and 3rd IMFs are orthogonal, an orthogonal constraint can be added to the loss function to ensure the orthogonality of the resulting components, i.e.,

\[\begin{split} L(\hat{\mathbf{Y}},\mathbf{Y})&=\sum_{ i\in\hat{\Omega}_{2}}\|\hat{Y}_{i}-Y_{i}\|_{2}^{2}+\sum_{(i,j)\in\Omega_{2}}\| \mathbf{W}^{o_{ij}}\hat{Y}_{i}-Y_{i}\|_{2}^{2}+\|\mathbf{W}^{o_{ij}}\hat{Y}_{j }-Y_{j}\|_{2}^{2},\\ &\text{s.t. }\mathbf{W}^{o_{ij}}\mathbf{W}^{o_{ij}\top}=\mathbf{I}, \end{split} \tag{6}\]

where \(\mathbf{W}^{o_{ij}}\in\mathbb{R}^{N\times N}\) stands for the orthogonal matrix to be determined by \(\min_{\{\mathbf{W},\mathbf{W}^{o_{ij}}\}}L(\hat{Y},Y)\), \(\Omega_{2}\) (here \(\Omega_{2}=\{(2,3)\}\)) denotes the subscript pairs of those orthogonal components, \(\hat{\Omega}_{2}\) (here \(\hat{\Omega}_{2}=\{2,3\}\)) represents the set consisting of all subscripts that appear in \(\Omega_{2}\), and \(\hat{\Omega}_{2}^{c}=\{1,2,\ldots,M\}-\hat{\Omega}_{2}\). In specific execution process, the orthogonal transformation \(\mathbf{W}^{o_{ij}}\) of \(\hat{Y}_{i}\) and \(\hat{Y}_{j}\) can be regarded as adding a fully connected layer after the outputs of \(\hat{Y}_{i}\) and \(\hat{Y}_{j}\). The two fully connected layers share weights, i.e., \(\mathbf{W}^{o_{ij}}\), and satisfy orthogonality, i.e., \(\mathbf{W}^{o_{ij}}\mathbf{W}^{o_{ij}\top}=\mathbf{I}\). In this case, the result of any IMF whose subscript meeting \(i\in\hat{\Omega}_{2}\) is updated from \(\hat{Y}_{i}\) to \(\mathbf{W}^{o_{ij}}\hat{Y}_{i}\), and the results of other components remain unchanged.

Compared with IF, RRCNN also contains two loops. The outer loop successively acquires \(M\) IMFs and one
residue. And the purpose of the inner loop, i.e., the RRCNN inner loop block, is mainly to compute the local average by several iterations, and finally get the IMF and local average of the current input signal. On the other hand, these two methods have important differences, which can be summed up as:

1. The number of outer loop of the IF method is data-driven, which might be different for different signals. While it is given in advance for the proposed RRCNN approach. Therefore, RRCNN is suitable for decomposing the non-stationary signals containing the same number of mono-components.
2. In the IF method, the filter length that is a key parameter, is determined only by the signal under decomposition. However, its filter weights lack of adaptability, and they are basically determined by the filter length. For the RRCNN approach, the filter length in each convolutional layer is adaptive, and can be selected by hyper-parameter optimization. Moreover, its filter weights are data-driven, which makes RRCNN flexible in the characterization of the local average.
3. In addition to more flexibly computing the local average of the signal, the proposed RRCNN method has all the advantages of deep learning, such as: the non-linear expression ability, the customized loss function according to specific decomposition tasks, etc., which are missing in the traditional signal decomposition methods.

Although the proposed RRCNN method is developed based on the IF, its ambition is not to replace the IF, or even to replace any existing signal decomposition method. The RRCNN is essentially a supervised deep learning model. This not only provides all the advantages that the existing signal decomposition methods are lacking, but also gives all the limitations of the supervision model itself. Such as, for any given signal, regardless of the decomposition performance, the existing signal decomposition methods can basically decompose it; however, the RRCNN model does not work without any training based on a set of input signals for which is known the ground truth. In addition, the training process of the RRCNN model can be regarded as the process of exploring the potential patterns of the input signals. If the patterns of the input signals are unique, their decomposition performance will also be greatly reduced.

### Training process and convergence analysis of RRCNN

In the training process, the gradient-based back-propagation method is used to learn the filter weights that appear in the convolutional layers by using the training data. For the hyper-parameters in the RRCNN, including the filter lengths and the numbers of neurons in the convolutional layers, Hyperopt1 is adopted to select their values by using the validation data, which is a Python library for serial and parallel optimization over awkward search spaces for hyper-parameters.

Footnote 1: Github website: [https://github.com/hyperopt/hyperopt](https://github.com/hyperopt/hyperopt)
Since when the parameter \(\eta\) in Eqn. (5) is equal to \(0\), it degenerates to the case of Eqn. (4), we only discuss the training processes and convergences of the models when their loss functions are given in Eqns. (5) and (6), respectively.

We first discuss the situation of the RRCNN model equipped with the loss function (5). For convenience, we consider that the output of the model only produces one IMF, and the number of recursion is also limited to \(1\), that is, \(M=1\) and \(S_{1}=1\) in Fig. 2. Suppose that \((X^{j},Y^{j})_{j=1}^{J}\subset\mathbb{R}^{N}\times\mathbb{R}^{N}\) is a given set of training samples, where \(J\) denotes the number of samples. The processe and the loss function of RRCNN are as follows:

\[\left\{\begin{array}{l}X_{C1}^{j}=\sigma(X^{j}*W_{1}),\\ \\ X_{C2}^{j}=X_{C1}^{j}*\tilde{W}_{2},\\ \\ \hat{Y}^{j}=X^{j}-X_{C2}^{j},\end{array}\right. \tag{7}\]

\(L(\hat{Y},Y)=\sum_{j=1}^{J}\|\hat{Y}^{j}-Y^{j}\|_{2}^{2}+\eta QTV(\hat{Y}^{j})\), where \(\tilde{W}_{2}=softmax(W_{2})\), \(\sigma=tanh\) in our model, \(\eta\) is a non-negative parameter, \(W_{1}\in\mathbb{R}^{K_{1}}\) and \(W_{2}\in\mathbb{R}^{K_{2}}\) are the undetermined convolution filters.

According to the gradient descent method and the chain derivation rule, \(W_{1}\) and \(W_{2}\) are learned following the back propagation method, that is, \(W_{h}^{(n+1)}=W_{h}^{(n)}-\lambda\nabla_{W_{h}}L,(h=1,2)\) where \(\lambda\) denotes the learning rate, \(\nabla_{W_{h}}L=(\frac{\partial}{\partial W_{h1}},\ldots,\frac{\partial}{ \partial W_{hK_{1}}})^{\top}L\), \(\frac{\partial L}{\partial W_{2i}}=-\sum_{j=1}^{J}\sum_{t,l=1}^{N}\frac{ \partial L}{\partial Y_{t}^{j}}\frac{\partial X_{C2i}^{j}}{\partial W_{2i}} \frac{\partial\tilde{W}_{2i}}{\partial W_{2i}}\), \(\frac{\partial L}{\partial W_{1i}}=-\sum_{j=1}^{J}\sum_{t,l=1}^{N}\frac{ \partial L}{\partial Y_{t}^{j}}\frac{\partial X_{C2i}^{j}}{\partial X_{C1i}^{ j}}\sigma^{\prime}(R_{l}^{j})\frac{\partial R_{l}^{j}}{\partial W_{1i}}\), \(\frac{\partial L}{\partial Y_{t}^{j}}=2(\tilde{Y}_{t}^{j}-Y_{t}^{j})\), \(\frac{\partial X_{C2i}^{j}}{\partial W_{2i}}=X_{C1(t-[K_{1}/2]+l)}^{j}\), \(R_{l}^{j}=(X^{j}*W_{1})_{l}\), \(\frac{\partial R_{l}^{j}}{\partial W_{1i}}=X^{j}{}_{l-[K_{2}/2]+i}\),

\[\frac{\partial\tilde{W}_{2l}}{\partial W_{2i}}=\left\{\begin{array}{l}\frac{ exp(W_{2i}+W_{2i})}{(\sum_{k}exp(W_{2k}))^{2}},\ \mbox{if}\ l\neq i;\\ \frac{exp(W_{2i})\sum_{k\neq i}exp(W_{2k})}{(\sum_{k}exp(W_{2k}))^{2}},\ \mbox{otherwise},\end{array}\right.\mbox{ and } \frac{\partial X_{C2i}^{j}}{\partial X_{C1i}^{j}}=\left\{\begin{array}{l} \tilde{W}_{2(l-t+[K_{1}/2])},\ \mbox{if}\ 1\leq l-t+[K_{1}/2]\leq K_{1};\\ \\ 0,\ \mbox{otherwise}.\end{array}\right.\]

Next, we discuss the convergence of the training process of the RRCNN model expressed in Eqn. (7).

**Theorem 3.1**.: _For the RRCNN model defined in Eqn. (7), the sequences \(\{W_{1}^{(n)}\}\) and \(\{W_{2}^{(n)}\}\) are generated by the gradient descent algorithm. Then, there exists a positive constant \(L\) independent of the input data, so that when the learning rate \(\lambda\) is less than \(L\), \(\{W_{1}^{(n)}\}\) and \(\{W_{2}^{(n)}\}\) converge to their corresponding stationary points, i.e., the sets \(\{W_{1}:\nabla_{W_{1}}L=\mathbf{0}\}\) and \(\{W_{2}:\nabla_{W_{2}}L=\mathbf{0}\}\), respectively._

Proof. From the Lagrange's mean value theorem, it is obvious to find that the composition function composed by two functions, the \(\beta_{1}\)- and \(\beta_{2}\)-smooth respectively, is still \(\beta\)-smooth (\(\beta\leq\beta_{1}\beta_{2}\)) under the condition that it can be composited.

According to Eqn. (3), the RRCNN model can be seen as a composition function composed of a series of functions. Then, combining with Theorem 2.1, we have that if the function of the \(i\)-layer in the RRCNN model is \(\beta_{i}\)-smooth, it can be proved that the weights sequences obtained by the gradient descent method, converge under the condition that the learning rate satisfies \(\lambda\leq 2/\Pi\beta_{i}\).
For the function \(\tilde{W}_{2}=softmax(W_{2})\) that is included in the RRCNN model, the second partial derivative of each \(\tilde{W}_{2l}\) (\(l=1,\ldots K_{2}\)) with respect to its independent variable vector \(W_{2}\) exists and satisfies:

\[\left|\frac{\partial^{2}\tilde{W}_{2l}}{\partial W_{2i}\partial W_{2j}}\right|= \frac{exp(W_{2i})}{(\sum_{k}exp(W_{2k}))^{3}}\begin{cases}\sum_{k\neq i}exp(W_{2 k})|\sum_{k}exp(W_{2k})-2exp(W_{2i})|,&\text{if }l=i=j,\\ exp(W_{2l})|2exp(W_{2i})-\sum_{k}exp(W_{2k})|,&\text{if }l\neq i=j,\,\leq 2.\\ exp(W_{2i}+W_{2l}+W_{2j}),&\text{if }l\neq i\neq j,\end{cases}\]

Hence, it states that each \(\tilde{W}_{2l}\) is \(\beta\)-smooth (and \(\beta\leq 2\)) according to the Lagrange's mean value theorem.

For the rest functions involved in RRCNN include quadratic function, tanh, and 1-D convolution operation, these functions can be easily proved to be \(\beta\)-smooth (\(\beta\) here is a general representation, and the \(\beta\) value of each function may not be the same.) by judging that their second derivative functions exist and bounded. Therefore, the conclusion is proved.

For the case of the model with an orthogonal constraint in the loss function \(L\) in Eqn. (6), the orthogonal constraint is a new obstacle compared to the previous model. However, in the field of optimization, the study of optimization problems with orthogonal constraints has become very common. And the gradient-based projection method can be used to find \(\mathbf{W}^{o}\) with convergence guarantees [36, 37]. Furthermore, under the idea of back propagation used in updating the weights of the neural networks, the solutions of problem (6) can be obtained according to Alg. 4.

```
1\(i=0\), given the learning rates \(\lambda_{1},\lambda_{2}\), and initialize the \(\mathbf{W}^{o},\mathbf{W}\);
2while\(i<Max\_Iter\)do
3\(\mathbf{W}^{o}\leftarrow\mathbf{W}^{o}-\lambda_{1}\nabla_{\mathbf{W}^{o}}L\);
4\(\mathbf{W}^{o}\leftarrow\mathcal{P}_{\mathcal{S}_{p,q}}(\mathbf{W}^{o})\), where \(\mathcal{P}_{\mathcal{S}_{p,q}}(\mathbf{W}^{o})\) denotes the projection of \(\mathbf{W}^{o}\in\mathbb{R}^{p\times q}\) to the Stiefel manifold \(\mathcal{S}_{p,q}\), \(\mathcal{P}_{\mathcal{S}_{p,q}}(\mathbf{W}^{o})=\tilde{\mathbf{U}}\mathbf{ \nabla}^{\top}\), and \(\tilde{\mathbf{U}}\tilde{\mathbf{x}}\tilde{\mathbf{V}}^{\top}\) is the reduced singular value decomposition of \(\mathbf{W}^{o}\);
5\(\mathbf{W}\leftarrow\mathbf{W}-\lambda_{2}\nabla_{\mathbf{W}}L\) according to the back propagation method;
6\(i\gets i+1\);
7 end while
```

**Algorithm 4**RRCNN with orthogonal constraint

Since the convergence analysis of \(\mathbf{W}\) calculated based on the gradient descent method is consistent with Theorem 3.1, and that of \(\mathbf{W}^{o}\) calculated based on the gradient projection method has been discussed in the literature [36, 37], so the convergence of Alg. 4 can be obtained intuitively.

**Remark 1**.: _Similar to the case with loss function in Eqn. (5), we assume that in the RRCNN model, there are only two RRCNN inner loop blocks, i,e., \(M=2\) in Fig. 2. Furthermore, the IMFs obtained by the two blocks satisfy orthogonality, i.e., \(\Omega_{2}=\{(1,2)\}\), \(\hat{\Omega}_{2}=\{1,2\}\) and \(\hat{\Omega}_{2}^{c}=\emptyset\) in Eqn. (6). In this case, the RRCNN model can be reduced to a simpler formula, which looks close to the expressions of some orthogonal constraint algorithms, which reads:_

\[\min_{\mathbf{W},\mathbf{W}^{o}}\|\mathbf{W}^{o}\hat{\mathbf{Y}}-\mathbf{Y}\| _{F}^{2},s.t.\ \mathbf{W}^{o}\mathbf{W}^{o\,\top}=\mathbf{I}, \tag{8}\]

_where \(\hat{\mathbf{Y}}\) depends on \(\mathbf{W}\)._
## 4 Experiments

To evaluate the performance of the proposed RRCNN inner loop block and RRCNN models, we test them against seven aspects, which are: 1) Can RRCNN inner loop block be used to find the local average of the non-stationary signal? 2) Is RRCNN inner loop block still effective on noisy signal? 3) Can RRCNN be used to decompose the non-stationary signals? 4) Can RRCNN be effective on noisy signals? 5) Can RRCNN be effective on the signals composed of orthogonal mono-components? 6) Can RRCNN be effective on solutions of differential equations? 7) Is RRCNN capable of processing real signals? Furthermore, we discuss the computational time of RRCNN, and compare it with other methods. In addition, limited by the length of the paper, we submit some experiments as supplementary materials.

In the experiments, we divide each of the constructed input data into the training and validation datasets with ratio \(7:3\). Moreover, our proposed signal average method, i.e., the RRCNN inner loop block, is compared with the existing signal average methods based on cubic spline envelopes (CSE) [1], optimization model (OP), a segment power-function based envelopes (SPFE) [38] and iterative filtering (IF) [9], respectively. For simplicity, we denote the averages obtained from the CSE, OP, SPFE, and IF as CSA, OPA, SPFA, and IF respectively. And the RRCNN2 method will be compared with the state-of-the-art methods, including EMD, IF3, VMD4[17], continuous wavelet transform based synchrosqueezing (SYNSQ_CWT) [20], short time Fourier transform based synchrosqueezing (SYNSQ_STFT5) [20], INCMD6[11], EWT7[21], FDM8[22], and its variant called DCT_GAS_FDM9[23]. In addition, in the experiments for verifying the robustness of noise interference, the proposed RRCNN method is compared with the ensemble EMD (called EEMD10) model [39]; and it is compared with the M-LFBF11[15, 16] model to verify the orthogonality of the decomposed components.

Footnote 2: Code of RRCNN is available at [https://github.com/zhoudafa08/RRCNN](https://github.com/zhoudafa08/RRCNN)

Footnote 3: Code of IF: [http://people.disim.univaq.it/](http://people.disim.univaq.it/)\(\sim\)antonio.cione/Software.html

Footnote 4: Code of VMD: [https://www.mathworks.com/help/wavelet/ref/vmd.html](https://www.mathworks.com/help/wavelet/ref/vmd.html)

Footnote 5: Codes of SYNSQ_CWT and SYNSQ_STFT: [https://github.com/ebrevd](https://github.com/ebrevd) o/synchrosqueezing

Footnote 6: Code of INCMD: [https://github.com/sheadan/IterativeNCMD](https://github.com/sheadan/IterativeNCMD)

Footnote 7: Code of EVT: [https://www2.mathworks.cs/help/wavelet/ug/empirical-wave](https://www2.mathworks.cs/help/wavelet/ug/empirical-wave) let-transform.html

Footnote 8: Code of FDM: [https://www.researchgate.net/publication/274570245_Matl](https://www.researchgate.net/publication/274570245_Matl) ab_Code_Of_The_Fourier_Decomposition_Method_FDM

Footnote 9: Code of DCT_GAS_FDM: [https://www.researchgate.net/publication/32629](https://www.researchgate.net/publication/32629) 4577_MATLABCodeOfFDM_DCT_DFT_FIR_FSASJuly2018

Footnote 10: Code of EEMD: [http://perso.ens-lyon.fr/patrick.flandrin/emd.html](http://perso.ens-lyon.fr/patrick.flandrin/emd.html)

Footnote 11: Code of M-LFBF: [http://perso.ens-lyon.fr/nelly.pustelnik/](http://perso.ens-lyon.fr/nelly.pustelnik/)

The results are measured by the metrics listed in Tab. 1, where Mean Absolute Error (MAE) and Root Mean Squared Error (RMSE) are used to measure the errors between the predicted results \(\hat{Y}\in\mathbb{R}^{N}\) and the ground truth \(Y\in\mathbb{R}^{N}\), and \(\rho(c_{1},c_{2})\) is used to evaluate the orthogonality between the resulting components \(c_{1}\in\mathbb{R}^{N}\) and \(c_{2}\in\mathbb{R}^{N}\).

\begin{table}
\begin{tabular}{c c c} \hline \hline Metric & MAE & RMSE & \(\rho(c_{1},c_{2})\) \\ \hline Expression & \(\frac{1}{2}\sum_{t=1}^{N}|\hat{Y}_{t}-Y_{t}|\) & \(\sqrt{\frac{1}{2}\sum_{t=1}^{N}(\hat{Y}_{t}-Y_{t})^{2}}\) & \(\frac{|\langle c_{1},c_{2}\rangle|}{|C_{1}|c_{2}|c_{2}|c_{2}}\) \\ \hline \hline \end{tabular}
\end{table}
Table 1: Evaluation indices
**4.1. Can RRCNN inner loop block be used to find the local average of the non-stationary signal?**

We first evaluate the performance of the proposed RRCNN inner loop block in solving the local average of the non-stationary signal. Several signals composed of the linear function and the mono-component function, or just the mono-component function, are constructed as the inputs, where the mono-component function can generally be expressed as \(a(t)\cos\theta(t)\), which meets \(a(t),\theta^{\prime}(t)>0\)\(\forall t\), and the changes in time of \(a(t)\) and \(\theta^{\prime}(t)\) are much slower than that of \(\theta(t)\). Ideally, the average of the mono-component signal is zero. The input signals and the corresponding labels we construct in this part are listed in Tab. 2. It should be pointed out that the label here represents the first IMF of the corresponding input signal, not the local average. The local average can be computed by subtracting the label from the input signal.

After the RRCNN inner loop block is trained with the inputs in Tab. 2, we select three mono-component signals with different instantaneous frequencies and instantaneous amplitudes, discussed in Examples 1-2, which are not included in the inputs, to test the performance of the RRCNN inner loop block, respectively.

**Example 1**: \(x(t)=(3+2\cos(2t))\cos(2t^{2})\), \(t\in[0,3]\).

**Example 2**: \(x(t)=(2t+\cos(2t^{2}))\cos(12t+t^{2}+2\cos(t))\), \(t\in[0,3]\).

**Example 3**: \(x(t)=(3+2\cos(3t))\cos(5t^{2})\), \(t\in[0,3]\).

[MISSING_PAGE_POST]

**Example 35**: \(x(t)=(2t+\cos(2t))\cos(12t+t^{2}+2\cos(t))\
The moving averages of the signals in Examples 1-3 obtained from different methods are shown in Fig. 3 (a)-(c), respectively, and the errors between the obtained moving averages and the true average are listed in Tabs. 3-5, respectively. According to the results, we can observe the following phenomena:

1. The existing methods are prone to boundary effects, which can be seen from the left boundaries of Fig. 3 (a)-(c). However, the RRCNN inner loop block method can avoid this problem to a certain extent.
2. When the signal is in a situation where the amplitude changes quickly and the frequency changes slowly, the RRCNN inner loop block performs best among all models according to the left parts of Fig. 3 (a)-(c). When the amplitude change is reduced and the frequency change is accelerated, its performance may even be inferior to other models, which can be seen from the right half of Fig. 3 (c).
3. Even though the RRCNN inner loop block has some dazzling and bleak performance compared with other local average methods, it can be seen from Tabs. 3-5 that the MAE and RMSE of RRCNN are significantly reduced compared to other models.

The reason for the phenomenon above can be attributed to:

1. The averages obtained by the comparison methods are basically determined by the local information of the signal, which makes the results reasonable when the information is sufficient (e.g., the part of the amplitude change is reduced and the frequency change is accelerated); and the results differ greatly when the information is insufficient (e.g., the part of the amplitude changes quickly and the frequency changes slowly).
2. The filter weights of each convolutional layer in the RRCNN are shared, they are determined by all the information contained in the whole signal. Therefore, the average obtained by the RRCNN is relatively

\begin{table}
\begin{tabular}{c c c c c} \hline \hline Metric & CSA & OPA & SPFA & IF & **RRCNN** \\ \hline MAE & 0.2067 & 0.2448 & 0.3646 & 0.8624 & **0.1418** \\ RMSE & 0.4185 & 0.5427 & 0.6388 & 1.5833 & **0.1712** \\ \hline \hline \end{tabular}
\end{table}
Table 5: Metrics of the moving averages of the input in Example 3.

Figure 3: Moving averages by different methods of Examples 1-3.
stable, and it is not easy to cause interference due to the large difference in the changing of the signal amplitude and frequency.

### **Rs RRCNN inner loop block still effective on noisy signal?**

In this part, we consider the robustness of the RRCNN inner loop block model to noise based on the constructed inputs and labels in Section 4.1. Specifically, each input signal is perturbed with additive Gaussian noise with the signal-to-noise ratio (SNR) set to \(15dB\), and the corresponding label remains unchanged, as detailed in Tab. 6. Similar to the section above, we select a noisy signal, which is essentially the signal in Examples 2-3 with additive Gaussian noise, and is described in Examples 4-5, to test the performance.

**Example 4**: \(x(t)=(2t+\cos(2t^{2}))\cos(20t+t^{2}+2\cos(t))+\varepsilon(t)\), where \(t\in[0,3]\) and \(SNR=15dB\).

**Example 5**: \(x(t)=(3+2\cos(3t))\cos(5t^{2})+\varepsilon(t)\), where \(t\in[0,3]\) and \(SNR=15dB\).

The results of Examples 4-5 are shown in Fig. 4 (a)-(b) and Tabs. 7-8. From the results, we can find that the RRCNN inner loop block performs the the most robust among all models for the signal with additive Gaussian noise. Specifically, in Example 4, the RRCNN inner loop block reduces the MAE and RMSE from \(0.1936,0.2420\), for the second the best method (i.e., SPFA), to \(0.1190,0.1500\), respectively; and in Example 5, the RRCNN inner loop block reduces the MAE and RMSE from \(0.2004,0.2520\), for the second the best method (i.e., SPFA), to \(0.1290,0.1623\), respectively.

\begin{table}
\begin{tabular}{c c c c c} \hline \hline \(x_{1}(t)\) & \(x_{2}(t)\) & Inputs & Labels & Notes \\ \hline \multirow{2}{*}{\(0.1kt\)} & \(\cos(3lt)\) & & & \(k=2,3,\ldots,9\) \\  & \(\cos(3lt+t+\cos(t))\) & & & \(k=2,4,6,8\) \\  & \(\cos(3lt+t+\cos(t))\) & & & \(l=2,4,6,8\) \\ \hline \multirow{2}{*}{\(0.1k\)} & \(\sin(3lt)\) & & & \(k=1,2,\ldots,10\) \\  & \(\sin(3klt)\) & \(x_{1}(t)+x_{2}(t)+\varepsilon(t)\) & \(x_{2}(t)\) & \\  & \(\sin(3klt+t^{2}+\cos(t))\) & & & \(l=2,4,\ldots,28\) \\ \hline \multirow{2}{*}{\(3+2\cos(0.5kt)\)} & \(\cos(0.5klt^{2})\) & & & \(k=2,3,\ldots,6\) \\  & \(\cos(0.5li^{2}+l\cos(t))\) & & & \(l=4,5,\ldots,9\) \\ \hline \hline \end{tabular}
\end{table}
Table 6: Inputs disturbed by the Gaussian noise with the SNR set to \(15dB\), and the labels used in Section 4.2, where \(t\in[0,3]\).

\begin{table}
\begin{tabular}{c c c c c c} \hline \hline Metric & CSA & OPA & SPFA & IF & **RRCNN** \\ \hline MAE & 0.2023 & 0.1998 & 0.1936 & 0.2183 & **0.1190** \\ RMSE & 0.2529 & 0.2501 & 0.2420 & 0.2762 & **0.1500** \\ \hline \hline \end{tabular}
\end{table}
Table 7: Metrics of the moving averages of the input in Example 4 obtained from different methods.

\begin{table}
\begin{tabular}{c c c c c} \hline \hline Metric & CSA & OPA & SPFA & IF & **RRCNN** \\ \hline MAE & 0.2141 & 0.2126 & 0.2004 & 0.2339 & **0.1290** \\ RMSE & 0.2676 & 0.2658 & 0.2520 & 0.2967 & **0.1623** \\ \hline \hline \end{tabular}
\end{table}
Table 8: Metrics of the moving averages of the input in Example 5.
### Can RRCNN be used to decompose the non-stationary signals?

After evaluating the proposed RRCNN inner loop block model in calculating the local average, the next task is to assess the decomposition performance of RRCNN for non-stationary signals. To demonstrate the problem, we only consider decomposing the signals consisting of two components. The input signals in the part can be divided into two categories: one is composed of a mono-component signal and a zero signal, and the other is of two mono-components with close frequencies. The former is to train RRCNN describe the zero local average of the mono-component signal; and the latter is to enable RRCNN to decompose signals with close frequencies, which is the main factor causing the mode mixing effect. The inputs and labels are constructed and shown in Tab. 9. To challenge the proposed model, we hereby choose the signal composed of two cosine signals with frequencies of \(2.5\)Hz and \(3.4\)Hz respectively, described in Example 4, that is more prone to introduce the mode mixing effect in the existing methods.

**Example 6**: \(x(t)=\cos(5\pi t)+\cos(6.8\pi t)\), \(t\in[0,6]\).

In Example 6, the components predicted by the trained RRCNN model are compared with those obtained from the state-of-the-art methods in signal decomposition. The metrics of the errors between the obtained components

\begin{table}
\begin{tabular}{c c c c} \hline \hline \(x_{1}(t)\) & \(x_{2}(t)\) & Inputs & Labels & Notes \\ \hline \(\cos(k\pi t)\) & \(\cos((k+1.5)\pi t+t^{2}+\cos(t))\) & \(x_{1}(t)+x_{2}(t)\) & \([x_{2}(t),x_{1}(t)]\) & \(k=5,6,\ldots,14\) \\ \(0\) & \(\cos((k+1.5)\pi t+t^{2}+\cos(t))\) & & \\ \hline \(\cos(k\pi t)\) & \(\cos(kl\pi t+t^{2}+\cos(t))\) & & \\ \(0\) & \(\cos(kl\pi t)\) & \(x_{1}(t)+x_{2}(t)\) & \([x_{2}(t),x_{1}(t)]\) & \\ \(0\) & \(\cos(kl\pi t+t^{2}+\cos(t))\) & & \(l=2,3,\ldots,19\) \\ \hline \hline \end{tabular}
\end{table}
Table 9: Inputs and labels used in Section 4.3, where \(t\in[0,6]\).

Figure 4: Moving averages by different methods of Examples 4-5.
and the labels, measured by MAE and RMSE, are shown in Tab. 10. In addition, to compare RRCNN and the existing methods more intuitively, we select the top three methods with comprehensive performance from Tab. 10, i.e., RRCNN, EMD, and INCMD, and plot their obtained components and the corresponding time-frequency-energy (TFE) representations in Fig. 5 (a), (b), respectively. It should be noted that the identification of an optimal TFE representation is a research topic on its own, and it is out of the scope of this work. Here, we set as TFE representation the Fourier quadrature transforms that was proposed in [23].

According to the results, we have the following conclusions:

1. The mode mixing problem is indeed a big challenge for some of the existing methods. For example, the maximum value of \(x(t)\) is 2, but the MAEs of the two components obtained by the SYNSQ_STFT method are as high as 0.6620 and 0.6778, respectively, which basically does not separate the cosine signals with frequencies of 2.5Hz and 3.4Hz from \(x(t)\).
2. Many methods achieve satisfactory decomposition for \(x(t)\). For example, it can be seen from the left plots in the 2nd and 3rd rows of Fig. 5 (a) that the components obtained by the EMD, INCMD, and RRCNN methods have relatively consistent oscillation modes with the ground truths. This viewpoint can also be drawn from Fig. 5 (b), although there are some obvious fluctuations, the TFE representations of the two components, obtained by EMD, INCMD, and RRCNN methods, are basically separated just like the those of the real components.
3. Nonetheless, a closer look at the right plots in the 2nd and 3rd rows of Fig. 5 (a) reveals the subtle advantage of the RRCNN model at the boundaries. Due to the incompleteness of the waveform at the boundary, many existing methods are deeply affected by it, as are EMD and INCMD. However, the weights of the convolution filters in the RRCNN model depend on the entire waveform of the whole training samples, which reduce consistently the boundary effects.

### Can RRCNN be effective on noisy signals?

Similar to the RRCNN inner loop block, we verify the robustness of RRCNN against additive Gaussian noise in this part. The constructed inputs and labels are listed in Tab. 11, where the inputs are generated by introducing additive Gaussian noise with the SNR set to \(25dB\) to the signals in Tab. 9. After the RRCNN model is trained,

\begin{table}
\begin{tabular}{c c c c c c c c c} \hline \multirow{2}{*}{Method} & \multicolumn{2}{c}{\(c_{1}\)} & \multicolumn{2}{c}{\(c_{2}\)} & \multirow{2}{*}{Method} & \multicolumn{2}{c}{\(c_{1}\)} & \multicolumn{2}{c}{\(c_{2}\)} \\ \cline{2-2} \cline{5-8}  & MAE & RMSE & MAE & RMSE & MAE & RMSE & MAE & RMSE \\ \hline EMD[1] & 0.1049 & 0.2110 & 0.1132 & 0.2058 & INCMD[11] & 0.1144 & 0.1667 & 0.1160 & 0.1509 \\ VMD[17] & 0.2193 & 0.2479 & 0.2193 & 0.2479 & SYNSQ\_CWT[20] & 0.2227 & 0.2496 & 0.2589 & 0.2965 \\ EWT[21] & 0.2145 & 0.2605 & 0.2145 & 0.2605 & SYNSQ\_STFT[20] & 0.6620 & 0.7402 & 0.6778 & 0.7608 \\ FDM[22] & 0.2438 & 0.2947 & 0.2429 & 0.2893 & DCT\_GAS\_FDM[23] & 0.1608 & 0.2005 & 0.1605 & 0.2003 \\ IF[9] & 0.2253 & 0.2756 & 0.2060 & 0.2489 & **RRCNN** & **0.1013** & **0.1242** & **0.0331** & **0.0422** \\ \hline \end{tabular}
\end{table}
Table 10: Metrics of the errors between the obtained components by different methods and the ground truth of Example 6.
we choose the signal consisting of two mono-components and additive Gaussian noise with a SNR of \(15dB\) as the test data, which is given in Example 7. Since the smaller SNR value, the greater the noise, the noise of \(x(t)\) is larger than that in the training data.

**Example 7**: \(x(t)=\cos(5\pi t)+\sin(8\pi t+2t^{2}+\cos(t))+\varepsilon(t)\), \(t\in[0,6]\), SNR=\(15dB\).

Figure 5: Results of Examples 6-8.
The errors between the ground truths and the components obtained by different methods, measured by MAE and RMSE, are reported in Tab. 12. Furthermore, the components, errors, and TFE representations of the three best performing methods, i.e., FDM, DCT_GAS_FDM, and RRCNN, are shown in Fig. 5 (d), (e), respectively.

According to the results, we can find that RRCNN works for the signals with additive Gaussian noise, although there is no overwhelming advantage, especially over FDM. Specifically, from the left plots in the 2nd-3rd rows of Fig. 5 (d), RRCNN basically separates the two mono-components, and the resulting components are consistent with the ground truths in the oscillation mode. Moreover, as shown in the right plots in the 2nd-3rd rows of Fig. 5 (d), the errors of RRCNN are relatively evenly dispersed in the entire time period, while those of the FDM and DCT_GAS_FDM methods are both small in the middle and large at the boundaries, which is consistent with the observations in Section 4.3.

At last, since RRCNN is a method designed with the help of CNN in the time domain, it can obtain an effect comparable to the existing methods in the time domain. Due to the lack of a prior information of RRCNN in the frequency domain, the effect of this method might be slightly reduced from the time-frequency domain. According to the results in Fig. 5 (e), the TFE distributions of the two mono-components obtained by the FDM, DCT_GAS_FDM, and RRCNN, are obviously spaced apart, but the TFE distribution of the component \(c_{1}\) obtained by RRCNN has a much more severe jitter than the true component \(\sin(8\pi t+2t^{2}+\cos(t))\) in the interval \(t\in[2,4]\). However, from the TFE representations, we can also see that the RRCNN method is able to reduce boundary effects compared to other methods.

### Can RRCNN be effective on the signals composed of orthogonal mono-components?

In this part, we test the proposed RRCNN model on the signal composed of the orthogonal components. As discussed in Section 3.2, the RRCNN model in this part should have been equipped with the loss function with

\begin{table}
\begin{tabular}{c c c c c c c} \hline \hline \multirow{2}{*}{Method} & \multicolumn{2}{c}{\(c_{1}\)} & \multicolumn{2}{c}{\(c_{2}\)} & \multirow{2}{*}{Method} & \multicolumn{2}{c}{\(c_{1}\)} & \multicolumn{2}{c}{\(c_{2}\)} \\ \cline{2-2} \cline{5-8}  & MAE & RMSE & MAE & RMSE & MAE & RMSE \\ \hline EMD[1] & 0.1693 & 0.2961 & 0.2029 & 0.3485 & INCMD[11] & 0.2436 & 0.3870 & 0.3343 & 0.4720 \\ VMD[17] & 0.1440 & 0.2371 & 0.1363 & 0.2286 & SYNSQ_CWT[20] & 0.5817 & 0.6477 & 0.0915 & 0.1522 \\ EWT[21] & 0.1840 & 0.2945 & 0.1754 & 0.2884 & SYNSQ_STFT[20] & 0.4353 & 0.4992 & 0.2009 & 0.2463 \\ FDM[22] & **0.1039** & **0.1765** & 0.0857 & 0.1327 & DCT_GAS_FDM[23] & 0.1288 & 0.1801 & 0.0826 & 0.1252 \\ IF[9] & 0.1559 & 0.2318 & 0.1601 & 0.2256 & EEMD[39] & 0.1870 & 0.2308 & 0.1243 & 0.1925 \\
**RRCNN** & 0.1556 & 0.1977 & **0.0805** & **0.1014** & & & & \\ \hline \hline \end{tabular}
\end{table}
Table 12: Metrics of the errors between the results obtained by different methods and the ground truths of Example 7.

\begin{table}
\begin{tabular}{c c c c} \hline \hline \(x_{1}(t)\) & \(x_{2}(t)\) & Inputs & Labels & Notes \\ \hline \(\cos(k\pi t)\) & \(\cos((k+1.5)\pi t)\) & & & \\ \(\cos((k+1.5)\pi t+t^{2}+\cos(t))\) & & & \\ \(0\) & \(\cos((k+1.5)\pi t+t^{2}+\cos(t))\) & & & \\ \(\cos((k+1.5)\pi t+t^{2}+\cos(t))\) & & & & \\ \hline \(\cos(k\pi t)\) & & & & \\ \(\cos(k\pi t+t^{2}+\cos(t))\) & & & & \\ \(0\) & \(\cos(kl\pi t+t^{2}+\cos(t))\) & & & \\ \hline \hline \end{tabular}
\end{table}
Table 11: Inputs disturbed by the Gaussian noise with the SNR set to \(25dB\), and the labels used in Section 4.4, where \(t\in[0,6]\).
an orthogonal constraint, we directly add an inner product term to the loss function to control the orthogonality, i.e., \(\gamma\frac{|c_{1},c_{2}>|}{\|c_{1}\|_{2}\|c_{2}\|_{2}}\), where \(\gamma\) denotes the positive penalty parameter, and \(c_{1}\) and \(c_{2}\) are two components that are orthogonal. Then we train the model by the back propagation method based on the loss function. Although there is no guarantee of convergence in this case, it is simple, computationally efficient, and basically meets the expectation of orthogonality from the experimental results.

The Fourier basis is adopted for the construction of the inputs of RRCNN because they are orthogonal. The constructed inputs and labels are given in Tab. 13. After the RRCNN is trained, we use it to decompose the signal given in Example 8, which is composed of two mono-components that are orthogonal and close in frequency.

**Example 8**: \(x(t)=\cos(7t)+\sin(9t)\), \(t\in[0,2\pi]\).

The orthogonality and the errors between the resulting components and the corresponding ground truths are reported in Tab. 14. According to the results, the EWT and DCT_GAS_FDM methods perform best in terms of orthogonality, which is attributed to the fact that the former is based on the wavelet transform and the latter is based on the discrete cosine transform, both of which have strong orthogonal constraint on the decomposition results. For RRCNN, its orthogonality is promoted by minimizing the loss function. Therefore, on one hand, the results of RRCNN tend to find a balance in each item of the loss function. On the other hand, the limited iterative process cannot ensure that the results are completely converges to the true solution. Combined with orthogonality and error metrics, the overall performance of the RRCNN model is still satisfactory. Specifically, it is not the best in terms of orthogonality, but it is also very close to orthogonality, outperforming the other models except EWT and DCT_GAS_FDM; in terms of error, its two components are also the closest to the true components.

Moreover, we select the top three methods in terms of the error metrics from Tab. 14, that is, FDM, DCT_GAS_FDM, and RRCNN, and plot their obtained components, errors, and TFE representations in Fig. 5 (c), (f), respectively. From the plots, we can draw a conclusion similar to Example 6, that is, the RRCNN model can obtain the performance that is comparable to the state-of-the-art methods in the time domain, especially at the boundary, which can weaken the the impact of incomplete waveforms at the boundary to a certain extent.

\begin{table}
\begin{tabular}{c c c c c c c c c c c} \hline \hline \multirow{2}{*}{Method} & \multirow{2}{*}{\(\rho(c_{1},c_{2})\)} & \multicolumn{2}{c}{\(c_{1}\)} & \multicolumn{2}{c}{\(c_{2}\)} & \multirow{2}{*}{Method} & \multirow{2}{*}{\(\rho(c_{1},c_{2})\)} & \multicolumn{2}{c}{\(c_{1}\)} & \multicolumn{2}{c}{\(c_{2}\)} \\ \cline{5-10}  & & MAE & RMSE & MAE & RMSE & & MAE & RMSE & MAE & RMSE \\ \hline EMD[1] & 0.1141 & 0.1796 & 0.2776 & 0.2615 & 0.3581 & INCMD[11] & 0.1785 & 0.1846 & 0.2485 & 0.6205 & 0.7004 \\ VMD[17] & 0.4351 & 0.5682 & 0.6303 & 0.5682 & 0.6303 & SYNSQ_CWT[20] & 0.6268 & 0.7697 & 0.9299 & 0.4094 & 0.4823 \\ EWT[21] & **0.0000** & 0.6364 & 0.7070 & 0.6364 & 0.7070 & SYNSQ_STFT[20] & 0.1561 & 0.3225 & 0.3880 & 0.3202 & 0.3835 \\ IF[9] & 0.6896 & 0.2262 & 0.2776 & 0.2519 & 0.3092 & M-LFBF[16] & 0.9543 & 0.3858 & 0.4721 & 0.3857 & 0.4720 \\ FDM[22] & 0.0939 & 0.1229 & 0.2047 & 0.1244 & 0.2035 & DCT_GAS_FDM[23] & **0.0000** & 0.1222 & 0.1898 & 0.1211 & 0.1824 \\
**RRCNN** & 0.0595 & **0.1195** & **0.1759** & **0.0639** & **0.0889** & & & & & \\ \hline \hline \end{tabular}
\end{table}
Table 14: Metrics of the orthogonality, and errors of the obtained components by different methods in Example 8.

\begin{table}
\begin{tabular}{c c c c c} \hline \hline \(x_{1}(t)\) & \(x_{2}(t)\) & Inputs & Labels & Notes \\ \hline \(\cos kt\) & \(\sin(k+l)t\) & \(x_{1}(t)+x_{2}(t)\) & \(\left[x_{2}(t),x_{1}(t)\right]\) & \(k=6,7,8,9\) \\  & & & \(l=3,5,\ldots,33\) \\ \hline \hline \end{tabular}
\end{table}
Table 13: Inputs and labels used in Section 4.5, where \(t\in[0,2\pi]\).
However, due to the lack of a priori information in the frequency domain in the design of RRCNN, its advantage in terms of TFE distribution might be reduced when compared with other methods especially in the middle of the signal. However its performance are overall better than existing methods at the boundaries, also in the time-frequency domain.

### Can RRCNN be effective on solutions of differential equations?

To address this question, and following what has been done in the seminal work by Huang et al. [1], we test RRCNN against the solutions of a Duffing and Lorenz equations. The labels are necessary in the training process of the RRCNN model, however, they are not known in this instance. Since the EMD method works well for these two types of signals, we decided to use the results of EMD as the labels to verify the learning and generalization capabilities of the RRCNN model.

**Example 9**: We consider the following Duffing equation: \(\ddot{x}+\alpha x+\beta x^{3}=\gamma\cos(\omega t).\) We focus our attention on the decomposition of \(\dot{x}(t)\). We first construct the inputs by solving the equation using the Euler method with the parameters set to \(\alpha\in[-1,0)\), \(\beta\in[0.8,1.2]\), \(\omega\in[1.0,1.4]\) with step size \(0.1\), and \(\gamma\in[0.04,0.06]\) with step size \(0.005\), respectively, where \(t\in[0,300]\) and the initial conditions: \(\{x(0)=1,\dot{x}(0)=1\}\). And then, the first two IMFs of each input obtained by the EMD are collected as the labels.

After the RRCNN is trained, we apply it to decompose the signal \(\dot{x}(t)\) with \(\alpha=-0.85,\beta=1.05,\gamma=0.045,w=1.3\), which was not included in the training data. The results are depicted in Fig. 6. Both the resulting components produced by RRCNN are very close to those of EMD.

**Example 10**: The Lorenz equation is mathematically expressed as: \(\dot{x}=-\sigma(x-y),\dot{y}=rx-y-xz,\dot{z}=-bz+xy.\) We take into account the decomposition of \(x(t)\). The inputs are the \(x(t)\) achieved by the ode45 (code: [https://www.mathwo](https://www.mathwo) rks.com/help/matlab/ref/ode45.html) method by changing the parameters \(\sigma\in[9,11]\) with step size \(0.2\), \(r\in[19,21]\) with step size \(0.2\), and \(b\in[2,5]\) with step size \(1\), where \(t\in[0,50]\) and initial conditions: \(\{x(0)=-10,y(0)=0,z(0)=0\}\). Similarly to Example 10, we treat the first two IMFs of each input produced

Figure 6: Components of \(\dot{x}(t)\) of Duffing equation with \(\alpha=-0.85,\beta=1.05,\gamma=0.045,w=1.3\) by EMD and RRCNN. First row: the original signal. The 2nd-3rd rows, left: component obtained by EMD (blue solid curve) and RRCNN (red dotted curve); right: the error between the obtained components by EMD and RRCNN.
by the EMD as the labels. The results for signal \(x(t)\) with parameters set to \(\sigma=10.5,r=20.5,b=3.5\), predicted by the trained RRCNN model are shown in Fig. 7. The results show that the RRCNN has good learning and generalization capabilities for the solution to Lorenz equation, and can basically achieve the performance of EMD.

### 4.7 Is RRCNN capable of processing real signals?

We hereby employ the RRCNN model to process the real data, i.e., the length of day (LOD, data source: [http://hpiers.obspm.fr/eoppc/eopc04/eopc04.62-now](http://hpiers.obspm.fr/eoppc/eopc04/eopc04.62-now), start date: Jan. 1,1962, end date: May 24, 2022), which is widely used in the verification of the signal decomposition methods.

To generate the signals from LOD data for use as the inputs to train the RRCNN model, we first split the LOD into a series of length 720 (about two years) segments with a stride of 180 (about half year). And then, similar to the situation in which we were dealing with the solutions of differential equations, the main challenge in

Figure 8: Comparisons of the components obtained by EWT and RRCNN for two LOD data given in Examples 11-12. Top panel: the original signal. The 2nd-6st panels: the obtained components by EWT (blue solid curve) and RRCNN (red dotted curve).

Figure 7: Components of \(x(t)\) of Lorenz equation with \(\sigma=10.5,r=20.5,b=3.5\) by EMD and RRCNN. Same as Fig. 6.
training the RRCNN model on real signals is the ground truth calibration. We use the EWT method to produce the labels for each segment, because it can produce an artificially pre-set number of components.

**Example 11**: After the RRCNN model is trained, we apply it to the input, which is the LOD data ranging from Dec. 12, 2018 to Nov. 30, 2020, and is excluded in the training dataset.

**Example 12**: We also apply the trained RRCNN model to another input, which is the LOD data ranging from Dec. 7, 2019 to Nov. 25, 2021, and is excluded in the training dataset.

The components obtained from EWT and RRCNN for the input are depicted in Fig. 8 (a)-(b). The plots prove that RRCNN can approximate EWT well.

### Computational time of different methods

To compare the computational time of the proposed method with other methods, we apply them to decompose the signal in Example 4, and record the corresponding runtime for comparison. Since in the prediction phase, RRCNN has the advantage of processing the decomposition of multiple signals in parallel, we also compare the running time of these methods when decomposing multiple signals. In this case, we no longer construct new signals to be decomposed, but just repeat the processing of the signal in Example 4, which does not affect the comparison of running time.

All the experiments are performed in Python 3.8.12 or Matlab R2020b on a Dell Precision 5820 tower with Intel(R) Xeon(R) W-2102 processor (2.90 GHz), 64G memory, and Ubuntu 18.04.3 operating system.

We depict the running time of different methods when decomposing different numbers of signals in Fig. 9. From it, we have the following findings: When decomposing only one signal, the computational efficiency of RRCNN ranks in the middle among the 11 methods compared in this work. However, as the number of signals increases to 10, the computational efficiency of RRCNN begins to improve, and its ranking rises to 4th, surpassing the IF and VMD methods; when the number increases to 100, its computational efficiency ranks 2nd, beyond the EMD and DCT_GAS_FDM methods. Although there is still a gap between RRCNN and the EWT method that ranks 1st in computational efficiency, it can be seen that as the number of signals increases, the gap between them shrinks significantly.

Figure 9: Comparison of computational time of different signal decomposition methods when decomposing different numbers of signals.
## 5 Conclusion

In the paper, we use deep learning techniques to tackle the issue of decomposing a non-stationary signal into oscillatory components. Specifically, we first construct the RRCNN inner loop block for obtaining the local average of a given signal, and then these blocks are cascaded into a deeper network, called RRCNN, which is used to decompose a given signal.

Since the proposed RRCNN model is based on the deep learning framework and is a supervised model, it has the advantages of these two types of models. First, the convolutional filter weights in the model are learnt according to the input signals, which makes the proposed method more adaptive. Second, some common tools in deep learning, like the residual structure and the nonlinear activation function, can be added to increase the expressive ability and flexibility of the proposed model. Third, the proposed model can be customized according to the application. For example, when processing signals composed of orthogonal components, an inner product term can be added to the loss function to enhance the orthogonality of the derived components. To verify the performance of the proposed model, we compare it with other existing models from seven aspects. And the artificial and real data are used in the experiments. All results show that the proposed method works better in handling the boundaries, mode mixing effects and the orthogonality of the decomposed components, and is more robust than the existing methods.

On the other hand, RRCNN has the classical limitations of supervised models. For example, the labels must be given in advance in the training phase. Therefore, the goal of this work is not to replace any existing method, but to propose a completely new kind of approach which is based on deep learning and supervised learning and to add it to the existing methods, to derive a more flexible and adaptive processing approach for signals. In the future, we plan to work on the extension of this work to multivariate signals. Furthermore, given that the main limitation of the RRCNN approach is that it requires other decomposition methods to calibrate the ground truths before the training phase, we plan to work in the future on the development of an unsupervised learning method to produce a new technique with similar performance to the RRCNN algorithm, but that does not require any training based on other techniques.

## Declaration of Competing Interest

No conflict of interest exits in the submission of this manuscript, and manuscript is approved by all authors for publication. We declare that the work described was original research that has not been published previously, and not under consideration for publication elsewhere, in whole or in part.
## Acknowledge

F. Zhou research was partially supported by the National Natural Science Foundation of China [grant number 11901113], the Guangdong Basic and Applied Basic Research Foundation [grant number 2020A1515110951], and the China Scholarship Council [grant number 202008440024]. A. Cicone is a member of the Italian GNCS of the INdAM. H. Zhou research was partially supported by the National Science Foundation of US [grant number DMS-1830225], and the Office of Naval Research [grant number N00014-18-1-2852].
