【信号处理】心电信号PQRST峰值检测工具箱.zip

function [qrs_amp_raw,qrs_i_raw,delay]=pan_tompkin2(ecg,fs) %#codegen %% function [qrs_amp_raw,qrs_i_raw,delay]=pan_tompkin(ecg,fs) % Complete implementation of Pan-Tompkins algorithm %% Inputs % ecg : raw ecg vector signal 1d signal % fs : sampling frequency e.g. 200Hz, 400Hz and etc % gr : flag to plot or not plot (set it 1 to have a plot or set it zero not % to see any plots %% 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 %% Method % See Ref and supporting documents on researchgate. % https://www.researchgate.net/publication/313673153_Matlab_Implementation_of_Pan_Tompkins_ECG_QRS_detector %% References : %[1] Sedghamiz. H, "Matlab Implementation of Pan Tompkins ECG QRS %detector.",2014. (See researchgate) %[2] PAN.J, TOMPKINS. W.J,"A Real-Time QRS Detection Algorithm" IEEE %TRANSACTIONS ON BIOMEDICAL ENGINEERING, VOL. BME-32, NO. 3, MARCH 1985. %% ============== Licensce ========================================== %% % THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS % "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT % LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS % FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT % OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, % SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED % TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR % PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF % LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING % NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS % SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. % Author : % Hooman Sedghamiz, Feb, 2018 % MSc. Biomedical Engineering, Linkoping University % Email : Hooman.sedghamiz@gmail.com %% ============ Update History ================== %% % Feb 2018 : % 1- Cleaned up the code and added more comments % 2- Added to BioSigKit Toolbox %% ================= Now Part of BioSigKit ==================== %% ecg = ecg(:); % vectorize %% ======================= Initialize =============================== % delay = 0; skip = 0; % becomes one when a T wave is detected m_selected_RR = 0; mean_RR = 0; ser_back = 0; %% ============ Noise cancelation(Filtering)( 5-15 Hz) =============== %% if fs == 200 % ------------------ remove the mean of Signal -----------------------% ecg = ecg - mean(ecg); %% ==== Low Pass Filter H(z) = ((1 - z^(-6))^2)/(1 - z^(-1))^2 ==== %% %%It has come to my attention the original filter doesnt achieve 12 Hz % b = [1 0 0 0 0 0 -2 0 0 0 0 0 1]; % a = [1 -2 1]; % ecg_l = filter(b,a,ecg); % delay = 6; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Wn = 12*2/fs; N = 3; % order of 3 less processing [a,b] = butter(N,Wn,'low'); % bandpass filtering ecg_l = filtfilt(a,b,ecg); ecg_l = ecg_l/ max(abs(ecg_l)); %% ======================= start figure ============================= %% %% ==== High Pass filter H(z) = (-1+32z^(-16)+z^(-32))/(1+z^(-1)) ==== %% %%It has come to my attention the original filter doesn achieve 5 Hz % b = zeros(1,33); % b(1) = -1; b(17) = 32; b(33) = 1; % a = [1 1]; % ecg_h = filter(b,a,ecg_l); % Without Delay % delay = delay + 16; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Wn = 5*2/fs; N = 3; % order of 3 less processing [a,b] = butter(N,Wn,'high'); % bandpass filtering ecg_h = filtfilt(a,b,ecg_l); ecg_h = ecg_h/ max(abs(ecg_h)); else %% bandpass filter for Noise cancelation of other sampling frequencies(Filtering) f1=5; % cuttoff low frequency to get rid of baseline wander f2=15; % cuttoff frequency to discard high frequency noise Wn=[f1 f2]*2/fs; % cutt off based on fs 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)); end %% ==================== derivative filter ========================== %% % ------ H(z) = (1/8T)(-z^(-2) - 2z^(-1) + 2z + z^(2)) --------- % if fs ~= 200 int_c = (5-1)/(fs*1/40); b = interp1(1:5,[1 2 0 -2 -1].*(1/8)*fs,1:int_c:5); else b = [1 2 0 -2 -1].*(1/8)*fs; end ecg_d = filtfilt(b,1,ecg_h); ecg_d = ecg_d/max(ecg_d); %% ========== Squaring nonlinearly enhance the dominant peaks ========== %% ecg_s = ecg_d.^2; %% ============ 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 + round(0.150*fs)/2; %% ===================== Fiducial Marks ============================== %% % 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 msec distance [pks,locs] = findpeaks(ecg_m,'MINPEAKDISTANCE',round(0.2*fs)); %% =================== Initialize Some Other Parameters =============== %% LLp = length(pks); % ---------------- Stores QRS wrt Sig and Filtered Sig ------------------% qrs_c = zeros(1,LLp); % amplitude of R qrs_i = zeros(1,LLp); % index qrs_i_raw = zeros(1,LLp); % amplitude of R qrs_amp_raw= zeros(1,LLp); % Index % ------------------- Noise Buffers ---------------------------------% nois_c = zeros(1,LLp); nois_i = zeros(1,LLp); % ------------------- Buffers for Signal and Noise ----------------- % SIGL_buf = zeros(1,LLp); NOISL_buf = zeros(1,LLp); SIGL_buf1 = zeros(1,LLp); NOISL_buf1 = zeros(1,LLp); THRS_buf1 = zeros(1,LLp); THRS_buf = zeros(1,LLp); %% initialize the training phase (2 seconds of the signal) to determine the THR_SIG and THR_NOISE THR_SIG = max(ecg_m(1:2*fs))*1/3; % 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 noise SIG_LEV= THR_SIG; NOISE_LEV = THR_NOISE; %% Initialize bandpath filter threshold(2 seconds of the bandpass signal) THR_SIG1 = max(ecg_h(1:2*fs))*1/3; % 0.25 of the max amplitude THR_NOISE1 = mean(ecg_h(1:2*fs))*1/2; SIG_LEV1 = THR_SIG1; % Signal level in Bandpassed filter NOISE_LEV1 = THR_NOISE1; % Noise level in Bandpassed filter %% ============ Thresholding and desicion rule ============= %% Beat_C = 0; % Raw Beats Beat_C1 = 0; % Filtered Beats Noise_Count = 0; % Noise Counter for i = 1 : LLp %% ===== 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))); else if i == 1 [y_i,x_i] = max(ecg_h(1:locs(i))); ser_back = 1; elseif locs(i)>= length(ecg_h) [y_i,x_i] = max(ecg_h(locs(i)-round(0.150*fs):end)); else [y_i,x_i] = max(ecg_h(locs(i)-round(0.150*fs):end)); e