Biomedical signal separation method based on fast multichannel non-negative matrix factorization

Abstract

The invention provides a biomedical signal separation method based on fast multi-channel non-negative matrix factorization, and relates to the technical field of biomedical signal processing. According to the method, useful signals and noise can be effectively separated based on the NMF separation method, and the noise and the useful signals can be separated by NMF through extracting base components of different sources from original signals, so that de-noising and source separation of the signals are realized, and the signal quality is improved. According to the method, the NMF decomposition principle, multichannel data processing and rapid optimization algorithm application are combined, so that efficient separation, denoising and source extraction of biomedical signals are achieved, the processing efficiency is improved, and innovation is made in the aspects of signal quality improvement and application precision. And due to the characteristic of top-speed calculation, the method can be applied to a large-scale data analysis scene and has a wide application prospect.

Description

Technical Field

[0001] The present invention relates to the technical field of biomedical signal processing, and in particular to a biomedical signal separation method based on fast multi-channel non-negative matrix decomposition. Background Art

[0002] In the field of biomedical signal processing, especially in the monitoring and analysis of electrocardiogram (ECG), electroencephalogram (EEG), electromyogram (EMG) and other signals, complex signal mixing problems are often faced. These signals are usually interfered by noise, artifacts and different physiological processes, making signal separation an important research task. Although traditional signal separation methods such as independent component analysis (ICA) and principal component analysis (PCA) are widely used, they often assume a linear relationship between signals and noise or between channel signals, and have many limitations, such as high computational complexity, poor robustness, and dependence on the number of signal sources. Especially in the face of complex situations of nonlinear, multi-noise, and multi-source signals, their effects are often unsatisfactory. In order to improve the efficiency and accuracy of biomedical signal separation, the non-negative matrix factorization (NMF) method has gradually attracted attention in recent years, especially in the application of multi-channel signal separation.

[0003] Signal separation methods based on NMF have many advantages, such as being able to handle nonlinear relationships, having good sparse representation capabilities, and having good suppression of noise and redundant information. With the advancement of computer technology, fast multi-channel NMF algorithms for biomedical signals have emerged. It can effectively extract useful information from complex signals in a short time, and has higher processing efficiency and lower computational burden than traditional methods. Separation methods based on emerging technologies such as NMF have become a hot topic in current research. These methods can make full use of the non-negative characteristics and sparsity of signals, and while improving computational efficiency, they can adapt to complex signal source interactions and noise environments, thus having greater potential in biomedical signal separation.

[0004] Although the traditional NMF method can effectively decompose signals, it has high computational complexity in the application of multi-channel signals, especially for large-scale data sets. In order to overcome this problem, the fast multi-channel NMF method was proposed, which solves the challenges in multi-channel signal separation by accelerating the optimization algorithm, reducing the amount of calculation and improving the decomposition efficiency. Summary of the invention

[0005] The technical problem to be solved by the present invention is to address the deficiencies of the above-mentioned prior art and to provide a biomedical signal separation method based on fast multi-channel non-negative matrix decomposition, which is more efficient, accurate and has strong robustness in separating signals, can improve the accuracy and efficiency of biomedical signal separation, solve the interference problems of noise and artifacts, and enhance the robustness of the algorithm.

[0006] In order to solve the above technical problems, the technical solution adopted by the present invention is:

[0007] A biomedical signal separation method based on fast multi-channel non-negative matrix decomposition includes data preprocessing and a formal separation process; the data preprocessing includes wavelet denoising on an input observed mixed biomedical signal; the formal separation process includes inputting the preprocessed mixed biomedical signal, converting the signal from a time domain image to a time-frequency domain image by short-time Fourier transform, fast multi-channel non-negative matrix decomposition (Fast-MNMF) algorithm separation, and Wiener filter separation; the Fast-MNMF algorithm separation includes assuming signal distribution, mixed spectrum covariance matrix modeling, power spectral density non-negative matrix decomposition, space covariance matrix joint diagonalization, mixed spectrum modeling, maximum likelihood estimation, and iterative optimization of parameters using the Fast-MNMF algorithm; the parameters obtained by iterative optimization of the Fast-MNMF algorithm are substituted into the formula of the Wiener filter to separate the biomedical signal.

[0008] Furthermore, the wavelet denoising adopts a wavelet transform method based on an improved threshold, combines a soft threshold function with a hard threshold function, and performs denoising processing; according to the change of noise with the wavelet transform, the mean square error σ of the wavelet system in each layer is estimated layer by layer i , calculate the threshold of each layer The threshold is trimmed by ln (1+i) to obtain the threshold at different scales; the specific formula of wavelet transform for improving the threshold is as follows:

[0009]

