A magnetocardiogram signal automatic noise reduction method based on independent component analysis

Abstract

The application provides a magnetocardiogram automatic noise reduction method based on independent component analysis, which comprises the following steps: collecting multi-channel resting-state magnetocardiogram signals by using a sensor, removing high-frequency interference and baseline drift, and performing timing segmentation and standardization processing on the magnetocardiogram signals; calculating each independent component of the magnetocardiogram signals by using an independent component analysis algorithm, and obtaining labeled information by visual inspection; performing sliding window clipping processing and spectrum analysis on each independent component; performing feature extraction based on the time-domain and frequency-domain independent component data sets; inputting the time-domain independent component data set and the extracted features into a time series classification network, training the labeled information as a label, and obtaining an automatic independent component analysis model. The method provides a standardized and automatic preprocessing procedure for magnetocardiogram signals, can be applied to complex use scenarios, has a simple processing procedure, high automation degree, good robustness, strong flexibility, and high medical application value.

Description

TECHNICAL FIELD

[0001] The application belongs to the field of biomedical signal processing, and particularly relates to a magnetocardiogram signal automatic denoising method based on independent component analysis. BACKGROUND

[0002] Magnetocardiography is a comprehensive use of extremely weak magnetic field measurement technology to measure the spatial magnetic field distribution of the heart due to the electrical physiological activity of the heart. Magnetocardiogram signals contain rich physiological and pathological information and have become a promising non-invasive diagnostic tool with important medical application value. However, in the actual clinical acquisition process, magnetocardiogram signals will be disturbed by many noises, such as other medical instrument interference, power frequency noise and other artifacts, etc. Therefore, in order to obtain high-quality magnetocardiogram signals with diagnostic value, perfect magnetocardiogram signal preprocessing technology is needed, and the robustness and automation degree of the preprocessing technology are also required.

[0003] The magnetocardiogram signal preprocessing method is mainly divided into single-channel processing method and multi-channel processing method. The single-channel processing method mainly includes frequency domain filtering, wavelet transform, empirical mode decomposition, etc.; the multi-channel processing method mainly includes independent component analysis, signal space projection, etc. In actual application, two or three preprocessing methods are usually combined to obtain magnetocardiogram signals with diagnostic value. For high-frequency noise and baseline drift in magnetocardiogram signals, traditional methods such as frequency domain filtering or wavelet transform can achieve good results. However, the respiratory artifacts and specific frequency noise that occur in the actual acquisition process usually have burstiness or are mixed in the frequency band range of the magnetocardiogram signals, and using traditional filtering methods will lose the original information of the magnetocardiogram signals and cause serious distortion. The signal space projection method projects the measured signal onto a subspace orthogonal to the interference space distribution to achieve noise reduction, but it often needs to measure the noise data of the acquisition environment synchronously, making it difficult to implement in actual application.

[0004] Independent component analysis is a blind source separation algorithm, which is based on the assumption that each independent component is linear and non-Gaussian. In the magnetocardiogram signal preprocessing process, it can effectively separate the signal source and the noise source, and has better robustness and practicality compared with other preprocessing methods. Although independent component analysis has obvious advantages in removing artifacts and specific frequency noise, visual inspection and manual removal of artifact signals are needed after each decomposition. Therefore, in order to achieve better preprocessing effect and improve the automation degree of the whole preprocessing process, it is necessary to develop an automatic independent component analysis algorithm to realize the automatic denoising of magnetocardiogram signals. SUMMARY

[0005] In order to solve the above technical problems, the application provides a magnetocardiogram automatic noise reduction method based on independent component analysis, which automatically identifies signal sources and noise sources for each independent component decomposed from the magnetocardiogram, so as to solve the automatic noise reduction demand of the magnetocardiogram collected in the clinic.

[0006] In order to achieve the above-mentioned purpose, the technical scheme adopted by the application is as follows,

[0007] A magnetocardiogram automatic noise reduction method based on independent component analysis comprises the following steps:

