A method for optimizing magnetocardiogram signal superposition averaging based on heartbeat classification

Abstract

The application provides a magnetocardiogram signal superposition average optimization method based on heart beat classification, and is suitable for the field of magnetocardiogram signal processing and analysis. The method comprises the following steps: obtaining a denoised magnetocardiogram signal; performing R wave detection on the denoised magnetocardiogram signal and performing heart beat segmentation; setting the segmented first heart beat as a heart beat classification, and calculating the correlation of subsequent heart beats with the first heart beat; the heart beat with a correlation greater than a threshold value is considered to belong to the existing classification, and the heart beat with a correlation less than the threshold value is taken as a new heart beat classification; repeating the above process until the correlation calculation of all heart beats is completed; superimposing and averaging all heart beats in the same classification; calculating the correlation between the superimposed and averaged results again, merging the classifications with high correlation, recalculating the superimposed average, and obtaining the result of the heart beat classification superimposed average. The present application improves the existing magnetocardiogram superposition average algorithm and provides a more accurate magnetocardiogram superposition average image.

Description

TECHNICAL FIELD

[0001] The application belongs to the field of processing and analysis of magnetocardiogram signals, and particularly relates to a magnetocardiogram signal superposition average optimization method based on heartbeat classification. BACKGROUND

[0002] With the periodic contraction and expansion of the heart, a periodic potential difference is generated on the surface of the heart, and thus a cardiac magnetic field is generated. Magnetocardiography has the advantages of non-invasiveness, non-contact, spatial resolution, and the like, contains rich physiological and pathological information, and has important clinical significance. In the processing of magnetocardiogram signals, due to the irregularity of heart activity and the interference of the external environment, the magnetocardiogram signals often contain a large amount of noise, which will affect the analysis and recognition of the signals. Therefore, it is necessary to perform signal processing to extract useful information.

[0003] At present, the superposition average method of magnetocardiogram signals is to collect synchronous magnetocardiogram signals, take the R-wave peak of the magnetocardiogram signals as a reference, segment the entire magnetocardiogram data into a single cycle, and then perform average processing to achieve the purpose of suppressing noise. This method can effectively improve the signal-to-noise ratio of the magnetocardiogram signals, so as to better detect the activity of the heart.

[0004] The superposition average method of magnetocardiogram signals may be affected by abnormal signals, which include external noise that cannot be processed, muscle artifacts, and heart rhythm disorders. Introducing these signals into the superposition average of magnetocardiogram signals will cause inaccurate results, so that the superposition average image cannot correctly display the P wave, QRS wave, and T wave, and the final result loses its meaning.

[0005] Therefore, it is of great significance to explore an algorithm that classifies heartbeats before superposition average. SUMMARY

[0006] To solve the above technical problems, the application provides a magnetocardiogram signal superposition average optimization method based on heartbeat classification, which identifies and classifies irregular heartbeats, calculates the superposition average of normal heartbeats and irregular heartbeats respectively, solves the problem of periodical average signal deformation caused by large differences between heartbeats, and improves the accuracy of the superposition average signal.

[0007] To achieve the above purpose, the application adopts the following technical solutions:

[0008] A magnetocardiogram signal superposition average optimization method based on heartbeat classification, comprising the following steps:

[0009] Step 1: performing denoising processing on the magnetocardiogram signals obtained by a multi-channel magnetocardiograph;

[0010] Step 2: performing R-wave detection and heartbeat segmentation on the magnetocardiogram signals after denoising processing;

[0011] Step three, set the first segmented heart beat as the initial heart beat classification, and calculate the correlation between the next heart beat and the initial heart beat classification, and compare it with the first threshold, the heart beat with the correlation greater than the first threshold is determined as the existing heart beat classification, the heart beat with the correlation less than the first threshold is determined as a new heart beat classification, repeat the above steps until the correlation calculation of the subsequent heart beat with all existing heart beat classifications is completed.

[0012] Step four: the heart beats in the same classification are subjected to first superposition average calculation to obtain the first superposition average result of multiple heart beat classifications.

[0013] Step five: calculate the correlation between the first superposition average results of multiple heart beat classifications, merge the heart beat classifications with high correlation and perform second superposition average calculation to obtain the second superposition average result of the heart beat classification.

