Портативная система сбора данных о здоровье человека

Скрипт программы для работы с ЭКГ сигналом

clc

clear all

close all

pkg load signal

maxSignals=csvread('7.csv');

%3,4 ppg

%5 not filtering ecg

ecgFullNF=real(maxSignals(2:end,5));

n1=8000;n2=15104;

Fs=200;%200 Гц

ecgFullNF_cutted=ecgFullNF(n1:n2)';

figure;

plot(ecgFullNF)

figure;

plot(ecgFullNF_cutted)

[b1,a1]=butter(5,1/100,'high');

[H,w]=freqz(b1,a1,1000);

figure;

plot(w/2/pi*1000,10*log10(abs(H)))

[ecgFullMalabFiltered20,d2]=filter(b1,a1,ecgFullNF_cutted(1:end));

figure;

plot((0:length(ecgFullNF)-1)*1/Fs,ecgFullNF*(12.247e-6)/8/9.5*1000)

xlabel('f,Гц'),ylabel('Амлитудагармоник,мB'),grid

title('ЭКГсигнал')

if 1

p=[1 0 0 0 -1];%Знаменатель

z=zeros(1,201);

z(1)=-1;z(201)=1;

z(101)=50;z(105)=-50;

k=1/50;

[ecgFullMalabFiltered]=filter(z*k,p,ecgFullNF_cutted);

ecgFullMalabFiltered=ecgFullMalabFiltered(101:end);

if 1

nOrder=4;

[b1,a1]=butter(nOrder,1/100,'high');

[b2,a2]=butter(nOrder,40/100,'low');

[b3,a3]=butter(nOrder,[49 51]/100,'stop');

[ecgFullMalabFiltered1,d1]=filter(b2,a2,ecgFullMalabFiltered);

figure

plot(ecgFullMalabFiltered1)

[ecgFullMalabFiltered,d2]=filter(b1,a1,ecgFullMalabFiltered1(3:end));

figure

plot(ecgFullMalabFiltered)

%[ecgFullMalabFiltered3]=filter(z*k,p,ecgFullMalabFiltered1);

%[ecgFullMalabFiltered2,d2]=filter(b2,a2,ecgFullMalabFiltered1(1:end));

%[ecgFullMalabFiltered3,d3]=filter(b3,a3,ecgFullMalabFiltered2(1:end));

%ecgFullMalabFiltered=ecgFullMalabFiltered3(3:end);

%ecgFullMalabFilteredFetched=ecgFullMalabFiltered(nStart:end);

end

if 1

figure;

hold on

plot((0:length(ecgFullNF_cutted)-1)*1/Fs,ecgFullNF_cutted*(12.247e-6)/8/5*1000,'-r')

plot((0:length(ecgFullMalabFiltered)-1)*1/Fs,ecgFullMalabFiltered*(12.247e-6)/8/5*1000,'-b')

xlabel('t,c'),ylabel('АмлитудаЭКГ,мB'),grid,legend('Безфильтра','Сфильтром')

title('Нефильтрованный и фильтрованный ЭКГ сигналы')

hold off

ecgFullMalabFilteredFetched=ecgFullMalabFiltered;

end

nStart=1;%

ecgV1=ecgFullNF_cutted(nStart:end)*(12.247e-6)/8/9.5*1000;

ecgV2=ecgFullMalabFilteredFetched*(12.247e-6)/8/9.5*1000;

fetchSig=ecgFullMalabFilteredFetched(1:Fs);%Находимпервыйпик

HRM_time=find(max(fetchSig)==ecgFullMalabFilteredFetched);%Находимеговремя

startInd=HRM_time+11;

while (length(ecgFullMalabFilteredFetched)-startInd>Fs)

fetchSig=ecgFullMalabFilteredFetched((startInd):(startInd+Fs));

tekMax=find(max(fetchSig)==ecgFullMalabFilteredFetched);

if ((tekMax-HRM_time(end))>Fs/2)

HRM_time=[HRM_timetekMax];

startInd=tekMax+10;

else startInd=tekMax+10;end

endwhile

%HRM_time=[2:length(HRM_time)];

HRM=[];

for m=1:length(HRM_time)-1

HRM=[HRM 60/(HRM_time(m+1)-HRM_time(m))*Fs];

endfor

HRM_mean=mean(HRM);

figure;

plot(0:length(HRM)-1,[HRM],'-ob'),grid,title('HRM')

xlabel('interval,n'),ylabel('HR, beats/min')

legend(['HRM mean=' num2str(HRM_mean)])

figure;

hold on

%tEcg=(0:length(ecgV1)-1)*1/Fs;

plot((0:length(ecgV1)-1)*1/Fs,ecgV1)%not filtered