[0008] S10: using a sensor to collect multi-channel resting-state magnetocardiogram, removing high-frequency interference and baseline drift, and performing time sequence segmentation and standardization processing on the magnetocardiogram;

[0009] S20: calculating each independent component of the magnetocardiogram by using an independent component analysis algorithm, and obtaining labeled information by visual inspection, performing sliding window clipping processing and spectrum analysis on each independent component, and obtaining an independent component dataset containing time domain and frequency domain information;

[0010] S30: performing feature extraction based on the independent component dataset containing time domain and frequency domain information;

[0011] S40: inputting the independent component dataset and the features extracted in step S30 into a time sequence classification network, using the labeled information as a label to train the time sequence classification network, finally obtaining an automatic independent component analysis model, and realizing magnetocardiogram automatic noise reduction based on independent component analysis; the time sequence classification network comprises an internal feature extractor based on a deformable convolutional neural network DCNN and a long short-term memory neural network LSTM, and an external feature fusioner based on a dendritic network DD.

[0012] Optionally, the S10 comprises:

[0013] S11: using an optically pumped magnetometer to collect multi-channel resting-state magnetocardiogram;

[0014] S12: using a wavelet threshold denoising method to remove high-frequency interference and baseline drift in the magnetocardiogram;

[0015] S13: performing time sequence segmentation and standardization processing on the magnetocardiogram to obtain a magnetocardiogram segment dataset.

[0016] Optionally, the S20 comprises:

[0017] S21: performing principal component analysis on each magnetocardiogram segment, calculating the cumulative variance contribution rate of the first N principal components, determining the number of independent components calculated by using the independent component analysis algorithm for each magnetocardiogram segment according to the variance contribution rate threshold 99.5%, manually labeling the noise source and the signal source through visual inspection to obtain a labeled independent component dataset.

[0018] S22, a sliding window clipping process is performed on the independent component data set, and a spectrum analysis is performed on each independent component to obtain an independent component data set containing time domain and frequency domain information.

[0019] Optionally, the S30 comprises:

[0020] The time domain and frequency domain features of the independent component data set are extracted respectively as the input part of the external feature fusioner in the next step; the time domain features include kurtosis, skewness, pulse factor, waveform factor, peak factor and margin factor; and the frequency domain features include barycentric frequency, root mean square frequency, frequency variance and 25Hz energy value.

[0021] Optionally, the S40 comprises:

[0022] S41, the DCNN-LSTM internal feature extractor serves as a backbone network, the input of the DCNN-LSTM internal feature extractor is the independent component data set decomposed from each magnetocardiogram signal segment, and the output of the DCNN-LSTM internal feature extractor is a prediction result of each independent component, the prediction result contains prediction labels of noise sources and signal sources, and the number of output results is adaptively adjusted according to the number of independent components of each magnetocardiogram signal segment.

[0023] S42, the external feature fusioner based on the dendritic network DD, the logical relationship between the features is calculated and linearly combined by introducing a logic extractor, feature engineering is realized, and the features are fused through Hadamard product, and after the external features are changed, the backbone network is taken as a pre-training model, and the first two modules of the network are fixed, so that the back propagation process is prevented.

[0024] The present application has the following advantages:

[0025] The present application realizes the magnetocardiogram signal automatic denoising method based on independent component analysis, which greatly improves the preprocessing efficiency of the magnetocardiogram signal without manual operation, and can better suppress various noises and artifacts in the clinical environment, and makes up for the shortcomings of previous algorithms in automation degree and preprocessing effect.

[0026] On the other hand, the present application overcomes the problem of unstable order of independent components by introducing DCNN, DCNN can more flexibly extract different independent component features of each magnetocardiogram signal segment, and through LSTM, the time sequence features are extracted, so that the expression ability of the whole model for the independent components of the magnetocardiogram signal is enhanced.

[0027] On the other hand, the present application realizes the fine tuning of the backbone network by introducing the external feature fusioner, when switching to other application scenarios, the internal network structure does not need to be adjusted, and the transfer learning can be realized by changing the external features.

