Contents

function [R_i,R_amp,S_i,S_amp,T_i,T_amp,Q_i,Q_amp,buffer_plot] = SimpleRST(ecg,fs,gr)

function [R_i,R_amp,S_i,S_amp,T_i,T_amp]=peakdetect(ecg,fs,view)

=================== Online Adaptive QRS detector ==================== %%

========================== Description ============================= %%

QRS detection Detects Q , R and S waves,T Waves Uses the state-machine logic to determine different peaks in an ECG signal. It has the ability to confront noise by canceling out the noise by high pass filtering and baseline wander by low pass. Besides, check out criterion to stop detection of spikes. The code is written in a way for future online implementation.

Inputs

ecg : raw ecg vector fs : sampling frequency view : display results? (0: no, 1: Yes)

Outputs

indexes and amplitudes of R_i, R_amp, etc heart_rate computed heart rate buffer_plot : processed signal

============== 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

Updates :

Feb, 2018 : Clean up and fixes.

============== Now Part of BioSigKit ======================== %%

if nargin < 3
    gr = 1;                                                                % on show Sig
    if nargin <2
       fs = 250;                                                           % default Sampling frequency
       if nargin < 1
           error('You need to provide a signal!');
       end
    end
end
Error using SimpleRST (line 45)
You need to provide a signal!

========================= initialize ============================ %%

R_i = zeros(1,length(ecg));                                                % save index of R wave
R_amp = zeros(1,length(ecg));                                              % save amp of R wave
S_i = zeros(1,length(ecg));                                                % save index of S wave
S_amp = zeros(1,length(ecg));                                              % save amp of S wave
T_i = zeros(1,length(ecg));                                                % save index of T wave
T_amp = zeros(1,length(ecg));                                              % save amp of T wave
Q_i = zeros(1,length(ecg));                                                % vectors to store Q wave
Q_amp = zeros(1,length(ecg));                                              % Vectors to store Q wave
S_amp1 = zeros(1,length(ecg));                                             % Buffer to set the adaptive T wave onset
thres_p =zeros(1,length(ecg));                                             % For plotting adaptive threshold
S_amp1_i = zeros(1,length(ecg));                                           % To save indices of S thres
buffer_plot = zeros(1,length(ecg));
thres2_p = zeros(1,length(ecg));                                           % T wave threshold indices
window = round(0.04*fs);                                                   % averaging window size
buffer_long= zeros(1,window);                                              % buffer for online processing
state = 0 ;                                                                % determines the state of the machine in the algorithm
c = 0;                                                                     % counter to determine that the state-machine doesnt get stock in T wave detection wave
T_on = 0;                                                                  % counter showing for how many samples the signal stayed above T wave threshold
T_on1=0;                                                                   % counter to make sure its the real onset of T wave
S_on = 0;                                                                  % counter to make sure its the real onset of S wave
sleep = 0;                                                                 % counter that avoids the detection of several R waves in a short time
buffer_base=zeros(1,2*fs);                                                 % buffer to determine online adaptive mean of the signal
dum = 0;                                                                   % counter for detecting the exact R wave
weight = 1.8;                                                              % initial value of the weigth
co = 0;                                                                    % T wave counter to come out of state after a certain time
thres_p_i = zeros(1,length(ecg));                                          % To save indices of main thres
thres2_p_i = zeros(1,length(ecg));                                         %to save indices of T threshold

========================= preprocess ================================ %%

ecg = ecg (:);                                                             % make sure its a vector
ecg_raw =ecg;                                                              % take the raw signal for plotting later

==================== Noise cancelation(Filtering) =================== %%

f1=0.5;                                                                    % cuttoff low frequency to get rid of baseline wander
f2=45;                                                                     % 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 = filtfilt(a,b,ecg);

============== define two buffers ================= %%

buffer_mean=mean(abs(ecg(1:2*fs)-mean(ecg(1:2*fs))));                      % adaptive threshold DC corrected (baseline removed)
buffer_T = mean(ecg(1:2*fs));                                              % second adaptive threshold to be used for T wave detection

================== Counters ============================ %%

B_Lcounter = 0;
B_counter = 0;
SP_counter = 0;
thres_p_C = 0;
R_C = 0;
S_C = 0;
T_C = 0;
Q_C = 0;
thres2_p_C = 0;

=start online inference (Assuming the signal is being acquired online) %%