[0010] Among them, d i,j is the wavelet coefficient of each layer, is the wavelet coefficient adjusted by threshold, σ is the mean square error of noise, σ i is the mean square error of each layer, i is the number of layers of wavelet decomposition, j is the length of the detail signal, μ is the universal threshold estimate, and μ i is the threshold estimation for each layer; L is the length of the original signal, a is the adjustment factor, when a takes different values, the wavelet transform formula of the improved threshold changes between soft and hard functions, and has the advantages of both soft and hard functions;

[0011] In the short-time Fourier transform, the parameters are set as a window length of 1024, a frame shift of 50% of the window length, and a FFT point number of 1024, and a time-frequency energy diagram of a three-channel mixture is generated.

[0012] Furthermore, in the hypothetical signal distribution, it is assumed that the mixed spectrum of the short-time Fourier transform obeys the centralized multivariate complex Gaussian distribution x~N c(∑), that is, the distribution of the signal has nothing to do with the mean, but only with the covariance. Assuming that the observed time-frequency energy diagram is the superposition of the three source signal spectra, since the spectrum of each source obeys the complex Gaussian distribution, the following formula is obtained:

[0013]

[0014] Where p(x) is the probability density function of the time-frequency point x in the energy map obtained after short-time Fourier transform, n is the signal source, N is the number of signal sources, f is the frequency, and t is the time frame; π is a mathematical constant, which is 3.14; M is the number of signal channels, and one channel is equivalent to one microphone; H represents the conjugate transpose, Σ is the covariance of the centered multivariate complex Gaussian distribution, and x nft For each source, the energy diagram, x ft is the energy diagram of the mixed signal.

[0015] Furthermore, in the mixed spectrum covariance matrix modeling, it is assumed that the covariance matrix of the probability density function of the mixed spectrum is equal to the product of the signal power spectrum density and the spatial covariance matrix. According to the additivity of the spectrum and the additivity of the Gaussian distribution, the expression of the mixed spectrum is obtained as follows:

[0016]

[0017] Among them, λ nft is the power spectral density of the low-rank source signal, G nf is the spatial covariance matrix of the source signal, Y ft is the covariance matrix of the multivariate complex Gaussian distribution, N c It is the abbreviation of multivariate complex Gaussian distribution.

[0018] Furthermore, in the power spectral density non-negative matrix decomposition, the power spectral density λ of the low-rank source signal nft Perform non-negative matrix decomposition, the expression is as follows:

[0019]

[0020] Where K is the number of bases, w nkf and h nkt represent the amplitude and activation degree of the basis, respectively.

[0021] Furthermore, in the joint diagonalization of the spatial covariance matrix, the weight-sharing joint diagonalization spatial covariance matrix G is adopted. nf ,Will Defined as the non-negative weight vector of the original signal, the expression is as follows:

[0022]

[0023] Among them, G nfCovariance matrix; Q f is the unmixing matrix, with a size of M×M×F, where M is the number of channels and F is the number of frequencies. Each F is a unitary matrix; H represents conjugate transpose; Q f Initialize to the identity matrix; q f1 ,…,q fM Indicates the value of different channels at the same frequency; g nfm is an element in the weight matrix g, the size of g is N×M×F, N is the number of signal sources, It represents the weights of different channels under the same source and the same frequency. Since the weights under different frequencies are shared, that is, the weights under different frequencies are the same, the weights are only related to the source and the number of channels, so g nfm Equivalent to g nm , g nm is an element in the weight matrix g, and the size of g matrix is ​​N×M.

[0024] Furthermore, in the mixed graph modeling, the mixed graph z is obtained ft The probability density function is a special form of non-negative tensor decomposition. Assume that z ft The elements of x are independent of each other. ft Then the relevant expression is as follows:

[0025]

[0026] Among them, Q f x ft The product of is defined as z ft , called the projection matrix; according to the formula derived above, we get z ft Variance of probability density

[0027] Furthermore, in the maximum likelihood estimation, for z ft The multivariate complex Gaussian probability density function is used for maximum likelihood estimation, and the cost function of the fast multi-channel non-negative matrix decomposition algorithm is derived. The formula is as follows:

[0028]

[0029] Among them, n is the signal source, k is the basis, t is the time frame, f is the frequency, m is a channel; W is the basis matrix, the size is N×K×F, w nkf It is equivalent to an element of the matrix W under a certain source n, a certain basis k, and a certain frequency f; H is the activation matrix, with a size of N×K×T, h nkf It is equivalent to an element of the matrix H under a certain source n, a certain basis k, and a certain time t; w nkf and h nkf The power spectral density λ is nkfIt is obtained by non-negative matrix decomposition; g is the weight matrix, the size is N×M, g nm It is equivalent to an element of matrix g under a certain source n and a certain channel m; Q is the unmixing matrix under all frequencies f, q fm It is equivalent to an element of the matrix Q at a certain frequency f and a certain channel m; is the energy of the actual signal source, To estimate the energy of the signal source, when In the iterative process, When the maximum likelihood estimate approaches a certain value, it means that the iteration has converged. ft It is equivalent to an element of the matrix Z at a certain frequency f and a certain time t. In the fast multi-channel non-negative matrix decomposition algorithm, the input is a multi-channel signal. It is equivalent to an element of the matrix Z at a certain frequency f, a certain time t, and a certain channel m;

