An Innovative Nonlinear Bounded Component Analysis Algorithm Based on Multivariate Nonlinear Chirp Mode Decomposition

Abstract

In complex and diverse practical application scenarios, the challenge of blind source separation under underdetermined and nonlinear conditions is often encountered. To address this challenge, this paper proposes an innovative underdetermined nonlinear bounded component analysis method. This method first employs Multivariate Nonlinear Chirp Mode Decomposition (MNCMD) to process and reconstruct the observed signals, transforming the original underdetermined problem into a positive definite problem. Subsequently, Gaussianization techniques are introduced as a means of nonlinear compensation, successfully converting the nonlinear model into an analyzable linear model, laying a solid foundation for subsequent signal separation. Finally, the signal is separated by the bounded component analysis method, which does not require the source signals to be independent of each other. To validate the effectiveness and superiority of the proposed algorithm, detailed simulation experiments were designed and implemented. The experimental results demonstrate that compared to traditional underdetermined blind source separation algorithms, the algorithm presented in this paper exhibits significant advantages in terms of universality, convergence speed, separation accuracy, and robustness. Furthermore, this paper successfully applies the algorithm to the blind extraction of fetal electrocardiogram (FECG) signals from real datasets. The experimental results show that the algorithm can rapidly and effectively extract clearer and more accurate FECG signals, demonstrating its great potential and value in practical applications.

1. Introduction

Blind source separation (BSS) is a signal processing method where, when the characteristics and quantity of the source signal cannot be confirmed and the mixing matrix of the transmission channel is also unknown, only the observed mixed signal can be used to obtain the desired source signal for processing [1]. Its purpose is to restore the source signal to the maximum extent through an estimation method. The underdetermined problem refers to the fact that the number of observable signals is lower than the number of source signals [2]. In practical applications, underdetermined BSS has become a hot research topic due to various objective factors that can easily cause underdetermined problems. For underdetermined BSS, research began in 1991, when J F. Cardoso proposed the use of the fourth order cumulant method [3], followed by the underdetermined system localization method [4] and sparse component analysis (SCA) [5], which became the mainstream. Domestic scholars followed suit and proposed algorithms such as the parallel factor method (PFA) [6], the embedded vector method [7], and the density space clustering method [8]. However, most of these algorithms are aimed at the simplest linear mixing and have limitations when facing nonlinear mixing problems [9]. A nonlinear mixing problem refers to the mixing model between source signals being nonlinear [10], which is usually caused by the observed signal undergoing nonlinear distortion, and often occurs in daily life. The research on nonlinear BSS began relatively early. In 1992, Borel analyzed deterministic nonlinear mixed models with unknown parameters in detail [11]. In recent years, popular algorithms began to include the Bayesian method [12], the radial basis function neural network method [13], the back propagation neural network algorithm [14], etc. However, most of the current nonlinear BSS algorithms are still based on the positive definite case. The common constraint of the above algorithms is that the source signals are required to be independent of each other.

Bounded component analysis (BCA) is a new BSS method that has attracted scholars’ attention in recent years. In BCA, the hypothetical conditions of traditional BSS methods are replaced with the boundedness of the signal, which utilizes the set compactness and Cartesian separability of the source signal for separation. This assumption is more lenient and allows for separation without the need for source signals to be independent of each other [15]. Many scholars have made improvements to the BCA algorithm. Cruces proposed the goal of minimizing the convex perimeter function to extract a single source signal from a mixed signal [16]. Erdogan proposed a BCA method with the maximizing volume ratio as the objective function [17]. Gong et al. [18] proposed a normalized boundary objective function that simplifies the objective function of BCA.

Currently, most FECG signal extraction algorithms are devised based on ideal linear instantaneous mixture models. However, these methods may encounter issues like often necessitating multi-channel ECG data, which adds complexity and resource consumption [19]. Considering that the maternal electrocardiogram (MECG) signal source is located far from the fetal position in the abdomen, as it travels through internal organs to reach the fetus, it inevitably undergoes nonlinear distortions. Hence, the measured mixture of maternal and FECG signals in the abdomens of pregnant women is better modeled as a nonlinear mixture model [20]. Therefore, FECG extraction can be used in applications involving underdetermined nonlinear BSS algorithms.

This article proposes a nonlinear BCA algorithm based on MNCMD to solve problems such as underdetermined observation signals, nonlinear mixing of source signals, and non-mutual independence. This algorithm can extract and separate independent and non-independent source signals from mixed signals without the need for prior information. The computer experiments and the results of simulations involving clinical electrocardiogram signal datasets and artificial synthesized electrocardiogram signal datasets have verified the correctness and effectiveness of the algorithm proposed in this paper. At the same time, it has been proven that this algorithm not only has stronger universality, but also outperforms traditional FECG signal extraction algorithms in terms of accuracy.