[0028] This invention improves the preprocessing technology for magnetic resonance imaging (MRI) signals. Due to its high degree of automation and excellent real-time performance, it can be directly applied to real-time clinical MRI acquisition systems. Compared with other preprocessing algorithms, it has better robustness, practicality, and promotional value. Attached Figure Description

[0029] Figure 1 This is a flowchart of an automatic noise reduction method for magnetic cardiomyography signals based on independent component analysis according to the present invention.

[0030] Figure 2 This is a diagram of the DCNN-LSTM network structure of the present invention;

[0031] Figure 3 This is a diagram of the convolutional layer structure of the present invention;

[0032] Figure 4 This is a structural diagram of the DCNN module of the present invention;

[0033] Figure 5 These are the accuracy and loss curves for the network training and validation sets of this invention. Detailed Implementation

[0034] The present invention will now be described in detail with reference to the accompanying drawings and embodiments. It should be understood that the specific embodiments described herein are only for explaining the present invention and are not intended to limit the present invention.

[0035] like Figure 1 As shown, the automatic noise reduction method for magnetic resonance imaging (MRI) signals based on independent component analysis of the present invention is implemented as follows:

[0036] S10: Multi-channel resting-state magnetocardiogram (MCC) signals were acquired using an optically pumped magnetometer. The number of acquisition channels was 36, and the sampling rate was 1000Hz. The optically pumped magnetometer is an extremely weak magnetic field measurement sensor based on spin-exchange without relaxation. Wavelet transform was used to remove high-frequency noise and baseline drift from the acquired MCC signals, and time-series segmentation and standardization were performed. The specific process is as follows:

[0037] (1) The Mallat pyramid algorithm was used to perform wavelet decomposition on the magnetocardiogram signal. The wavelet basis function was selected as "Daubechies". The number of decomposition layers was 8. The threshold was determined by unbiased risk estimation and the threshold was processed using a soft threshold function.

[0038] (2) Use a 30-second sliding window with an overlap rate of 20% to perform temporal segmentation of the magnetic heart signal and perform Z-Score normalization on each magnetic heart signal segment.

[0039] S20: Determine the number of independent components decomposed by principal component analysis on each magnetocardiogram segment, calculate each independent component of the magnetocardiogram signal by independent component analysis algorithm. Visually inspect each independent component and manually label noise sources and signal sources to obtain an independent component dataset. The specific process is as follows:

[0040] (1) Perform principal component analysis on each 36-channel magnetocardiogram segment as observation signal X, calculate the covariance matrix C of the 36-channel magnetocardiogram segment and calculate its eigenvalue matrix S and eigenvector matrix E; the eigenvalues in the eigenvalue matrix S are arranged in descending order: , wherein the eigenvalue represents the principal component variance, and the eigenvector represents the principal component direction. Generally, the principal component with a larger eigenvalue contains more information, so the number of independent components generated in the independent component analysis algorithm is determined by the cumulative variance contribution rate; in this invention, the cumulative variance contribution rate threshold is 99.5%, that is:

[0041]

[0042] , wherein M is the number of calculated eigenvalues, and N represents the number of cumulative variance contribution rates of the first N principal component eigenvalues greater than 99.5%.

[0043] (2) Generate H independent components using the independent component analysis algorithm.

[0044] a. First, center the magnetocardiogram segment: , wherein represents the original magnetocardiogram segment, represents the expectation of the original magnetocardiogram segment, represents the magnetocardiogram segment after centering.

[0045] b. Whitening, the centered magnetocardiogram segment is whitened to obtain data , , wherein the whitening matrix is represented as , the matrix is the eigenvalue matrix in the above principal component analysis step, and the matrix E is the eigenvector matrix in the above principal component analysis step.

[0046] c. Initialize the separation matrix W, , randomly initialize all , wherein represents each column of the separation matrix W, and satisfies ;

[0047] d. For each , the update rule is:

[0048]

[0049] where g(y) represents the cumulative distribution function of the signal source s, where y represents the input signal, is a constant, usually taken in the range of 1 to 2, and in the present application tanh(y) is the hyperbolic tangent function. E() represents the expectation of a matrix. e. orthogonalizing the matrix

