CN110141244A - Electrocardiogram personal identification method - Google Patents

Abstract

The invention discloses an electrocardiogram identification method, which comprises: 1) importing data; 2) removing noise of electrocardiogram signals; 3) extracting R waves; 4) constructing features; 5) automatically identifying samples; 6) outputting classification results a . Compared with the prior art, the advantages of the present invention are: the heartbeat signal is used for identification, which effectively reduces the signal acquisition time required for identification and reduces the calculation cost. In this method, k-medoids clustering is used for the training set, The robustness of the method for identification of noisy signals is better, and the effect on small samples is more obvious.

Description

ECG identification method

technical field

The invention relates to the technical field of medical signal processing, more specifically an electrocardiogram identification method.

Background technique

Biometrics refers to the close combination of computer and high-tech means such as optics, acoustics, biosensors and biostatistics principles, using the inherent physiological characteristics of the human body (such as fingerprints, face, iris, brain waves, pulse, etc.). ) or behavioral characteristics (such as handwriting, voice, gait, etc.) for personal identity authentication. Biometric identification technology has the advantages of not being forgotten, not easily forged or stolen, portable and available anytime, anywhere, and is more secure, confidential and convenient than traditional identity authentication methods.

In recent years, the method of identification by means of ECG (Electrocardiogram) inherent in the human body has attracted wide attention. ECG is a potential signal collected from the human body surface that reflects the heartbeat. The differences in the physiological conditions of the human body make ECG have many individual characteristics. Compared with fingerprints, voices and palms, ECG, as a living biological signal, is easy to detect and difficult to replicate. The definition of using ECG signal for identification is: given an ECG signal, determine the identity of the person to which the signal belongs. A system or device for automatic identification using ECG signals is usually implemented based on pattern recognition technology.

In this case, the auxiliary equipment related to ECG is developing rapidly. With the advancement of science and technology in the field of information, especially with the progress of pattern recognition technology, the effect of using ECG equipment for identification in practical applications is getting better and better. .

The core of the ECG equipment is the ECG identification system. Most of the published identification technologies based on ECG signals adopt multi-feature point extraction schemes, which are complicated in operation and huge in computation. Therefore, the identification methods of heartbeat signals have begun to emerge. The main advantage of such methods is that some independent heartbeat signals are identified from the ECG signals, which improves the practical application effect of ECG signal identification.

SUMMARY OF THE INVENTION

The purpose of the present invention is to provide an electrocardiogram identification method in order to solve the problems of long time and high cost for identification using a long period of electrocardiogram.

ECG identification method, which includes:

1) Import data

Collect the single-lead ECG signals of N people from the outside, and collect two pieces of each person, N is greater than 50, and store it as an array variable S, with a total of 2N pieces. 2N labels of 2N ECG signals are stored as variable label_L;

2) Remove the noise of the ECG signal

Execute the wavelet decomposition and reconstruction algorithm for each array in the array variable S, and store the result as the array variable S1;

3) Extract the R wave

Find m R wave peak point positions in each array in the array variable S1, a total of 2N*m positions, take 99 points before each position, take 100 points after each position, and take a total of 200 points including the current position. point as an R wave. 2N*m positions are extracted to 2N*m R waves, which are sequentially placed in a matrix of 2N*m rows and 200 columns, and this variable is recorded as Y, as the input signal of the feature construction method; 2N*m waves in the corresponding variable Y The 2N*m labels of the R wave are stored as the variable label;

4) Build features

a. Divide the R waves in Y into training set and test set, take N*m R waves extracted from the first signal of N people as the training set, and N*m R waves extracted from the second signal as the training set. The test set is stored in a matrix of N*m rows and 200 columns, denoted as train and test; the labels corresponding to the training set and the test set are divided from the variable label, and stored in an array of length N*m respectively, The training set label is recorded as train_label, and the test set label is recorded as test_label; the training set train and test set test are divided into N*m*20 segments, each segment is 10 in length, and stored in N*m*20 rows and 10 columns respectively In the matrix of , denoted as train_f and test_f;