2. Materials and Methods

2.1. Underdetermined Nonlinear Mixed Model

Assuming that the transmission of the signal is instantaneous, the mixed signals can be expressed as follows:

$$x_i(t)=f_i({\sum }_{j=1}^n{a}_{ij}{s}_j\left(t\right)),i=1,\dots ,m$$

(1)

where $x_i\left(t\right)$ denotes the mixed signals received by sensors, $s_j(t)$ is the source signals, and the above equation can be expressed as Equation (2) using a matrix.

$$x=f\left(As\left(t\right)\right)$$

(2)

where $s(t)={\left[s_1\right(t),s_2(t),s_3(t)\dots s_m(t\left)\right]}^T$. A is an m $\times$ n dimensional mixed matrix [21], where m < n means that the number of observed signals is less than that of the source signal, $f$ is an unknown reversible nonlinear function, and $a_{ij}$ (i, j = 1,2,… n) is an element of an unknown nonsingular matrix A. The entire underdetermined nonlinear mixing process is shown in Figure 1, and the nonlinear part occurs post nonlinear mixing [22].

Separation under this condition should first solve the underdetermined problem of signal elevation and then nonlinear BSS, as shown in Figure 2.

The components of the output $y(t)={[y_1(t),y_2(t)\dots y_n(t)]}^T$ in the figure are as follows:

$$y_i(t)={\sum }_{j=1}^n w_{ij}{g}_j(h_j(t);{\theta }_j),i=1,\dots ,n.$$

(3)

A nonlinear function $g_j$(j = 1, 2…, n) with undetermined parameters ${\mathsf{\theta }}_{\mathrm{i}}$ is a parameter adaptation function, which is used to correct the distortion caused by $f_i$(i = 1, 2…, n); $w_{ij}$(i, j = 1, 2…, n) is an element of an unknown n $\times$ n-order nonsingular matrix W. The complete mixed separation system can be represented as follows:

$$\begin{array}{c}e\left(t\right)=A^{*}s\left(t\right)\\ x\left(t\right)=f\left(e\left(t\right)\right)\\ x\left(t\right)\to \left(t\right)\\ \begin{array}{c}z\left(t\right)=g\left(\left(t\right)\right)\\ y\left(t\right)=W^{*}z\left(t\right)\end{array}\\ \end{array}.$$

(4)

If $g=f^{-1}$ and W = PDA⁻¹ (P and D are permutation matrices and diagonal scaling matrices, respectively), then $y\left(t\right)=P \ast D \ast s\left(t\right)$, i.e., separation is achieved, and the output y is equal to the source signal after rearranging and scaling the components [23].

2.2. Nonlinear Bounded Component Analysis (NLBCA)

The biggest advantage of the BCA algorithm is that it can separate non-independent signals, which brings stronger universality to blind source separation algorithms. Its basic principle requires some necessary assumptions to be made [24,25]:

(1) Assume that the mixing matrix A has rank n, i.e., column full rank, i.e., the number of sensors should be greater than or equal to the number of source signals.

(2) Assume that the distribution of source signals is bounded.

(3) Assume that the branch set S, consisting of the joint distribution of all source signals, can be expressed by the Cartesian product of the branch set Si of each source signal:

$$\mathrm{S}={\mathrm{S}}_1\times {\mathrm{S}}_2\times \cdots \times {\mathrm{S}}_{\mathrm{n}}$$

(5)

where n is the number of source signals and × denotes the Cartesian product.

Assumptions (2) and (3) are weaker and more loosely constrained. Assumption (3) is automatically satisfied when the signals are independent of each other, so the BSS algorithm based on BCA has stronger generalizability.

With the help of the analysis above, it is clear that the BCA algorithm’s steps are to find the separation matrix W when the mixed signal X is known, and then multiply W by X to determine Y, which is the estimated value of the source signal S. Then, the process of the BCA algorithm can be simply summarized as solving the separation matrix to approximate the independent source signal, and the specific approximation method is actually an optimization search of the established objective function.

In this paper, the normalized boundary objective function is selected as follows:

$$J(y)={\displaystyle \frac{\hat{R}(y)}{\sqrt{E\left\{y^2\right\}}}}$$

(6)

where $\mathrm{y}$ denotes the separation signal, $\mathrm{E}\left\{{\mathrm{y}}^2\right\}$ denotes the variance of $\mathrm{y}$, and $\hat{\mathrm{R}}\left(\mathrm{y}\right)$ is the boundary operator.

The idea of the NLBCA algorithm is to introduce a nonlinear unmixing function to compensate for the observed signal, then select a normalized boundary objective function, and solve the unmixing matrix W through an optimization algorithm, ultimately achieving signal separation.

Assuming that v is the observed signal after nonlinear compensation,