for i = 1 : length(ecg)
 B_Lcounter = B_Lcounter + 1;
 buffer_long(B_Lcounter) = ecg(i);                                         % save the upcoming new samples
 if B_Lcounter > window
    B_Lcounter = 0;
 end

 B_counter = B_counter + 1;
 buffer_base(B_counter) = ecg(i);                                          % save the baseline samples

============================= Renew Mean ======================= %%

 if B_counter >= 2*fs
    buffer_mean = mean(abs(buffer_base - mean(buffer_base)));
    buffer_T = mean(buffer_base);
    B_counter = 0;
 end

========= Smooth 15 samples and add the new upcoming samples ======== %%

  if i >= window                                                  % take a window with length 15 samples for averaging
      mean_online = mean(buffer_long);                       % take the mean
      SP_counter = SP_counter + 1;
      buffer_plot(SP_counter) = mean_online;                               % save the processed signal

============== Enter state 1(putative R wave) ================ %%

    if state == 0
     if SP_counter >= 3                                                                         % added to handle bugg for now
      if (mean_online > buffer_mean*weight) && (buffer_plot(i-1-window) > buffer_plot(i-window))    % 2.4*buffer_mean
          state = 1;                                                                            % entered R peak detection mode
          currentmax = buffer_plot(i-1-window);
          ind = i-1-window;
          thres_p_C = thres_p_C + 1;
          thres_p(thres_p_C) = buffer_mean*weight;
          thres_p_i(thres_p_C) = ind;
      else
          state = 0;
      end
     end
    end

============= Locate R by finding highest Peak =================== %%

      if state == 1                                                        % look for the highest peak
            if  currentmax > buffer_plot(i-window)
                dum = dum + 1;
                if dum > 4
                   R_C = R_C + 1;
                   R_i(R_C) = ind;                                          % save index
                   R_amp(R_C) = buffer_plot(ind);                          % save index
                %-------------- Locate Q wave --------------------%
                  [Q_tamp,Q_ti] = min(buffer_plot(ind-round(0.040*fs):(ind)));
                  Q_ti = ind-round(0.040*fs) + Q_ti -1;
                  Q_C = Q_C + 1;
                  Q_i(Q_C) = Q_ti;
                  Q_amp(Q_C) = Q_tamp;
                  if R_C > 8
                     weight = 0.30*mean(R_amp(R_C-7:R_C));                  % calculate the 35% of the last 8 R waves
                     weight = weight/buffer_mean;
                  end
                  state = 2;                                                % enter S detection mode state 2
                  dum = 0;
                end
            else
                dum = 0;
                state = 0;
            end

      end

=== check if Sig drops below the threshold to look for S wave === %%

      if state == 2
        if  mean_online <= buffer_mean                                     % check the threshold
             state = 3;                                                    % enter S detection
        end
      end

============ Enter S wave detection state3 (S detection) =========== %%

         if state == 3
            co = co + 1;
          if co < round(0.200*fs)
            if buffer_plot(i-window-1) <= buffer_plot(i-window)            % see when the slope changes
             S_on = S_on + 1;                                              % set a counter to see if its a real change or just noise
             if S_on >= round(0.0120*fs)
                S_C = S_C + 1;
                S_i(S_C) = i-window-4;                                     % save index of S wave
                S_amp(S_C) = buffer_plot(i-window-4);                      % save index
                S_amp1(S_C) = buffer_plot(i-window-4);                     % ecg(i-4)
                S_amp1_i(S_C) = ind;                                       % index of S_amp1_i
                state = 4;                                                 % enter T detection mode
                S_on = 0;
                co = 0;
             end
            end
          else
                state = 4;
                co = 0;
          end
          end

======= enter state 4 possible T wave detection ============ %%

       if state == 4
         if mean_online < buffer_mean                                      % See if the signal drops below mean
           state = 6;                                                      % Confirm
         end
       end