b. Perform k-medoids clustering on train_f, set the number of clusters to 50, get a matrix of size 50*10 after clustering, denoted as C1, transpose C1 to get a matrix of 10*50, denoted as C, Calculate the Euclidean distance for train_f and C, and store the calculation result in a matrix of N*m*20 rows and 50 columns, denoted as train_sample_f; calculate the Euclidean distance for test_f and C, and store the calculation results in N*m*20 rows and 50 columns In the matrix of , denoted as test_sample_f;

c. Transform the training set sample train_sample_f and the test set sample test_sample_f into a matrix of N*m rows and 1000 columns according to the principle of column priority, and normalize them to [0,1], denoted as train_sample and test_sample;

train_sample is the constructed test set feature, and test_sample is the constructed test set feature;

5) Automatically identify samples

Sparse the training set feature train_sample, and input all the sparsed training set features and the corresponding label train_label into the LR classifier for training. After the LR classifier is trained, keep the parameters, and then sparse the test set feature test_sample. The sparse test set features and the corresponding label test_label are input to the trained LR classifier for testing;

The recognition result is the judgment accuracy rate a for judging the individual to which the R wave belongs;

6) Output the classification result a.

The wavelet decomposition and reconstruction algorithm is:

a. Wavelet decomposition: choose Harr function as the wavelet basis. Each array in the array variable S is decomposed to obtain approximate coefficient vector cA1, detail coefficient vector cD1; decompose cA1 to obtain approximate coefficient vector cA2, detail coefficient vector cD2; decompose cA2 to obtain approximate coefficient vector cA3, detail coefficient vector cD3; decompose cA3 Obtain approximate coefficient vector cA4, detail coefficient vector cD4; decompose cA4 to obtain approximate coefficient vector cA5, detail coefficient vector cD5; decompose cA5 to obtain approximate coefficient vector cA6, detail coefficient vector cD6; decompose cA6 to obtain approximate coefficient vector cA7, detail coefficient Vector cD7; decompose cA7 to obtain approximate coefficient vector cA8, detail coefficient vector cD8;

b. Remove noise: For the signal cD2 that needs to be removed, use the soft threshold method to remove noise. The thresholds required for noise reduction of each layer coefficient are different, and STDC is used as the noise intensity estimation;

STDC formula: stdc2 = median(|cD2|)/0.6745, where median(M) is the mean of calculating M;

Calculate the threshold for each layer:

thr2=(2*log(length(cD2)))^(1/2)*stdc2, where length(cD2) is the length of the calculated detail vector cD2;

For the signal cD2 that needs to be de-noised, set the coefficients whose modulus value is less than 3* thr2 in cD2 to zero, and process the coefficients whose modulus value is greater than 3* thr2: uniformly subtract 3* thr2 from the coefficients greater than 3* thr2, and less than -3 * The coefficient of thr2 is uniformly added 3* thr2

After the calculation, the denoised signal cd2 is obtained;

The above soft threshold method is also used to remove noise for cD3, cD4, cD5, cD6, cD7, and cD8 to obtain cd3, cd4, cd5, cd6, cd7, and cd8.

c. Wavelet reconstruction: set the coefficient vector cA8 obtained by decomposition to zero and assign it to ca8, complete wavelet reconstruction for cd8 to obtain signal ca7; complete wavelet reconstruction for cd7 to obtain signal ca6; complete wavelet reconstruction for cd6 to obtain signal ca5; complete wavelet reconstruction for cd6 to obtain signal ca5; Complete wavelet reconstruction to obtain signal ca4; complete wavelet reconstruction for ca4 to obtain signal ca3; complete wavelet reconstruction for cd3 to obtain signal ca2; complete wavelet reconstruction for cd2 to obtain signal ca1; ca1 is an array in S to perform wavelet decomposition and reconstruction the result of the algorithm.

Compared with the prior art, the present invention has the advantages that the heartbeat signal is used for identification, which effectively reduces the signal acquisition time required for identification and reduces the calculation cost. In this method, k-medoids clustering is used for the training set, The robustness of the method for identification of noisy signals is better, and the effect on small samples is more obvious.

Detailed ways

