Blind source separation method and system for single-channel extremely-low-frequency signals

Abstract

The invention provides a blind source separation method and system for a single-channel extremely-low-frequency signal, and the method comprises the steps: firstly carrying out the preprocessing of a collected signal through a Hanpur filter, and eliminating the pulse noise, aiming at a large amount of pulse noise generated by a natural environment in an extremely-low-frequency signal channel; secondly, decomposing the single-channel signal sequence into a plurality of independent components through singular spectrum analysis and a kmeans clustering algorithm to construct virtual multiple channels; and finally, separating a source signal with a high interference-to-signal ratio by using a feature matrix-based joint approximate diagonalization algorithm. According to the method, the accuracy of source signal estimation is improved, any probability distribution does not need to be assumed, parameters do not need to be manually selected, and good robustness and universality are achieved.

Description

Technical Field

[0001] The present invention relates to the technical field of blind source separation of extremely low frequency signals, and specifically, to a blind source separation method and system for single-channel extremely low frequency signals. Background Art

[0002] Blind source separation evolved from the "cocktail party problem", that is, in a noisy cocktail party, there are musical instruments and various noises. When many participants talk together, people need to exclude the unwanted sounds in these mixed sound sources and extract the speech of the target speaker. Mapping this problem to the field of signal processing. According to the relationship between the number of source signals and the number of observed signals, blind source separation can be divided into three categories: overdetermined blind source separation, well-determined blind source separation, and underdetermined blind source separation. In the extremely low frequency communication environment, the wavelength of the signal can reach 10 Mm to 100 Mm. Therefore, a single long antenna is usually required to transmit and receive single-channel signals. These signals can propagate very long distances. Due to their extremely low frequency, these signals can penetrate the atmosphere and the earth's crust during propagation, so they are widely used in fields such as deep-sea communication and geophysical exploration. Traditional blind source separation algorithms need to fully explore different potential characteristics of the signals and use the time-frequency characteristics and statistical characteristics of the signals to estimate the source signals. However, this algorithm is no longer applicable in the single-channel case of extremely low frequency communication.

[0003] At present, the development of multi-channel blind source separation algorithms has been relatively perfect and certain research results have been achieved. However, most of these algorithms are only applicable to the underdetermined blind source separation scenario and have obvious defects under single-channel reception conditions, and the separation effect is poor. Because there is only one received signal in single-channel blind source separation and the available effective information is less, there is a fundamental change in the processing method. Single-channel blind source separation requires estimating a large amount with very little amount, and the separation process needs to rely on the potential characteristics of the source signals. Therefore, there is still no general processing algorithm.

[0004] Patent application document CN118503659A discloses a single-channel signal separation method based on wavelet denoising and modal decomposition. However, the modal decomposition of this method is sensitive to impulse noise and may lead to mode mixing. And the parameter selection of this method is relatively complex and the calculation amount is large.

[0005] Patent application document CN118964917B discloses an underdetermined blind source separation method, device and electronic equipment based on continuous wavelet transform. This patented method is for multi-channel signals, and it may not be able to effectively separate for single-channel signals or cases with a small number of signal channels.

[0006] Patent application document CN118227981A discloses a single-channel vibration signal blind source separation method in multiple frequency bands. The effectiveness of this method depends on the accurate estimation of parameters such as the mixing matrix, the number of models, and signal separation. If these parameters are not well estimated, it may affect the performance of the model. In addition, this method involves multiple stages and requires the selection of parameters such as the number of layers and model size, which increases the complexity of implementation.

[0007] Patent application document CN116052707A discloses a single-channel blind source separation method and system. This method has relatively high requirements for the input signal: to ensure the separation effect, this method requires a relatively high quality of the input signal, otherwise it may affect the final separation effect, especially in the case of relatively high environmental noise.

[0008] Patent application document CN116415138A discloses a blind source separation method for single-channel ultrasonic signals, including obtaining ultrasonic signals in the target environment, identifying the number of local power sources, and quickly performing virtual separation of multi-channel signals, etc. The performance of this method depends on the effective acquisition of ultrasonic signals in the target environment. If the signal acquisition is inaccurate or the noise is too high, it may affect the separation effect. In addition, it is necessary to accurately identify the number of local power sources, which may become difficult in some complex environments, especially in the case of mixed multiple signal sources.

[0009] The model proposed in the literature "Electrical Drive Single-Channel Noise Source Separation and Identification Based on EEMD-AESSAICA" is sensitive to factors such as the sampling frequency of data, noise intensity, or sensor arrangement, which may affect its scope of application. Although good results have been achieved in electrical drive noise separation, it is not clear whether it is applicable to other types of noise sources. If the noise type and environment change greatly, it may be necessary to readjust or develop corresponding strategies.

[0010] The literature "Single-Channel Blind Source Separation Algorithm for Energy Difference Mixed Signals". This algorithm assumes that the mixed signal comes from two main signal sources, which may not hold in the case of more mixed signal sources, resulting in a decrease in the separation effect. This method relies on frequency thresholds and correlations to distinguish source signals, but when the frequency overlap of source signals is large, signal discrimination will be affected. Summary of the Invention

[0011] Aiming at the defects in the prior art, the purpose of the present invention is to provide a blind source separation method and system for single-channel extremely low-frequency signals.

[0012] The blind source separation method for single-channel extremely low-frequency signals provided by the present invention includes:

[0013] Step 1: Collect the single-channel mixed signal and set the parameters for separating the mixed signal, including the sample length and the number of mixed signals; preprocess the collected mixed signal to remove strong impulse noise and enhance the separability and stability of the signal;

[0014] Step 2: Perform singular spectrum analysis on the preprocessed mixed signal, construct the Hankel matrix of the signal, and perform singular spectrum decomposition on it to decompose the mixed signal into multiple independent signal components; record the eigenvectors and eigenvalues obtained from the decomposition, and construct the eigenmatrix of each signal component;

[0015] Step 3: Use the K-means clustering method to perform clustering analysis on the eigenvalues extracted from the multiple independent components obtained by singular spectrum decomposition to form virtual multi-channels;

