Estimating heart rate and SpO2 levels using PPG (photoplethysmography) (Matlab code implementation)

Welcome to this blog

Advantages of bloggers:Blog content should be as thoughtful and logical as possible for the convenience of readers.

Motto:He who travels a hundred miles is half as good as ninety.

The directory of this article is as follows:

Table of Contents

1 Overview

2 Operation results

3 References

4 Matlab code implementation


1 Overview

Heart rate and SpO2 levels can be accurately estimated using PPG (photoplethysmography) technology. PPG is a non-invasive biometric technology that obtains heart rate and oxygen saturation information by measuring changes in light absorption of blood on the skin surface.

During the measurement process, the PPG device captures the reflection and absorption of light by shining light at wavelengths of 940 nanometers and 660 nanometers onto the skin. The absorption properties of hemoglobin cause the amount of light absorption in the blood to change during heartbeat, which can be used to estimate heart rate. In addition, oxyhemoglobin and deoxygenated hemoglobin absorb different wavelengths of light in different proportions, so the SpO2 level can be estimated based on this difference.

To improve accuracy, stable signal quality and accurate data processing methods are required. Common methods include filtering, peak detection, time domain and frequency domain analysis, etc. Through these processing steps, accurate and reliable heart rate and SpO2 level estimation results can be obtained.

Additionally, there are some key considerations and technical details to consider. For example, keeping the measured object still to avoid motion artifacts, ensuring that the PPG sensor is positioned correctly against the skin, and resisting ambient light interference, etc. These factors can affect the measurement results.

Overall, PPG is a very useful technique to estimate heart rate and SpO2 levels in a non-invasive and convenient way. It has broad application prospects in medical care, sports health monitoring, sleep quality assessment and other fields, and provides people with better biological monitoring and health management methods.

2 Operation results

Part of the code:

%% Data frame
for n=1:fix((length(X)/(2*fs))-2)
y1=X(n*fs:(n*fs + FFT_size-1),1); %RED
y2=X(n*fs:(n*fs + FFT_size-1),2); %IR

%% FFT Transform RED
NFFT = FFT_size; % Next power of 2 from length of y
Y1 = fft(y1,NFFT);
f1 = fs/2*linspace(0,1,NFFT/2 + 1);

figure(1)
plot(f1,abs(Y1(1:NFFT/2 + 1)));
axis([0.5 2.5 0 3e5])

%% FFT Transform IR
NFFT = FFT_size; % Next power of 2 from length of y
Y2 = fft(y2,NFFT);
f2 = fs/2*linspace(0,1,NFFT/2 + 1);

figure(2)
plot(f2,abs(Y2(1:NFFT/2 + 1)));
axis([0.5 2.5 0 2e5])
hold on

%% Find local maximum in RED spectrum

YY=abs(Y1(6:12));
local_max_i=1;
local_max=YY(1);
for i=2:(length(YY)-1)
    if local_max<(YY(i))
        local_max_i=i;
        local_max=YY(i);
    end
end
pk_RED_i=6-1 + local_max_i;

%% Find local maximum in IR spectrum

YY=abs(Y2(6:12));
local_max_i=1;
local_max=YY(1);
for i=2:(length(YY)-1)
    if local_max<(YY(i))
        local_max_i=i;
        local_max=YY(i);
    end
end
pk_IR_i=6-1 + local_max_i;

%%Heart rate
HEART_RATE(n) = f2(pk_IR_i)*60 %%In fact, using FFT limits the accuracy of heart rate estimation. See the points on f1/ f2 arrays and you know why. I wrote a peak detection algorithm for heart rate only. See the second .m file.

%%SpO2
R_RED = abs(Y1(pk_RED_i)/abs(Y1(1)));
R_IR = abs(Y2(pk_IR_i)/abs(Y2(1)));
R=R_RED/R_IR;
SpO2(n) = 104 - 28*R

end

%% Take average value of heart rate and SpO2
HR=sum(HEART_RATE(2:(length(HEART_RATE)-1)))/(length(HEART_RATE)-2);
S=sum(SpO2(2:(length(SpO2)-1)))/(length(SpO2)-2);
Heart_Rate=round(HR)
SpO2_Level=round(S)

%% Plot result
y1=X(:,1);
y2=X(:,2);

% Denoising
for i=1:(length(y1)-1)
    if ((y1(i + 1)-y1(i))>50000)
        y1(i + 1)=y1(i + 1) + 65100;
    elseif ((y1(i)-y1(i + 1))>50000)
        y1(i + 1)=y1(i + 1) + 65100;
    end
end

for i=1:(length(y2)-1)
    if ((y2(i + 1)-y2(i))>50000)
        y2(i + 1)=y2(i + 1) + 65100;
    elseif ((y2(i)-y2(i + 1))>50000)
        y2(i + 1)=y2(i + 1) + 65100;
    end