[0030] The fast multi-channel non-negative matrix factorization (Fast-MNMF) algorithm is used to iteratively optimize the parameters W, H, g, and Q.

[0031] Furthermore, the specific method of the Fast-MNMF algorithm for iteratively optimizing the parameters W, H, g, and Q is as follows:

[0032] Step 8.1: Initialization; Initialize W and H to random Gaussian distribution, set K = 2, that is, the basis of non-negative matrix decomposition is 2, Q f is the identity matrix, Initialize according to cyclic redundancy. The cyclic redundancy formula is as follows:

[0033] ∈=10 -2 ;

[0034] Since W, H, g, and Q have all been initialized, the parameters The initialization formula is as follows:

[0035]

[0036] Step 8.2: Normalization; Normalize the W, H, g, and Q parameters, and substitute the normalized matrix into the following iterative formula. The normalization update order is as follows:

[0037]

[0038] Step 8.3: Parameter iteration; update the values ​​of W, H, and g respectively, and adjust them according to the cost function. and Adjust and update the iteration formula as follows:

[0039]

[0040] Step 8.4: Update the unmixing matrix Q; adjust Q by optimizing the diagonalization operation f , ensure Q f The stability and validity of the value, e m is a one-hot vector, and the relevant formula is as follows:

[0041]

[0042] q fm =(Q f V fm ) -1 e m ,

[0043] Among them, X ft is x ft The covariance matrix of fm A row for Q.

[0044] Furthermore, in the Wiener filter separation, the finally obtained W, H, g, Q parameters are substituted into the Wiener filter formula to separate the biomedical signal. The specific method is:

[0045] Step 9.1: Calculate the Wiener filter matrix W as follows nft , used to estimate the time-frequency spectrum of each source signal:

[0046]

[0047] Among them, λ nft is the power spectral density of the source signal, is the directional weight of the source signal, represents the weighted sum of all source signals;

[0048] Step 9.2: Signal separation; using the Wiener filter matrix W nft For the input mixed graph X ft Separate and obtain the time-frequency spectrum X of each source signal nft :

[0049] X nft =W nft X ft ;

[0050] Step 9.3: Time domain signal reconstruction: for each separated source signal X nft Perform short-time inverse Fourier transform to reconstruct the corresponding time domain signal, thereby obtaining a pure biomedical signal.

[0051] The beneficial effect of adopting the above technical solution is that the biomedical signal separation method based on fast multi-channel non-negative matrix decomposition provided by the present invention can effectively separate useful signals from noise based on the NMF separation method. NMF can separate noise and useful signals by extracting basis components from different sources from the original signal, thereby realizing signal denoising and source separation, and improving signal quality. The key to the fast multi-channel NMF method is to improve the computational efficiency, and accelerate the iterative process by jointly diagonalizing the spatial covariance matrix. The present invention combines the decomposition principle of NMF, the processing of multi-channel data, and the application of fast optimization algorithms to achieve efficient separation, denoising and source extraction of biomedical signals, which not only improves the processing efficiency, but also makes innovations in signal quality improvement and application accuracy. The method of the present invention performs well in terms of signal separation accuracy, computational efficiency, noise suppression, application flexibility, etc. It not only improves the accuracy of signal separation, but also can make full use of spatial and temporal information in multi-channel signal processing, significantly improving the separation effect, and the characteristics of extremely fast calculation enable the method to be applied to large-scale data analysis scenarios, with broad application prospects, especially in the fields of biomedical detection, management, disease diagnosis, etc. BRIEF DESCRIPTION OF THE DRAWINGS

[0052] Figure 1 A flow chart of generating and preprocessing a mixed signal provided by an embodiment of the present invention;

[0053] Figure 2 A flow chart of biomedical signal separation based on fast multi-channel non-negative matrix decomposition provided in an embodiment of the present invention. DETAILED DESCRIPTION

[0054] The specific implementation of the present invention is further described in detail below in conjunction with the accompanying drawings and examples. The following examples are used to illustrate the present invention, but are not intended to limit the scope of the present invention.

[0055] A biomedical signal separation method based on fast multi-channel non-negative matrix factorization aims to improve the separation accuracy and efficiency of biomedical signals by optimizing the traditional non-negative matrix factorization algorithm, while enhancing the ability to suppress noise interference. It includes data preprocessing steps and formal separation process.