[0014] Further, step one uses median filtering, low-pass filtering and band-stop filtering to eliminate the baseline drift, high-frequency noise and power frequency noise of the magnetocardiogram signal, uses ICA method to separate the noise component and signal component of the magnetocardiogram signal, sets the noise component to zero and reconstructs to obtain the denoised magnetocardiogram signal.

[0015] Further, step two uses the method combining signal double difference and RR interval to locate the R wave vertex of the heart beat, and obtains a complete heart beat from each R wave vertex forward and backward by taking a plurality of sampling points.

[0016] Further, for each located R wave vertex, N alternative R wave vertices are set, and each alternative R wave vertex corresponds to an alternative complete heart beat.

[0017] Further, step three uses Pearson correlation to calculate the correlation between the next heart beat and the initial heart beat classification, including the correlation calculation between the R wave complete heart beat of the next heart beat and the N alternative complete heart beats and the initial heart beat classification, to obtain N+1 correlation values between the next heart beat and the initial heart beat classification, select the maximum correlation value in the N+1 correlation values and compare it with the set first threshold, if it is greater than the first threshold, it is considered that the next heart beat belongs to the initial heart beat classification, if it is less than the first threshold, it is considered that the next heart beat belongs to a new heart beat classification; repeat the above steps until the correlation calculation of the subsequent heart beat with all existing heart beat classifications is completed.

[0018] Further, the first superposition average calculation in step four refers to adding all heart beat values and then dividing by the number of heart beats in the classification to obtain multiple first superposition average results corresponding to multiple heart beat classifications.

[0019] Further, the step five comprises: calculating a Pearson correlation matrix for the plurality of first superimposed average results, merging the categories with a correlation greater than a second threshold value into one category, and finally performing a second superimposed average calculation on the new category to obtain a second superimposed average result.

[0020] The present application has the following beneficial effects:

[0021] In view of the problem of periodic average signal distortion caused by large differences between continuous heartbeats, a magnetocardiogram superimposed average optimization method based on heartbeat classification is provided. The method uses a digital filter and independent component analysis method to pre-process the magnetocardiogram signal, uses the signal interval double difference method to detect the R wave position, uses the alternative R wave vertex method to compensate for the offset of the R wave positioning, and compares the Pearson correlation between heartbeats to realize the classification of abnormal heartbeats before superimposed average, thereby improving the problem that the superimposed average result is affected by abnormal heartbeats. BRIEF DESCRIPTION OF DRAWINGS

[0022] Figure 1 The present application is a magnetocardiogram superimposed average optimization method based on heartbeat classification flowchart.

[0023] Figure 2 The present application is a normal and abnormal magnetocardiogram signal comparison chart, wherein Figure 2 (a) is a normal magnetocardiogram signal schematic diagram, Figure 2 (b) is a magnetocardiogram signal schematic diagram of ventricular premature beat;

[0024] Figure 3 The present application is a normal and abnormal magnetocardiogram signal classification superimposed average result schematic diagram, wherein Figure 3 (a) is a normal heartbeat classification schematic diagram, Figure 3 (b) is an abnormal heartbeat classification schematic diagram;

[0025] Figure 4 The present application is a classification superimposed average and non-classification superimposed average result comparison chart. DETAILED DESCRIPTION

[0026] In order to make the present application easy to understand, the present application will be further described below in conjunction with the drawings and examples.

[0027] As Figure 1 shown, the present application is a magnetocardiogram superimposed average algorithm based on heartbeat classification, which comprises the following steps:

[0028] Step one: denoising the magnetocardiogram signal obtained by the multi-channel magnetocardiograph. The signal measured by the multi-channel magnetocardiograph contains a large amount of noise. A median filter with a window length of 501 is designed to first remove the baseline drift of the signal measured by the magnetocardiograph, then a 4th order Butterworth low-pass filter with a cutoff frequency of 50 Hz is connected in series to filter out high-frequency noise above 50 Hz, and then a 4th order Butterworth notch filter with a cutoff frequency of 50 Hz is connected in series to filter out 50 Hz power frequency noise.

[0029] The noise-filtered signal is then subjected to independent component analysis, which includes the following steps:

[0030] (1) Centering processing