plot((0:length(ecgV2)-1)*1/Fs,ecgV2)%filtered

legend('recorded ECG','filtered ECG')

grid,xlabel('t,c'),ylabel('ЭКГсигнал,мB')

hold off

%FFT

fftEcg1=abs(fft(ecgV1));fftEcg0=fftEcg1(1);fftEcg1(1)=0;

fftEcg2=abs(fft(ecgV2));

L=length(fftEcg1)/2;

figure;

subplot(2,1,1)

plot((0:L)/(2*L)*Fs,fftEcg1(1:L+1)/(2*L))%not filtered

xlabel('f,Гц'),ylabel('Амлитудагармоник,мB'),grid

title(['Нефильтрованный ЭКГ сигнал без постоянной составляющей(' num2str(fftEcg0/(2*L)) ' B)'])

subplot(2,1,2)

plot((0:L)/(2*L)*Fs,fftEcg2(1:L+1)/(2*L))%filtered

xlabel('f,Гц'),ylabel('Амлитудагармоник,мB'),grid

title('Фильтрованный ЭКГ сигнал')

%Записываем в файл данные

save -mat7-binary 'ecg_data.mat' 'ecgV2';

end

Скрипт программы для работы с ФПГ сигналом

clc

clear all

close all

pkg load signal

ppgAll=csvread('7.csv');

%3,4 ppg

%5 not filtering ecg

Fs=200;%Частота дискретизации

ppgRed=ppgAll(2:end,4);%4 столбец красный

ppgIRed=ppgAll(2:end,3);%3 инфракрасный

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if 1 % комментарий fft

%Plot 2000 samples

figure;

hold on

plot(ppgRed,'-b')

ppgRed1000=ppgRed(8000:10000)';

ppgIRed1000=ppgIRed(8000:10000)';

fftRed=abs(fft(ppgRed1000))/length(ppgRed1000);fftRed0=fftRed(1);fftRed(1)=0;

fftIRed=abs(fft(ppgIRed1000))/length(ppgIRed1000);fftIRed0=fftIRed(1);fftIRed(1)=0;

L=length(fftRed)/2;

figure;

subplot(2,1,1)

plot((0:L/5)/2/L*Fs,fftRed(1:(L/5+1))),legend(['Ampl(0)=' num2str(fftRed0)]),grid

title('Red light AC')

xlabel('t'),ylabel('Amplitude')

subplot(2,1,2)

plot((0:L/5)/2/L*Fs,fftIRed(1:(L/5+1))),legend(['Ampl(0)=' num2str(fftIRed0)]),grid

title('IRed light AC')

xlabel('t'),ylabel('Amplitude')

%

order=6;

wCutoff=10;

Wn=wCutoff/(Fs/2);

[b a] = butter(order,Wn,'low');

%f_filt=0:(Fs/2)/1000:Fs/2;

%[H w]=freqz(b,a,1000);

%figure;

%plot(w,10*log10(abs(H)),'-b')

[ppgRed1000Filt,delay]=filter(b,a,ppgRed1000);

figure;

hold on

plot(ppgRed1000,'-b')

plot(ppgRed1000Filt(14:end),'-r')

grid

hold off

endif

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if 1

%Фильтр

order=6;%порядок

wCutoff=10;%Частота среза

Wn=wCutoff/(Fs/2);

[b a] = butter(order,Wn,'low');

%Algorithm

nStart=8000;

nEnd=15000;

ppgRedFiltered=filter(b,a,ppgRed(nStart:nEnd+order*2));

ppgIRedFiltered=filter(b,a,ppgIRed(nStart:nEnd+order*2));

%не сдвинутые во времени, фильтрованные отсчеты с 8000 до 15000

ppgRed30s=ppgRedFiltered(2*order:end)';

ppgIRed30s=ppgIRedFiltered(2*order:end)';

%Compute SpO2

%Firstly: search for minimum value - dc value, start algorithm

%v1Red - minimum 1

%v2Red - maximum 1

%v3Red - minimum 2

%v4Red - maximum 2

v1Red=0;v2Red=0;v3Red=0;v4Red=0;

v1IRed=0;v2IRed=0;v3IRed=0;v4IRed=0;

vRed=0;errorVRed=0;

vIRed=0;errorVIRed=0;

flag=0;

RedAc=[];RedDc=[];IRedAc=[];IRedDc=[];

nakR=[];nakIR=[];

tK=[];tK1=0;tK2=0;tK3=0;tK4=0;

counter=0;

stepCheck=10;

for k=1:length(ppgIRed30s)

%Red

vNewRed=(ppgRed30s(k));

vNewIRed=(ppgIRed30s(k));

switch (flag)

case 0

if (vNewIRed>vIRed)

if(counter<stepCheck)

