Homec4science
Diffusion SPARE code (master) Browse

R10935/SPARE.m
8a5fd2c8d253master

Last Change
Raw File

SPARE.m
View Options

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_idx<numel(p) && p(max_idx+1)<=p(max_idx))
max_idx = max_idx+1;
end
min_idx = idx;
while (min_idx-1>0 && 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_idx<numel(p) && p(max_idx+1)<=p(max_idx) )
max_idx = max_idx+1;
end
min_idx = idx;
while (min_idx-1>0 && 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_idx<numel(p) && p(max_idx+1)<=p(max_idx))
max_idx = max_idx+1;
end
min_idx = idx;
while (min_idx-1>0 && 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<signal_t) & (signal_t<ending));
data_segment = signal(index);
end
c4science · Help