Embodiment 1 A kind of electrocardiogram identification method

The present invention will be further described below in conjunction with specific embodiments.

In the specific example, the international electrocardiogram database ECG-ID is used, and the data and usage instructions of the database are published on the well-known physionet.org website in the industry. The ECG-ID dataset contains 310 recordings obtained from 90 subjects (44 males and 46 females). Each recording contains the signal of one lead, recorded over a 30 s duration, and digitized at 500 Hz with 12-bit resolution over the nominal ±10 mV range. In this example, a total of 55 people's single-lead ECG signals are used, and each person collects two signals through a software system working on a computer and a well-known Matlab simulation environment in the industry.

The detailed steps of this embodiment are as follows:

one. Import Data

The single-lead ECG signals of 55 people are collected from the outside, and each person collects two pieces of ECG signals, which are stored as an array variable S, and the array variable S stores a total of 110 ECG signals. The 110 labels corresponding to the 110 ECG signals in the array variable S are stored as variable label_L.

two. Remove noise from ECG signals

The wavelet decomposition and reconstruction algorithm is performed on each array in the array variable S, and the result is stored as an array variable S1, and each array in S1 corresponds to the result of the wavelet decomposition and reconstruction algorithm performed on an array in S. The wavelet decomposition and reconstruction algorithm is:

(1) Wavelet decomposition

Select the Harr function as the wavelet basis. Decompose each array in the array variable S to obtain approximate coefficient vector cA1, detail coefficient vector cD1; decompose cA1 to obtain approximate coefficient vector cA2, detail coefficient vector cD2; decompose cA2 to obtain approximate coefficient vector cA3, detail coefficient vector cD3; Decompose to obtain approximate coefficient vector cA4, detail coefficient vector cD4; decompose cA4 to obtain approximate coefficient vector cA5, detail coefficient vector cD5; decompose cA5 to obtain approximate coefficient vector cA6, detail coefficient vector cD6; decompose cA6 to obtain approximate coefficient vector cA7, detail coefficient vector Coefficient vector cD7; decompose cA7 to obtain approximate coefficient vector cA8 and detail coefficient vector cD8.

(2) remove noise

For the signal cD2 that needs to be de-noised, a soft threshold method is used to de-noise. The thresholds required for the noise reduction of the coefficients of each layer are different, and STDC is used as the noise intensity estimation

STDC formula: stdc2 = median(|cD2|)/0.6745, where median(M) is the mean of calculating M

Calculate the threshold for each layer:

thr2=(2*log(length(cD2)))^(1/2)*stdc2, where length(cD2) is the length of the calculated detail vector cD2

For the signal cD2 that needs to be de-noised, set all coefficients in cD2 whose modulus is less than 3* thr2 to zero, and do a special treatment for those whose modulus is greater than 3* thr2. The coefficients greater than 3* thr2 are uniformly subtracted by 3* thr2, Coefficients less than -3* thr2 are uniformly added to 3* thr2

After the calculation, the denoised signal cd2 is obtained.

The above soft threshold method is also used to remove noise for cD3, cD4, cD5, cD6, cD7, and cD8. Get cd3, cd4, cd5, cd6, cd7, cd8.

(3) Wavelet reconstruction

The coefficient vector cA8 obtained by decomposing is directly set to zero and assigned to ca8, and the wavelet reconstruction of ca8 and cd8 is completed to obtain the signal ca7; the wavelet reconstruction of ca7 and cd7 is completed to obtain the signal ca6; the wavelet reconstruction of ca6 and cd6 is completed to obtain the signal ca5 ; Complete wavelet reconstruction for ca5, cd5 to obtain signal ca4; complete wavelet reconstruction for ca4, cd4 to obtain signal ca3; complete wavelet reconstruction for ca3, cd3 to obtain signal ca2; complete wavelet reconstruction for ca2, cd2 to obtain signal ca1. ca1 is the result of performing the wavelet decomposition and reconstruction algorithm for an array in S.

three. Extract R waves

