CN108113665B - A method for automatic noise reduction of ECG signals - Google Patents

Abstract

The invention discloses an automatic noise reduction algorithm for ECG signals. The processing process is: A) establishing an echo state network and initializing; B) acquiring human body ECG signals and constructing a training set on this basis; C) using recursion-based The online training algorithm of the least squares method trains the echo state network and obtains the trained echo state network; D) constructs a test set, and inputs the test set into the trained echo state network to obtain a clean ECG signal. After being processed by the method of the present invention, the clean ECG signal after denoising not only effectively filters out noise, but also restores the low-frequency characteristic wave of the ECG signal and retains the effective information of the ECG signal.

Description

A method for automatic noise reduction of ECG signals

technical field

The invention relates to a method for automatic detection and analysis of electrocardiographic signals, in particular to an automatic noise reduction method for electrocardiographic signals.

Background technique

Cardiovascular disease has become the number one killer of people's health, with high morbidity, disability and mortality. Among them, the low-frequency components of the ECG signal are very important for the intelligent diagnosis of cardiovascular diseases. For example, the ST segment, one of the low-frequency components, represents the potential change in a period of time after the end of ventricular depolarization and the beginning of ventricular repolarization. An important indicator for the diagnosis of myocardial ischemia. However, at the same time, the ECG signal itself has the common characteristics of the bioelectric signal: weak amplitude, low frequency, large impedance, randomness, etc. These characteristics make the ECG signal extremely susceptible to various noise interference during the acquisition process. Especially in the current telemedicine background, the noise during ECG signal acquisition presents more complex characteristics, which can easily cause the low-frequency part of the ECG signal to be overwhelmed by complex noise, and the phenomenon of loss of effective information occurs.

There are three commonly used methods: the first, empirical mode decomposition (EMD) is to extract the intrinsic mode function (IMF) to denoise the noisy ECG signal, and the result can remove 95% of the Gaussian white noise. However, EMD is not able to separate out signals with frequencies similar to ECG signals. Therefore, this noise reduction method has the potential to filter the P and T waves from the signal, leading to misdiagnosis. Second, Principal Component Analysis (PCA) and Independent Component Analysis (ICA) are both blind source separation techniques used to find independent source signals from non-Gaussian noise signals, which have obvious advantages for in-band filtering. However, its model is very sensitive to the subtle changes of the noisy ECG signal, and is not suitable for the situation that the ECG signal contains more complex noise. The third, wavelet denoising method has achieved good results in removing high-frequency noise, but the selection of the wavelet base and the setting of the threshold need to be determined through a large number of experiments, and have a direct impact on the final filtering effect. It can be seen that the existing algorithms tend to lose the effective low-frequency information of the ECG signal during the noise reduction process, which affects the doctor's judgment.

SUMMARY OF THE INVENTION

The purpose of the present invention is to provide an automatic noise reduction algorithm for electrical signals, so as to retain the low-frequency characteristic wave while realizing the noise reduction of the ECG signal, and overcome the problem that the effective low-frequency information of the ECG signal is easily lost in the noise reduction process of the existing algorithm. .

The object of the present invention is achieved in this way:

An automatic noise reduction algorithm for ECG signals, including the following steps:

A) establish the echo state network and initialize;

B) Obtain human ECG signals and build a training set on this basis;

C) utilize the online training algorithm based on recursive least squares to train the echo state network, obtain and save the trained echo state network;

D) Construct a test set with the ECG signal immediately after the training set, and input the test set into the trained echo state network to obtain a clean ECG signal.

The automatic noise reduction algorithm of described electrocardiogram, step A) comprises initialization weight matrix and network parameter initialization setting, wherein:

① The process of initializing the weight matrix is:

Randomly initialize the weight matrix W _in that connects the input unit to the inside of the reserve pool and the weight matrix W of the neuron connection inside the reserve pool, so that they all obey a uniform distribution;

Initialize the connection weight matrix W _back of the output unit and the reserve pool and the connection weight matrix W _out of the reserve pool and the output unit to a zero matrix;

②The network parameter initialization setting process is as follows:

The number of input neurons is set to 4, the number of output neurons is set to 1, and the number of neurons in the reserve pool is set to N=1000±300;