[0056] In this embodiment, appropriate biomedical signal data is first selected, such as electroencephalogram (EEG), electrocardiogram (ECG), and electromyogram (EMG), and then the respective signals are mixed to generate a mixed signal, and the wavelet with improved threshold is used for denoising and the Fast-MNMF algorithm is used for separation, so as to extract pure biomedical signals for subsequent analysis and diagnosis. The preprocessing process of this embodiment is as follows Figure 1 shown.

[0057] First, generate biomedical signals, download clean electroencephalogram (EEG), electrocardiogram (ECG), and electromyogram (EMG) from the PhysioNet website and import them into MATLAB for processing. These biomedical signals are usually affected by electromyographic signals (EMG), baseline drift, and power frequency interference. Electromyographic signals are low-frequency to medium-frequency signals caused by muscle contraction, and the frequency range is usually between 50-100Hz; baseline drift is a low-frequency fluctuation caused by factors such as breathing and body movement, with a frequency below 1Hz, which affects the stability of the signal; power frequency interference is caused by a 50Hz power supply frequency, which is common in ECG signals and will significantly reduce the clarity of the signal.

[0058] Then, a spatial structure is randomly generated. In order to simulate the above interference, a spatial structure is randomly generated and simulated by the following method: the electromyographic signal is generated using the Hodgkin-Huxley model, which is based on biophysical mechanisms and can generate real electromyographic signals with frequencies ranging from low to medium frequencies and small amplitudes; the baseline drift can be represented by a low-frequency sine wave or polynomial function to restore the baseline fluctuation in the actual signal as much as possible; the power frequency interference is generated by Gaussian white noise and directly superimposed on the biomedical signal, electromyographic signal and baseline drift. In this randomly generated spatial structure, the number of sensors is set to 3 and the number of sound sources is also set to 3, representing the biomedical signal, electromyographic signal and baseline drift respectively.

[0059] Finally, a mixed signal is generated. Using the generated spatial structure, three signals with power frequency interference are recorded through sensors to obtain three mixed signals.

[0060] For the above mixed signal, data preprocessing is performed, and denoising is performed using a wavelet transform method based on an improved threshold. This embodiment combines the advantages of the soft threshold function and the hard threshold function by adjusting between the two functions to more efficiently filter out power frequency interference and other noises while retaining the key features of the signal.

[0061] The essence of wavelet denoising is to transform the noisy signal into the wavelet domain and process the corresponding wavelet coefficients. The general threshold wavelet denoising algorithm can be divided into the following steps:

[0062] Step 1: Determine the wavelet transform, select the appropriate wavelet basis and decomposition layer number, perform wavelet transform on the noisy signal to obtain the wavelet coefficients d of each layer i,j ;

[0063] Step 2: Estimate the threshold μ according to the corresponding threshold rule, and use the threshold function to calculate the wavelet coefficient d i,j Perform threshold processing to obtain the wavelet coefficients after threshold processing

[0064] Step 3: Using the processed wavelet coefficients The signal is reconstructed with the approximate coefficients to obtain the denoised signal.

[0065] Common threshold functions include hard threshold function and soft threshold function. The formula is as follows:

[0066] The hard threshold function is

[0067] The soft threshold function is

[0068] The hard threshold function simply removes and retains the wavelet coefficients, and it is easy to filter out the mutation point information while denoising. The soft threshold function smoothes the wavelet coefficients whose absolute values ​​are greater than the threshold according to a certain rule, and achieves the denoising effect by weakening the wavelet coefficients. However, this method reduces the wavelet coefficients with large absolute values, resulting in the loss of high-frequency information.

[0069] In general, in the input noisy signal, the wavelet coefficients of the noise are concentrated in the detail coefficients of the first layer, and they decrease as the wavelet scale increases. In conventional methods, the threshold μ of the general wavelet denoising uses the noise variance of the first layer of wavelet coefficients as the overall noise variance. The threshold calculated from this method is constant and does not conform to the law of noise distribution with wavelet transform. Therefore, the denoising algorithm in this paper estimates the mean square error σ of the wavelet system in each layer layer by layer according to the change of noise with wavelet transform. i , calculate the threshold of each layer The threshold is trimmed by ln(1+i) to obtain the threshold at different scales. The specific formula of wavelet transform for improving the threshold is as follows:

[0070]

[0071] Among them, d i,j is the wavelet coefficient of each layer, is the wavelet coefficient adjusted by threshold, σ is the mean square error of noise, σ i is the mean square error of each layer, i is the number of layers of wavelet decomposition, j is the length of the detail signal, μ is the universal threshold estimate, and μ i is the threshold estimation for each layer, L is the length of the original signal, a is the adjustment factor, and when a takes different values, the function changes between soft and hard functions and has the advantages of both soft and hard functions.