(1) Find 13 R wave peak points in each array in the array variable S1 (the ECG database ECG-ID has only 13 R wave peak points for individual single-lead ECG signals), a total of 1430 positions. Take 99 points before each position and 100 points after each position, including the current position, take a total of 200 points as an R wave, extract 1430 R waves from 1430 positions, and put these R waves into 1430 lines in turn , 200-column matrix, name the variable Y.

(2) Y will be used as the input signal to construct the feature method. The 1430 labels corresponding to the 1430 R waves in variable Y are stored as variable label.

Four. build feature

The input signal of the feature building method is Y, and Y stores 1430 R waves extracted from 110 ECG signals.

(1) Divide the R waves in Y into the training set and the test set, and use the 715 R waves extracted from each person's first signal as the training set, which are stored in a matrix of 715 rows and 200 columns and named as train, The 715 R waves extracted from each person's second signal are stored in a matrix of 715 rows and 200 columns and named test. The labels corresponding to the training set and the test set are divided from the variable label, and the training set labels are stored in an array of length 715, named train_label. The obtained test set label is stored in an array of length 715 named test_label. Divide the training set train into 14300 segments, each segment is 10 in length, stored in a matrix of 14300 rows and 10 columns, named train_f; the test set test is divided into 14300 segments, each segment is 10 in length, stored in 14300 rows, 10-column matrix, named test_f.

(2) Perform k-medoids clustering on train_f, the number of clusters is 50, and a matrix of size 50*10 is obtained after clustering, named C1. Transpose C1 to get a 10*50 matrix, named C, and calculate the Euclidean distance for train_f and C. The calculation result is stored in a matrix of 14300 rows and 50 columns, named train_sample_f; calculate the Euclidean distance for test_f and C. The calculation result Stored in a matrix of 14300 rows and 50 columns, named test_sample_f.

The usage of kmedoids in matlab is:

[C1,~] = kmedoids(train_f,50)

The dist function in matlab calculates the Euclidean distance:

train_sample_f=dist(train_f,C);

(3) Transform the training set sample train_sample_f into a matrix of 715 rows and 1000 columns according to the principle of column priority and normalize it (implemented with the mapminmax function in matlab, the same below) to [0,1], named as train_sample; transform the test set sample train_sample_f into a matrix with 715 rows and 1000 columns according to the principle of column priority and normalize it to [0,1], named test_sample.

train_sample is the constructed test set feature, and test_sample is the constructed test set feature.

five. Automatic identification of samples

Sparse the training set feature train_sample, and input all the sparsed training set features and the corresponding label train_label into the LR classifier for training. After the LR classifier is trained, keep the parameters, and then sparse the test set feature test_sample. The sparse test set features and the corresponding label test_label are input to the trained LR classifier for testing.

After the test, the network will automatically give the recognition result. The recognition result is the judgment accuracy rate a for judging the individual to which the R wave belongs.

6. Output the classification result a.

The results are shown in Table 1

.

Claims (2)