The sparsity of the reserve pool SD=m/N, where m is the number of interconnected neurons in the reserve pool, and the SD value ranges from 1 to 5%;

The value of the spectral radius SR of the reserve pool is debugged according to the step size of 0.1 in the range of 0 to 1, and the SR value corresponding to the optimal result is the final value, SR=max{abs(eigenvalue of W)} ;

The value of the input unit scale IS of the reserve pool is in the range of -1 to 1 for debugging with a step size of 0.01, and the IS value corresponding to the optimal result is the final value.

The automatic noise reduction algorithm of described electrocardiogram, step B) concrete process is as follows:

①Intercept several consecutive heartbeats x(n), n=1, 2,..., M, M=55±2 from the acquired human ECG signal data, and then use the noise-containing ECG signal x(n) at the current moment , the noisy ECG signal x(n-1) at the previous moment, the first derivative x′(n) of the noisy ECG signal at the current moment, and the second derivative x″ of the noisy ECG signal at the current moment (n) The input signal x _r (n) that constitutes the echo state network, expressed as

Define the clean heartbeat data corresponding to the noisy ECG signal x(n) at the current moment as d(n);

A training set consisting of _xr (n) and d(n) ( _xr (n), d(n), n=1, 2,...,M);

其中g(·)为输出神经元激活函数，选为恒等函数，即输出神经元为一个线性神经元，该函数中，矩阵v(n)定义为

其中r(n)为储备池状态变量；②Definition

Among them, g( ) is the activation function of the output neuron, which is selected as the identity function, that is, the output neuron is a linear neuron. In this function, the matrix v(n) is defined as

where r(n) is the state variable of the reserve pool;

③Define error e(n)=d(n)-y(n).

The automatic noise reduction algorithm of described electrocardiogram, step C) concrete process is:

① The initial value of the reserve pool state variable of the selected network is 0, that is, r(0)=0, the sample x in the training set (x _r (n), d(n), n=1, 2, ..., M) _r (n) passes through W _in , d(n) passes through W _back and is added to the reserve pool respectively, using the formula r(n)=(1-α)[f(W·r(n-1)+W _in ·u (n)+W _back d(n))]+α r(n-1) to update the state variable r(n) of the reserve pool;

where f( ) is the activation function of the neurons in the reserve pool, which is selected as the nonlinear tanh( ) function;

来进行更新，其中，

λ为遗忘因子，其取值在0.95～1范围以0.0001为步长进行调试，以结果最优时λ值为最终取值，ψ(0)＝δ·I，I为单位矩阵，δ为一个极小的正数，其取值在0～0.001范围内以0.00001为步长调试，以结果最优时δ值为最终取值；u(n)=ψ(n-1)·v(n), ψ(n) is the correlation coefficient matrix of v(n) and is obtained by formula

to update, where,

λ is the forgetting factor, and its value is in the range of 0.95 to 1 with 0.0001 as the step for debugging. When the result is optimal, the λ value is the final value, ψ(0)=δ·I, I is the identity matrix, and δ is a A very small positive number, its value is in the range of 0 to 0.001, with 0.00001 as the step for debugging, and the final value of δ when the result is optimal;

α is the forgetting rate, and its value is in the range of 0 to 1 for debugging with a step size of 0.1, and the α value when the result is optimal is the final value;

② Update the output weight W _out according to the formula W _out (n)=W _out (n-1)+k(n)·e(n);

③ By updating W _out (n) and r (n), y(n) is approximated to the clean ECG signal d(n), so that the error e(n) can be approximated to 0, and the trained echo state network is obtained.

The automatic noise reduction algorithm of described electrocardiogram, step D) concrete process is:

即其包含当前时刻的心电信号X(n)、其上一个心电信号X(n-1)、其一阶导数X′(n)和二阶导数X″(n)；The test set is constructed by using the T heart beats immediately after the M heart beats intercepted when constructing the training set, the test set (X _r (n), n=1, 2,..., T, T=200±10), wherein,

That is, it contains the current ECG signal X(n), the previous ECG signal X(n-1), its first derivative X'(n) and second derivative X"(n);