[0072] This method can dynamically adjust the threshold according to the signal characteristics, improve the adaptability and effectiveness of signal processing, and extract purer biomedical signals through wavelet denoising, providing a reliable data basis for subsequent analysis and application.

[0073] This is followed by the formal separation step, e.g. Figure 2 As shown, the details are as follows.

[0074] Step 1: Input signal. Input the mixed signal observed by the three channels (biomedical signal, electromyographic signal, baseline drift signal) into the system as the basic data for subsequent processing.

[0075] Step 2: Short-time Fourier transform. Perform short-time Fourier transform (STFT) on the signals of the three channels respectively, with the parameters set to a window length of 1024, a frame shift of 50% of the window length, and an FFT point number of 1024 to generate a time-frequency energy map of the three channel mixture. Since the sensor signal is more in line with the convolution model (time domain convolution is equal to frequency domain product), the signal is first converted to a time-frequency domain image for subsequent operations. The time-frequency domain image is an energy map with time as the horizontal axis and frequency as the vertical axis.

[0076] Step 3: Assume signal distribution. Assume that the mixed spectrum of the short-time Fourier transform obeys the centralized multivariate complex Gaussian distribution x~N c (Σ), that is, the distribution of the signal has nothing to do with the mean, but only with the covariance. Assuming that the observed time-frequency energy graph is the superposition of the three source signal spectra, since the spectrum of each source obeys the complex Gaussian distribution, the following formula can be obtained:

[0077]

[0078] Where p(x) is the probability density function of the time-frequency point x in the energy map obtained after short-time Fourier transform, n is the signal source, N is the number of signal sources, f is the frequency, and t is the time frame; π is a mathematical constant, which is 3.14; M is the number of signal channels (one channel is equivalent to one microphone); H represents conjugate transposition. Since each point on the energy map obtained after the signal is short-time Fourier transform is a complex number, the method of transposition followed by conjugation is adopted instead of direct transposition; ∑ is the covariance of the centralized multivariate complex Gaussian distribution, x nft For each source, the energy diagram, x ft is the energy diagram of the mixed signal.

[0079] Step 4: Modeling the covariance matrix of the mixed spectrum. Assuming that the covariance matrix of the probability density function of the mixed spectrum is equal to the product of the signal power spectrum density (PSD) and the spatial covariance matrix (SCM), according to the additivity of the spectrum, the expression of the mixed spectrum can be obtained:

[0080]

[0081] Among them, λ nft is the power spectral density of the low-rank source signal, G nf is the spatial covariance matrix of the source signal, Y ft is the covariance matrix of the multivariate complex Gaussian distribution; N c It is the abbreviation of multivariate complex Gaussian distribution.

[0082] Step 5: Power spectral density non-negative matrix decomposition. The power spectral density λ of the low-rank source signal nft Perform non-negative matrix factorization (NMF), the expression is as follows:

[0083]

[0084] Where K is the number of bases, w nkf and h nkt represent the amplitude and activation degree of the basis, respectively.

[0085] Step 6: Joint diagonalization of the spatial covariance matrix. The joint diagonalization of the spatial covariance matrix G with weight sharing is adopted nf ,Will Defined as the non-negative weight vector of the original signal, the expression is as follows:

[0086]

[0087] Among them, G nf Covariance matrix; Q f is the unmixing matrix, with a size of M×M×F, where M is the number of channels and F is the number of frequencies. Each F is a unitary matrix; H represents conjugate transpose; Q f Initialize to the identity matrix; q f1 ,…,q fM Indicates the value of different channels at the same frequency; g nfm is an element in the weight matrix g, the size of g is N×M×F, N is the number of signal sources, It represents the weights of different channels under the same source and the same frequency. Since the weights under different frequencies are shared, that is, the weights under different frequencies are the same, the weights are only related to the source and the number of channels, so g nfm Equivalent to g nm , g nm is an element in the weight matrix g, and the size of g matrix is ​​N×M.

[0088] Step 7: Mixed graph modeling. Get the mixed graph z ft The probability density function is a special form of non-negative tensor decomposition. Assume that z ft The elements of x are independent of each other. ft The expression is as follows:

[0089]

[0090] Among them, x ft To observe the mixed energy diagram of the signal, Q f x ft The product of is defined as z ft, which can be called the projection matrix. According to the formula derived above, we can get z ft Variance of probability density

[0091] Step 8: Maximum likelihood estimation. ft The multivariate complex Gaussian probability density function is used for maximum likelihood estimation, and the cost function of the fast multi-channel non-negative matrix factorization algorithm (Fast-MNMF) is derived. The formula is as follows:

[0092]