1. An electrocardiogram identification method, which includes: 1) Import data Collect the single-lead ECG signals of N people from the outside, collect two pieces per person, N is greater than 50, and store them as an array variable S, a total of 2N pieces; 2N labels of 2N ECG signals are stored as variable label_L; 2) Remove the noise of the ECG signal Execute the wavelet decomposition and reconstruction algorithm for each array in the array variable S, and store the result as the array variable S1; 3) Extract the R wave Find m R wave peak point positions in each array in the array variable S1, a total of 2N*m positions, take 99 points before each position, take 100 points after each position, and take a total of 200 points including the current position. point as an R wave; 2N*m positions are extracted to 2N*m R waves, which are sequentially placed in a matrix of 2N*m rows and 200 columns, and this variable is recorded as Y, as the input signal of the feature construction method; 2N*m waves in the corresponding variable Y The 2N*m labels of the R wave are stored as the variable label; 4) Build features a. Divide the R waves in Y into training set and test set, take N*m R waves extracted from the first signal of N people as the training set, and N*m R waves extracted from the second signal as the training set. The test set is stored in a matrix of N*m rows and 200 columns, denoted as train and test; the labels corresponding to the training set and the test set are divided from the variable label, and stored in an array of length N*m respectively, The training set label is recorded as train_label, and the test set label is recorded as test_label; the training set train and test set test are divided into N*m*20 segments, each segment is 10 in length, and stored in N*m*20 rows and 10 columns respectively In the matrix of , denoted as train_f and test_f; b. Perform k-medoids clustering on train_f, set the number of clusters to 50, get a matrix of size 50*10 after clustering, denoted as C1, transpose C1 to get a matrix of 10*50, denoted as C, Calculate the Euclidean distance for train_f and C, and store the calculation results in a matrix of N*m*20 rows and 50 columns, denoted as train_sample_f; calculate the Euclidean distance for test_f and C, and store the calculation results in N*m*20 rows and 50 columns In the matrix of , denoted as test_sample_f; c. Transform the training set sample train_sample_f and the test set sample test_sample_f into a matrix of N*m rows and 1000 columns according to the principle of column priority, and normalize them to [0,1], denoted as train_sample and test_sample; train_sample is the constructed test set feature, and test_sample is the constructed test set feature; 5) Automatically identify samples Sparse the training set feature train_sample, and input all the sparsed training set features and the corresponding label train_label into the LR classifier for training. After the LR classifier is trained, keep the parameters, and then sparse the test set feature test_sample. The sparse test set features and the corresponding label test_label are input to the trained LR classifier for testing; The recognition result is the judgment accuracy rate a for judging the individual to which the R wave belongs; 6) Output the classification result a. 2. a kind of electrocardiogram identification method according to claim 1, is characterized in that, described wavelet decomposition reconstruction algorithm is: a. Wavelet decomposition: choose Harr function as the wavelet basis; Each array in the array variable S is decomposed to obtain approximate coefficient vector cA1, detail coefficient vector cD1; decompose cA1 to obtain approximate coefficient vector cA2, detail coefficient vector cD2; decompose cA2 to obtain approximate coefficient vector cA3, detail coefficient vector cD3; decompose cA3 Obtain approximate coefficient vector cA4, detail coefficient vector cD4; decompose cA4 to obtain approximate coefficient vector cA5, detail coefficient vector cD5; decompose cA5 to obtain approximate coefficient vector cA6, detail coefficient vector cD6; decompose cA6 to obtain approximate coefficient vector cA7, detail coefficient Vector cD7; decompose cA7 to obtain approximate coefficient vector cA8, detail coefficient vector cD8; b. Noise removal: For the signal cD2 that needs to be denoised, use the soft threshold method to denoise; The thresholds required for noise reduction of each layer coefficient are different, and STDC is used as the noise intensity estimation; STDC formula: stdc2 = median(|cD2|)/0.6745, where median(M) is the mean of calculating M; Calculate the threshold for each layer: thr2=(2*log(length(cD2)))^(1/2)*stdc2, where length(cD2) is the length of the calculated detail vector cD2; For the signal cD2 that needs to be de-noised, set the coefficients whose modulus value is less than 3* thr2 in cD2 to zero, and process the coefficients whose modulus value is greater than 3* thr2: uniformly subtract 3* thr2 from the coefficients greater than 3* thr2, and less than -3 * The coefficient of thr2 is uniformly added 3* thr2 After the calculation, the denoised signal cd2 is obtained; For cD3, cD4, cD5, cD6, cD7, and cD8, the above soft threshold method is also used to remove noise, and cd3, cd4, cd5, cd6, cd7, cd8 are obtained; c. Wavelet reconstruction: set the coefficient vector cA8 obtained by decomposition to zero and assign it to ca8, complete wavelet reconstruction for cd8 to obtain signal ca7; complete wavelet reconstruction for cd7 to obtain signal ca6; complete wavelet reconstruction for cd6 to obtain signal ca5; complete wavelet reconstruction for cd6 to obtain signal ca5; Complete wavelet reconstruction to obtain signal ca4; complete wavelet reconstruction for ca4 to obtain signal ca3; complete wavelet reconstruction for cd3 to obtain signal ca2; complete wavelet reconstruction for cd2 to obtain signal ca1; ca1 is an array in S to perform wavelet decomposition and reconstruction the result of the algorithm.