[0050] according to the convergence, repeating steps d and e until the above-mentioned number of independent component sets are decomposed. (3) visually inspecting each signal and making a judgment in combination with kurtosis threshold processing, first determining the kurtosis threshold of the target source signal to be 0.05 according to prior information, and the kurtosis calculation formula is:

[0051]

[0052] E() represents the expectation of a matrix.

[0053] According to the visual inspection and the kurtosis threshold, the independent component categories are preliminarily determined, and the actual magnetocardiogram segment processing result is combined to make a final interpretation of each independent component category, determine the signal source and the noise source, and label to obtain a labeled independent component dataset. The visual inspection process is discussed together with clinicians and experts, and mainly identifies cardiac source signals, impact noise, various artifacts and specific frequency noise.

[0054] (4) Data augmentation is performed on the dataset, the labeled independent component dataset is processed by sliding window cropping, the cropping window size is 90% of the length of the independent component data, and the cropping position is randomly divided. Each independent component obtains 20 pieces of augmented data after sliding window cropping processing, and an augmented independent component dataset is generated. The augmented independent component dataset is subjected to spectral analysis, and the power spectral density of each independent component is calculated.

[0055] S30: six time domain features and four frequency domain features are extracted based on the independent component dataset. The time domain features include kurtosis, skewness, pulse factor, waveform factor, peak factor and margin factor. The frequency domain features include center of gravity frequency, root mean square frequency, frequency variance and 25Hz energy value. The time domain features are explained as follows:

[0056] ​​​​​​​The kurtosis formula is described above, which is used to measure whether a signal conforms to a normal distribution;

[0057] The skewness S represents the degree to which the signal distribution deviates from the normal distribution to the left or right:

[0058] Assuming that the magnetocardiogram signal segment is N represents the length of the magnetocardiogram signal segment. The skewness S is represented as:

[0059]

[0060] wherein, represents the mean, represents the standard deviation.

[0061] The pulse factor I is represented as:

[0062]

[0063] wherein, represents the maximum value in the magnetocardiogram signal segment.

[0064] The waveform factor W is:

[0065]

[0066] The peak factor C is:

[0067]

[0068] The margin factor L is:

[0069]

[0070] wherein, the kurtosis and the skewness are measurement indexes of the difference in signal form between the signal source and the noise source, the pulse factor, the waveform factor, the peak factor, and the margin factor have similar physical meanings, and all can be used as indexes for detecting impact signals or alternating current components in the signal. Therefore, the above six features can be selected to distinguish specific frequency noise sources or impact noise or artifacts at a certain moment from the time domain signal. However, in actual application scenarios, part of the specific frequency noise will be mixed in the effective magnetocardiogram signal frequency band range, and therefore it is necessary to extract features from the frequency domain. The extracted frequency domain features are as follows:

[0071] The center of gravity frequency is:

[0072]

[0073] The root mean square frequency is:

[0074]

[0075] The frequency variance is:

[0076]

[0077] wherein, represents the corresponding frequency, represents the energy value of the corresponding frequency. In the actual acquisition process, other frequency domain features such as the energy value of a specific frequency can also be introduced. In the application scenario, the energy value of 25 Hz frequency is used to identify 25 Hz power frequency noise.

[0078] S40: input the independent component data set and the extracted features into a time series classification network, the input signal dimension is (N x L), N represents the number of independent components, L represents the length of the independent components, the learning rate is 0.001, the input batch is 8, the cross-entropy loss function is used as the loss function of the network, and the SGD optimizer is used for training, and the network is trained for 100 epochs on the training set and the validation set. The main framework of the time series classification network is as shown in Figure 2 The main framework of the time series classification network is further explained as follows:

[0079] (1) DCNN-LSTM internal feature extractor:

