clear all; close all; clc;
%%
addpath('Data');
load 'TEST_S01_T01';
%load 'BPM_S04_T01';
% load 'DATA_10_TYPE02.mat'
% load 'DATA_10_TYPE02_BPMtrace'
%bpm=BPM0';
ppg1=sig(1,:); ppg2=sig(2,:);
acx1=sig(3,:); acy1=sig(4,:); acz1=sig(5,:);
accelaration1=delayseq(acx1',85)';
accelaration2=delayseq(acy1',85)';
accelaration3=delayseq(acz1',85)';
accelaration11=delayseq(acx1',0)'; %54
accelaration12=delayseq(acy1',0)'; %12
accelaration13=delayseq(acz1',0)'; %13
ac1=sqrt((accelaration11).^2+(accelaration12).^2+(accelaration13).^2);
steps=floor(((length(ppg1))-1000)/250+1);
%%
W1 = 1/(62.5); W2 = (5/62.5);
[b,a]=butter(4,[W1,W2]);
%ppg1 = filtfilt(b,a,ppg1);
%%
BPMest=zeros(1,steps);
figure('units','normalized','outerposition',[0 0 1 1])
k=1;
flag=0;BPM_flag=0;step=0;
for ss = 0:2:(steps-1)*2% ss = start time in seconds
clf
N = 4096*16;
start = ss*125+1;
stop = start + 8*125-1;
delay = 85;
%
ppg = ppg1(start:stop);
ppg = filtfilt(b,a,ppg);
ppg = detrend(ppg);
ppg=detrendnonlin(ppg,5);
ppg=detrend(ppg);
accelarationx=accelaration1(start:stop);accelarationy=accelaration2(start:stop);accelarationz=accelaration3(start:stop);
ac=ac1(start:stop);ac=ac-mean(ac);
acx=accelaration11(start:stop);acx=acx-mean(acx);
acy=accelaration12(start:stop);acy=acy-mean(acy);
acz=accelaration13(start:stop);acz=acz-mean(acz);
ac=acx+acy+acz;
if k~=1
[c,e]=butter(4,[(Nprev-15)/(60*62.5),(Nprev+15)/(60*62.5)],'stop') ;
acx=filtfilt(c,e,acx);acy=filtfilt(c,e,acy);acz=filtfilt(c,e,acz);ac=sqrt(acx.^2+acy.^2+acz.^2);ac=filtfilt(c,e,ac);
display('butter');
end
ac=ac-mean(ac);[Pac, fac]=periodogram(ac,hamming(length(ac)),8192,125);
acx=acx-mean(acx);[Pacx, facx]=periodogram(acx,hamming(length(acx)),8192,125);
acy=acy-mean(acy);[Pacy, facy]=periodogram(acy,hamming(length(acx)),8192,125);
acz=acz-mean(acz);[Pacz, facz]=periodogram(acz,hamming(length(acx)),8192,125);
bpm_index = ss/2 + 1;
%result = bpm(bpm_index);
%% delaying sequence
[ppg_ac,lag_ac]=xcorr(ppg,ac);[ppg_acx,lag_acx]=xcorr(ppg,acx);[ppg_acy,lag_acy]=xcorr(ppg,acy);[ppg_acz,lag_acz]=xcorr(ppg,acz);
[max1,index1]=max(abs(ppg_acx(1000:1125)));[max2,index2]=max(abs(ppg_acy(1000:1125)));[max3,index3]=max(abs(ppg_acz(1000:1125)));
acx=delayseq(acx',lag_acx(1000+index1))';acy=delayseq(acy',lag_acy(1000+index2))';acz=delayseq(acz',lag_acz(1000+index3))';
%% assigning ppg to ppg2
%ppg2=ppg;
% %% adaptive filtering
%
% filterorder=40;
% initialcoefficients=rand(filterorder+1,1)*2*3-3;
% S=struct('step',1e-2,'filterOrderNo',filterorder,'initialCoefficients',initialcoefficients,'gamma',1e-12 );
%
% u = acx;
% d = ppg;
% delta = 1;
% ppg2=ppg;
% if k~=1
% if (rms(acx))>0.28
% [outputVector,errorVector,thetaVector]=NLMS(d,u,S);
% ppg2 = errorVector';
% display('runx')
% else
% ppg2=ppg;
% display('not run')
% end
% u=acy;
% d=ppg2;
% if (rms(acy))>=0.28
% [outputVector,errorVector,thetaVector]=NLMS(d,u,S);
% ppg2=errorVector';
% display('runy');
% else
% ppg2=ppg2;
% display('not runy');
% end
% u=acz;
% d=ppg2;
% if (rms(acz))>=0.28
% [outputVector,errorVector,thetaVector]=NLMS(d,u,S);
% ppg2=errorVector';
% display('runz');
% else
% ppg2=ppg2;
% display('not runz');
% end
% end
%% Adaptive Filtering RLS
ac_bound=0.9;
% RLS
l = 600; % order/length of filter
lambda = 1;
coeffs = zeros(1,l);
states = zeros(1, l-1);
x = ppg;
f1 = 0; f2 = 0; f3 = 0;
if k~=1
if rms(acx)>ac_bound
invcov = (1/(cov(ppg)))*eye(l);
ha = adaptfilt.rls(l,lambda,invcov,coeffs,states);
[~,x] = filter(ha, acx, x);f1=1;
% x = filter(ha, acx, x);f1=1;
end
if rms(acy)>ac_bound
invcov = (1/(cov(ppg)))*eye(l);
ha = adaptfilt.rls(l,lambda,invcov,coeffs,states);
[~,x] = filter(ha, acy, x);f2=1;
% x = filter(ha, acy, x);f2=1;
end
if rms(acz)>ac_bound
invcov = (1/(cov(ppg)))*eye(l);
ha = adaptfilt.rls(l,lambda,invcov,coeffs,states);
[~,x] = filter(ha, acz, x);f3=1;
% x = filter(ha, acz, x);f3=1;
end
end
ppg2 = x;
ppg2 = filtfilt(b,a,ppg2);
%% eemd
PPG2=eemd(ppg2,1,10)';
ppg_eemd1=PPG2(1,:);ppg_eemd2=PPG2(6,:);ppg_eemd3=PPG2(1,:);
if k~=1
[numerator,denominator]=butter(4,[(Nprev-30)/(60*62.5),(Nprev+30)/(60*62.5)]);
ppg_eemd1=filtfilt(numerator,denominator,ppg_eemd1);ppg_eemd1=filtfilt(numerator,denominator,ppg_eemd1);ppg_eemd1=filtfilt(numerator,denominator,ppg_eemd1);
end
%% plotting EEMD
for i=1:10
[Pi,Fi]=periodogram(PPG2(i,:),rectwin(1000),4096,125);
figure(6);subplot(10,1,i);plot(Fi*60,Pi/max(Pi));axis([0 300 0 1]);
end
%% For EEMD1
[Pxx_eemd1,f_eemd1]=periodogram(ppg_eemd1,rectwin(length(ppg_eemd1)),8192,125);nPxx_eemd1=Pxx_eemd1/max(Pxx_eemd1);[peaks,locs]=findpeaks(nPxx_eemd1);
[Pppg,fppg]=periodogram(ppg,rectwin(length(ppg)),2048,125);Pppg=Pppg/max(Pppg);[peaks1,locs1]=findpeaks(Pppg);
figure(1);subplot(211);plot(f_eemd1*60,nPxx_eemd1,f_eemd1(locs)*60,nPxx_eemd1(locs),'k^','markerfacecolor',[1 0 0]);axis([0 250 0 1]);subplot(212);plot(fppg*60,Pppg,fppg(locs1)*60,Pppg(locs1),'k^','markerfacecolor',[1 0 0]);axis([0 250 0 1]);legend(strcat('Start Time: ', num2str(ss), ' seconds'));
% [ppg_ac,lag_ac]=xcorr(ppg,ac);[ppg_acx,lag_acx]=xcorr(ppg,acx);[ppg_acy,lag_acy]=xcorr(ppg,acy);[ppg_acz,lag_acz]=xcorr(ppg,acz);
figure(2);subplot(411);plot(lag_ac,abs(ppg_ac));subplot(412);plot(lag_acx,abs(ppg_acx));subplot(413);plot(lag_acy,abs(ppg_acy));subplot(414);plot(lag_acz,abs(ppg_acz));
t = linspace(ss, ss+8, length(ppg));
figure(3);subplot(211);plot(t,ppg);subplot(212);plot(t,ppg2);
figure(4);subplot(411);plot(t,ac);subplot(412);plot(t,acx);subplot(413);plot(t,acy);subplot(414);plot(t,acz);
figure(5);subplot(411);plot(fac*60,Pac/max(Pac));axis([0 250 0 1]);subplot(412);plot(facx*60,Pacx/max(Pacx));axis([0 250 0 1]);subplot(413);plot(facy*60,Pacy/max(Pacy));axis([0 250 0 1]);subplot(414);plot(facz*60,Pacz/max(Pacz));axis([0 250 0 1]);
if k~=1
nPxx_eemd1=nPxx_eemd1-Pacx/max(Pacx)-Pacy/max(Pacy)-Pacz/max(Pacz);
end
%% initializing the first BPM EEMD1
if k==1;
[maximum_max,index_max]=findpeaks(nPxx_eemd1);
p=0;
for i=1:length(index_max)
if f_eemd1(index_max(i))*60>=55 && f_eemd1(index_max(i))*60<=130
p=p+1;
new_index_max(p)=index_max(i);
new_maximum_max(p)=maximum_max(i);
end
end
[A,B]=max(new_maximum_max);
Nprev1=f_eemd1(new_index_max(B))*60;
end
%% For EEMD2
[Pxx_eemd2,f_eemd2]=periodogram(ppg_eemd2,hanning(length(ppg_eemd2)),8192,125);nPxx_eemd2=Pxx_eemd2/max(Pxx_eemd2);[peaks,locs]=findpeaks(nPxx_eemd2);
if k~=1
nPxx_eemd2=nPxx_eemd2-Pacx/max(Pacx)-Pacy/max(Pacy)-Pacz/max(Pacz);
end
%% initializing the first BPM EEMD2
if k==1;
clearvars maximum_max; clearvars index_max;clearvars new_index_max;clearvars new_maximum_max;clearvars peaks; clearvars locs;
[maximum_max,index_max]=findpeaks(nPxx_eemd2);
p=0;
for i=1:length(index_max)
if f_eemd2(index_max(i))
