function [new_noNorm, old, gt] = SPARE(ground_truth, raw_PPG, time_PPG, winLen, stride, fs_PPG, True_BPMs, plotting, log) %len = round(winLen*0.3086*fs_PPG); len = stride * fs_PPG; % Initialize old = zeros(1, len); gt = zeros(1, len); new_noNorm = zeros(1, len); times = time_PPG(1:len); count = 0; % Filter the data [dataFull, ~] = PPG_filter(raw_PPG, fs_PPG); signal_end = time_PPG(end); segment_start = time_PPG(1); segment_end = time_PPG(1) + winLen; while(segment_end < signal_end) count = count+1; %disp(count) if mod(count,1) == 0 disp((1-(signal_end-segment_end)/(signal_end-time_PPG(1)))*100) end data = getSignalSegmentTEST(dataFull, time_PPG, segment_start, segment_end); data_t = getSignalSegmentTEST(time_PPG, time_PPG, segment_start, segment_end); data_ref = getSignalSegmentTEST(ground_truth, time_PPG, segment_start, segment_end); data_noise = getSignalSegmentTEST(raw_PPG, time_PPG, segment_start, segment_end); %% if (~all(isnan(data))) RR_mean = 1/True_BPMs(count)*60; if (log) disp('BPM from ECG:') disp((1/RR_mean*60)) end PPG = data'; % remove NaNs from signal PPG(isnan(PPG)) = 0; % devide the PPG in components comps = SSA(PPG, 10, floor(length(PPG)*0.75/2) -1); N = 60*fs_PPG; b = zeros(N,1); % define frequency axis %f = (0:fs_PPG/(N-1):fs_PPG)*60; % do FFT per components for i=1:length(comps(:,1)) [BPM, s] = SSR(comps(i,:)', fs_PPG, N); b = b + s; end b = b(40:500); BPM = BPM(40:500); p = abs(b).^2; BPMs = BPM'; % find the index of all peaks [~,index_P] = findpeaks(p); index_Full = index_P; % BPM of all peaks BPM_peaks_full = BPMs(index_Full); % compute second derivative second = diff(diff(comps'))'; % FFT per component of the second derivative a = zeros(N,1); for i=1:length(second(:,1)) [~, Comps] = SSR(second(i,:)', fs_PPG, N); a = a + Comps; end a = a(40:500); a = abs(a).^2; % find all peaks of the spectrum of the second derivative [pks,index_P] = findpeaks(a); % BPM of all the peaks BPM_peaks = BPMs(index_P); % compute distance between the biggests peaks and the HR from ECG distance = abs(BPM_peaks-((1/RR_mean)*60)); % remove the index that are less than 10% of the biggest one among % the candidates index_cand = find(distance<7); pks_cands = pks(index_cand); maximum = max(pks(index_cand)); index_P_first = index_cand(pks_cands>=0.2*maximum); % candidates for being the representative of the fundamental of the signal BPM_cand = BPM_peaks(index_P_first); % same but for the second harmonic distance = abs(BPM_peaks-((1/RR_mean)*60*2)); BPM_cand2 = BPM_peaks(distance<14); % same but for the third harmonic distance = abs(BPM_peaks-((1/RR_mean)*60*3)); BPM_cand3 = BPM_peaks(distance<20); % from the three sets of peaks (BPM_cand, BPM_cand2, BPM_cand3), % find the triplet that staisfy the condition: % (HR_cand in BPM_cand, HR_cand2 in BPM_cand2, HR_cand3 in BPM_cand3) s.t. HR_cand2/HR_cand = 2, HR_cand3/HR_cand = 3 valMin = +inf; valMin2 = +inf; for i=1:numel(BPM_cand) for j=1:numel(BPM_cand2) for k=1:numel(BPM_cand3) val = abs(BPM_cand2(j)/BPM_cand(i)-2); val2 = abs(BPM_cand3(k)/BPM_cand(i)-3); if (val <= valMin && val2 <= valMin2) valMin = val; valMin2 = val2; HR_cand = BPM_cand(i); HR_cand2 = BPM_cand2(j); HR_cand3 = BPM_cand3(k); end end end end % sanity check: HR_cand/HR_cand is at most 30% away from 2 and % HR_cand do not differ much from HR from ECG if (valMin<0.3 && abs(HR_cand-(1/RR_mean)*60) < 15) HR_p = HR_cand; HR_p2 = HR_cand2; HR_p3 = HR_cand3; else % otherwise, use HR from ECG (extrema ratio) HR_p = (1/RR_mean)*60; HR_p2 = (1/RR_mean)*60*2; HR_p3 = (1/RR_mean)*60*3; end % find index of the peaks computed from the spectrum of the orignal % signal [~, idx] = min(abs(BPM_peaks_full-(HR_p))); [~, idx] = min(abs(BPMs-(BPM_peaks_full(idx)))); if (log) disp('HR peak from PPG:') disp(BPMs(idx)) end HR_2 = BPMs(idx); BPM_fund = BPMs(idx); p_fund = p(idx); max_idx = idx; while (max_idx0 && p(min_idx-1)<=p(min_idx)) min_idx = min_idx-1; end f_min = BPMs(min_idx)/60; f_max = BPMs(max_idx)/60; distance2HR = abs( ((f_max+f_min)/2 - (1/RR_mean))*60 ); if ((f_max-f_min)<0.02 || (f_max-f_min>1) || distance2HR>10) f_min = ((1/RR_mean)*60 -15)/60; f_max = ((1/RR_mean)*60 +15)/60; end x = bandpass(PPG, [f_min f_max], fs_PPG); HR = ((f_max+f_min)/2)*60; [~, idx] = min(abs(BPM_peaks_full-(HR_p2))); [~, idx] = min(abs(BPMs-(BPM_peaks_full(idx)))); if (log) disp('2th HR peak from PPG:') disp(BPMs(idx)) end BPM_second = BPMs(idx); p_second = p(idx); max_idx = idx; while (max_idx0 && p(min_idx-1)<=p(min_idx)) min_idx = min_idx-1; end f_min = BPMs(min_idx)/60; f_max = BPMs(max_idx)/60; distance2HR = abs( ((f_max+f_min)/2 - (1/RR_mean*2))*60 ); if ((f_max-f_min)<0.02 || (f_max-f_min>1) || distance2HR>15) f_min = ((1/RR_mean)*60*2 -15)/60; f_max = ((1/RR_mean)*60*2 +15)/60; end y = bandpass(PPG, [f_min f_max], fs_PPG); [~, idx] = min(abs(BPM_peaks_full-(HR_p3))); [~, idx] = min(abs(BPMs-(BPM_peaks_full(idx)))); BPM_third = BPMs(idx); p_third = p(idx); max_idx = idx; while (max_idx0 && p(min_idx-1)<=p(min_idx)) min_idx = min_idx-1; end f_min = BPMs(min_idx)/60; f_max = BPMs(max_idx)/60; distance2HR = abs( ((f_max+f_min)/2 - (1/RR_mean*3))*60 ); if ((f_max-f_min)<0.02 || (f_max-f_min>1) || distance2HR>25) f_max = 1/RR_mean*3+30/60; f_min = 1/RR_mean*3-30/60; end z = bandpass(PPG, [f_min f_max], fs_PPG); segment_start = segment_start+len/fs_PPG; segment_end = segment_start + winLen; Hal = x+y+z; %disp(kurtosis(raw_PPG)) % if (kurtosis(raw_PPG)<2.5) % Hal = data_noise.'; % end if (len*2 <= numel(Hal)) if (max(abs(data_noise(len:len*2))) > 2.5) new_noNorm = [new_noNorm, Hal(len:len*2-1)]; else new_noNorm = [new_noNorm, data_noise(len:len*2-1)]; end old = [old, data_noise(len:len*2-1).']; gt = [gt, data_ref(len:len*2-1)]; times = [times, data_t(len:len*2-1)]; else new_noNorm = [new_noNorm, Hal(len:end)]; old = [old, data_noise(len:end).']; gt = [gt, data_ref(len:end)]; times = [times, data_t(len:end)]; end if (plotting) try figure plot(subplot(2,1,1), Hal(len:len*2-1)); hold on; plot(subplot(2,1,1), PPG(len:len*2-1)); legend('filtered', 'original') %[~, p_ref] = mySpectrum(data_ref, fs_PPG); plot(subplot(2,1,2), BPMs, p); hold on; plot(subplot(2,1,2), BPMs, p); plot(subplot(2,1,2), BPM_fund,p_fund,'r*'); plot(subplot(2,1,2), BPM_second,p_second,'r*'); plot(subplot(2,1,2), BPM_third,p_third,'r*'); xlim([40, 350]) catch e disp(e) end end oldRR = RR_mean; else aiaia = 0; segment_start = segment_start+len/fs_PPG; segment_end = segment_start + winLen; end end end function data_segment = getSignalSegmentTEST(signal, signal_t, start, ending) index = find((start