[0031] The independent component analysis algorithm assumes that all mixed variables and independent components have zero mean, which can greatly simplify the algorithm. In actual situations, it is impossible to guarantee that the obtained signal meets this condition, so the original observation data needs to be centered first, as follows:

[0032] ,

[0033] The above formula represents the input signal, represents the mean of the input signal, represents the result of the centering processing of the input signal. After such preprocessing, the separated independent components will also have zero mean, and the mixing matrix calculated will not change. After obtaining the separated independent components, only the final result needs to be added to recover the mean value subtracted.

[0034] (2) Whitening processing

[0035] The whitening processing is performed on the centered data to obtain data , which is represented as:

[0036] ,

[0037] In the formula, the whitening matrix , the matrix is a diagonal matrix obtained by performing eigenvalue decomposition on the variance matrix of , and E represents the matrix composed of eigenvectors obtained by performing eigenvalue decomposition on the variance matrix . represents taking the square root of each element in D, and the superscript T represents the transpose of the matrix.

[0038] (3) Select the number m of independent components to be estimated

[0039] m represents the number of independent components desired, for magnetocardiogram noise reduction, the principal component analysis method can be selected, and the threshold value is set according to the proportion of different eigenvalues of different principal components to determine the number of independent components to be estimated. According to actual testing, 6 components can contain 99% of the magnetocardiogram signal information, and therefore the number of estimated independent components is 6.

[0040] (4) Initialize the demixing matrix , randomly select data to initialize all , , and satisfy .

[0041] (5) For each , update each element of the matrix, and the rule is , wherein represents the cumulative distribution function of the signal source , represents the transpose of , and represents the mathematical expectation. There are three functions as follows: ,

[0042] ,

[0043] ,

[0044] ,

[0045] In the formula, the constant is suitable to be 1, and is generally 1. is suitable for the recovery of super-Gaussian signal sources, is suitable for the recovery of sub-Gaussian signal sources, and is more suitable if the signal contains super-Gaussian signal and sub-Gaussian signal.

[0046] (6) Symmetric orthogonalization is performed on the matrix :

[0047] ,

[0048] In the formula, is the inverse square root of , which can be obtained by eigenvalue decomposition.

[0049] (7) If the obtained result does not converge, return to step (5).

[0050] ​​The filtered signal is processed by ICA to obtain 6 independent components, and the signal components are reserved by visual observation, the noise components are set to 0, and the independent components are reconstructed to obtain the denoising magnetocardiogram.

[0051] Step two: segmenting each heart beat of the magnetocardiogram. R-wave detection is performed on the preprocessed heart beat sequence, and the R-wave detection method is as follows:

[0052] (1) Detecting the position of QRS wave

[0053] Let the magnetocardiogram data be , , where N is the total number of samples, and the first-order difference operation is performed on the data:

[0054] ,

[0055] Then, the first-order difference operation is performed again on to obtain :

[0056] ,

[0057] Finally, square the to obtain the square double difference array of the signal :

[0058] ,

[0059] Sort the square double difference array in descending order of amplitude, and take 3% of the maximum value as the threshold. Mark the points in greater than the threshold as 1, and mark the points less than the threshold as 0. Perform connectivity analysis on the result, and consider consecutive 1s as a QRS region, and record the start and end points of each region. Since the maximum duration of the QRS region is 150ms, in order to eliminate the possibility of detecting multiple peaks in the same QRS region, all difference peaks within a ±75ms interval of each selected difference peak will be eliminated. The QRS region is identified as a ±75ms window around each selected peak on the ECG data array.

[0060] (2) Detecting R-wave

[0061] The R-wave is a positive peak in the QRS wave segment, and is detected by comparing the relative amplitude of the QRS wave. For each QRS wave region, iterate through , compare the amplitude of each point with the amplitudes of the previous and next two points, and if the amplitude of the current point is greater than the amplitudes of the previous and next two points, consider the current point as a candidate R-wave and record its position and amplitude. Do not use the absolute value of the QRS wave segment to calculate, so as not to calculate the S wave by mistake.

[0062] (3) Optimization according to RR interval

[0063] The R peaks thus obtained can not be accurate. A peak value can be missed or detected incorrectly. To ensure detection accuracy, the RR intervals are processed according to certain criteria. The time interval of the two adjacent candidate peaks obtained in step (2) is calculated, and the average and standard deviation of all time intervals are calculated, denoted as and .