Input X _r (n) in the test set (X _r (n), n=1, 2, . . . , T) into the trained echo state network, and the network outputs a clean ECG signal Y (n).

After being processed by the method of the present invention, the clean ECG signal after denoising not only effectively filters out noise, but also restores the low-frequency characteristic wave of the ECG signal and retains the effective information of the ECG signal.

Description of drawings

Fig. 1 is a flow chart of the implementation process of the method of the present invention.

Figure 2 is a schematic diagram of an echo state network structure.

FIG. 3 is a schematic diagram of the waveform structure of the human ECG signal.

Figure 4 is a graph of a noisy ECG signal before filtering.

Figure 5 is a graph of the clean ECG signal after filtering.

Detailed ways

This embodiment is implemented in a computer with Intel Xeon CPU E5-2697@2.70GHz, memory of 128.00GB, Win7, 64-bit operating system, and the entire ECG signal automatic classification algorithm is implemented in Matlab language.

In conjunction with Fig. 1, the implementation process of the present invention is as follows:

A) Establish an Echo State Network (ESN):

①, initialize the weight matrix:

Randomly initialize the weight matrix W _in that the input unit is connected to the interior of the reserve pool and the weight matrix W that connect the neurons inside the reserve pool, and make them all obey a uniform distribution;

The connection weight matrix W _back of the output unit and the reserve pool, and the connection weight matrix W _out of the reserve pool and the output unit are initialized to a zero matrix;

②, parameter setting:

The number of input neurons of ESN is set to 4, the number of output neurons is set to 1, and the number of neurons in the reserve pool is the size of the reserve pool N=1000;

The sparse degree SD of the ESN’s reserve pool is calculated according to the formula SD=m/N, where m is the number of interconnected neurons in the reserve pool. In this example, m is taken as m=20, then SD=2%;

The spectral radius SR of the reserve pool is defined as SR=max{abs(eigenvalue of W)}. Since the ESN has the property of echo state only when SR<1, the value of SR in this example is based on the range of 0 to 1 in steps of 0.1 Debugging is carried out for a long time, and the optimal result is taken as SR=0.4;

Define W _S as a randomly generated sparse matrix with a uniform distribution and a sparse degree of SD, then the neuron connection weight matrix in the reserve pool W = SR · W _S to ensure that the spectral radius of W is SR;

The input unit scale IS of the reserve pool is a scale factor that needs to be multiplied before the input signal is connected to the reserve pool. By multiplying the input signal with it, the input signal is transformed into the input range of the neuron activation function of the reserve pool. The neuron activation function The input range is [-1, 1]. In this example, the value of IS is in the range of -1 to 1 and is debugged with a step size of 0.01. When the result is optimal, IS=[0.01; 0.01; 0.9; 0.8];

B) Collect the original ECG signal of the human body with the acquisition frequency of 250Hz, and its waveform structure is shown in Figure 3, and store it in the data form of the TXT file, and then use Matlab software to read the original ECG signal data stored in the TXT file into In the computer, as the original ECG signal, as shown in Figure 4, a training set is constructed:

为一个N×4的矩阵；定义与当前心拍相应的干净心拍数据为d(n)；由x_r(n)与d(n)组成训练集(x_r(n),d(n),n＝1,2,…,55)。Define x _r (n) as the input signal of ESN, which is to intercept M=55 consecutive heartbeats x(n), n=1, 2, . . . , 55 from the original ECG signal. The signal x(n), the noisy ECG signal x(n-1) at the previous moment, the first derivative x′(n) of the noisy ECG signal at the current moment, the The second derivative x″(n) is formed, which is expressed as

is an N×4 matrix; the clean heartbeat data corresponding to the current heartbeat is defined as d(n); the training set is composed of _xr (n) and d(n) ( _xr (n),d(n),n =1,2,...,55).

其中g(·)为输出神经元激活函数，选为恒等函数，即输出神经元为一个线性神经元，该函数中，矩阵v(n)定义为

其中r(n)为储备池的状态变量，diag(IS)表示以IS为对角线生成一个对角矩阵；definition

Among them, g( ) is the activation function of the output neuron, which is selected as the identity function, that is, the output neuron is a linear neuron. In this function, the matrix v(n) is defined as

