Datasets
Models
Downloads: 1
Diff of
/preprocessOfApneaECG/BioSigKit/Algorithms/SimpleRST.m
[000000]
..
[409112]
Switch to side-by-side view
--- a +++ b/preprocessOfApneaECG/BioSigKit/Algorithms/SimpleRST.m @@ -0,0 +1,310 @@ +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 + +%% ========================= 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 \ No newline at end of file