心电图预处理、分割和获取平均 P-QRS-T 片段附matlab代码

classdef MBassiouni properties filePath; frequency; loadedSamples; plotFlag; signal; time;preproc;orgsign; baseLineCorrectedSignal; waveletBaseLineCorrectedSignal; bandStopWaveletBaseLineCorrectedSignal; lowPassBandStopWaveletBaseLineCorrectedSignal; smoothedLowPassBandStopWaveletBaseLineCorrectedSignal; left; right; R; PLocation; QLocation; RLocation; SLocation; TLocation; PValue; QValue; RValue; SValue; TValue; RQValue; RSValue; PRLocation; RTLocation; QSLocation; RMean; RQValueMean; RSValueMean; PRLocationMean; RTLocationMean; QSLocationMean; RRLocation; RRLocationMean; RRLocationMedian; RRLocationThreshold; PQRSTFragments; PQRSTFragmentMean; feature; end methods %% function [ECG] = MBassiouni(filePath, frequency, loadedSamples, plotFlag) ECG.filePath = filePath; ECG.frequency = frequency; ECG.loadedSamples = loadedSamples; ECG.plotFlag = plotFlag; % [ECG] = loadSignaldata(ECG); [ECG] = loadSignal(ECG); ECG.orgsign = ECG.signal; [ECG] = preProcessing(ECG); %[ECG.left, ECG.right, ECG.time] = FarukUYSAL(ECG.signal, ECG.frequency); [ECG] = FarukUYSAL(ECG); [ECG] = detectPeaks(ECG); [ECG] = allign(ECG); [ECG] = checkR(ECG); if(plotFlag == 1) plotPeaks(ECG); end [ECG] = plotFragments(ECG); end %% function [ECG] = loadSignal(ECG) load(ECG.filePath,'val'); ECG.signal = val(1,ECG.loadedSamples); ECG.orgsign = val(1,ECG.loadedSamples); figure; plot(ECG.signal,'b'); end %% function [ECG]=preProcessing(ECG) % plot(ECG.signal); % Base Line Correction [C,L]=wavedec(ECG.signal,9,'db8'); cA=appcoef(C,L,'db8',9); l=length(cA); C(1:l)=0; ECG.baseLineCorrectedSignal=waverec(C,L,'db8'); sorh = 's'; % Specified soft or hard thresholding thrSettings = 4.294226749492056; % wavelet filter ECG.waveletBaseLineCorrectedSignal = cmddenoise(ECG.baseLineCorrectedSignal,'db8',9,sorh,NaN,thrSettings); ECG.preproc = ECG.waveletBaseLineCorrectedSignal; % band stop filter % Design a filter with a Q-factor of Q=35 to remove a 60 Hz tone from % system running at 300 Hz. Wo = 50/(500/2); Bw = Wo/3; [b,a] = iirnotch(Wo,Bw); ECG.bandStopWaveletBaseLineCorrectedSignal=filter(b,a,ECG.waveletBaseLineCorrectedSignal); % low pass filter hFs=500/2; Wp=40/hFs; Ws=60/hFs; [Lp_n,Lp_Wn] = buttord(Wp,Ws,0.1,30); [b,a] = butter(Lp_n,Lp_Wn); % freqz(b,a); ECG.lowPassBandStopWaveletBaseLineCorrectedSignal=filter(b,a,ECG.bandStopWaveletBaseLineCorrectedSignal); % smooth filter ECG.smoothedLowPassBandStopWaveletBaseLineCorrectedSignal = smooth(ECG.lowPassBandStopWaveletBaseLineCorrectedSignal,5); % Transposed ECG.signal = ECG.smoothedLowPassBandStopWaveletBaseLineCorrectedSignal'; ECG.preproc = ECG.signal; figure; plot(ECG.signal,'b'); end %% function [ECG] = FarukUYSAL(ECG) x1 = ECG.signal; fs = ECG.frequency; %{ QRS Detection Example shows the effect of each filter according to Pan-Tompkins algorithm. Note that, the decision algorithm is different than the mentioned algorithm. by Faruk UYSAL %} N = length(x1); % Signal length t = (0:(N - 1)) / fs; % time index % CANCELLATION DC DRIFT AND NORMALIZATION x1 = x1 - mean(x1); % cancel DC conponents x1 = x1 / max(abs(x1)); % normalize to one % LOW PASS FILTERING % LPF (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_LP = filter(b, a, [1 zeros(1, 12)]); % transfer function of LPF x2 = conv (x1 ,h_LP); x2 = x2 (6 + (1:N)); % cancle delay x2 = x2 / max(abs(x2)); % normalize , for convenience . % HIGH PASS FILTERING % HPF = Allpass-(Lowpass) = z^-16-[(1-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_HP=filter(b, a, [1 zeros(1, 32)]); % impulse response iof HPF x3 = conv (x2, h_HP); x3 = x3 (16 + (1:N)); %cancle delay x3 = x3 / max(abs(x3)); % DERIVATIVE FILTER % MECG impulse response h = [-1 -2 0 2 1]/8; % Apply filter x4 = conv (x3 ,h); x4 = x4 (2+(1:N)); x4 = x4/ max( abs(x4 )); % SQUARING x5 = x4 .^2; x5 = x5 / max(abs(x5)); % MOVING WINDOW INTEGRATION % MECG impulse response h = ones(1, 31)/31; % Delay = 15; % Delay in samples % Apply filter x6 = conv(x5, h); x6 = x6 (15 + (1:N)); x6 = x6 / max(abs(x6)); % FIND QRS POINTS WHICH IT IS DIFFERENT THAN PAN-TOMPKINS ALGORITHM maxHeight = max(x6); thresh = mean(x6); poss_reg = (x6 > thresh * maxHeight)'; left = find(diff([0 poss_reg']) == 1); right = find(diff([poss_reg' 0]) == -1); ECG.left = left; ECG.right = right; ECG.time = t; end %% function [ECG]=detectPeaks(ECG) ECG.PLocation = zeros(1, length(ECG.left)); ECG.QLocation = zeros(1, length(ECG.left)); ECG.RLocation = zeros(1, length(ECG.left)); ECG.SLocation = zeros(1, length(ECG.left)); ECG.TLocation = zeros(1, length(ECG.left)); ECG.PValue = zeros(1, length(ECG.left)); ECG.QValue = zeros(1, length(ECG.left)); ECG.RValue = zeros(1, length(ECG.left)); ECG.SValue = zeros(1, length(ECG.left)); ECG.TValue = zeros(1, length(ECG.left)); ECG.R = zeros(1, length(ECG.left)); for i=1:length(ECG.left) try window = ECG.signal((ECG.left(i) - 60):ECG.left(i)); [ECG.PValue(i), ECG.PLocation(i)] = max(window); ECG.PLocation(i) = ECG.PLocation(i) + ECG.left(i) - length(window); catch indexExceedsMatrixDimensions end window = ECG.signal(ECG.left(i):ECG.right(i)); [ECG.RValue(i), ECG.RLocation(i)] = max(window); ECG.RLocation(i) = ECG.RLocation(i) + ECG.left(i) - 1; ECG.R(i) = length(window); window = ECG.signal(ECG.left(i):ECG.RLocation(i)); [ECG.QValue(i), ECG.QLocation(i)] = min(window); ECG.QLocation(i) = ECG.QLocation(i) + ECG.left(i) - 1; window = ECG.signal(ECG.RLocation(i):ECG.right(i)); [ECG.SValue(i