Calculate SpO₂ value from raw fingertip plethysmography (PPG)

The question seems rather simple, however I have yet to figure out the right method for calculating the SpO₂ value.

The formula for SpO₂ is: $$R = {\mathrm{AC_{red} / DC_{red}} \over \mathrm{AC_{ir} / DC_{Ir}}}$$ and $$\mathrm{SpO_{2}} = A - B \cdot R$$

So, the question is, how to obtain the correct value for $AC$ and $DC$?

So far, all articles I red seem to agree on what the AC has to do with. The AC is the fluctuating part of the signal. e.g. $$x = A + \sin(B \cdot x)$$ The AC component has to do with the $\sin(B \cdot x)$ and not the $A$. But, what is the correct value for $AC$? I've seen people take the maximum points and the minimum points of the fluctuation, and subtract the max from the min. But, "How to Design Peripheral Oxygen Saturation (SpO₂) and Optical Heart Rate Monitoring (OHRM) Systems Using the AFE4403" document from TI, uses the RMS (root mean square) of "$AC$".

Furthermore, people don't seem to agree on what the $DC$ component is. Most of the articles say the $DC$ value is the bottom line of the fluctuation: they claim the $DC$ component is the static absorption by bone, tissue and muscle. The above document implies that the $DC$ value is the "offset" of the signal (thus the $A$ in the above formula). Masimo (Pleth Variability Index: A Dynamic Measurement to Help Assess Physiology and Fluid Responsiveness) (page 4) even has a document that implies the $DC$ value is from zero till the top of the $AC$ envelope, however, they use that variable to calculate the Perfusion index, so i guess that's just an unfortunate naming of variable.

The above image displays the raw data (50Hz) recorded with an SpO2 device, recorded from a "ProSim SPOT Light SpO₂ Pulse Oximeter Analyzer" which is set to 80BMP, 100% SpO₂, 10% PI. Using Matlab, I calculated the "envelope" of the fluctuating signal. (top and bottom horizontal line in the image).

The question: How do I calculate the SpO₂ from here?

• Do I use $DC_{1}$ or $DC_{2}$? in which $DC_{2}$ is the mean.
• In case of $DC_{2}$, should I take the mean of the signal (over a window of x seconds) or the mean of the envelope?
• Is the $AC$ value: the max envelope value - the min envelope value or should (instead of envelope) I use the RMS of the signal (after applying a 0.5-5Hz bandpass filter over the signal).

Note: I only showed the signal of the red sensor. Obviously the same process needs to be done for the IR sensor.

------------ Edit ------------

The question is which AC value are you after? That is indeed the question...

I've red that it is possible to get the SpO₂ using FFT as well. In this method, you'll get a window of X seconds, then obtain the value for 0Hz (this is DC) as well as the max value in 0.5 - 3Hz (this is AC). The advantage of this method is that only the Heart Rate frequenties are used, and thus the method is less sensitive to noise.

So, I've made a simple test in Matlab:

Offset = 15;
Fs = 50;                      % Sampling frequency (samples per second)
dt = 1/Fs;                    % seconds per sample
Duration = 10;                % seconds
t = (0:dt:Duration-dt)';      % seconds
F = 80/60;                    % Sine wave frequency (hertz) 80 BMP
data = Offset + sin(2*pi*F*t);

%% Apply low pass filter to obtain DC
DCl = Offset - 1;
DC = Offset;

%%
disp( ['Lim Sup/Inf: ' sprintf('%f', 2/DC) ] )
disp( ['Lim Sup/Inf / DCl: ' sprintf('%f', 2/DCl) ] )

%% Apply bandpass filter to obtain AC
fcl = 0.25; % Hz
fch = 3.25; % Hz
[b,a] = butter( 6, [fcl/(Fs/2) fch/(Fs/2)], 'bandpass' );
AC = filtfilt( b, a, data );
AC = rms( AC );
disp( ['Bandpass: ' sprintf('%f', AC/DC) ] )
disp( ['Bandpass / DCl: ' sprintf('%f', AC/DCl) ] )

% or perhaps highpass
fcl = 0.25; % Hz
[b,a] = butter( 6, fcl/(Fs/2), 'high' );
AC = filtfilt( b, a, data );
AC = rms( AC );
disp( ['High: ' sprintf('%f', AC/DC) ] )
disp( ['High / DCl: ' sprintf('%f', AC/DCl) ] )

%% FFT method to get AC/DC
L = length( data );
NFFT = 2^nextpow2( L );
Y = abs( fft( data, NFFT ) );

%Find local maximum in Heart Rate spectrum
DC = Y(1);
[AC,~] = max( Y( 6:31 ) ); % 0.5-3Hz

disp( ['FFT: ' sprintf('%f', AC/DC) ] )


Results