Among them, r(n) is the state variable of the reserve pool, and diag(IS) means that a diagonal matrix is generated with IS as the diagonal;

Definition error e(n)=d(n)-y(n);

C) Using the online training algorithm based on recursive least squares to train the echo state network:

①The initial value of the state variable of the reserve pool of the selected network is 0, that is, r(0)=0, the samples in the training set (x _r (n), d(n), n=1, 2, ..., 55) x _r (n) passes through W _in , and d(n) passes through W _back to be added to the reserve pool respectively, using the formula r(n)=(1-α)[f(W·r(n-1)+W _in · u(n)+W _back d(n))]+α r(n-1) to update the reserve pool state r(n):

Among them, f( ) is the activation function of the neurons in the reserve pool. Since the network needs strong nonlinearity for modeling, it is selected as the nonlinear tanh( ) function;

In this function, u(n)=ψ ^-1 (n-1)·v(n), ψ(n) is the inverse matrix of the correlation coefficient matrix of v(n) and is updated by the following formula:

其中，

λ为遗忘因子，定义为一个范围为0.95＜λ＜1的常数，其取值以0.0001为步长，进行调试，以结果最优时λ值为最终取值，本例中取λ＝0.9999，ψ(0)＝δ^-1·I，I为单位矩阵，δ为一个极小的正数，其取值范围为0～0.001，以0.00001为步长调试，以结果最优(即训练集中的输入信号再放进网络中，得到的y(n)与d(n)的误差近似于0时)时δ值为最终取值，本例中取δ＝0.00172；

in,

λ is the forgetting factor, which is defined as a constant in the range of 0.95 < λ < 1. Its value takes 0.0001 as the step size, and is debugged. When the result is optimal, the λ value is the final value. In this example, λ = 0.9999, ψ(0)=δ ^-1 ·I, I is the identity matrix, δ is a very small positive number, its value range is 0~0.001, and the step size is 0.00001 for debugging, and the result is the best (that is, the When the input signal is put into the network again, and the obtained error between y(n) and d(n) is approximately 0), the δ value is the final value, in this example, δ=0.00172;

In this function, α is the forgetting rate, which is defined as a positive number less than 1. The value ranges from 0 to 1. The debugging is performed with a step size of 0.1. The final value is the α value when the result is optimal. In this example The value is α=0.8.

② The output weight W _out is updated according to the formula W _out (n)=W _out (n-1)+k(n)·e(n).

③ By updating W _out (n) and r (n), y(n) is approximated to the clean ECG signal d(n), so that the error e(n) can be approximated to 0, and the well-trained echo state network is obtained and saved.

When y(n) is close to d(n), that is, the error e(n) is close to 0, the network is in a stable state, and the trained back-to-state network can be used for ECG signal denoising.

D) Verify:

即其包含当前时刻的心电信号X(n)、其上一个心电信号X(n-1)、其一阶导数X′(n)和二阶导数X″(n)；The test set is constructed with T=200 heartbeats immediately after the 55 heartbeats intercepted when constructing the training set, and the test set is ( _Xr (n), n=1,2,...,200), wherein,

That is, it contains the current ECG signal X(n), the previous ECG signal X(n-1), its first derivative X'(n) and second derivative X"(n);

Input X _r (n) in the test set (X _r (n), n=1, 2,..., 200) into the echo state network after training, and the network outputs a clean ECG signal Y (n), as shown in Figure 5 Show.

It can be seen from the above processing results that the noise-containing ECG signal (as shown in Figure 4) is compared with the normal ECG signal waveform (as shown in Figure 2). Distortion, processed by the method of the present invention, the clean ECG signal after denoising (as shown in Figure 5) not only effectively filters out the noise, but also restores the low-frequency characteristic wave of the ECG signal and retains the effective information of the ECG signal.

Claims (1)