======= Enter state 6 which is T wave possible detection ======%%

       if state ==6
         c = c + 1;                                                        % set a counter to exit the state if no T wave detected after 0.3 second
         if c <= 0.7*fs
             %------------------------------------------------------------%
             % set a double threshold based on the last detected S wave and
             % baseline of the signal and look for T wave in between these
             % two threshold
             %------------------------------------------------------------%
             thres2 = ((abs(abs(buffer_T)-abs(S_amp1(S_C))))*3/4 + S_amp1(S_C));
             thres2_p_C = thres2_p_C + 1;
             thres2_p(thres2_p_C) = thres2;
             thres2_p_i(thres2_p_C) = ind;
             if mean_online > thres2
              T_on = T_on +1;                                              % make sure it stays on for at least 3 samples
              if T_on >= round(0.0120*fs)
               if buffer_plot(i-window-1)>= buffer_plot(i-window)
                   T_on1 = T_on1+1;                                        % make sure its a real slope change
                  if T_on1 > round(0.0320*fs)
                   T_C = T_C + 1;
                   T_i(T_C) = i-window-11;                                 % save index of T wave
                   T_amp(T_C) = buffer_plot(i-window-11);                  % save index
                   state = 5;                                              % enter sleep mode
                   T_on = 0;
                   T_on1 = 0;
                  end
               end
              end
             end
         else
             state= 5;                                                     % enter Sleep mode
         end

       end

==== Sleep To avoid multiple detections ================== %%

       if state==5
           sleep =sleep+c+1;
           c = 0;
           if sleep/fs >= 0.400
               state = 0;
               sleep = 0;
           end
       end
  end
end

============== Adjust Length of Signals ===================== %%

R_i = R_i(1:R_C);
S_i = S_i(1:S_C);
S_amp1 = S_amp1(1:S_C);
S_amp1_i = S_amp1_i(1:S_C);
T_i = T_i(1:T_C);
Q_i = Q_i(1:Q_C);
thres_p_i = thres_p_i(1:thres_p_C);
thres_p = thres_p(1:thres_p_C);
buffer_plot = buffer_plot(1:SP_counter);
thres2_p = thres2_p(1:thres2_p_C);
thres2_p_i = thres2_p_i(1:thres2_p_C);

conditions

%heart_rate=R_C/(time_scale/60); % calculate heart rate
%msgbox(strcat('Heart-rate is = ',mat2str(heart_rate)));

plottings

if gr
   view = length(ecg)/fs;
   time = 1/fs:1/fs:view;
   R = find(R_i <= view*fs);                                                  % determine the length for plotting vectors
   S = find(S_i <= view*fs);                                                  % determine the length for plotting vectors
   T = find(T_i <= view*fs);                                                  % determine the length for plotting vectors
   Q = find(Q_i <= view*fs);                                                  % determine the length for plotting vectors
   L1 = find(thres_p_i <= view*fs);
   L2 = find(S_amp1_i <= view*fs);
   L3 = find(thres2_p_i <= view*fs);
   if view*fs > length(buffer_plot)
      ax(1) = subplot(211);plot(time(1:length(buffer_plot)),buffer_plot(1:end));
   else
      ax(1) = subplot(211);plot(time,buffer_plot(1:(view*fs)));
   end
   axis tight;
   hold on,scatter(R_i(1:R(end))./fs,R_amp(1:R(end)),'r');
   hold on,scatter(S_i(1:S(end))./fs,S_amp(1:S(end)),'g');
   hold on,scatter(T_i(1:T(end))./fs,T_amp(1:T(end)),'k');
   hold on,scatter(Q_i(1:Q(end))./fs,Q_amp(1:Q(end)),'m');
   hold on,plot(thres_p_i(1:L1(end))./fs,thres_p(1:L1(end)),'LineStyle','-.','color','r',...
    'LineWidth',2.5);
   hold on,plot(S_amp1_i(1:L2(end))./fs,S_amp1(1:L2(end)),'LineStyle','--','color','c',...
    'LineWidth',2.5);
   hold on,plot(thres2_p_i(1:L3(end))./fs,thres2_p(1:L3(end)),'-k','LineWidth',2);
   legend('Raw ECG Signal','R wave','S wave','T wave','R adaptive thres','Latest S wave','T wave adaptive threshold threshold','Location','NorthOutside','Orientation','horizontal');
   xlabel('Time(sec)'),ylabel('V');
   axis tight;
   title('Zoom in to see both signal details overlaied');
   title('Filtered, smoothed and processed signal');
   ax(2) =subplot(212);
   plot(time,ecg_raw(1:(round(view*fs))));
   title('Raw ECG')
   xlabel('Time(sec)'),ylabel('V');
   legend();
   linkaxes(ax,'x');
   zoom on;
   axis tight;
end