$$v=F_v^{-1}\left[F_x\right(x\left)\right].$$

(7)

The cumulative distribution function, $F_x\left(x\right)$, of the observed signal x can be directly obtained through the kernel estimation method. Based on the previous model, assuming that Q is the covariance matrix of x, the normalized boundary objective function of nonlinear BCA can be written as follows:

$$J(W)={\displaystyle \frac{\hat{R}\left(Wv\right)}{\sqrt{WQ{W}^T}}}.$$

(8)

Since Equation (8) is non-differentiable and its gradient does not exist, a subgradient algorithm is used to solve it. The update rules are as follows:

$$W\left(k\right)=W\left(k-1\right)+\delta {\left(h\left(k-1\right)\right)}^{-{\displaystyle \frac{1}{2}}}({\displaystyle \frac{\hat{R}\left(W\left(k-1\right)v\right)}{h\left(k-1\right)}}W\left(k-1\right)Q-b_m^{\left(k-1\right)T})$$

(9)

where $\mathrm{h}\left(\mathrm{k}-1\right)=\mathrm{W}\left(\mathrm{k}-1\right)\mathrm{Q}{\mathrm{W}\left(\mathrm{k}-1\right)}^{\mathrm{T}}$, δ is the iteration step size. If the step size is fixed, it will cause the algorithm to always miss the best advantage. It is necessary to gradually reduce the step size value to ensure that the algorithm ultimately converges to the best advantage. Here, the adaptive weight is used as the attenuation factor, and the specific method is to handle the step size as follows:

$${\delta }_k=[\mathit{cos}\left({\displaystyle \frac{\pi \times k}{2\times M}}+{\displaystyle \frac{\pi }{2}}\right)+1]\times {\delta }_0$$

(10)

where M is the maximum number of iterations, ${\mathsf{\delta }}_0$ is the initial step size, and k represents the k-th iteration. Let the set of moments at which the output value of the m-th path is maximized and minimized during the (k − 1)-th iteration be ${\mathsf{\gamma }}_{\mathrm{m}+}$ and ${\mathsf{\gamma }}_{\mathrm{m}-}$, then:

$$b_m^{\left(k-1\right)}=\sum _{{\left(k-1\right)}_{m+}\in {\gamma }_{m+}}{\tau }_{m+}^{\left(k-1\right)}{(k-1)}_{m+}y{(k-1)}_{m+}-\sum _{{\left(k-1\right)}_{m-}\in {\gamma }_{m-}}{\tau }_{m-}^{\left(k-1\right)}{(k-1)}_{m-}y{(k-1)}_{m-}$$

(11)

where $\mathrm{y}{(\mathrm{k}-1)}_{\mathrm{m}+}$, $\mathrm{y}{(\mathrm{k}-1)}_{\mathrm{m}-}$ are the estimated vectors of the source signal at time ${(\mathrm{k}-1)}_{\mathrm{m}+}$ and time ${(\mathrm{k}-1)}_{\mathrm{m}-}$, respectively, and ${\mathsf{\tau }}_{\mathrm{m}+}^{\left(\mathrm{k}-1\right)}{(\mathrm{k}-1)}_{\mathrm{m}+}$ and ${\mathsf{\tau }}_{\mathrm{m}-}^{\left(\mathrm{k}-1\right)}{(\mathrm{k}-1)}_{\mathrm{m}-}$ are joint coefficients that satisfy the following:

$${\tau }_{m+}^{\left(k-1\right)}{(k-1)}_{m+}\ge 0,{\left(k-1\right)}_{m+}\in {\gamma }_{m+}, {\sum }_{{\left(k-1\right)}_{m+}\in {\gamma }_{m+}}{\tau }_{m+}^{\left(k-1\right)}{\left(k-1\right)}_{m+}=1$$

(12)

$${\tau }_{m-}^{\left(k-1\right)}{(k-1)}_{m-}\ge 0,{\left(k-1\right)}_{m-}\in {\gamma }_{m-}, {\sum }_{{\left(k-1\right)}_{m-}\in {\gamma }_{m-}}{\tau }_{m-}^{\left(k-1\right)}{\left(k-1\right)}_{m-}=1.$$

(13)

As the iteration progresses, the output scale may become larger or smaller, so it is necessary to form constraints through normalization. The pseudocode of NLBCA is given in Algorithm 1.

$$W\left(k\right)={\displaystyle \frac{W\left(k\right)}{\sqrt{W\left(k\right){W\left(k\right)}^T}}}$$

(14)

Algorithm 1 NLBCA

