matlabpan_tompkin算法
function [qrs_amp_raw,qrs_i_raw,delay]=pan_tompkin(ecg,fs)
Complete implementation of Pan-Tompkins algorithm
完成PT算法的实现
function [qrs_amp_raw,qrs_i_raw,delay]=pan_tompkin(ecg,fs,gr)
%% function [qrs_amp_raw,qrs_i_raw,delay]=pan_tompkin(ecg,fs)
% Complete implementation of Pan-Tompkins algorithm
%完成PT算法的实现
%% Inputs
% ecg : raw ecg vector signal 1d signal
% fs : sampling 200Hz, 400Hz and etc
% gr : flag to plot or not plot (t it 1 to have a plot or t it zero not
% to e any plots
%ecg：原始⼼电图向量
%fs：采样频率如200Hz、400Hz等
%gr：标记绘图与否（设置为1绘图，设置为0不查看任何绘图）
%% Outputs
% qrs_amp_raw : amplitude of R waves amplitudes
% qrs_i_raw : index of R waves
% delay : number of samples which the signal is delayed due to the
% filtering
% qrs_amp_raw :R波幅度
% qrs_i_raw : R波位置
% delay ：由于滤波⽽使信号延迟的样本数
%% Method :
%% PreProcessing
% 1) Signal is preprocesd , if the sampling frequency is higher then it is downsampled
% and if it is lower upsampled to make the sampling frequency 200 Hz
% with the same filtering tups introduced in Pan
% tompkins paper (a combination of low pass and high pass filter 5-15 Hz)
% to get rid of the baline wander and muscle noi.
%信号预处理，如果采样频率较⾼，则向下采样，如果采样频率较低，则向上采样，使采样频率为200Hz，与PT算法论⽂中介绍的滤波设置相同（低通和⾼通滤波器的组合，5-15Hz），以消除基线漂移和肌电⼲扰% 2) The filtered signal is derivated using a derivating filter to high light the QRS compl
ex.
%滤波后信号由⼀个导数滤波器导出，⽤以⾼亮QRS波(微分)
% 3) Signal is squared.
%信号平⽅
% 4)Signal is averaged with a moving window to get rid of noi (0.150 conds length).
%信号与⼀个移动窗⼝进⾏平均，以消除噪声（0.15秒的长度）
% 5) depending on the sampling frequency of your signal the filtering
%options are changed to best match the characteristics of your ecg signal
%根据信号的采样频率，滤波选项更改为和⼼电信号最佳匹配的特征
% 6) Unlike the other implementations in this implementation the desicion
% rule of the Pan tompkins is implemented completely.
%% Decision Rule
% At this point in the algorithm, the preceding stages have produced a roughly pul-shaped
% waveform at the output of the MWI . The determination as to whether this pul
% corresponds to a QRS complex (as oppod to a high-sloped T-wave or a noi artefact) is
% performed with an adaptive thresholding operation and other decision
% rules outlined below;
% a) FIDUCIAL MARK - The waveform is first procesd to produce a t of weighted unit
% samples at the location of the MWI maxima. This is done in order to localize the QRS
% complex to a single instant of time. The w[k] weighting is the maxima value.
% b) THRESHOLDING - When analyzing the amplitude of the MWI output, the algorithm us
% two threshold values (THR_SIG and THR_NOISE, appropriately initialized during a brief
% 2 cond training pha) that continuously adapt to changing ECG signal quality. The
% first pass through y[n] us the thresholds to classify the each non-zero sample
% (CURRENTPEAK) as either signal or noi:
% If CURRENTPEAK > THR_SIG, that location is identified as a QRS complex
% candidate?and the signal level (SIG_LEV) is updated:
% SIG _ LEV = 0.125 CURRENTPEAK + 0.875?SIG _ LEV
% If THR_NOISE < CURRENTPEAK < THR_SIG, then that location is identified as a
% Noi peak?and the noi level (NOISE_LEV) is updated:
% NOISE _ LEV = 0.125CURRENTPEAK + 0.875?NOISE _ LEV
% Bad on new estimates of the signal and noi levels (SIG_LEV and NOISE_LEV,
% respectively) at that point in the ECG, the thresholds are adjusted as follows:
% THR _ SIG = NOISE _ LEV + 0.25 ?(SIG _ LEV-NOISE _ LEV )
% THR _ NOISE = 0.5?(THR _ SIG)
% The adjustments lower the threshold gradually in signal gments that are deemed to
% be of poorer quality.
% c) SEARCHBACK FOR MISSED QRS COMPLEXES - In the thresholding step above, if
% CURRENTPEAK < THR_SIG, the peak is deemed not to have resulted from a QRS
% complex. If however, an unreasonably long period has expired without an abovethreshold
% peak, the algorithm will assume a QRS has been misd and perform a
% archback. This limits the number of fal negatives. The minimum time ud to trigger
% a archback is 1.66 times the current R peak to R peak time period (called the RR
% interval). This value has a physiological origin - the time value between adjacent
% heartbeats cannot change more quickly than this. The misd QRS complex is assumed
% to occur at the location of the highest peak in the interval that lies between THR_SIG and
邓超演的电影% THR_NOISE. In this algorithm, two average RR intervals are stored,the first RR interval is
% calculated as an average of the last eight QRS locations in order to adapt to changing heart
% rate and the cond RR interval mean is the mean
% of the most regular RR intervals . The threshold is lowered if the heart rate is not regular
% to improve detection.
% d) ELIMINATION OF MULTIPLE DETECTIONS WITHIN REFRACTORY PERIOD - It is
% impossible for a legitimate QRS complex to occur if it lies within 200ms after a previously
% detected one. This constraint is a physiological one ?due to the refractory period during
% which ventricular depolarization cannot occur despite a stimulus[1]. As QRS complex
% candidates are generated, the algorithm eliminates such physically impossible events,
% thereby reducing fal positives.
% e) T WAVE DISCRIMINATION - Finally, if a QRS candidate occurs after the 200ms
% refractory period but within 360ms of the previous QRS, the algorithm determines
% whether this is a genuine QRS complex of the next heartbeat or an abnormally prominent
% T wave. This decision is bad on the mean slope of the waveform at that position. A slope of % less than one half that of the previous QRS complex is consistent with the slower
% changing behaviour of a T wave ?otherwi, it becomes a QRS detection.
% Extra concept : beside the points mentioned in the paper, this code also
% checks if the occured peak which is less than 360 mc latency has also a
% latency less than 0,5*mean_RR if yes this is counted as noi
% f) In the final stage , the output of R waves detected in smoothed signal is analyzed and double % checked with the help of the output of the bandpass signal to improve
% detection and find the original index of the real R waves on the raw ecg
% signal
%% References :
%[1]PAN.J, TOMPKINS. W.J,"A Real-Time QRS Detection Algorithm" IEEE
%TRANSACTIONS ON BIOMEDICAL ENGINEERING, VOL. BME-32, NO. 3, MARCH 1985. %% Author : Hooman Sedghamiz
% Linkoping university
% email : hoo792@student.liu.
% hooman.
% Any direct or indirect u of this code should be referenced
% Copyright march 2014
%%
if ~isvector(ecg)
error('ecg must be a row or column vector');
end
丹山赤水if nargin < 3
gr = 1; % on default the function always plots
end
ecg = ecg(:); % vectorize
%% Initialize
qrs_c =[]; %amplitude of R
qrs_i =[]; %index
SIG_LEV = 0;
nois_c =[];
nois_i =[];
delay = 0;
skip = 0; % becomes one when a T wave is detected
not_nois = 0; % it is not noi when not_nois = 1
lected_RR =[]; % Selected RR intervals
m_lected_RR = 0;
mean_RR = 0;
qrs_i_raw =[];
qrs_amp_raw=[];
r_back = 0;
test_m = 0;
SIGL_buf = [];
NOISL_buf = [];
THRS_buf = [];
SIGL_buf1 = [];
NOISL_buf1 = [];
THRS_buf1 = [];
%% Plot differently bad on filtering ttings
if gr
if fs == 200
figure, ax(1)=subplot(321);plot(ecg);axis tight;title('Raw ECG Signal');
el
figure, ax(1)=subplot(3,2,[1 2]);plot(ecg);axis tight;title('Raw ECG Signal');
end
end
%% Noi cancelation(Filtering) % Filters (Filter in between 5-15 Hz)
if fs == 200
%% Low Pass Filter H(z) = ((1 - z^(-6))^2)/(1 - z^(-1))^2
b = [1 0 0 0 0 0 -2 0 0 0 0 0 1];
a = [1 -2 1];
h_l = filter(b,a,[1 zeros(1,12)]);
ecg_l = conv (ecg ,h_l);
ecg_l = ecg_l/ max( abs(ecg_l));
delay = 6; %bad on the paper
微波炉可以热汤吗
if gr
ax(2)=subplot(322);plot(ecg_l);axis tight;title('Low pass filtered');
end
%% High Pass filter H(z) = (-1+32z^(-16)+z^(-32))/(1+z^(-1))
b = [-1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 32 -32 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1];
a = [1 -1];
h_h = filter(b,a,[1 zeros(1,32)]);
ecg_h = conv (ecg_l ,h_h);
咀嚼造句ecg_h = ecg_h/ max( abs(ecg_h));
delay = delay + 16; % 16 samples for highpass filtering
if gr
ax(3)=subplot(323);plot(ecg_h);axis tight;title('High Pass Filtered');
end
el
%% bandpass filter for Noi cancelation of other sampling frequencies(Filtering)
硕士几年制
f1=5; %cuttoff low frequency to get rid of baline wander
f2=15; %cuttoff frequency to discard high frequency noi
Wn=[f1 f2]*2/fs; % cutt off bad on fsv
N = 3; % order of 3 less processing
[a,b] = butter(N,Wn); %bandpass filtering
ecg_h = filtfilt(a,b,ecg);
ecg_h = ecg_h/ max( abs(ecg_h));
if gr
ax(3)=subplot(323);plot(ecg_h);axis tight;title('Band Pass Filtered');
end
end
%% derivative filter H(z) = (1/8T)(-z^(-2) - 2z^(-1) + 2z + z^(2))
h_d = [-1 -2 0 2 1]*(1/8);%1/8*fs
ecg_d = conv (ecg_h ,h_d);
ecg_d = ecg_d/max(ecg_d);
delay = delay + 2; % delay of derivative filter 2 samples
if gr
ax(4)=subplot(324);plot(ecg_d);axis tight;title('Filtered with the derivative filter');
end
%% Squaring nonlinearly enhance the dominant peaks
ecg_s = ecg_d.^2;
if gr
ax(5)=subplot(325);plot(ecg_s);axis tight;title('Squared');
end
%% Moving average Y(nt) = (1/N)[x(nT-(N - 1)T)+ x(nT - (N - 2)T)+...+x(nT)]
ecg_m = conv(ecg_s ,ones(1 ,round(0.150*fs))/round(0.150*fs));
delay = delay + 15;
if gr
ax(6)=subplot(326);plot(ecg_m);axis tight;title('Averaged with 30 samples length,Black noi,Green Adaptive Threshold,RED Sig Level,Red circles QRS adaptive threshold'); axis tight;
描写雪景的优美句子
end
%% Fiducial Mark
% Note : a minimum distance of 40 samples is considered between each R wave
% since in physiological point of view no RR wave can occur in less than
% 200 mc distance
[pks,locs] = findpeaks(ecg_m,'MINPEAKDISTANCE',round(0.2*fs));
%% initialize the training pha (2 conds of the signal) to determine the THR_SIG and THR_NOISE
THR_SIG = max(ecg_m(1:2*fs))*1/4; % 0.25 of the max amplitude
THR_NOISE = mean(ecg_m(1:2*fs))*1/2; % 0.5 of the mean signal is considered to be noi
SIG_LEV= THR_SIG;
NOISE_LEV = THR_NOISE;
%% Initialize bandpath filter threshold(2 conds of the bandpass signal)
THR_SIG1 = max(ecg_h(1:2*fs))*1/4; % 0.25 of the max amplitude
THR_NOISE1 = mean(ecg_h(1:2*fs))*1/2; %
SIG_LEV1 = THR_SIG1; % Signal level in Bandpasd filter
NOISE_LEV1 = THR_NOISE1; % Noi level in Bandpasd filter
%% Thresholding and online desicion rule
for i = 1 : length(pks)
%% locate the corresponding peak in the filtered signal
if locs(i)-round(0.150*fs)>= 1 && locs(i)<= length(ecg_h)
[y_i, x_i] = max(ecg_h(locs(i)-round(0.150*fs):locs(i)));
el
if i == 1
[y_i, x_i] = max(ecg_h(1:locs(i)));
r_back = 1;
elif locs(i)>= length(ecg_h)
[y_i, x_i] = max(ecg_h(locs(i)-round(0.150*fs):end));
end
end
%% update the heart_rate (Two heart rate means one the moste recent and the other lected)
if length(qrs_c) >= 9
diffRR = diff(qrs_i(end-8:end)); %calculate RR interval
mean_RR = mean(diffRR); % calculate the mean of 8 previous R waves interval
comp =qrs_i(end)-qrs_i(end-1); %latest RR
if comp <= 0.92*mean_RR || comp >= 1.16*mean_RR
% lower down thresholds to detect better in MVI
THR_SIG = 0.5*(THR_SIG);
%THR_NOISE = 0.5*(THR_SIG);
% lower down thresholds to detect better in Bandpass filtered
THR_SIG1 = 0.5*(THR_SIG1);
%THR_NOISE1 = 0.5*(THR_SIG1);
el
m_lected_RR = mean_RR; %the latest regular beats mean
end
end
%% calculate the mean of the last 8 R waves to make sure that QRS is not
% missing(If no R detected , trigger a arch back) 1.66*mean
if m_lected_RR
test_m = m_lected_RR; %if the regular RR availabe u it
elif mean_RR && m_lected_RR == 0
test_m = mean_RR;
el
test_m = 0;
end
if test_m
if (locs(i) - qrs_i(end)) >= round(1.66*test_m)% it shows a QRS is misd
[pks_temp,locs_temp] = max(ecg_m(qrs_i(end)+ round(0.200*fs):locs(i)-round(0.200*fs))); % arch back and locate the max in this interval
locs_temp = qrs_i(end)+ round(0.200*fs) + locs_temp -1; %location
if pks_temp > THR_NOISE
qrs_c = [qrs_c pks_temp];
qrs_i = [qrs_i locs_temp];
% find the location in filtered sig
if locs_temp <= length(ecg_h)
[y_i_t x_i_t] = max(ecg_h(locs_temp-round(0.150*fs):locs_temp));
el
[y_i_t x_i_t] = max(ecg_h(locs_temp-round(0.150*fs):end));
end
% take care of bandpass signal threshold
if y_i_t > THR_NOISE1
qrs_i_raw = [qrs_i_raw locs_temp-round(0.150*fs)+ (x_i_t - 1)];% save index of bandpass
qrs_amp_raw =[qrs_amp_raw y_i_t]; %save amplitude of bandpass
SIG_LEV1 = 0.25*y_i_t + 0.75*SIG_LEV1; %when found with the cond thres
end
not_nois = 1;
SIG_LEV = 0.25*pks_temp + 0.75*SIG_LEV ; %when found with the cond threshold
end
el
not_nois = 0;
end
end
%% find noi and QRS peaks
if pks(i) >= THR_SIG
% if a QRS candidate occurs within 360ms of the previous QRS
% ,the algorithm determines if its T wave or QRS
if length(qrs_c) >= 3
if (locs(i)-qrs_i(end)) <= round(0.3600*fs)
Slope1 = mean(diff(ecg_m(locs(i)-round(0.075*fs):locs(i)))); %mean slope of the waveform at that position Slope2 = mean(diff(ecg_m(qrs_i(end)-round(0.075*fs):qrs_i(end)))); %mean slope of previous R wave
if abs(Slope1) <= abs(0.5*(Slope2)) % slope less then 0.5 of previous R
nois_c = [nois_c pks(i)];
nois_i = [nois_i locs(i)];
skip = 1; % T wave identification
% adjust noi level in both filtered and
% MVI
NOISE_LEV1 = 0.125*y_i + 0.875*NOISE_LEV1;
NOISE_LEV = 0.125*pks(i) + 0.875*NOISE_LEV;
el
skip = 0;
end
end
end
if skip == 0 % skip is 1 when a T wave is detected
qrs_c = [qrs_c pks(i)];
qrs_i = [qrs_i locs(i)];
% bandpass filter check threshold
if y_i >= THR_SIG1
if r_back
qrs_i_raw = [qrs_i_raw x_i]; % save index of bandpass
el
qrs_i_raw = [qrs_i_raw locs(i)-round(0.150*fs)+ (x_i - 1)];% save index of bandpass
end
qrs_amp_raw =[qrs_amp_raw y_i];% save amplitude of bandpass
SIG_LEV1 = 0.125*y_i + 0.875*SIG_LEV1;% adjust threshold for bandpass filtered sig
end
% adjust Signal level
SIG_LEV = 0.125*pks(i) + 0.875*SIG_LEV ;
end
elif THR_NOISE <= pks(i) && pks(i)<THR_SIG
%adjust Noi level in filtered sig
NOISE_LEV1 = 0.125*y_i + 0.875*NOISE_LEV1;
%adjust Noi level in MVI
NOISE_LEV = 0.125*pks(i) + 0.875*NOISE_LEV;
elif pks(i) < THR_NOISE
nois_c = [nois_c pks(i)];
nois_i = [nois_i locs(i)];
% noi level in filtered signal
NOISE_LEV1 = 0.125*y_i + 0.875*NOISE_LEV1;
%end
%adjust Noi level in MVI
NOISE_LEV = 0.125*pks(i) + 0.875*NOISE_LEV;
end
%% adjust the threshold with SNR
if NOISE_LEV ~= 0 || SIG_LEV ~= 0
THR_SIG = NOISE_LEV + 0.25*(abs(SIG_LEV - NOISE_LEV));
THR_NOISE = 0.5*(THR_SIG);
end
% adjust the threshold with SNR for bandpasd signal
if NOISE_LEV1 ~= 0 || SIG_LEV1 ~= 0
THR_SIG1 = NOISE_LEV1 + 0.25*(abs(SIG_LEV1 - NOISE_LEV1));
THR_NOISE1 = 0.5*(THR_SIG1);
end
% take a track of thresholds of smoothed signal
SIGL_buf = [SIGL_buf SIG_LEV];
NOISL_buf = [NOISL_buf NOISE_LEV];
THRS_buf = [THRS_buf THR_SIG];
卡佛连% take a track of thresholds of filtered signal
SIGL_buf1 = [SIGL_buf1 SIG_LEV1];
NOISL_buf1 = [NOISL_buf1 NOISE_LEV1];
THRS_buf1 = [THRS_buf1 THR_SIG1];
skip = 0; %ret parameters
not_nois = 0; %ret parameters
r_back = 0; %ret bandpass param
end
if gr
hold on,scatter(qrs_i,qrs_c,'m');
hold on,plot(locs,NOISL_buf,'--k','LineWidth',2);
hold on,plot(locs,SIGL_buf,'--r','LineWidth',2);
hold on,plot(locs,THRS_buf,'--g','LineWidth',2);
if ax(:)
宽带连接错误651
linkaxes(ax,'x');
zoom on;
end
end
%% overlay on the signals
if gr
figure,az(1)=subplot(311);plot(ecg_h);title('QRS on Filtered Signal');axis tight;
hold on,scatter(qrs_i_raw,qrs_amp_raw,'m');
hold on,plot(locs,NOISL_buf1,'LineWidth',2,'Linestyle','--','color','k');
hold on,plot(locs,SIGL_buf1,'LineWidth',2,'Linestyle','-.','color','r');
hold on,plot(locs,THRS_buf1,'LineWidth',2,'Linestyle','-.','color','g');
az(2)=subplot(312);plot(ecg_m);title('QRS on MWI signal and Noi level(black),Signal Level (red) and Adaptive Threshold(green)');axis tight; hold on,scatter(qrs_i,qrs_c,'m');
hold on,plot(locs,NOISL_buf,'LineWidth',2,'Linestyle','--','color','k');
hold on,plot(locs,SIGL_buf,'LineWidth',2,'Linestyle','-.','color','r');
hold on,plot(locs,THRS_buf,'LineWidth',2,'Linestyle','-.','color','g');
az(3)=subplot(312);plot(ecg-mean(ecg));title('Pul train of the found QRS on ECG signal');axis tight;
line(repmat(qrs_i_raw+2,[2 1]),repmat([min(ecg-mean(ecg))/2; max(ecg-mean(ecg))/2],size(qrs_i_raw+2)),'LineWidth',2.5,'LineStyle','-.','Color','r'); linkaxes(az,'x');
zoom on;
end
end