1. an automatic noise reduction algorithm of electrocardiogram, is characterized in that, comprises the following steps: A) establish the echo state network and initialize; B) Obtain human ECG signals and build a training set on this basis; C) utilize the online training algorithm based on recursive least squares to train the echo state network, obtain and save the trained echo state network; D) construct a test set with the ECG signal immediately after the training set, and input the test set into the trained echo state network to obtain a clean ECG signal; Described step A) includes initialization weight matrix and network parameter initialization setting, wherein: ① The process of initializing the weight matrix is: Randomly initialize the weight matrix W _in that connects the input unit to the inside of the reserve pool and the weight matrix W of the neuron connection inside the reserve pool, so that they all obey a uniform distribution; Initialize the connection weight matrix W _back of the output unit and the reserve pool and the connection weight matrix W _out of the reserve pool and the output unit to a zero matrix; ②The network parameter initialization setting process is as follows: The number of input neurons is set to 4, the number of output neurons is set to 1, and the number of neurons in the reserve pool is set to N=1000±300; The sparsity of the reserve pool SD=m/N, where m is the number of interconnected neurons in the reserve pool, and the SD value ranges from 1 to 5%; The value of the spectral radius SR of the reserve pool is debugged according to the step size of 0.1 in the range of 0 to 1, and the SR value corresponding to the optimal result is the final value, SR=max{abs(eigenvalue of W)} ; The value of the input unit scale IS of the reserve pool is in the range of -1 to 1, with a step size of 0.01 for debugging, and the IS value corresponding to the optimal result is the final value; Described step B) concrete process is as follows: ①Intercept several consecutive heartbeats x(n), n=1, 2,..., M, M=55±2 from the acquired human ECG signal data, and then use the noise-containing ECG signal x(n) at the current moment , the noisy ECG signal x(n-1) at the previous moment, the first derivative x′(n) of the noisy ECG signal at the current moment, and the second derivative x″ of the noisy ECG signal at the current moment (n) The input signal x _r (n) that constitutes the echo state network, expressed as

Define the clean heartbeat data corresponding to the noisy ECG signal x(n) at the current moment as d(n); A training set consisting of _xr (n) and d(n) ( _xr (n), d(n), n=1, 2,...,M); ②Definition

Among them, g( ) is the activation function of the output neuron, which is selected as the identity function, that is, the output neuron is a linear neuron. In this function, the matrix v(n) is defined as

Among them, r(n) is the state variable of the reserve pool, and diag(IS) means that a diagonal matrix is generated with IS as the diagonal; ③Define error e(n)=d(n)-y(n); Described step C) concrete process is: ① The initial value of the reserve pool state variable of the selected network is 0, that is, r(0)=0, the sample x in the training set (x _r (n), d(n), n=1, 2, ..., M) _r (n) passes through W _in , d(n) passes through W _back and is added to the reserve pool respectively, using the formula r(n)=(1-α)[f(W·r(n-1)+W _in ·u (n)+W _back d(n))]+α r(n-1) to update the state variable r(n) of the reserve pool; where f( ) is the activation function of the neurons in the reserve pool, which is selected as the nonlinear tanh( ) function; u(n)=ψ(n-1)·v(n), ψ(n) is the correlation coefficient matrix of v(n) and is obtained by formula

to update, where,

λ is the forgetting factor, and its value is in the range of 0.95 to 1 with 0.0001 as the step for debugging. When the result is optimal, the λ value is the final value, ψ(0)=δ·I, I is the identity matrix, and δ is a A very small positive number, its value is in the range of 0 to 0.001, with 0.00001 as the step for debugging, and the final value of δ when the result is optimal; α is the forgetting rate, and its value is in the range of 0 to 1 for debugging with a step size of 0.1, and the α value when the result is optimal is the final value; ② Update the output weight W _out according to the formula W _out (n)=W _out (n-1)+k(n)·e(n); ③ By updating W _out (n) and r (n), y(n) is approximated to the clean ECG signal d(n), so that the error e(n) can be approximated to 0, and the trained echo state network is obtained; Described step D) concrete process is: The test set is constructed by using the T heart beats immediately after the M heart beats intercepted when constructing the training set, the test set (X _r (n), n=1, 2,..., T, T=200±10), wherein,

That is, it contains the current ECG signal X(n), the previous ECG signal X(n-1), its first derivative X'(n) and second derivative X"(n); Input X _r (n) in the test set (X _r (n), n=1, 2, . . . , T) into the trained echo state network, and the network outputs a clean ECG signal Y (n).

Images