Lim Sup/Inf:       0.133333
Lim Sup/Inf / DCl: 0.142857
Bandpass:          0.047864
Bandpass / DCl:    0.051283
High:              0.047018
High / DCl:        0.050377
FFT:               0.045220


If I assume that the FFT method is valid. Than, it follows that using the limit superior/inferior results in "significant" different results.

However...

• R is the ratio of ratios.
• R-table needs to be calibrated for each method individually (which converts and R value to a SpO₂ value)
• AC and DC are calculated the same way for the Red and Infra Red signal

so, does it even matter what I use for AC and DC?

• I think that the guiding principle here should be the accuracy of the SPO2 estimate and the risk an over-estimation or under-estimation brings with it. A "Consumer" device is different than a "medical" device. This is what gives meaning to "Significant".
– A_A
Jan 25, 2018 at 9:34
• I'm using a custom designed device. However, using its output, I get weird, non-consistent, results. This leads me to think that either I have made a mistake in my calculation method, or there is a mistake in the firmware. The "results" in the question show x3 values for the lim sup/inf. But if Cedron Dawg is correct, the ratio of ratios should not matter. The results I've so far, for same settings with a Fluke device, give different R values for measurement in Dark and Light enviroment. This makes me believe stronger that the firmware is the problem, assuming that my code is correct. Jan 25, 2018 at 11:17
• electronics.stackexchange.com/questions/252898/pulse-oximeter-spo2-calculation <== I did a little research on the topic and found this. Seemed sensible to me. Jan 26, 2018 at 4:43
• @CedronDawg I did came across that post a long time ago. I re-red it, with all the articles and did come across some useful information. Furthermore, I received a firmware update from the manufacturer: turns out they normalise the signal internally.... sigh. I'll implement a beat-to-beat method (as you describe) and use that to obtain the R = (AC{_r} / AC{_i}) in my case... Your short answer kept me sane the past days: thank you for that! Jan 26, 2018 at 15:27
• @Jeffrey, Thank you for the bounty. Good luck on your project. It would be nice to know how it turns out. My email address is cedron at protonmail dot com. Feel free to also ask further questions there. Jan 26, 2018 at 16:31

There are a couple of important considerations you have left out.

What granularity is required? It seems to me the most granular reading you expect are from a single cycle. Yet you mention frame sizes, so what is acceptable?

What precision is needed?

Disclaimer, I have no experinece with SpO2 processing. Having said that I can make some observations which may be helpful.

You state: However... R is the ratio of ratios.
R-table needs to be calibrated for each method individually (which converts and R value to a SpO₂ value) AC and DC are calculated the same way for the Red and Infra Red signal so, does it even matter what I use for AC and DC?

Short answer: Roughly not, but some.

The ratio of ratios observations makes it possible to use your "RMS" (which is really a standard deviation of the signal as you define it) or a difference of extremes as they should be proportional.

The trending of the signal is going to throw off any frame based method. With that in mind, this is how I would do it, looking at your graph as a guide.

For each peak:
Find the value call it P
Find the value of the previous trough, call it Tp
Find the value of the next trough, call it Tn
Use the average of the trough values (Tp+Tn)/2 call it T
Calculate AC as P-T
Calculate DC as (P+T)/2


Process the signals simultaneously to calculate your ratio of ratios to generate a new sequence of values.

Use a short term smoothing filter, maybe five points wide to smooth out your sequence. Maybe include a rejection criteria for readings that are "out of whack".

No, this answer doesn't meet your bounty requirements, but I hope it helps.

Ced

• Do I use $DC_1$ or $DC_2?$ in which $DC_2$ is the mean.

From a DSP point of view, direct current is a constant component around which there may be additive variation. The closest to this is what you are calling $DC_2$. What you are calling $DC_1$ is the limit inferior.

The estimation of an average value implies that you derive the mean over a time interval. For example, in deriving the expected value of a process, the interest is what is the value of the mean as time tends to infinity.

So, to estimate $SpO_2$, you are going to have to integrate that waveform over some time interval. Think of this as maintaining a running sum of the value of the waveform over a finite length window.

From this point of view, $DC_2$ will emerge naturally as the values that "drop-out" of the moving buffer do not contribute to the sum anymore and the new values oscillate around some average. Therefore, this running sum will "hover" around a value. This however would be a "value over a specific amount of time". This amount of time is the time constant of the integrator.

• In case of $DC_2$, should I take the mean of the signal (over a window of x seconds) or the mean of the envelope?

• Is the $AC$ value: the max envelope value - the min envelope value or should (instead of envelope) I use the RMS of the signal (after applying a 0.5-5Hz bandpass filter over the signal).
The question is which $AC$ value are you after? Peak Amplitude? Peak-To-Peak Amplitude? or RMS Amplitude?.
Again, if you are after an estimation over a window of time, you can subtract your running mean (from above) and then process what is left (it will be $\approx\pm 0.15$).