if (vIRed<errorVIRed)

counter=0;

errorVRed=vRed;

errorVIRed=vIRed;

tK1=k-1;

endif

else

v1IRed=errorVIRed;

v1Red=errorVRed;

flag=1;

counter=0;

endif

counter=counter+1;

endif%search minimum

%

case 1

if (vNewIRed<vIRed)

if(counter<stepCheck)

if (vIRed>errorVIRed)

counter=0;

errorVRed=vRed;

errorVIRed=vIRed;

tK2=k-1;

endif

else

v2Red=errorVRed;

v2IRed=errorVIRed;

flag=2;

counter=0;

endif

counter=counter+1;

endif%search minimum

%

case 2

if (vNewIRed>vIRed)

if(counter<stepCheck)

if (vIRed<errorVIRed)

counter=0;

errorVRed=vRed;

errorVIRed=vIRed;

tK3=k-1;

endif

else

if (v1IRed>errorVIRed)

flag=1;

tK1=tK3;

v1Red=errorVRed;

v1IRed=errorVIRed;

else

v3Red=errorVRed;

v3IRed=errorVIRed;

flag=3;

endif

counter=0;

endif

counter=counter+1;

endif

case 3

if (vNewIRed<vIRed)

if(counter<stepCheck)

if (vIRed>errorVIRed)

counter=0;

errorVRed=vRed;

errorVIRed=vIRed;

tK4=k-1;

endif

else

if (v2IRed<errorVIRed)

v4Red=errorVRed;

v4IRed=errorVIRed;

flag=0;

nakR=[nakR v1Red v2Red v3Red v4Red];

nakIR=[nakIR v1IRed v2IRed v3IRed v4IRed];

tK=[tK tK1 tK2 tK3 tK4];

else

v1Red=v3Red;

v1IRed=v3IRed;

v2Red=errorVRed;

v2IRed=errorVIRed;

tK1=tK3;tK2=tK4;flag=2;

endif

counter=0;

endif

counter=counter+1;

endif

endswitch

vRed=vNewRed;

vIRed=vNewIRed;

endfor

%Удаляемпервыйвычесленныйотчет

nakR=nakR(5:end-8);nakIR=nakIR(5:end-8);tK=tK(5:end-8);

RedDc=nakR(1:4:end);

RedAc=nakR(4:4:end)-RedDc;

IRedDc=nakIR(1:4:end);

IRedAc=nakIR(4:4:end)-IRedDc;

Ratio=(RedAc./RedDc)./(IRedAc./IRedDc);

SpO2=104-17*Ratio;

mean_Spo2=sum(SpO2)/length(SpO2);

%Heart rate measurement

HRM_time=tK(8:4:(end-8));

HRM=[];

for m=1:length(HRM_time)-1

HRM=[HRM 60/(HRM_time(m+1)-HRM_time(m))*Fs];

endfor

HRM_mean=mean(HRM);

figure;

plot(0:length(HRM)-1,[HRM],'-ob'),grid,title('HRM')

xlabel('interval,n'),ylabel('HR, beats/min')

legend(['HRM mean=' num2str(HRM_mean)])

figure;

hold on

plot(1:length(ppgIRed30s),fliplr(ppgIRed30s),'-r')

plot(1:length(ppgRed30s),fliplr(ppgRed30s),'-b'),grid

%legend('Red light','IR light')

legend('IR light')

xlabel('t'),ylabel('Light reflection coefficient')

hold off

figure;

hold on

plot(1:length(ppgIRed30s),fliplr(ppgIRed30s),'-r')

plot(length(ppgRed30s)-tK+1,(nakIR),'-ok')

plot(1:length(ppgRed30s),fliplr(ppgRed30s),'-b'),grid

plot(length(ppgRed30s)-tK+1,(nakR),'-ok')

xlabel('t'),ylabel('Light reflection coefficient')

%legend('Red light PPG','Max and Min of func ','IR light PPG')

legend('IR light PPG','Max and Min of func ')

hold off

if 1

figure;

plot(Ratio),grid

endif

if 0

figure;

hold on

plot(nakR),grid

plot(nakIR),grid

title('nakr')

hold off

endif

if 1

figure;

plot(0:length(SpO2)-1,SpO2,'-ob'),grid

xlabel('T,n'),ylabel('Spo2,%')

legend(['mean value Spo2 = ' num2str(mean_Spo2) '%'])

%title('SpO2')

endif

endif

ppgData.ppgRed30s=ppgRed30s;

ppgData.ppgIRed30s=ppgIRed30s;

ppgData.nakR=nakR;

ppgData.nakIR=nakIR;

save -mat7-binary 'ppg_data.mat' 'ppgData';

Размещено на Allbest.ru