[0093] Among them, n is the signal source, k is the basis, t is the time frame, f is the frequency, m is a channel; W is the basis matrix, the size is N×K×F, w nkf It is equivalent to an element of the matrix W under a certain source n, a certain basis k, and a certain frequency f; H is the activation matrix, with a size of N×K×T, h nkt It is equivalent to an element of the matrix H under a certain source n, a certain basis k, and a certain time t; w nkf and h nkt The power spectral density λ is nft It is obtained by non-negative matrix decomposition; g is the weight matrix, the size is N×M, g nm It is equivalent to an element of matrix g under a certain source n and a certain channel m; Q is the unmixing matrix under all frequencies, q fm It is equivalent to an element of the matrix Q at a certain frequency f and a certain channel m; is the energy of the actual signal source, To estimate the energy of the signal source, when In the iterative process, When the maximum likelihood estimate approaches a certain value, it means that the iteration has converged. ft It is equivalent to an element of the matrix Z at a certain frequency f and a certain time t. In the fast multi-channel non-negative matrix factorization algorithm, the input is a multi-channel (multiple microphones) signal. It is equivalent to an element of the matrix Z at a certain frequency f, a certain time t, and a certain channel m.

[0094] The Fast-MNMF algorithm is used to iteratively optimize the parameters W, H, g, and Q. The steps are as follows:

[0095] Step 8.1: Initialization; Randomly initialize W and H to Gaussian distribution, set K = 2, that is, the basis of the non-negative matrix is ​​2, Q f is the unit matrix, initialized according to cyclic redundancy, the cyclic redundancy formula is as follows:

[0096] ∈=10 -2 ;

[0097] Since W, H, g, and Q have all been initialized, the parameters The initialization formula is as follows, which is equivalent to the formula defined in the cost function above:

[0098]

[0099] Step 8.2: Normalization: Normalize the W, H, G, and Q parameters. The update order is as follows:

[0100]

[0101]

[0102] Step 8.3: Parameter iteration; update the values ​​of W, H, and G respectively, and adjust them according to the cost function. and Make adjustments and the iterative formula is as follows:

[0103]

[0104] The iterative formula is derived from the above maximum likelihood estimation and non-negative matrix decomposition related formulas. This iteration adopts the multiplication update rule. The multiplication update rule is a commonly used weight parameter adjustment method. Its basic idea is to minimize the loss function by continuously adjusting the weight parameters. nkf is an element of the basis matrix W at the fth frequency of the kth basis of the nth source. The basis matrix W and the activation matrix H are derived from the power spectrum density λ through non-negative matrix decomposition. nm is the element in the weight matrix g.

[0105] Step 8.4: Update the unmixing matrix Q; adjust Q by optimizing the diagonalization operation f , ensure Q f The stability and validity of the value, e m is a one-hot vector, and the relevant formula is as follows:

[0106]

[0107] q fm =(Q f V fm ) -1 e m ,

[0108] The algorithm is an iterative projection algorithm (IP), in which: is the actual signal source energy, is the estimated energy of the signal source, T is the total number of time frames, x ft Observed mixed energy spectrum, X ft is x ft The covariance matrix of fm A row for Q.

[0109] Step 9: After completing the iterative optimization, substitute the final W, H, g, Q parameters into the Wiener filter formula to separate the biomedical signal. The specific steps are as follows:

[0110] Step 9.1: Calculate the Wiener filter matrix W as follows nft , used to estimate the time-frequency spectrum of each source signal:

[0111]

[0112] Among them, λ nft is the power spectral density of the source signal, is the directional weight of the source signal, Represents the weighted sum of all source signals.

[0113] Step 9.2: Signal separation. Use the Wiener filter matrix W nft For the input mixed graph X ft Separate and obtain the time-frequency spectrum X of each source signal nft :

[0114] X nft =W nft X ft ;

[0115] Step 9.3: Time domain signal reconstruction. For each separated source signal X nft Perform inverse short-time Fourier transform (ISTFT) to reconstruct the corresponding time domain signal, thereby obtaining a pure biomedical signal.

[0116] First, in biomedical signals, especially in electroencephalogram (EEG) and electrocardiogram (ECG) signals, noise and interference signals such as artifacts often affect diagnosis and analysis. The NMF separation method can effectively separate useful signals from noise. NMF can separate noise and useful signals by extracting basis components from different sources from the original signal, thereby achieving signal denoising and source separation and improving signal quality. Secondly, since biomedical signals usually involve a large amount of data, standard methods may be very time-consuming, especially in the case of multi-channel signals. The key to the fast multi-channel NMF method is to improve computational efficiency and accelerate the iteration process by jointly diagonalizing the spatial covariance matrix. Finally, the key point of biomedical signal separation based on fast multi-channel non-negative matrix factorization is to combine the decomposition principle of NMF, the processing of multi-channel data, and the application of fast optimization algorithms to achieve efficient separation, denoising and source extraction of biomedical signals. This method not only improves processing efficiency, but also makes innovations in signal quality improvement and application accuracy, and has broad application potential.