Inputs: maximum number of iterations M, line progress length ${\delta }_0$, initial matrix $W\left(0\right)$, iterations k = 0   1. Perform non-linear transformation on the observed signal x to obtain the nonlinear compensated signal v   2. Normalize $W\left(0\right)$ using Equation (14)   3.   While k < M   4.   k = k + 1   5.   Calculate ${\tau }_{m+}^{\left(k\right)}{\left(k\right)}_{m+}$ using Equation (12)   6.   Calculate ${\tau }_{m-}^{\left(k\right)}{\left(k\right)}_{m-}$ using Equation (13)   7.   Update the line progress length ${\delta }_0$ using Equation (10)   8.    Update the $W\left(k\right)$ using Equation (9)   9.   Normalize $W\left(k\right)$ using Equation (14) 10. end while 11.   $y=Wx$ 12. output results

2.3. Nonlinear Bounded Component Analysis Based on Multivariate Nonlinear Chirp Mode Decomposition (MNCMD-NLBCA)

MNCMD is a new signal feature extraction and decomposition denoising method proposed by Chen [26] that can effectively handle time-varying, nonlinear, and non-stationary multivariate signals. The principle is to first define a multivariate nonlinear frequency modulation mode based on the instantaneous frequency information between all input data channels, and then establish an objective function based on minimizing the sum of mode bandwidth on all signal channels to achieve accurate and stable signal decomposition. The core idea of MNCMD is as follows:

The signal $x\left(t\right)$ can be written as follows in nonlinear chirp mode [27]:

$$x(t)={\sum }_{i=1}^k g_i(t)$$

(15)

where k represents the number of nonlinear chirp mode components, and $g_i(t)$ is defined in MNCM as follows:

$${g_i\left(t\right)=a}_i(t)cos(2\pi {\int }_0^t f_i(e)de+{\gamma }_i)$$

(16)

where $a_i(t)$ is the instantaneous amplitude, $f_i(e)$ is the instantaneous frequency, and ${\mathsf{\gamma }}_{\mathrm{i}}$ denotes the initial phase. According to reference [28], Equation (16) can be rewritten as follows:

$$g_i\left(t\right)={\mu }_i(t)cos(2\pi {\int }_0^t {\tilde{f}}_i(e)de)+{\rho }_i(t)sin(2\pi {\int }_0^t {\tilde{f}}_i(e)de)$$

(17)

where ${\mu }_i(t)$ and ${\rho }_i(t)$ are written as follows:

$$\begin{array}{cc}& {\mu }_i(t)=a_i(t)cos(2\pi ({\int }_0^t f_i(e)de-{\int }_0^t {\tilde{f}}_i(e)de)+{\gamma }_i)\\ & {\rho }_i(t)=-a_i(t)sin(2\pi ({\int }_0^t f_i(e)de-{\int }_0^t {\tilde{f}}_i(u)du)+{\gamma }_i)\end{array}.$$

(18)

Assuming the number of signal channels is n, perform k-order decomposition and construct the following problem:

$$\begin{array}{cc}& \underset{\left\{{\mu }_{i,m}\left(t\right)\right\},\left\{{\rho }_{i,m}\left(t\right)\right\},\left\{{\tilde{f}}_i\left(t\right)\right\}}{min}{\displaystyle \sum _{i=1}^k} {\displaystyle \sum _{m=1}^n}\left\{{‖{\mu }_{i,m}^{″}\left(t\right)‖}_2^2+{‖{\rho }_{i,m}^{″}\left(t\right)‖}_2^2\right\}\\ & \mathrm{s}.\mathrm{t}. {‖x_m\left(t\right)-{\displaystyle \sum _{i=1}^k} {\mu }_{i,m}\left(t\right)\mathit{cos}\left(2\pi {\int }_0^t\widetilde{f_i}\left(e\right)de\right)-{\displaystyle \sum _{i=1}^k} {\rho }_{i,m}\left(t\right)\mathit{sin}\left(2\pi {\int }_0^t\widetilde{f_i}\left(e\right)de\right)‖}_2\le {\epsilon }_m,\\ & m=\mathrm{1,2},\dots n.\end{array}$$

(19)

Due to the fact that the collected signals are mostly in discrete form, Equation (19) is rewritten as follows:

$$\begin{array}{cc}& \underset{\left\{{\mu }_{i,m}\right\},\left\{{\rho }_{i,m}\right\},\left\{{\tilde{f}}_i\right\}}{min}\sum _{i=1}^k \sum _{m=1}^n\left\{{‖\Psi {\mu }_{i,m} ‖}_2^2+{‖\Psi {\rho }_{i,m} ‖}_2^2\right\}\\ & \mathrm{s}.\mathrm{t}. {‖x_m-\sum _{i=1}^k{L}_i {\mu }_{i,m}-\sum _{i=1}^k V_i{\rho }_{i,m}‖}_2\le {\epsilon }_m, \\ & m=\mathrm{1,2},\dots n\end{array}$$

(20)