[0016] Step 4: For the virtual multi-channels obtained by clustering, use the method of joint matrix diagonalization to restore the signal one by one; for low-frequency signals, perform diagonal averaging operations using the eigenmatrix reconstruction method to achieve blind source separation;

[0017] Step 5: Evaluate the separation quality of the separated signal, calculate the similarity index of signal separation. If the signal separation effect does not meet the expectation, adjust the number of components decomposed by singular spectrum or the clustering parameters, and return to Step 2; otherwise, output the multi-channel signal separation result.

[0018] Preferably, in Step 1, a Hampel filter is used to preprocess the mixed signal, including:

[0019] Step 1.1: For each point in the data sequence, calculate the median and standard deviation of this point within a sliding window of a fixed size;

[0020] Step 1.2: For each data point, calculate its deviation from the window median to reflect the difference between the data point and its neighborhood;

[0021] Step 1.3: According to the standard deviation and median, set a threshold. If the deviation of the data point from the median is greater than a preset threshold, then this data point is considered an outlier;

[0022] Step 1.4: If the data point is considered an outlier, replace it with the median of this window.

[0023] Preferably, the singular spectrum decomposition process in Step 2 is as follows:

[0024] Step 2.1: Construct a trajectory matrix. Cut the original time series into multiple vectors according to a preset window length, and then form a matrix with these vectors arranged in rows, which is called a trajectory matrix. Let the time series be [x1, x2, x3, …, x N , then the obtained trajectory matrix is:

[0025]

[0026] Step 2.2: Singular value decomposition. Perform singular value decomposition on the trajectory matrix to obtain its singular values and singular vectors. Arrange the singular values in descending order as, where d is the rank of matrix X. The singular value decomposition result is:

[0027] X = UΛV T

[0028] where Λ = diag{λ1, λ2, …, λ d , 0, …, 0}, U is the eigenvector matrix of XX T , and U i is the i-th column vector of U, called the left singular vector; V is the eigenvector matrix of X T X, and V i is the i-th column vector of V, called the right singular vector. Then X can be written as:

[0029] X = X1 + X2 + … + X d

[0030]

[0031] where λ i is the i-th singular value of the singular spectrum decomposition, and X i is the grouped component matrix;

[0032] Step 2.3: Grouping. Divide the index set {1, …, d} into m non-overlapping subsets I1, …, I m , and let I = {i1, …, i p}, where p is the number of matrices in a group. Then the composite matrix corresponding to I There is:

[0033]

[0034] Step 2.4: Reconstruction. For each of the divided matrices , corresponding to different components of the signal. Then, by the method of diagonal averaging, each matrix is converted into a sequence of length N. Let Y be an M×K matrix with its elements being y ij . Let M * = min(M, K), K * = max(M, K), N = M + K - 1. If M < K, y * ij = y ij , otherwise y * ij = y ji; where M is the window length, K is the number of subsequences generated after window sliding, and N is the sequence length.

[0035] Preferably, the process of the K-means clustering method in step 2.3 is as follows:

[0036] Step 3.1: Randomly select K initial cluster centers, select k sample points in the sample data set D, and assign the k sample point values to the initial cluster centers respectively

[0037] Step 3.2: For each point p in the data set t (t = 1,..., n), calculate its Euclidean distance d(t, i) from the k cluster centers and assign this point to the cluster center closest to it;

[0038]

[0039] Step 3.3: For each cluster, calculate the new cluster center;

[0040]

[0041] Step 3.4: Calculate the squared error E of all points in the data set i , and compare it with the previous error E i-1 ;

[0042]

[0043] If |E i+1 - E i | < δ, the algorithm ends, where δ is a preset threshold; otherwise, go to step 3.2 for another iteration;

[0044] Step 3.5: Output the final k clusters and their cluster centers.

[0045] Preferably, the process of using the method of joint matrix diagonalization in step 4 is as follows:

[0046] Step 4.1: Whiten the received data, convert the received data into a random variable with zero mean and variance of 1, so that the signal satisfies E[s(t)s(t) H = I, restoring the second-order independence between signals, E[] is the expectation operation; H is the conjugate transpose, and I is the identity matrix;

[0047] Assume that the covariance matrix of the observed data is R xx , then R xx = E[x(t)x(t) H = FΛF H , where Λ is a diagonal matrix, and the diagonal elements are the diagonal elements of Rxx The eigenvalues, F is an orthogonal matrix, and each column corresponds to the orthonormal eigenvector of the eigenvalue. Let the whitening matrix be Q, and the whitening matrix should satisfy:

[0048] R zz = E[z(t)z(t) H

[0049] = QE[x(t)x(t) H H

[0050] = QFΛF H Q H

[0051] = I

[0052] where s(t) is the original signal source to be separated, z(t) is the signal after whitening processing, and x(t) is the observed signal; R zz is the fourth-order cumulant matrix of the whitened signal;

[0053] Solving gives Q = Λ -1 / 2 F H , so after whitening, z(t) = Qx(t) = QAs(t), and the estimation of the mixing matrix A is transformed into the estimation problem of the unitary matrix U, and the unitary matrix U is an orthogonal matrix;

[0054] Step 4.2: Establish the fourth-order cumulant matrix;

[0055] Use the objective function to measure the independence between vectors. The JADE algorithm uses the fourth-order cumulant matrix as the objective function, and its definition is as follows:

[0056]

[0057] In the formula, cum(z p ,z q ,z k ,z l ) represents the fourth-order cumulant of the vector z, z p ,z q ,z k ,z l are signal components, and m kl is the (k, l)th element of any N×N-dimensional matrix M;

[0058] Step 4.3: Joint diagonalization of the eigenmatrix;

[0059] Obtain a unitary matrix U to get the estimation of the source signal y(t) = U H ​​The estimation of \(z(t)\) and the unitary matrix \(U\) is obtained by jointly diagonalizing the fourth-order cumulant matrix. The fourth-order cumulant matrix is eigen-decomposed to obtain:

[0060]

[0061] wherein, is the estimated value of the unitary matrix \(U\), \(\Lambda(M)\) represents the diagonal matrix of the fourth-order cumulant matrix, and at this time, the matrix \(M\) is the eigen-matrix of the fourth-order cumulant matrix;

[0062] Step 4.4: Source signal separation;

[0063] Multiply the estimated unitary matrix by the whitened data matrix to obtain the estimated source signal

[0064] According to the blind source separation system for single-channel extremely low-frequency signals provided by the present invention, it includes:

[0065] Module M1: Collect the single-channel mixed signal and set the parameters for separating the mixed signal, including the sample length and the number of mixed signals; preprocess the collected mixed signal to remove strong impulse noise and enhance the separability and stability of the signal;

[0066] Module M2: Perform singular spectrum analysis on the preprocessed mixed signal, construct the Hankel matrix of the signal, and perform singular spectrum decomposition on it to decompose the mixed signal into multiple independent signal components; record the eigenvectors and eigenvalues obtained by the decomposition, and construct the eigen-matrix of each signal component;

[0067] Module M3: Use the K-means clustering method to perform clustering analysis on the eigenvalues extracted from the multiple independent components obtained by singular spectrum decomposition to form virtual multi-channels;

[0068] Module M4: For the virtual multi-channels obtained by clustering, use the method of joint matrix diagonalization to restore the signals one by one; for low-frequency signals, use the eigen-matrix reconstruction method to perform diagonal averaging operation to achieve blind source separation;

[0069] Module M5: Evaluate the separation quality of the separated signals, calculate the similarity index of signal separation. If the signal separation effect does not meet the expectation, adjust the number of components decomposed by singular spectrum or the clustering parameters, and trigger Module M2; otherwise, output the multi-channel signal separation result.

[0070] Preferably, in the module M1, a Hampel filter is used to preprocess the mixed signal, including:

[0071] Module M1.1: For each point in the data sequence, calculate the median and standard deviation of this point within a sliding window of a fixed size;

[0072] Module M1.2: For each data point, calculate its deviation from the window median, which reflects the difference between the data point and its neighborhood;

[0073] Module M1.3: Set a threshold according to the standard deviation and the median. If the deviation of a data point from the median is greater than a preset threshold, then the data point is considered an outlier;

[0074] Module M1.4: If a data point is considered an outlier, replace it with the median of the window.

[0075] Preferably, the singular spectrum decomposition process in the module M2 is as follows:

[0076] Module M2.1: Construct a trajectory matrix. Cut the original time series into multiple vectors according to a preset window length, and then form these vectors into a matrix by rows, which is called a trajectory matrix. Let the time series be [x1, x2, x3, …, x N , then the obtained trajectory matrix is:

[0077]

[0078] Module M2.2: Singular value decomposition. Perform singular value decomposition on the trajectory matrix to obtain its singular values and singular vectors. Arrange the singular values in descending order as, d is the rank of the matrix X, and the singular value decomposition result is:

[0079] X = UΛV T

[0080] where, Λ = diag{λ1, λ2, …, λ d , 0, …, 0}, U is the eigenvector matrix of XX T and U i is the i-th column vector of U, called the left singular vector; V is the eigenvector matrix of X T X, and V i is the i-th column vector of V, called the right singular vector. Then X can be written as:

[0081] X = X1 + X2 + … + X d

[0082]

[0083] where, λ i is the i-th singular value of the singular spectrum decomposition, and X i is the grouped component matrix;

[0084] Module M2.3: Grouping. Divide the subscript set {1, …, d} into m non-overlapping subsets I1, …, I m, let \(I = \{i_1, \ldots, i p \}\), \(p\) be the number of matrices in a group, then the composite matrix corresponding to \(I\) There is:

[0085]

[0086] Module M2.4: Reconstruction, for each matrix divided corresponding to different components of the signal, and then by the method of diagonal averaging, each matrix is converted into a sequence of length \(N\). Let \(Y\) be an \(M\times K\) matrix, whose elements are \(y\) ij , let \(M\) * =\(\min(M, K)\), \(K\) * =\(\max(M, K)\), \(N = M + K - 1\). If \(M < K\), \(y\) * ij =\(y\) ij , otherwise \(y\) * ij =\(y\) ji ; where \(M\) is the window length, \(K\) is the number of subsequences generated after window sliding, and \(N\) is the sequence length.

[0087] Preferably, the process of the K-means clustering method in the module M2.3 is:

[0088] Module M3.1: Randomly select \(K\) initial cluster centers, select \(k\) sample points in the sample data set \(D\), and assign the values of the \(k\) sample points to the initial cluster centers

[0089] Module M3.2: For each point \(p\) in the data set t \((t = 1, \ldots, n)\), calculate its Euclidean distance \(d(t, i)\) from the \(k\) cluster centers , and assign this point to the cluster center closest to it;

[0090]

[0091] Module M3.3: For each cluster, calculate the new cluster center;

[0092]

[0093] Module M3.4: Calculate the sum of squared errors \(E\) of all points in the data set i , and compare it with the previous error \(E\) i-1 ;

[0094]

[0095] If \(|E\) i+1 - E\) iIf <δ, the algorithm ends, where δ is a preset threshold; otherwise, trigger module M3.2 for another iteration;

[0096] Module M3.5: Output the final k clusters and their cluster centers.

[0097] Preferably, the process of using the method of joint matrix diagonalization in the module M4 is as follows:

[0098] Module M4.1: Whiten the received data, convert the received data into a random variable with zero mean and variance of 1, so that the signal satisfies E[s(t)s(t) H = I, restoring the second-order independence between signals, E[] is the expectation operation; H is the conjugate transpose, and I is the identity matrix;

[0099] Assume that the covariance matrix of the observed data is R xx , then R xx = E[x(t)x(t) H = FΛF H , where Λ is a diagonal matrix, the diagonal elements are the eigenvalues of R xx , F is an orthogonal matrix, and each column corresponds to the orthonormal eigenvector of the eigenvalue. Let the whitening matrix be Q, and the whitening matrix should satisfy:

[0100] R zz = E[z(t)z(t) H

[0101] = QE[x(t)x(t) H Q H

[0102] = QFΛF H Q H

[0103] = I

[0104] where s(t) is the original signal source to be separated, z(t) is the signal after whitening, and x(t) is the observed signal; R zz is the fourth-order cumulant matrix of the whitened signal;

[0105] Solve to get Q = Λ -1 / 2 F H , so after whitening, z(t) = Qx(t) = QAs(t), and the estimation of the mixing matrix A is transformed into the estimation problem of the unitary matrix U, where the unitary matrix U is an orthogonal matrix;

[0106] Module M4.2: Establish the fourth-order cumulant matrix;

[0107] ​The independence between vectors is measured using the objective function. The JADE algorithm uses the fourth-order cumulant matrix as the objective function, and its definition is as follows:

[0108]

[0109] where cum(z p ,z q ,z k ,z l ) represents the fourth-order cumulant of the vector z, and z p ,z q ,z k ,z l are signal components, and m kl is the (k, l)-th element of an arbitrary N×N-dimensional matrix M;

[0110] Module M4.3: Joint diagonalization of the characteristic matrix;

[0111] An orthogonal matrix U is obtained, and the estimated source signal is y(t) = U H z(t). The estimation of the orthogonal matrix U is obtained by jointly diagonalizing the fourth-order cumulant matrix. The fourth-order cumulant matrix is eigen-decomposed to obtain:

[0112]

[0113] where is the estimated value of the orthogonal matrix U, and Λ(M) represents the diagonal matrix of the fourth-order cumulant matrix. At this time, the matrix M is the characteristic matrix of the fourth-order cumulant matrix;

[0114] Module M4.4: Source signal separation;

[0115] Multiply the estimated orthogonal matrix by the whitened data matrix to obtain the estimated source signal

[0116] Compared with the prior art, the present invention has the following beneficial effects:

[0117] The present invention decomposes a single-channel complex mixed signal into multiple independent components by introducing singular spectrum decomposition, providing a basis for subsequent signal separation. Secondly, before signal separation, the preprocessing step effectively removes strong impulse noise and other interfering signals, thus ensuring the robustness and reliability of subsequent separation. This preprocessing operation can effectively improve the quality of the separation result. In addition, the present invention combines the K-means clustering technique to reasonably cluster the decomposed signal components to construct a virtual multi-channel signal. This innovative introduction overcomes the limitations of traditional single-channel signal separation techniques and provides a more flexible and effective means for signal separation. At the same time, by using the joint approximate diagonalization method of the feature matrix, efficient source signal separation of the multi-channel signal is achieved, greatly improving the efficiency and accuracy of separation. Finally, the present invention has excellent recovery ability for low-frequency signals and is particularly suitable for application scenarios that require the use of low-frequency signals. BRIEF DESCRIPTION OF THE DRAWINGS

[0118] Other features, objects, and advantages of the present invention will become more apparent by reading the following detailed description of non-limiting embodiments with reference to the accompanying drawings:

[0119] Figure 1 It is a system block diagram of the blind source separation method for single-channel extremely low-frequency signals in the present invention;

[0120] Figure 2 It is a schematic diagram of the extremely low-frequency signal communication system in Embodiment 1 and Embodiment 2 of the present invention;

[0121] Figure 3 It is a time-domain diagram of the signal and interference in Embodiment 1 of the present invention;

[0122] Figure 4 It is the mixed signal with added impulse noise and the preprocessed mixed signal in Embodiment 1 of the present invention;

[0123] Figure 5 It is a time-domain diagram of the separated source signal and interference signal in Embodiment 1 of the present invention. DETAILED DESCRIPTION OF THE EMBODIMENTS

[0124] The present invention will be described in detail below with reference to specific embodiments. The following embodiments will help those skilled in the art to further understand the present invention, but do not limit the present invention in any form. It should be noted that those of ordinary skill in the art can make several changes and improvements without departing from the concept of the present invention. These all belong to the protection scope of the present invention.

[0125] Embodiment 1

[0126] As Figure 1, the present invention discloses a blind source separation method for single-channel extremely low frequency signals. Aiming at the problems of existing technologies, first, a Hampel filter is used to preprocess the collected signals to eliminate impulse noise. Secondly, the single-channel signal sequence is decomposed into several independent components by singular spectrum analysis and K-means clustering algorithm to construct a virtual multi-channel. Finally, the source signals with high signal-to-interference ratio are separated using the joint approximate diagonalization algorithm based on the feature matrix. The present invention improves the accuracy of source signal estimation while not assuming any probability distribution and not requiring manual parameter selection, and has good robustness and universality.

[0127] Refer to Figure 2 , the main application scenario of the present invention is in a long-distance wireless communication system. At the transmitting end, the source signal and the interference signal are linearly mixed and then enter the channel for transmission. Let the source signal vector be s(t) = [s1(t), s2(t),.... s N (t)] T , and the noise vector be n(t) = [n1(t), n2(t),... n M (t)] T , and the noise signal follows the Alpha-stable distribution. The observed signal can be obtained as x(t) = As(t) + n(t), where A is a mixing coefficient vector of size 1×N.

[0128] The method proposed by the present invention first uses a Hampel filter to preprocess the extremely low frequency collected signals to eliminate impulse noise, specifically including the following steps:

[0129] Step 1.1: For each point in the data sequence, calculate the median and standard deviation of this point within a fixed-size sliding window.

[0130] Step 1.2: For each data point, calculate its deviation from the window median. This deviation reflects the difference between the data point and its neighborhood.

[0131] Step 1.3: Set a threshold according to the standard deviation and the median. If the deviation of the data point from the median is greater than a preset threshold (usually 2 times the standard deviation), then this data point is considered an outlier.

[0132] Step 1.4: If the data point is considered an outlier, replace it with the median of this window. The replaced data will not be affected by a single outlier, thus protecting the overall trend of the data.

[0133] The preprocessed observed signal is decomposed into multiple independent components using singular spectrum decomposition to construct a virtual multi-channel. Each component represents an inherent characteristic of the data. These components can be trends, cycles, noises, etc., and they are sorted according to their importance. It includes the following steps:

[0134] Step 2.1: Construct the trajectory matrix. Cut the original time series into several vectors according to a certain window length, and then form these vectors into a matrix by rows, which is called the trajectory matrix. Suppose the preprocessed time series is [x1, x2, x3, …, x N , take the window length as M, and take K = N - M + 1, then the obtained trajectory matrix is:

[0135]

[0136] Step 2.2: Singular value decomposition. Perform singular value decomposition on the trajectory matrix to obtain its singular values and singular vectors. Arrange the singular values in descending order. Let d be the rank of matrix X, and the singular value decomposition result is:

[0137] X = UΛV T

[0138] where Λ = diag{λ1, λ2, …, λ d , 0, …, 0}, U is the eigenvector matrix of XX T , U i is the i-th column vector of U, called the left singular vector, V is the eigenvector matrix of X T X, V i is the i-th column vector of V, called the right singular vector, then X can be written as:

[0139] X = X1 + X2 + … + X d

[0140]

[0141] Step 2.3: Grouping. Divide the set {1, …, d} into m non-overlapping subsets I1, …, I m , let I = {i1, …, i p}, then the synthesis matrix corresponding to I then there is:

[0142]

[0143] In order to further construct virtual multi-channels, it is generally manually set by the user according to prior knowledge or specific rules (such as uniform division). In order to automatically and accurately divide components into different groups, the present invention proposes to automatically assign components through the K-means clustering algorithm, calculate the statistical characteristics of each component, and determine the grouping of components according to these characteristics. The specific K-means clustering includes the following steps:

[0144] Step 3.1: Randomly select k initial cluster centers. Select k sample points in the sample data set D, and assign the k sample point values to the initial cluster centers respectively

[0145] Step 3.2: For each point p in the dataset t (t = 1, ..., n), calculate its Euclidean distance d(t, i) from the k cluster centers and assign the point to the cluster center that is closest to it.

[0146]

[0147] Step 3.3: For each cluster, calculate the new cluster center (i.e., the mean of all points within the cluster).

[0148]

[0149] Step 3.4: Calculate the squared error E of all points in the dataset i , and compare it with the previous error E i-1 ;

[0150]

[0151] If |E i+1 - E i | < δ, the algorithm ends; otherwise, go back to Step 3.2 for another iteration.

[0152] Step 3.5: Output the final k clusters and their cluster centers.

[0153] Step 3.6: Reconstruction. For each matrix X divided in Step 2 IJ , corresponding to different components of the signal. Then, through the method of diagonal averaging, each matrix is converted into a sequence of length N. Let Y be an M×K matrix with its elements being y ij , let M * = min(M, K), K * = max(M, K), N = M + K - 1. If M < K, y * ij = y ij , otherwise y * ij = y ji .

[0154]

[0155] That is, through anti - diagonal averaging, each signal component matrix is converted into the corresponding time series y k . Then, according to the grouping rule, these time series are accumulated to obtain m decomposed time series components, each with a length of N.

[0156] The JADE algorithm mainly consists of several steps: preprocessing the received data, establishing the fourth-order cumulant matrix, jointly diagonalizing the feature matrices, and separating the source signals. The JADE algorithm is proposed based on the diagonalization of the fourth-order cumulant matrix and has good robustness. It can also well separate the source signals when the frequency differences of the signals are small.

[0157] Step 4.1: Mainly whiten the received data to transform the received data into random variables with zero mean and variance of 1, so that the signals satisfy E[s(t)s(t) H = I, and restore the second-order independence between the signals. Assume that the covariance matrix of the observed data is R xx , then R xx = E[x(t)x(t) H = FΛF H , where Λ is a diagonal matrix, and the diagonal elements are the eigenvalues of R xx , and F is an orthogonal matrix, and each column corresponds to the orthonormal eigenvector of the eigenvalue. Let the whitening matrix be Q, and the whitening matrix should satisfy:

[0158] R zz = E[z(t)z(t) H

[0159] = QE[x(t)x(t) H Q H

[0160] = QFΛF H Q H

[0161] = I

[0162] Solving gives Q = Λ -1 / 2 F H , so after whitening, z(t) = Qx(t) = QAs(t), and the estimation of the mixing matrix A can be transformed into the estimation problem of the unitary matrix U, and the unitary matrix U is an orthogonal matrix.

[0163] Step 4.2: Establish the fourth-order cumulant matrix;

[0164] Use the objective function to measure the independence between vectors. The JADE algorithm uses the fourth-order cumulant matrix as the objective function, and its definition is as follows:

[0165]

[0166] In the formula, cum(z p ,z q ,z k ,z l ) represents the fourth-order cumulant of the vector z, m​kl is the (k, l)-th element of an arbitrary N×N dimensional matrix M.

[0167] Step 4.3: Joint diagonalization of the characteristic matrix;

[0168] The purpose of the JADE algorithm is to obtain a unitary matrix U, and the estimated source signal is se(t) = U H The estimation of the unitary matrix U is obtained by jointly diagonalizing the fourth-order cumulant matrix. The eigenvalue decomposition of the fourth-order cumulant matrix gives:

[0169]

[0170] where is the estimated value of the unitary matrix U, Λ(M) represents the diagonal matrix of the fourth-order cumulant matrix, and at this time, the matrix M is the characteristic matrix of the fourth-order cumulant matrix.

[0171] Step 4.4: Source signal separation;

[0172] Multiply the estimated unitary matrix by the whitened data matrix, and the estimated source signal can be obtained

[0173] Example 2

[0174] The experimental conditions of the present invention are as follows: an extremely low frequency communication signal with a modulation method of MSK, a frequency of 10 Hz, a propagation rate of 1 bps, and a sampling frequency of 1000 Hz. Comb spectrum interference and broadband blocking interference are respectively added, and impulse noise with an Alpha stable distribution is added after linear mixing. Blind source separation is performed on the noisy mixed signal by the method provided by this invention. The source signal and interference signal before mixing are shown by Figure 3 , the time-domain waveform of the mixed signal after passing through the channel with impulse noise and Hampel filtering is shown by Figure 4 , and the time-domain waveforms of the separated source signal and interference signal are shown by Figure 5 . The specific numerical values of each of the above steps are shown below:

[0175] The observed signal in this embodiment is x(t) = As(t) + n(t), where A is a mixing coefficient vector of size 1×N

[0176] A = [0.164 0.469 0.389]

[0177] The three rows of the source signal s(t) respectively represent the specific data of the extremely low frequency communication signal, comb spectrum interference, and broadband blocking interference of MSK as follows:

[0178]

[0179] The specific data of the mixed observation signal part is as follows:

[0180] x(t) = [-0.666 0.547 -0.447 0.0723 0.119 -0.338 -0.405 -0.301 -1.031 -0.249 0.395 -3.787 -0.505....-0.499 0.747 0.167 0.127 -0.443 -0.895] 1×1000

[0181] The specific data of the observation signal preprocessed by the Hampel filter in Step 1 is as follows:

[0182] x Hampel (t) = [-0.4477 0.0723 0.1193 -0.3385 -0.4050 -0.3015 -0.3015.....0.7474 0.1671 0.1276 -0.4431 -0.8953] 1×1000

[0183] The specific data of constructing the trajectory matrix in Step 2.1 Singular Spectrum Analysis is as follows:

[0184]

[0185] The eigenvalues after singular value decomposition of the trajectory matrix in Step 2.2 are arranged from largest to smallest as:

[0186] λ i = [6999.5 6484.1 2892.1 2880.4 2381.5 2373.1 1096.5 1089.9...........7.4 7.3] 1×400

[0187] The clustering results obtained by using the K-means clustering algorithm to cluster the eigenvalues obtained in Step 2.2 in Step 3.1 are as follows:

[0188] idx = [2 2 3 3 3 3 1 1 1.............1 1 1] 1×400

[0189] Each eigenvalue will be labeled idx by the K-means clustering, and the same idx represents the same class.

[0190] In Step 3.4, for each matrix divided in Step 2 by the method of diagonal averaging The specific data of the reconstructed virtual multi-channel signal part is as follows:

[0191]

[0192] In step 4.1, the specific data of the received virtual multi-channel signal after whitening processing with the whitening matrix Q are as follows:

[0193]

[0194] In step 4.2, cum(z p ,z q ,z k ,z l ) represents the fourth-order cumulant of the vector z. The specific data are as follows: cum(z p ,z q ,z k ,z l ) = [-1.486 7.623e-05 -0.0221 7.623e-05 0.01813 -0.0026 -0.02219......0.5555 5.3704] 1×81

[0195] The estimation of the source signal in step 4.3 is For the estimated value of the unitary matrix U, its specific data are as follows:

[0196]

[0197] The specific data of the signal after final blind source separation are as follows:

[0198]

[0199] To better describe the similarity between the separated signal and the source signal, the similarity coefficient analysis of the separated signal and the source signal can be used as a measure of the separation performance, and its definition is as follows:

[0200]

[0201] This paper believes that when the similarity coefficient is greater than 0.8, the signal separation is successful. It can be seen from the similarity coefficient matrix that the separated signal has a high similarity with the source signal, and there is only one value close to 1 in each row and each column, indicating good separation performance.

[0202] From this, it is calculated that Figure 5 The similarity coefficient matrix of the separated signal and Figure 3 the source signal is:

[0203]

[0204] It can be seen from the similarity coefficient matrix that the similarity coefficients in the first row of the matrix are close to 1. The matrix indicates that the first signal of the separated signals is similar to the first source signal, which is the separated MSK signal. The separated signals are highly similar to the source signals, and the separation effect is good.

[0205] Those skilled in the art know that in addition to implementing the systems, devices, and their respective modules provided by the present invention in the form of pure computer-readable program codes, the method steps can be logically programmed to enable the systems, devices, and their respective modules provided by the present invention to be implemented in the form of logic gates, switches, application-specific integrated circuits, programmable logic controllers, embedded microcontrollers, etc. Therefore, the systems, devices, and their respective modules provided by the present invention can be regarded as a kind of hardware components, and the modules included therein for implementing various programs can also be regarded as the structures within the hardware components; the modules for implementing various functions can also be regarded as either software programs for implementing the methods or the structures within the hardware components.

[0206] The specific embodiments of the present invention have been described above. It should be understood that the present invention is not limited to the above specific implementation manners, and those skilled in the art can make various changes or modifications within the scope of the claims, which do not affect the essence of the present invention. Without conflict, the embodiments of the present application and the features in the embodiments can be combined with each other arbitrarily.

Claims

1. A blind source separation method for single-channel extremely low frequency signals, characterized in that Including: Step 1: Collect the single-channel mixed signal and set the parameters for mixed signal separation, including the sample length and the number of mixed signals; Preprocess the collected mixed signal to remove strong impulse noise and enhance the separability and stability of the signal; Step 2: Perform singular spectrum analysis on the preprocessed mixed signal, construct the Hankel matrix of the signal, and perform singular spectrum decomposition on it to decompose the mixed signal into multiple independent signal components; Record the eigenvectors and eigenvalues obtained from the decomposition, and construct the eigenmatrix of each signal component; Step 3: Use the K-means clustering method to perform clustering analysis on the eigenvalues extracted from the multiple independent components obtained by singular spectrum decomposition to form virtual multi-channels; Step 4: For the virtual multi-channels obtained by clustering, use the method of joint matrix diagonalization to restore the signals one by one; For low-frequency signals, use the eigenmatrix reconstruction method to perform diagonal averaging operation to achieve blind source separation; Step 5: Evaluate the separation quality of the separated signals, calculate the similarity index of signal separation. If the signal separation effect does not meet the expectation, adjust the number of components decomposed by singular spectrum or the clustering parameters, and return to Step 2; Otherwise, output the multi-channel signal separation result.

2. The blind source separation method for single-channel extremely low frequency signals according to claim 1, wherein, In Step 1, a Hampel filter is used to preprocess the mixed signal, including: Step 1.1: For each point in the data sequence, calculate the median and standard deviation of this point within a sliding window of a fixed size; Step 1.2: For each data point, calculate its deviation from the window median to reflect the difference between the data point and its neighborhood; Step 1.3: According to the standard deviation and median, set a threshold. If the deviation of the data point from the median is greater than a preset threshold, then this data point is considered an outlier; Step 1.4: If the data point is considered an outlier, then replace it with the median of this window.

3. The blind source separation method for single-channel extremely low frequency signals according to claim 2, characterized in that The process of singular spectrum decomposition in Step 2 is as follows: Step 2.1: Construct a trajectory matrix. Cut the original time series into multiple vectors according to a preset window length, and then form a matrix with these vectors in rows, which is called the trajectory matrix. Let the time series be [x1, x2, x3, …, x N , then the obtained trajectory matrix is: Step 2.2: Singular value decomposition. Perform singular value decomposition on the trajectory matrix to obtain its singular values and singular vectors. Arrange the singular values in descending order. d is the rank of matrix X, and the singular value decomposition result is: X = UΛV T where, Λ = diag{λ1, λ2, …, λ d , 0, …, 0}, U is the eigenvector matrix of XX T , Ui i is the i-th column vector of U, called the left singular vector; V is the eigenvector matrix of XX T , Vi i is the i-th column vector of V, called the right singular vector, then X can be written as: X = X1 + X2 + … + X d Among them, λ i is the i-th singular value of the singular spectrum decomposition, and X i is the component matrix after grouping; Step 2.3: Grouping. Divide the subscript set {1, …, d} into m non-overlapping subsets I1, …, I m , and let I = {i1, …, i p}, where p is the number of matrices in a group. Then the composite matrix corresponding to I There is: Step 2.4: Reconstruction. For each divided matrix corresponding to different components of the signal, and then by the method of diagonal averaging, each matrix is converted into a sequence of length N. Let Y be an M×K matrix, whose elements are y ij . Let M * = min(M, K), K * = max(M, K), N = M + K - 1. If M < K, y * ij = y ij , otherwise y * ij = y ji ; where M is the window length, K is the number of subsequences generated after window sliding, and N is the sequence length.

4. The blind source separation method for single-channel extremely low frequency signals according to claim 3, characterized in that The process of the K-means clustering method in Step 2.3 is as follows: Step 3.1: Randomly select K initial cluster centers, select k sample points from the sample data set D, and assign the k sample point values to the initial cluster centers respectively Step 3.2: For each point p in the dataset t (t = 1,..., n), calculate its Euclidean distance d(t, i) from the k cluster centers and assign the point to the cluster center that is closest to it; Step 3.3: For each cluster, calculate the new cluster center; Step 3.4: Calculate the squared error E of all points in the dataset i , and compare it with the previous error E i-1 . If |E i+1 - E i | < δ, the algorithm ends, where δ is a preset threshold value; Otherwise, transfer to Step 3.2 for another iteration; Step 3.5: Output the final k clusters and their cluster centers.

5. The blind source separation method for single-channel extremely low frequency signals according to claim 4, characterized in that The process of using the method of joint matrix diagonalization in Step 4 is as follows: Step 4.1: Whiten the received data to transform the received data into a random variable with zero mean and variance of 1, so that the signal satisfies E[s(t)s(t) H = I, restoring the second-order independence between signals, where E[] is the expected value operation; H is the conjugate transpose, and I is the identity matrix; Assume that the covariance matrix of the observed data is R xx , then R xx = E[x(t)x(t) H = FΛF H , where Λ is a diagonal matrix, and the diagonal elements are the eigenvalues of R xx , and F is an orthogonal matrix, and each column corresponds to the orthonormal eigenvector of the eigenvalue. Let the whitening matrix be Q, and the whitening matrix should satisfy: R zz = E[z(t)z(t) H ​ = QE[x(t)x(t) H Q H = QFΛF H Q H ＝I Among them, s(t) is the original signal source to be separated, z(t) is the signal after whitening processing, and x(t) is the observed signal; R zz is the fourth-order cumulant matrix of the whitened signal; The solution gives Q = Λ -1 / 2 F H , so after whitening, we get z(t) = Qx(t) = QAs(t), which transforms the estimation of the mixing matrix A into the problem of estimating the unitary matrix U. The unitary matrix U is an orthogonal matrix; Step 4.2: Establish a fourth-order cumulant matrix; Use the objective function to measure the independence between vectors. The JADE algorithm uses the fourth-order cumulant matrix as the objective function, and its definition is as follows: where cum(z p ,z q ,z k ,z l ) represents the fourth-order cumulant of the vector z, z p ,z q ,z k ,z l are signal components, and m kl is the (k, l)-th element of an arbitrary N×N-dimensional matrix M; Step 4.3: Joint diagonalization of the eigenmatrix; An unitary matrix U is obtained, and the estimated source signal y(t) = U H z(t). The estimation of the unitary matrix U is obtained by jointly diagonalizing the fourth-order cumulant matrix. By performing eigenvalue decomposition on the fourth-order cumulant matrix, we get: In the formula, is the estimated value of the unitary matrix U, and Λ(M) represents the diagonal matrix of the fourth-order cumulant matrix. At this time, the matrix M is the eigenmatrix of the fourth-order cumulant matrix; Step 4.4: Source signal separation; Multiply the estimated unitary matrix by the whitened data matrix to obtain the source signal estimate 6. A blind source separation system for single-channel extremely low frequency signals, characterized in that, Including: Module M1: Collect the single-channel mixed signal and set the parameters for mixed signal separation, including the sample length and the number of mixed signals; Preprocess the collected mixed signal to remove strong impulse noise and enhance the separability and stability of the signal; Module M2: Perform singular spectrum analysis on the preprocessed mixed signal, construct the Hankel matrix of the signal, and perform singular spectrum decomposition on it to decompose the mixed signal into multiple independent signal components; record the eigenvectors and eigenvalues obtained from the decomposition, and construct the eigenmatrix of each signal component; Module M3: Use the K-means clustering method to perform clustering analysis on the eigenvalues extracted from multiple independent components of the singular spectrum decomposition to form virtual multi-channels; Module M4: For the virtual multi-channels obtained by clustering, use the method of joint matrix diagonalization to restore the signals one by one; For low-frequency signals, use the eigenmatrix reconstruction method to perform diagonal averaging operation to achieve blind source separation; Module M5: Evaluate the separation quality of the separated signals, calculate the similarity index of signal separation. If the signal separation effect does not meet the expectation, adjust the number of components decomposed by the singular spectrum or the clustering parameters, and trigger Module M2; Otherwise, output the multi-channel signal separation result.

7. The blind source separation system for single-channel extremely low frequency signals according to claim 6, characterized in that, In the said Module M1, a Hampel filter is used to preprocess the mixed signal, including: Module M1.1: For each point in the data sequence, calculate the median and standard deviation of this point within a sliding window of a fixed size; Module M1.2: For each data point, calculate its deviation from the window median to reflect the difference between the data point and its neighborhood; Module M1.3: According to the standard deviation and median, set a threshold. If the deviation of the data point from the median is greater than a preset threshold, then this data point is considered an outlier; Module M1.4: If the data point is considered an outlier, then replace it with the median of this window.

8. The blind source separation system for single-channel extremely low frequency signals according to claim 7, characterized in that, The singular spectrum decomposition process in the said Module M2 is as follows: Module M2.1: Construct a trajectory matrix. Cut the original time series into multiple vectors according to a preset window length, and then form a matrix with these vectors as rows, which is called the trajectory matrix. Let the time series be [x1, x2, x3, …, x N , then the obtained trajectory matrix is: Module M2.2: Singular value decomposition. Perform singular value decomposition on the trajectory matrix to obtain its singular values and singular vectors. Arrange the singular values in descending order. Let \(d\) be the rank of matrix \(X\), and the singular value decomposition result is: X = UΛV T where, Λ = diag{λ1, λ2, …, λ d , 0, …, 0}, U is the eigenvector matrix of XX T , Ui i is the i-th column vector of U, called the left singular vector; V is the eigenvector matrix of XX T , Vi i is the i-th column vector of V, called the right singular vector, then X can be written as: X = X1 + X2 + … + X d Among them, λ i is the i-th singular value of the singular spectrum decomposition, and X i is the component matrix after grouping; Module M2.3: Grouping. Partition the subscript set {1, …, d} into m non - overlapping subsets I1, …, I m , and let I = {i1, …, i p}, where p is the number of matrices in a group. Then the composite matrix corresponding to I There is: Module M2.4: Reconstruction, each matrix Corresponding to different components of the signal, each matrix is averaged by diagonal line. Convert to a sequence of length N, let Y be an M×K matrix whose elements are y ij , let M * =min(M,K),K * =max(M,K),N=M+K-1, if M <K,y * ij =y ij , otherwise y * ij =y ji ; Where M is the window length, K is the number of subsequences generated after the window slides, and N is the sequence length.

9. The blind source separation system for single-channel extremely low frequency signals according to claim 8, wherein The process of the K-means clustering method in the said Module M2.3 is as follows: Module M3.1: Randomly select K initial cluster centers, select k sample points in the sample data set D, and assign the values of the k sample points to the initial cluster centers respectively Module M3.2: For each point p in the dataset t (t = 1,..., n), calculate its Euclidean distance d(t, i) from k cluster centers and assign the point to the cluster center that is closest to it; Module M3.3: For each cluster, calculate the new cluster center; Module M3.4: Calculate the squared error E of all points in the dataset i , and compare it with the previous error E i-1 ; If |E i+1 - E i | < δ, the algorithm ends, where δ is a preset threshold; otherwise, trigger module M3.2 for another iteration; Module M3.5: Output the final \(k\) clusters and their cluster centers.

10. The blind source separation system for single-channel extremely low frequency signals according to claim 9, wherein The process of using the method of joint matrix diagonalization in the said Module M4 is as follows: Module M4.1: Whiten the received data, transform the received data into a random variable with zero mean and variance of 1, so that the signal satisfies E[s(t)s(t) H = I, restore the second-order independence between signals, E[] is the expected value operation; H is the conjugate transpose, and I is the identity matrix; Assume that the covariance matrix of the observed data is R xx , then R xx = E[x(t)x(t) H = FΛF H , where Λ is a diagonal matrix, and the diagonal elements are the eigenvalues of R xx , F is an orthogonal matrix, and each column corresponds to the orthonormal eigenvector of the eigenvalue. Let the whitening matrix be Q, and the whitening matrix should satisfy: R zz = E[z(t)z(t) H ​ = QE[x(t)x(t) H Q H = QFΛF H Q H ＝I Among them, s(t) is the original signal source to be separated, z(t) is the signal after whitening processing, and x(t) is the observed signal; R zz is the fourth-order cumulant matrix of the whitened signal; The solution gives Q = Λ -1 / 2 F H , so after whitening, we get z(t) = Qx(t) = QAs(t). The estimation of the mixing matrix A is transformed into the problem of estimating the unitary matrix U, and the unitary matrix U is an orthogonal matrix; Module M4.2: Establish a fourth-order cumulant matrix; Use the objective function to measure the independence between vectors. The JADE algorithm uses the fourth-order cumulant matrix as the objective function, and its definition is as follows: where cum(z p ,z q ,z k ,z l ) represents the fourth-order cumulant of the vector z, z p ,z q ,z k ,z l are signal components, and m kl is the (k, l)-th element of an arbitrary N×N-dimensional matrix M; Module M4.3: Joint diagonalization of the eigenmatrix; Find a unitary matrix U to obtain the estimated source signal y(t) = U H z(t). The estimation of the unitary matrix U is obtained by jointly diagonalizing the fourth-order cumulant matrix. Perform eigenvalue decomposition on the fourth-order cumulant matrix to get: In the formula, is the estimated value of the unitary matrix U, Λ(M) represents the diagonal matrix of the fourth-order cumulant matrix, and the matrix M is the characteristic matrix of the fourth-order cumulant matrix; Module M4.4: Source signal separation; Multiply the estimated unitary matrix by the whitened data matrix to obtain the source signal estimate

Images