end

x=0:1:length(y1)-1;
plot(x,y1 + 3e5,'r');
hold on
plot(x,y2,'b');
legend('RED','IR');
legend('boxoff');
str1 = ['Heart rate = ',num2str(HR)];
text(length(y2)/2,max(y2)*4/5,str1);

%%Data frame
for n=1:fix((length(X)/(2*fs))-2)
y1=X(n*fs:(n*fs + FFT_size-1),1); %RED
y2=X(n*fs:(n*fs + FFT_size-1),2); %IR

%% FFT Transform RED
NFFT = FFT_size; % Next power of 2 from length of y
Y1 = fft(y1,NFFT);
f1 = fs/2*linspace(0,1,NFFT/2 + 1);

figure(1)
plot(f1,abs(Y1(1:NFFT/2 + 1)));
axis([0.5 2.5 0 3e5])

%% FFT Transform IR
NFFT = FFT_size; % Next power of 2 from length of y
Y2 = fft(y2,NFFT);
f2 = fs/2*linspace(0,1,NFFT/2 + 1);

figure(2)
plot(f2,abs(Y2(1:NFFT/2 + 1)));
axis([0.5 2.5 0 2e5])
hold on

%% Find local maximum in RED spectrum

YY=abs(Y1(6:12));
local_max_i=1;
local_max=YY(1);
for i=2:(length(YY)-1)
if local_max<(YY(i))
local_max_i=i;
local_max=YY(i);
end
end
pk_RED_i=6-1 + local_max_i;

%% Find local maximum in IR spectrum

YY=abs(Y2(6:12));
local_max_i=1;
local_max=YY(1);
for i=2:(length(YY)-1)
if local_max<(YY(i))
local_max_i=i;
local_max=YY(i);
end
end
pk_IR_i=6-1 + local_max_i;

%%Heart rate
HEART_RATE(n) = f2(pk_IR_i)*60 %%In fact, using FFT limits the accuracy of heart rate estimation. See the points on f1/ f2 arrays and you know why. I wrote a peak detection algorithm for heart rate only. See the second .m file.

%%SpO2
R_RED = abs(Y1(pk_RED_i)/abs(Y1(1)));
R_IR = abs(Y2(pk_IR_i)/abs(Y2(1)));
R=R_RED/R_IR;
SpO2(n) = 104 – 28*R

end

%% Take average value of heart rate and SpO2
HR=sum(HEART_RATE(2:(length(HEART_RATE)-1)))/(length(HEART_RATE)-2);
S=sum(SpO2(2:(length(SpO2)-1)))/(length(SpO2)-2);
Heart_Rate=round(HR)
SpO2_Level=round(S)

%% Plot result
y1=X(:,1);
y2=X(:,2);

% Denoising
for i=1:(length(y1)-1)
if ((y1(i + 1)-y1(i))>50000)
y1(i + 1)=y1(i + 1) + 65100;
elseif ((y1(i)-y1(i + 1))>50000)
y1(i + 1)=y1(i + 1) + 65100;
end
end

for i=1:(length(y2)-1)
if ((y2(i + 1)-y2(i))>50000)
y2(i + 1)=y2(i + 1) + 65100;
elseif ((y2(i)-y2(i + 1))>50000)
y2(i + 1)=y2(i + 1) + 65100;
end
end

x=0:1:length(y1)-1;
plot(x,y1 + 3e5,’r’);
hold on
plot(x,y2,’b’);
legend(‘RED’,’IR’);
legend(‘boxoff’);
str1 = [‘Heart rate = ‘,num2str(HR)];
text(length(y2)/2,max(y2)*4/5,str1);

3 References

Some of the content in the article is quoted from the Internet. The source will be indicated or cited as a reference. It is inevitable that there will be some unfinished information. If there is anything inappropriate, please feel free to contact us to delete it.

[1] Zhang Lieliang, Zhu Juan, Xu Lei. Research progress on clinical application of photoplethysm wave[J]. Journal of Clinical Anesthesiology, 2013, 29(11):3.DOI:CNKI:SUN:LCMZ.0.2013-11-040 .

[2] Wang Nan, Wang Xia, Fu Xiaojing, et al. Blind signal estimation method for photoplethysmographic wave extraction from video [J]. Electronic Design Engineering, 2017(3):5.DOI:CNKI:SUN:GWDZ.0.2017 -03-043.

[3] Yu Siyuan, Wei Wei. Research on the feasibility of using photoplethysmography (PPG) to evaluate postoperative pain [C]//National Academic Forum of Young Anesthesiologists. 2012.

4 Matlab code implementation

The knowledge points of the article match the official knowledge files, and you can further learn related knowledge. Algorithm skill tree Home page Overview 55,296 people are learning the system