[0064] The candidate R waves are screened and corrected. If the time interval of a candidate peak from the previous candidate peak is less than 0.3 seconds or greater than 1.5 seconds, the candidate peak is considered abnormal and is deleted from the candidate peaks. If the time interval of a candidate peak from the previous candidate peak is greater than 1.66 times the , it is considered that an R wave peak can have been missed before the candidate peak, and the candidate peak is deleted from the candidate peaks and a possible R wave peak is found before it and added to the candidate peaks. If the amplitude of a candidate peak is less than 0.3 times the average of the overall amplitude, the candidate peak is considered to be caused by noise and is deleted from the candidate peaks. If the amplitude of a candidate peak is greater than 1.2 times the average of the overall amplitude, the candidate peak is considered abnormal and is deleted from the candidate peaks and a possible R wave peak is found near it and added to the candidate peaks.

[0065] The last candidate peak is returned as the located R wave position, but the R wave vertex located by the algorithm can have a deviation. For the accuracy of subsequent correlation calculation, for each R wave vertex located by the algorithm, 10 sampling points are taken forward and backward as alternative R wave vertices, and 500 sampling points are taken forward and backward for each alternative R wave vertex to obtain the corresponding complete heartbeats. That is, for each R wave vertex located by the algorithm, there should be an actually located R wave vertex and 20 alternative R wave vertices, corresponding to an actually located heartbeat and 20 alternative heartbeats.

[0066] Step three: set the first segmented heartbeat as an initial heartbeat classification, calculate the correlation of the next heartbeat with the first heartbeat, and the correlation is calculated as Pearson correlation, the formula is:

[0067] ,

[0068] represents the size of each sampling point value of the first heartbeat, represents the average value of the first heartbeat sampling points, represents the size of each sampling point value of the next heartbeat, The average value representing the sample point of the latter heart beat. The correlation calculation between the latter heart beat in this step and the first heart beat means that the R-positioned heart beat of the latter heart beat and the 20 candidate heart beats are all calculated with the first heart beat, and finally 21 values of the correlation with the initial classification are obtained, the maximum value is selected, and compared with the first threshold value, if it is greater than the first threshold value, it is considered that the heart beat belongs to the initial heart beat classification, if the maximum value of the correlation is less than the first threshold value, it is considered that the heart beat belongs to a new classification, and the heart beat is added to the subsequent calculation as a new classification.

[0069] The above operation is repeated, and the correlation of each subsequent heart beat with the existing heart beat classification is calculated, the R-positioned heart beat of the subsequent heart beat and the 20 candidate heart beats are calculated with all the classified heart beats, the maximum correlation is retained for each classification, and after the correlation calculation of all classifications is completed, the maximum correlation of all classifications is selected and compared with the first threshold value, if it is greater than the first threshold value, it is considered that the heart beat belongs to the classification with the maximum correlation, if it is less than the first threshold value, the heart beat is set as a new classification for subsequent correlation calculation.

[0070] Step four: after the correlation calculation of all heart beats is completed, the first superposition average calculation of all heart beats in the same classification is performed, that is, the values of all heart beats in each classification are added and then divided by the number of heart beats in the classification to obtain an average value representing the first superposition average result of each classification.

[0071] Step five: the different first superposition average results obtained in step four are set with a new correlation second threshold value, the correlation between the first superposition average results of step four is calculated, the classifications higher than the correlation second threshold value are merged, and the new classification is re-calculated for the second superposition average calculation, the similar heart beat classification is merged, and the final heart beat classification superposition average result is obtained.