[0117] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention, rather than to limit it. Although the present invention has been described in detail with reference to the above embodiments, those skilled in the art should understand that they can still modify the technical solutions described in the above embodiments, or replace some or all of the technical features therein with equivalents. However, these modifications or replacements do not cause the essence of the corresponding technical solutions to deviate from the scope defined by the claims of the present invention.

Claims

1. A biomedical signal separation method based on fast multi-channel non-negative matrix factorization, characterized in that: The method includes data preprocessing and formal separation processes; the data preprocessing includes wavelet denoising of the input observed mixed biomedical signal; The formal separation process includes inputting the pre-processed mixed biomedical signal, converting the signal from the time domain image to the time-frequency domain image by short-time Fourier transform, fast multi-channel non-negative matrix decomposition algorithm separation, and Wiener filter separation; the fast multi-channel non-negative matrix decomposition algorithm separation includes assuming signal distribution, mixed spectrum covariance matrix modeling, power spectrum density non-negative matrix decomposition, spatial covariance matrix joint diagonalization, mixed spectrum modeling, maximum likelihood estimation, and iterative optimization of parameters using the fast multi-channel non-negative matrix decomposition algorithm; The parameters obtained by iterative optimization of the fast multi-channel non-negative matrix factorization algorithm are substituted into the formula of Wiener filtering to separate the biomedical signals.

2. The method for separating biomedical signals based on fast multi-channel non-negative matrix factorization according to claim 1, characterized in that: The wavelet denoising adopts a wavelet transform method based on an improved threshold, combines a soft threshold function with a hard threshold function, and performs denoising processing; according to the change of noise with the wavelet transform, the mean square error σ of the wavelet system in each layer is estimated layer by layer i , calculate the threshold of each layer The threshold is trimmed by ln (1+i) to obtain the threshold at different scales; the specific formula of wavelet transform for improving the threshold is as follows: Among them, d i,j is the wavelet coefficient of each layer, is the wavelet coefficient adjusted by threshold, σ is the mean square error of noise, σ i is the mean square error of each layer, i is the number of layers of wavelet decomposition, j is the length of the detail signal, μ is the universal threshold estimate, and μ i is the threshold estimation for each layer; L is the length of the original signal, a is the adjustment factor, when a takes different values, the wavelet transform formula of the improved threshold changes between soft and hard functions, and has the advantages of both soft and hard functions; In the short-time Fourier transform, the parameters are set as a window length of 1024, a frame shift of 50% of the window length, and a FFT point number of 1024, and a time-frequency energy diagram of a three-channel mixture is generated.

3. The method for separating biomedical signals based on fast multi-channel non-negative matrix factorization according to claim 2, characterized in that: In the hypothetical signal distribution, it is assumed that the mixed spectrum of the short-time Fourier transform obeys the centralized multivariate complex Gaussian distribution x~N c (∑), that is, the distribution of the signal has nothing to do with the mean, but only with the covariance. Assuming that the observed time-frequency energy diagram is the superposition of the three source signal spectra, since the spectrum of each source obeys the complex Gaussian distribution, the following formula is obtained: Where p(x) is the probability density function of the time-frequency point x in the energy map obtained after short-time Fourier transform, n is the signal source, N is the number of signal sources, f is the frequency, and t is the time frame; π is a mathematical constant, which is 3.14; M is the number of signal channels, and one channel is equivalent to one microphone; H represents the conjugate transpose, ∑ is the covariance of the centered multivariate complex Gaussian distribution, and x nft For each source, the energy diagram, x ft is the energy diagram of the mixed signal.

4. The method for separating biomedical signals based on fast multi-channel non-negative matrix factorization according to claim 3, characterized in that: In the mixed spectrum covariance matrix modeling, it is assumed that the covariance matrix of the probability density function of the mixed spectrum is equal to the product of the signal power spectrum density and the spatial covariance matrix. According to the additivity of the spectrum and the additivity of the Gaussian distribution, the expression of the mixed spectrum is obtained as follows: Among them, λ nft is the power spectral density of the low-rank source signal, G nf is the spatial covariance matrix of the source signal, Y ft is the covariance matrix of the multivariate complex Gaussian distribution, N c It is the abbreviation of multivariate complex Gaussian distribution.

5. The method for separating biomedical signals based on fast multi-channel non-negative matrix factorization according to claim 4, characterized in that: In the power spectral density non-negative matrix decomposition, the power spectral density λ of the low-rank source signal nft Perform non-negative matrix decomposition, the expression is as follows: Where K is the number of bases, w nkf and h nkt represent the amplitude and activation degree of the basis, respectively.