where $L_i=diag[\mathit{cos}\left({\alpha }_i\left(t_1\right)\right),\mathit{cos}\left({\alpha }_i\left(t_2\right)\right),\dots \mathit{cos}\left({\alpha }_i\left(t_N\right)\right)]$, $V_i=diag[\mathit{sin}\left({\alpha }_i\left(t_1\right)\right),\mathit{sin}\left({\alpha }_i\left(t_2\right)\right),\dots \mathit{sin}\left({\alpha }_i\left(t_N\right)\right)]$, ${\alpha }_i\left(t\right)=2\pi {\int }_0^t\widetilde{f_i}\left(e\right)de$, $\Psi$ is a second-order difference operator.

$$\Psi =\left[\begin{array}{ccccc}-1& 1& 0& \cdots & 0\\ 1& -2& 1& \cdots & 0\\ ⋮& \ddots & \ddots & \ddots & ⋮\\ 0& \cdots & 1& -2& 1\\ 0& \cdots & 0& 1& -1\end{array} \right]$$

Introduce the Lagrange function to solve the optimal solution of the aforementioned variational problem [29].

$$\begin{array}{cc}& L(\left\{{\mu }_{i,m}\right\},\left\{{\rho }_{i,m}\right\},\left\{{\tilde{f}}_i\right\},\left\{{\theta }_m\right\},\left\{{\lambda }_m\right\})={\displaystyle \sum _{m=1}^n}{\zeta }_{{\mathcal{C}}_{{\epsilon }_m}}\left({\theta }_m\right)+{\displaystyle \sum _{i=1}^k} {\displaystyle \sum _{m=1}^n}\left\{{‖\Psi {\mu }_{i,m} ‖}_2^2+{‖\Psi {\rho }_{i,m} ‖}_2^2\right\}+\\ & {\displaystyle \sum _{m=1}^n} \left({\displaystyle \frac{\psi }{2}}{\Vert {\theta }_m+{\displaystyle \sum _{i=1}^Q}(L_i{\mu }_{i,m}+V_i{\rho }_{i,m})-x_m+{\displaystyle \frac{1}{\psi }}{\lambda }_m\Vert }_2^2-{\displaystyle \frac{1}{2\psi }}{\Vert {\lambda }_m\Vert }_2^2\right)\end{array}$$

(21)