[0080] This part serves as the backbone network, which extracts deep features of the input N independent components through four convolutional layers. The first convolutional layer structure is as shown in Figure 3 In the figure, Conv represents the convolution kernel in the convolutional layer, 1 x 3 represents the convolution kernel size, 16 represents the number of convolution kernels, and Stride represents the step length during convolution operation. The number of convolution kernels of the four convolutional layers is (16, 32, 64, 128), the convolution kernel size is 3, all convolutional layers use the Relu activation function, and batch normalization (BN) operation is performed during training.

[0081] The deformable convolutional layer (DCNN) is used to further extract deep features in the feature map generated in the previous step. The basic structure of the deformable convolutional layer is as shown in Figure 4 In the figure, Convolutional Layer represents the bias layer in the deformable convolutional layer, offset field represents the bias layer output, and offsets represent the trained bias. The DCNN module-1 includes four deformable convolutional layers, and the number of convolution kernels of each deformable convolutional layer is (128, 32, 16, 6), and the convolution kernel size is 3. The offset refers to the input feature position actually participating in the calculation of the convolution kernel weight. Assuming that the input feature expression is The receptive field size of the deformable convolution is 3 x 3, and the receptive field position domain is wherein the center position of the receptive field corresponds to the coordinates The feature map calculation formula of the corresponding position is as follows:

[0082]

[0083] The offset is reflected in , is the offset, and since there will be a decimal, the three spline interpolation method is used in the present application to obtain . Through the offset of deformable convolution, the problem of generating independent components with an unstable order can be indirectly compensated, so that the convolution kernel can focus more on expressing single or independent components of the same category. After passing through the DCNN module-1, the number and size of the convolution kernel are reduced to keep consistent with the input dimension to facilitate subsequent external feature fusion.

[0084] After fusing external features, continue to use DCNN module-2 for deep feature extraction. The basic structure of DCNN module-2 is similar to that of DCNN module-1, and the number of convolution kernels is (16, 32, 64, 128), and the size of the convolution kernel is 3.

[0085] Then the final interpretation information is obtained through the LSTM and the fully connected layer. The LSTM network is composed of a series of LSTM units, and each unit is composed of a forgetting gate , an input gate , an output gate , an output and a cell state , wherein t represents the current time, and t-1 represents the previous time. The output state at the current time can be represented as , represents the cell state update value, , represents the state update weight vector parameter, represents the activation function. The matrix splicing of the output at the previous time and the input features at the current time is multiplied by the weight vector parameter, and the bias term is introduced, so that the cell state update value at the current time can be obtained. The forgetting gate at the current time can be represented as , wherein represents the weight vector parameter of the forgetting gate, represents the sigmoid activation function, and the bias term is . The input gate is represented as , represents the weight vector parameter of the input gate, and the bias term is . which is used to control which input features are used to update . The output gate is represented as , denotes the weight vector parameter of the input gate, and the bias term is which controls the output , and finally the output of each node .

[0086] (2) External feature fusioner based on dendritic network

[0087] The external feature fusioner based on dendritic network mainly realizes the function of feature engineering, and the dendritic network model can be expressed as:

[0088]

[0089] Among them, and denote the input and output, denote the weight matrix, denote the Hadamard product, and X denotes the input feature. The dendritic network is similar to a logical extractor, which can capture the logical relationship between features. The multi-layer dendritic network architecture is:

[0090] Y= W L, L-1 [⋯W l,l-1 (⋯W 21 ( W 10 X° X)° X⋯)° X⋯]

[0091] Among them, L is the number of modules, and 3-layer dendritic network is used in the present application to fuse features. The input of the external feature fusioner is 6 time domain features and 4 frequency domain features proposed by the S30, and the output is 6 combined feature values, which are weighted in the form of Hadamard product to the features output by the DCNN module-1. Since the external feature fusioner based on dendritic network only contains matrix multiplication and Hadamard product in calculation, the calculation complexity is low, which ensures the training efficiency of the feature fusioner and the flexibility of the whole model in the case of changing the use scene.

[0092] Figure 5 The training result of the time series classification network based on independent component data set of the present application, and the related evaluation indicators also include accuracy, recall rate, F1 score and confusion matrix.