6. The method for separating biomedical signals based on fast multi-channel non-negative matrix factorization according to claim 5, characterized in that: In the joint diagonalization of the spatial covariance matrix, the weight-sharing joint diagonalization spatial covariance matrix G is adopted. nf ,Will Defined as the non-negative weight vector of the original signal, the expression is as follows: Among them, G nf Covariance matrix; Q f is the unmixing matrix, with a size of M×M×F, where M is the number of channels and F is the number of frequencies. Each F is a unitary matrix; H represents conjugate transpose; Q f Initialize to the identity matrix; q f1 ,…,q fM Indicates the value of different channels at the same frequency; g nfm is an element in the weight matrix g, the size of g is N×M×F, N is the number of signal sources, It represents the weights of different channels under the same source and the same frequency. Since the weights under different frequencies are shared, that is, the weights under different frequencies are the same, the weights are only related to the source and the number of channels, so g nfm Equivalent to g nm , g nm is an element in the weight matrix g, and the size of g matrix is ​​N×M.

7. The method for separating biomedical signals based on fast multi-channel non-negative matrix factorization according to claim 6, characterized in that: In the mixed spectrum modeling, the mixed spectrum z is obtained ft The probability density function is a special form of non-negative tensor decomposition. Assume that z ft The elements of x are independent of each other, and ft Then the relevant expression is as follows: Among them, Q f x ft The product of is defined as z ft , called the projection matrix; according to the formula derived above, we get z ft Variance of probability density 8. The method for separating biomedical signals based on fast multi-channel non-negative matrix factorization according to claim 7, characterized in that: In the maximum likelihood estimation, for z ft The multivariate complex Gaussian probability density function is used for maximum likelihood estimation, and the cost function of the fast multi-channel non-negative matrix decomposition algorithm is derived. The formula is as follows: Among them, n is the signal source, k is the basis, t is the time frame, f is the frequency, m is a channel; W is the basis matrix, the size is N×K×F, w nkf It is equivalent to an element of the matrix W under a certain source n, a certain basis k, and a certain frequency f; H is the activation matrix, with a size of N×K×T, h nkt It is equivalent to an element of the matrix H under a certain source n, a certain basis k, and a certain time t; w nkf and h nkt The power spectral density λ is nft It is obtained by non-negative matrix decomposition; g is the weight matrix, the size is N×M, g nm It is equivalent to an element of matrix g under a certain source n and a certain channel m; Q is the unmixing matrix under all frequencies f, q fm It is equivalent to an element of the matrix Q at a certain frequency f and a certain channel m; is the energy of the actual signal source, To estimate the energy of the signal source, when In the iterative process, When the maximum likelihood estimate approaches a certain value, it means that the iteration has converged. ft It is equivalent to an element of the matrix Z at a certain frequency f and a certain time t. In the fast multi-channel non-negative matrix decomposition algorithm, the input is a multi-channel signal. It is equivalent to an element of the matrix Z at a certain frequency f, a certain time t, and a certain channel m; The fast multi-channel non-negative matrix factorization (Fast-MNMF) algorithm is used to iteratively optimize the parameters W, H, g, and Q. The specific method is as follows: Step 8.1: Initialization; Initialize W and H to random Gaussian distribution, set K = 2, that is, the basis of non-negative matrix decomposition is 2, Q f is the identity matrix, Initialize according to cyclic redundancy. The cyclic redundancy formula is as follows: Since W, H, g, and Q have all been initialized, the parameters The initialization formula is as follows: Step 8.2: Normalization; Normalize the W, H, g, and Q parameters, and substitute the normalized matrix into the following iterative formula. The normalization update order is as follows: Step 8.3: Parameter iteration; update the values ​​of W, H, and g respectively, and adjust them according to the cost function. and Make adjustments and update the iteration formula as follows: Step 8.4: Update the unmixing matrix Q; adjust Q by optimizing the diagonalization operation f , ensure Q f The stability and validity of the value, e m is a one-hot vector, and the relevant formula is as follows: Among them, X ft is x ft The covariance matrix of fm A row for Q.

9. The method for separating biomedical signals based on fast multi-channel non-negative matrix factorization according to claim 8, characterized in that: In the Wiener filter separation, the finally obtained W, H, g, Q parameters are substituted into the Wiener filter formula to separate the biomedical signal. The specific method is: Step 9.1: Calculate the Wiener filter matrix W as follows nft , used to estimate the time-frequency spectrum of each source signal: Among them, λ nft is the power spectral density of the source signal, is the directional weight of the source signal, represents the weighted sum of all source signals; Step 9.2: Signal separation; using the Wiener filter matrix W nft For the input mixed graph X ft Separate and obtain the time-frequency spectrum X of each source signal nft : X nft =W nft X ft ; Step 9.3: Time domain signal reconstruction: for each separated source signal X nft Perform short-time inverse Fourier transform to reconstruct the corresponding time domain signal, thereby obtaining a pure biomedical signal.

Images