where $\psi$ is the penalty parameter, ${\lambda }_m$ is the Lagrangian multiplier, and ${\theta }_m=-\sum _{i=1}^Q(L_i{\mu }_{i,m}+V_i{\rho }_{i,m})+x_m$, ${\zeta }_{{\mathcal{C}}_{{\epsilon }_m}}(\theta )\triangleq \{\begin{array}{cc}0,& \theta \in C_{\epsilon }\\ +\mathrm{\infty }& \theta \notin C_{\epsilon }\end{array}$, $C_{\epsilon }$ is a Euclidean ball whose radius is $\mathsf{\epsilon }$ and center is 0.

Finally, the alternating direction multiplier iterative algorithm is used to solve the variational problem. The expressions for ${\mu }_{i,m}$, ${\rho }_{i,m}{\widetilde{,f}}_i,{\theta }_m$, and ${\lambda }_m$ after alternating optimization iterations are as follows:

$${{\mu }_{i,m}}^{\omega +1}=\underset{\left\{{\mu }_{i,m}\right\}}{argmin}\left\{L(\left\{{\mu }_{i,m}\right\},\left\{{\rho }_{i,m}\right\},\left\{{\tilde{f}}_i\right\},\left\{{\theta }_m\right\},\left\{{\lambda }_m\right\})\right\}$$

(22)

$${{\rho }_{i,m}}^{\omega +1}=\underset{\left\{{\rho }_{i,m}\right\}}{argmin}\left\{L(\left\{{\mu }_{i,m}\right\},\left\{{\rho }_{i,m}\right\},\left\{{\tilde{f}}_i\right\},\left\{{\theta }_m\right\},\left\{{\lambda }_m\right\})\right\}$$

(23)

$${{\tilde{f}}_i}^{\omega +1}={\tilde{f}}_i+{\displaystyle \frac{\delta }{2\pi }}\times {\left({\displaystyle \frac{2}{\gamma }}{\Psi }^T\Psi +I\right)}^{-1}\left({\displaystyle \frac{d}{dt}}\left(arctan{\displaystyle \frac{{{\rho }_{i,m}}^{\omega +1}}{{{\mu }_{i,m}}^{\omega +1}}}\right)\right)$$

(24)

$${{\theta }_m}^{\omega +1}=Q_{{\mathcal{C}}_{{\epsilon }_m}}\left(x_m-{\sum }_{i=1}^k{L}_i {\mu }_{i,m}-{\sum }_{i=1}^k V_i{\rho }_{i,m}-{\displaystyle \frac{{\lambda }_m}{\psi }}\right)$$

(25)

$${{\lambda }_m}^{\omega +1}={\lambda }_m+\psi \left({\theta }_m+{\sum }_{i=1}^Q(L_i{\mu }_{i,m}+V_i{\rho }_{i,m})-x_m\right)$$

(26)

where $\delta$ is the learning rate, γ is the noise tolerance, and $Q_{{\mathcal{C}}_{{\epsilon }_m}}(·)$ is a proximity operator [30], $Q_{{\mathcal{C}}_{{\epsilon }_m}}(z)\triangleq \{\begin{array}{cc}{\displaystyle \frac{{\epsilon }_m}{\Vert z{\Vert }_2}}\cdot z,& \Vert z{\Vert }_2>{\epsilon }_m\\ z,& \Vert z{\Vert }_2\le {\epsilon }_m\end{array}$.

Based on the introductions of previous hybrid models, this paper proposes an MNCMD-NLBCA algorithm to solve the problem of underdetermined nonlinearity. The main idea is to perform MNCMD processing on the mixed signal, which is decomposed into several NCM components, and then these NCM components are reconstructed to obtain a new signal sequence. This new signal sequence can be merged with the original mixed signal to form a new observation signal, ultimately forming a positive definite BSS model. Users can perform nonlinear transformation on each observed signal to obtain nonlinear compensation, then select a normalized boundary objective function, and then solve the unmixing matrix W to recover the source signal.

The main process of the algorithm is as follows:

(1)

Pre-process the received mixed signal $X\left(t\right)$, including de-averaging and pre-whitening.

(2)

Perform MNCMD processing on the pre-processed observation signal to obtain k NCM components, $u_i\left(t\right)$.

(3)

Perform one-dimensional reconstruction on these k NCM components by assigning different random weights between (0, 1) and adding them together to obtain a new signal $X^{\prime }\left(t\right)$. Then, make $X\left(t\right)=[X\left(t\right);X^{\prime }\left(t\right)]$ to obtain a new observation signal.

(4)

Carry out nonlinear transformation of the new observation signals.

(5)

Select the normalized boundary objective function and use the subgradient descent algorithm to solve the mixed matrix W.

(6)

Complete signal separation.

Figure 3 shows a flow chart of the MNCMD-NLBCA algorithm.

3. Experiment and Results

3.1. Simulation Dataset

Using MATLAB, we generated three signals: two independent signals labeled s1 and s2, and a noise signal denoted as s3. s1 is a sinusoidal wave with a frequency of 10 Hz. s2 is a square wave featuring a 50% duty cycle. s3 represents the noise component.

All signals are sampled at a frequency Fs of 1000 Hz. Their respective waveforms and frequency spectra are presented in Figure 4, providing a clear visualization of their unique characteristics.

We generated a mixed signal through a mixed matrix and nonlinear transformation, with a signal-to-noise ratio (SNR) of 20 dB. Their waveforms and spectra are shown in Figure 5.

The waveforms and spectra of the observed signals after separation with the MNCMD-NLBCA algorithm are shown in Figure 6.

In order to see the separation effect more objectively and accurately, the correlation coefficient $\mathsf{\rho }$ between the extracted signal and the source signal was chosen as the evaluation criterion in this paper [31], as shown in (27), and the results are shown in

Table 1. Correlation coefficients.

	y1	y2	y3
s1	0.99782904	0.01841561	0.00165612
s2	0.00404252	0.98984297	0.04701977
s3	0.00741474	0.04752748	0.99377657

.

$$\rho ={\displaystyle \frac{\mathit{cov}(s_i,y_i)}{\sqrt{\mathit{cov}(s_i,s_i)\mathit{cov}\left(y_i,y_i\right)}}}$$

(27)

where $\mathit{cov}( )$ denotes the variance, $y_i$ is the separated signal, and $s_i$ is the source signal.

From Table 1, it can be seen that the correlation coefficients of s1 and y1, s2 and y2, and s3 and y3 are all close to 1, which corresponds to Figure 4 and Figure 6. From Figure 4 and Figure 6, it can be seen that independent source signals were successfully separated with the MNCMD-NLBCA algorithm.

To evaluate the robustness of this algorithm, we devised five distinct operating scenarios with varying SNRs by modulating the noise signals and subsequently performing BSS experiments. The signal interference ratio (SIR) [32] was used as the performance index. The equation of SIR is shown in Equation (28).

$$SIR=-10lg{\displaystyle \frac{{\Vert e_i-s_i\Vert }^2}{{\Vert s_i\Vert }^2}}$$

(28)

where $e_i$ is the separated signal and $s_i$ is the source signal. Signal processing-related knowledge tells us that the separation is improved by increasing the SIR value.

As depicted in Figure 7, it is evident that as the SNR increases, there is no significant decline in the SIR values. Furthermore, all SIR values remain above 24, indicating that the signals have been effectively separated. This demonstrates the strong robustness and resistance to interference exhibited by the algorithm.

To further validate the universality of this algorithm, the subsequent section of this paper covers three correlated source signals that follow a Copula-t distribution. Again, 1000 sample points are utilized for these non-independent sources. The waveform plots and spectrum diagrams of these generated correlated signals are presented in Figure 8.

With an SNR maintained at 20 dB, mixed signals are generated following the same methodology as previously described. The resulting waveform plots and spectrum diagrams of these mixed signals are presented in Figure 9.

Finally, the MNCMD-NLBCA algorithm is used to separate the observed signals. In order to generate a comparative effect, three different algorithms are used for blind source separation of the mixed signal: SCA, VMD-ICA, and EMD-ICA. The results are shown in Figure 10.

The correlation coefficients between the signals separated by the MNCMD-NLBCA algorithm and the source signals are shown in

Table 2. Correlation coefficients.

	y1	y2	y3
s1	0.97778428	0.52447726	0.21264333
s2	0.55792012	0.98634734	0.60709731
s3	0.2556597	0.54286605	0.96087454

.

From Table 2, it is evident that the correlation coefficients between s1 and y1, s2 and y2, and s3 and y3 are all nearly 1, aligning with Figure 8 and Figure 10. This indicates that the algorithm has successfully separated three highly correlated source signals from underdetermined nonlinear mixed signals at an SNR of 20 dB.

To facilitate an objective and quantitative assessment of the efficacy of these algorithms, we not only employed similarity coefficients and SIR, but also incorporated the mean squared error (MSE) [32] and the response time into our research methodology for a comprehensive comparison.

MSE is a fit metric commonly used to evaluate the difference between model predictions and actual observations. A smaller MSE value indicates a better fit of the model. Signal processing-related knowledge has shown that the separation is improved by increasing the SIR value.

The equation of MSE and SIR is shown in Equation (29).

$$\mathrm{M}\mathrm{S}\mathrm{E}={\displaystyle \frac{1}{\mathrm{N}}}{\sum }_{\mathrm{i}=1}^{\mathrm{N}}{\left({\mathrm{s}}_{\mathrm{i}}-{\mathrm{e}}_{\mathrm{i}}\right)}^2$$

(29)

It should be noted that due to the inability of blind source separation to accurately recover the amplitude of the source signal, in order to avoid the influence of amplitude when calculating evaluation indicators, it is necessary to normalize signals. The results are shown in

Table 3. Extraction performance evaluation index data.

	VMD-ICA	EMD-ICA	SCA	MNCMD-NLBCA
Similarity Coefficient	0.785	0.732	0.611	0.964
MSE	2.0204	3.1207	3.0299	0.0452
SIR	12.39	10.29	9.23	21.26
Response time	16.32	17.28	22.54	20.16

. It can be seen that the signal extracted by the MNCMD-NLBCA algorithm has the highest correlation coefficient with the source signal, and the SIR value is also the largest and produces the smallest MSE value, so the separation effect of the MNCMD-NLBCA algorithm is the best.

3.2. ADFECGD Dataset

ADFECGD is a small database of maternal abdominal ECG signals, created for the evaluation of detection algorithms. During the recording of the maternal abdominal electric signals, a repeatable arrangement of electrodes was used on the surface of a maternal abdomen. The configuration consisted of three measuring electrodes aligned in a horizontal line along the navel, and one measuring electrode that recorded the abdominal signal just over the navel. A reference electrode with “zero” potential for measuring was located above the pubic symphysis [33]. We selected the records “r01” and “r08”for the extraction of the FECG signal. They both consist of five-channel signals (5 min long) recorded during labor, with one channel of a direct fetal electrocardiogram. A more detailed description of the database can be found in [34].

Both records “r01” and “r08” exhibited high-quality signals with clear FECG components. These records were chosen to represent a diverse range of possible scenarios in the dataset and have been previously validated or used in other studies.

We use an ECG signal from one of five channels to investigate the accuracy of the extracted FECG morphology. The MNCMD-NLBCA method is compared with SCA, VMD-ICA, and EMD-ICA. Figure 11 and Figure 12 show visualizations of the extracted FECG to intuitively illustrate the effects of these methods.

Shown in the first row of the above figures are desired FECG signals that are direct fetal electrocardiograms in the records “r01” and “r08”. The second to fifth lines show the FECG signals extracted using the VMD-ICA, SCA, EMD-ICA, and MNCMD-NLBCA algorithms, respectively. From the two figures, it can be seen that the signal extracted by the SCA algorithm shows significant distortion compared to the source signal in the first row, while the distortion of the signal processed by the EMD-ICA algorithm has improved, but still exists. In Figure 11, the signals in the second and fifth lines indicate that the VMD-ICA and MNCMD-NLBCA algorithms have basically extracted the morphological features of FECG without significant distortion. In Figure 12, the FECG signal extracted by the VMD-ICA algorithm is morphologically close to the desired signal, but is still accompanied by slight distortion. From the waveform in the last line, it can be seen that the MNCMD-NLBCA algorithm effectively suppresses distortion and ensures the accuracy of signal extraction.

We normalized these signals and calculated the evaluation indicators, as shown in

Table 4. Extraction performance evaluation index data for simulation experiment.

Methods	Record “r01”			Record “r08”
	Similarity Coefficient	MSE	SIR	Similarity Coefficient	MSE	SIR
VMD-ICA	0.947	0.0312	11.23	0.914	0.0204	9.95
EMD-ICA	0.894	0.0391	6.02	0.901	0.0237	4.56
SCA	0.462	0.0263	8.01	0.611	0.0299	7.92
MNCMD-NLBCA	0.981	0.0306	14.36	0.964	0.0154	11.20

. From the perspective of the similarity coefficient, the MNCMD-NLBCA algorithm has the highest similarity between the extracted FECG signal and the desired signal in both experiments, with values of 0.981 and 0.964, respectively. It can be considered that this algorithm has extracted a relatively clean FECG signal, while the similarity coefficients of the other three algorithms’ extracted signals from high to low are VMD-ICA, EMD-ICA, and SCA.

In both experiments, the MSE values of all algorithms were relatively small, mostly around 0.03. The SIR value of MNCMD-NLBCA is significantly higher than that of the other algorithms, reaching 14.36 and 11.20, respectively. Therefore, based on the waveform comparisons in Figure 11 and Figure 12, as well as the evaluation index data in Table 4, it can be seen that the signals extracted by the MNCMD-NLBCA algorithm not only fully retain the morphological features of FECG, but are also the most similar to the desired signal, and have better performance compared to the other three mainstream algorithms.

Compared to the experiment reported in this article, the results of different recorded experiments may not be completely consistent. However, compared with the other three algorithms, the performance of our algorithm is still the best, so this does not affect the experimental conclusion.

4. Conclusions

In this paper, we propose a blind source separation algorithm based on MNCMD-NLBCA. This algorithm first processes the observed signals with MNCMD, then reconstructs the signals to increase their dimensionality, and finally combines with the NLBCA algorithm, thereby enhancing the universality and extraction accuracy of the algorithm. By introducing a nonlinear demixing function into BCA, this algorithm performs nonlinear compensation on the observed signals. Meanwhile, a step decay factor is employed during the iteration process to avoid falling into local optima and improve convergence performance. When combined with the MNCMD algorithm, the entire algorithm becomes applicable to underdetermined nonlinear scenarios, ensuring operational efficiency while achieving better optimization accuracy. The prior assumptions of this algorithm are more relaxed, as it does not require the source signals to be mutually independent, and fewer observed signals are needed. Therefore, it is highly suitable for FECG signal extraction.

To verify the effectiveness and superiority of this algorithm, simulation analysis and experimental verification were conducted. Firstly, BSS experiments were performed on three independent source signals that underwent underdetermined nonlinear mixing. The similarity coefficients between the separated signals and the source signals were 0.99782904, 0.98984297, and 0.99377657, respectively, which are all close to 1, thus proving the effectiveness of this algorithm. Additionally, robustness tests were conducted on this algorithm. As evident from Figure 7, varying SNRs had minimal impact on the experimental results of this algorithm, indicating its strong robustness. Then, under identical conditions, blind source separation simulations were conducted on cross-correlation signal sources, and the separation effects were compared with those of SCA, VMD-ICA, and EMD-ICA. The similarity coefficient, MSE, and SIR of the separated signals objectively and quantitatively demonstrated that this algorithm is also suitable for the separation of non-independent source signals, and it exhibits better separation accuracy compared to the other three algorithms. Finally, this algorithm was used to extract FECG signals from the ADFECGD real dataset and compared with three other algorithms. The morphological features presented in Figure 11 and Figure 12, along with the evaluation index data in Table 4, clearly show the effectiveness and superiority of this algorithm in extracting FECG signals.

Although this algorithm has strong universality, it still has certain limitations, such as being an offline algorithm that cannot be effectively applied in real-time signal processing. Moreover, the response time of the algorithm is still relatively long, and there is still room for improvement in the separation accuracy of the algorithm. In future research, we will further investigate online blind source separation algorithms to achieve blind source separation of real-time signals, and the application of this algorithm will be more extensive. We will also aim to further shorten the response time of the algorithm without affecting the separation accuracy.