[0072] Figure 2 Fig. 2 is a comparison chart of the normal and premature ventricular contraction magnetocardiogram signals after step one of the present application, Figure 2 (a) is a measured healthy person magnetocardiogram signal, Figure 2 (b) is a measured magnetocardiogram signal of a person with premature ventricular contraction, the horizontal axis is the measurement time, the unit is s, and the vertical axis is the signal intensity, the unit is pT. Figure 3 Fig. 3 is the classification superposition average result obtained by the present application after all processes, Figure 3 (a) shows the classification superposition average result of the normal magnetocardiogram signal, Figure 3 (b) shows the classification superposition average result of the premature ventricular contraction magnetocardiogram signal, Figure 2 Fig. 3b shows the classification superposition average result of the premature ventricular contraction magnetocardiogram signal, the horizontal axis is the measurement time, the unit is s, and the vertical axis is the signal intensity, the unit is pT.Figure 4 In order to compare the results of classified superimposed averaging and non-classified superimposed averaging, the horizontal axis is the measurement time length, in seconds, and the vertical axis is the signal intensity, in pT. The classified heartbeat superimposed average image can clearly show the P wave, QRS wave and T wave. However, the superimposed average image without classification processing is affected by abnormal heartbeats, and the characteristics such as P wave, T wave and PR segment are significantly different from the electrocardiogram of a normal person, proving that the present application has the ability to improve the superimposed average image of the magnetocardiogram containing abnormal heartbeat components.

[0073] The above specific embodiments further illustrate the purpose, technical solutions and beneficial effects of the present application. It should be understood that the above description is only a specific embodiment of the present application and is not intended to limit the present application. Any modification, equivalent replacement, improvement, etc. made within the spirit and principles of the present application should be included in the protection scope of the present application.

Claims

1. A method for optimizing the stacked averaging of magnetocardiographic signals based on heart beat classification, characterized by, The method comprises the following steps: Step one, denoising the magnetocardiogram signal obtained by the multi-channel magnetocardiograph; Step two, detecting the R wave of the denoised magnetocardiogram signal and segmenting the heartbeats; including: using the method of combining the square double difference and the RR interval to locate the R wave vertex of the heartbeat, taking a plurality of sampling points from each R wave vertex to obtain a complete heartbeat corresponding to the R wave vertex, setting N candidate R wave vertices for each located R wave vertex, and each candidate R wave vertex corresponding to a candidate complete heartbeat; Step three, setting the first segmented heartbeat as an initial heartbeat classification, calculating the correlation of the subsequent heartbeat with the initial heartbeat classification and comparing it with the first threshold value, and determining the heartbeat with the correlation greater than the first threshold value as an existing heartbeat classification, and determining the heartbeat with the correlation less than the first threshold value as a new heartbeat classification; including: calculating the correlation between the subsequent heartbeat and the initial heartbeat classification using the Pearson correlation, including calculating the correlation between the R wave complete heartbeat of the subsequent heartbeat and the N candidate complete heartbeats and the initial heartbeat classification, obtaining N+1 correlation values of the subsequent heartbeat and the initial heartbeat classification, selecting the maximum correlation value in the N+1 correlation values, and comparing the maximum correlation value with the set first threshold value, if the maximum correlation value is greater than the first threshold value, it is considered that the subsequent heartbeat belongs to the initial heartbeat classification, if the maximum correlation value is less than the first threshold value, it is considered that the subsequent heartbeat belongs to a new heartbeat classification; repeating step three until the correlation calculation of the subsequent heartbeat with all existing heartbeat classifications is completed; Step four, performing first superposition average calculation on the heartbeats in the same classification to obtain first superposition average results of multiple heartbeat classifications; Step five, calculating the correlation between the first superposition average results of multiple heartbeat classifications, merging the heartbeat classifications with high correlation and performing second superposition average calculation to obtain second superposition average results of the heartbeat classifications.

2. The method for optimizing the magnetocardiographic signal stack average based on heart beat classification according to claim 1, wherein, The step one comprises: eliminating the baseline drift, high frequency noise and power frequency noise of the magnetocardiogram signal by using median filtering, low pass filtering and band stop filtering, separating the noise component and signal component of the magnetocardiogram signal by using the ICA method, setting the noise component to zero and reconstructing to obtain the denoised magnetocardiogram signal.

3. The method for optimizing the magnetocardiography signal stack average based on heart beat classification according to claim 1, wherein, In step four, the first superposition average calculation refers to adding all the heartbeat values in each classification and then dividing by the number of heartbeats in the classification to obtain multiple first superposition average results corresponding to multiple heartbeat classifications.

4. The method for optimizing the magnetocardiography signal stack average based on heart beat classification according to claim 1, wherein, The step five comprises: calculating the Pearson correlation matrix for the multiple first superposition average results, merging the classifications with correlation greater than the second threshold value into one classification, and finally performing second superposition average calculation on the new classification to obtain the second superposition average result.

Images