[0093] In summary, the magnetocardiogram automatic noise reduction method based on independent component analysis of the present application provides a standardized and automated preprocessing process for magnetocardiogram, which can be applied to more complex use scenarios. Its advantages are high automation degree, good robustness and strong flexibility, and has high medical application value.

[0094] The above-described specific embodiments further illustrate the purpose, technical solutions and beneficial effects of the present application, and it should be understood that the above-described specific embodiments are merely specific embodiments of the present application and are not intended to limit the present application, and any modifications, equivalent replacements, improvements, etc. made within the spirit and principle of the present application shall be included in the protection scope of the present application.

Claims

1. A method for automatic noise reduction of magnetocardiographic signals based on independent component analysis, characterized in that, For each independent component decomposed from the magnetocardiogram signal, the signal source and the noise source are automatically identified, and the steps include the following: S10: using a sensor to collect multi-channel resting-state magnetocardiogram signals, removing high-frequency interference and baseline drift, and performing time sequence segmentation and standardization processing on the magnetocardiogram signals; S20: calculating each independent component of the magnetocardiogram signal by an independent component analysis algorithm, and obtaining labeled information by visual inspection, performing sliding window clipping processing and spectral analysis on each independent component to obtain an independent component dataset containing time domain and frequency domain information; S30: feature extraction based on the independent component dataset containing time domain and frequency domain information; S40: inputting the independent component dataset and the features extracted in step S30 into a time sequence classification network, training the time sequence classification network with the labeled information as labels, and finally obtaining an automatic independent component analysis model to realize automatic noise reduction of the magnetocardiogram signal based on independent component analysis; the time sequence classification network includes an internal feature extractor based on a deformable convolutional neural network (DCNN) and a long short-term memory neural network (LSTM), and an external feature fusioner based on a dendritic network (DD).

2. The independent component analysis-based automatic noise reduction method for magnetocardiographic signals according to claim 1, characterized in that, The step S10 includes: S11, using an optically pumped magnetometer to collect multi-channel resting-state magnetocardiogram signals; S12, using a wavelet threshold denoising method to remove high-frequency interference and baseline drift in the magnetocardiogram signal; S13, performing time sequence segmentation and standardization processing on the magnetocardiogram signal to obtain a magnetocardiogram signal segment dataset.

3. The independent component analysis-based automatic noise reduction method for magnetocardiographic signals according to claim 1, characterized in that, The S20 includes: S21, performing principal component analysis on each magnetocardiogram signal segment, calculating the cumulative variance contribution rate of the first N principal components, determining the number of independent components calculated by the independent component analysis algorithm for each magnetocardiogram signal segment according to the variance contribution rate threshold of 99.5%, manually labeling the noise source and the signal source by visual inspection to obtain a labeled independent component dataset; S22, performing sliding window clipping processing on the independent component dataset, and performing spectral analysis on each independent component to obtain an independent component dataset containing time domain and frequency domain information.

4. The independent component analysis-based automatic noise reduction method for magnetocardiographic signals according to claim 1, characterized in that, The S30 includes: Respectively extracting time domain and frequency domain features of the independent component dataset as input parts of the external feature fusioner in the next step; the time domain features include kurtosis, skewness, pulse factor, waveform factor, peak factor and margin factor; the frequency domain features include barycentric frequency, root mean square frequency, frequency variance and 25Hz energy value.

5. The independent component analysis-based automatic noise reduction method for magnetocardiographic signals according to claim 1, characterized in that, The S40 includes: S41, the DCNN-LSTM internal feature extractor as the backbone network, the input of which is the independent component dataset decomposed from each magnetocardiogram signal segment, and the output of which is the prediction result of each independent component, the prediction result containing the prediction labels of the noise source and the signal source, and the number of output results being adaptively adjusted according to the number of independent components of each magnetocardiogram signal segment; S42, the external feature fusioner based on the dendritic network DD, realizes feature engineering and feature fusion with the internal features through Hadamard product by introducing a logic extractor to calculate the logical relationship between features and linear combination; after the external features change, the backbone network is used as a pre-training model and the first two modules of the network are fixed to prevent the process of back propagation.

Images