3 views (last 30 days)

Hey!

My qeustion is: I need to remove all the R peaks. I have all the Q and S junctions points and I thought about connecting them with interp1 but I doesn't seem to work.

How can I connect all of the Q and S points at once?

Thank you in advance!

This is my code so far: (of course not the best, but it works)

- I am using this Pan–Tompkins algorithm implemention

(I have added one ECG file as an example for what I use)

SE = strel('line',2*Fs,0); % 2 seconds ("fast timescale")

SE_1 = strel('line',3*Fs,0); % 3 seconds ("slow timescale")

signal_open = imopen(signal, SE); % Morphologically opens "image"

signal_close = imclose(signal_open, SE_1); % Morphologically closes "opened image"

signal_filtered = signal - signal_close; % Subtracts computed reference value

%% Step 2: R peak detection

%

% Pan-Tompkins Algorithm to detect R peaks and to calculate the RR Intervals

%

[R_Peak_amp,R_Peak_ind] = pan_tompkin(signal_filtered,Fs,0);

if(debug_flag)

%TEST

N = numel(R_Peak_ind);

Max_QRS_duration = 0.12;

Min_QRS_duration = 0.08;

QRS_Duration = (Min_QRS_duration-Max_QRS_duration).*rand(1,N) + Max_QRS_duration;

return;

end

Q_j_ind = zeros(0,length(R_Peak_ind));

S_j_ind = zeros(0,length(R_Peak_ind));

%% Step 3: search for the junctions of Q wave and S wave

first_peak = R_Peak_ind(1);

last_peak = R_Peak_ind(end-1);

q = 1;

s = 1;

while s < param_limit_length_to_N_samples

if(s == first_peak)

R_Peak_ind = R_Peak_ind(2:end);

end

if(s == R_Peak_ind(2))

d = 0;

Left_Peak = R_Peak_ind(1);

Current_Peak = R_Peak_ind(2);

Right_Peak = R_Peak_ind(3);

Y_CurrentPeak = signal_filtered(Current_Peak);

Y_LeftPeak = signal_filtered(Left_Peak);

Y_RightPeak = signal_filtered(Right_Peak);

delta_Y = signal_filtered(Current_Peak) - signal_filtered(Left_Peak);

delta_X = Current_Peak - Left_Peak;

for Pi = Current_Peak:-1:Left_Peak

d = d + 1;

% delta_Y = Y_R - Y_L;

D_Pi = abs(delta_Y * Pi - delta_X * signal_filtered(Pi) + Current_Peak * Y_LeftPeak - Y_CurrentPeak * Left_Peak)/sqrt(delta_Y^2 + delta_X^2);

W_Pi = cos(2*pi*(Pi-Current_Peak)/(Current_Peak-Left_Peak));

WD_Pi_q = D_Pi * W_Pi;

WD_Pi_amp_q(d) = WD_Pi_q;

end

half = floor((Current_Peak-Left_Peak)/2);

[~,Q_j] = max(WD_Pi_amp_q(1:half));

Q_j_ind(q) = Current_Peak - Q_j + 1;

delta_Y = signal_filtered(Right_Peak) - signal_filtered(Current_Peak);

delta_X = Right_Peak - Current_Peak;

d = 0;

for Pi = Current_Peak:Right_Peak

d = d + 1;

% delta_Y = Y_R - Y_L;

D_Pi = abs(delta_Y * Pi - delta_X * signal_filtered(Pi) + Right_Peak * Y_CurrentPeak - Y_RightPeak * Current_Peak)/sqrt(delta_Y^2 + delta_X^2);

W_Pi = cos(2*pi*(Pi-Current_Peak)/(Right_Peak-Current_Peak));

WD_Pi_s = D_Pi * W_Pi;

WD_Pi_amp_s(d) = WD_Pi_s;

end

half = floor((Right_Peak-Current_Peak)/2);

[~,S_j] = max(WD_Pi_amp_s(1:half));

S_j_ind(q) = Current_Peak + S_j - 1;

q = q + 1;

R_Peak_ind = R_Peak_ind(2:end);

end

if(s == last_peak)

break;

end

s = s+1;

end

subplot(2,1,1)

plot(t,signal_filtered,'-o','MarkerEdgeColor','b','MarkerIndices',Q_j_ind);

subplot(2,1,2)

plot(t,signal_filtered,'-o','MarkerEdgeColor','r','MarkerIndices',S_j_ind);

Shubh Sahu
on 4 May 2020

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!
## 0 Comments

Sign in to comment.