11#include " BrainflowSpO2Algorithm.h"
22
3- #define MAX_FILTER_ORDER 8
4-
53/* *
64 * @brief Calculates the oxygen saturation level (SpO2) from PPG IR and RED signals.
75 *
1513 */
1614void get_oxygen_level (float *ppg_ir, float *ppg_red, int data_size, float *oxygen_level)
1715{
18- float *red_raw = new float [data_size];
19- float *ir_raw = new float [data_size];
20- int filter_shift = 25 ; // to get rif of filtereing artifact, dont use first elements
21- int new_size = data_size - filter_shift;
22- float *new_red_raw = red_raw + filter_shift;
23- float *new_ir_raw = ir_raw + filter_shift;
24- memcpy (red_raw, ppg_red, data_size * sizeof (float ));
25- memcpy (ir_raw, ppg_ir, data_size * sizeof (float ));
16+ float *red_raw = new float [data_size];
17+ float *ir_raw = new float [data_size];
18+ int filter_shift = 25 ; // to get rif of filtereing artifact, dont use first elements
19+ int new_size = data_size - filter_shift;
20+ float *new_red_raw = red_raw + filter_shift;
21+ float *new_ir_raw = ir_raw + filter_shift;
22+ memcpy (red_raw, ppg_red, data_size * sizeof (float ));
23+ memcpy (ir_raw, ppg_ir, data_size * sizeof (float ));
2624
27- // need prefiltered mean of red and ir for dc
28- float mean_red = mean (new_red_raw, new_size);
29- float mean_ir = mean (new_ir_raw, new_size);
25+ // need prefiltered mean of red and ir for dc
26+ float mean_red = mean (new_red_raw, new_size);
27+ float mean_ir = mean (new_ir_raw, new_size);
3028
31- // filtering(full size)
32- detrend (red_raw, data_size, (int )DetrendOperations::CONSTANT );
33- detrend (ir_raw, data_size, (int )DetrendOperations::CONSTANT );
34- perform_bandpass (red_raw, data_size, FILTER_SAMPLING_RATE , 0.7 , 1.5 , 4 , (int )FilterTypes::BUTTERWORTH , 0.0 );
35- perform_bandpass (ir_raw, data_size, FILTER_SAMPLING_RATE , 0.7 , 1.5 , 4 , (int )FilterTypes::BUTTERWORTH , 0.0 );
29+ // filtering(full size)
30+ detrend (red_raw, data_size, (int )DetrendOperations::CONSTANT );
31+ detrend (ir_raw, data_size, (int )DetrendOperations::CONSTANT );
3632
37- // calculate AC & DC components using mean & rms:
38- float redac = rms (new_red_raw, new_size);
39- float irac = rms (new_ir_raw, new_size);
40- float reddc = mean_red;
41- float irdc = mean_ir;
33+ const double start_f = 0.7 ;
34+ const double stop_f = 1.5 ;
4235
43- // https://www.maximintegrated.com/en/design/technical-documents/app-notes/6/6845.html
44- float r = (redac / reddc) / (irac / irdc);
45- float spo2 = (CALIB_COEFF_1 * r * r) + (CALIB_COEFF_2 * r) + (CALIB_COEFF_3 );
46- if (spo2 > 100.0 )
47- {
48- spo2 = 100.0 ;
49- }
50- if (spo2 < 0 )
51- {
52- spo2 = 0.0 ;
53- }
54- *oxygen_level = spo2;
36+ auto start_filter = butter<FILTER_ORDER >(2 * start_f / FILTER_SAMPLING_RATE );
37+ auto stop_filter = butter<FILTER_ORDER >(2 * stop_f / FILTER_SAMPLING_RATE );
38+
39+ for (int i = 0 ; i < data_size; i++) {
40+ red_raw[i] = red_raw[i] - start_filter (red_raw[i]);
41+ red_raw[i] = stop_filter (red_raw[i]);
42+ }
5543
56- delete[] red_raw;
57- delete[] ir_raw;
44+ // TODO better way to reset internal state of filters?
45+ start_filter = butter<FILTER_ORDER >(2 * start_f / FILTER_SAMPLING_RATE );
46+ stop_filter = butter<FILTER_ORDER >(2 * stop_f / FILTER_SAMPLING_RATE );
47+
48+ for (int i = 0 ; i < data_size; i++) {
49+ ir_raw[i] = ir_raw[i] - start_filter (ir_raw[i]);
50+ ir_raw[i] = stop_filter (ir_raw[i]);
51+ }
52+
53+ // calculate AC & DC components using mean & rms:
54+ float redac = rms (new_red_raw, new_size);
55+ float irac = rms (new_ir_raw, new_size);
56+ float reddc = mean_red;
57+ float irdc = mean_ir;
58+
59+ // https://www.maximintegrated.com/en/design/technical-documents/app-notes/6/6845.html
60+ float r = (redac / reddc) / (irac / irdc);
61+ float spo2 = (CALIB_COEFF_1 * r * r) + (CALIB_COEFF_2 * r) + (CALIB_COEFF_3 );
62+ if (spo2 > 100.0 )
63+ {
64+ spo2 = 100.0 ;
65+ }
66+ if (spo2 < 0 )
67+ {
68+ spo2 = 0.0 ;
69+ }
70+ *oxygen_level = spo2;
71+
72+ delete[] red_raw;
73+ delete[] ir_raw;
5874}
5975
6076void detrend (float *data, int data_len, int detrend_operation)
@@ -99,57 +115,4 @@ void detrend (float *data, int data_len, int detrend_operation)
99115 data[i] = data[i] - (grad * i + y_int);
100116 }
101117 }
102- }
103-
104- void perform_bandpass (float *data, int data_len, int sampling_rate, float start_freq, float stop_freq, int order, int filter_type, float ripple)
105- {
106- float center_freq = (start_freq + stop_freq) / 2.0 ;
107- float band_width = stop_freq - start_freq;
108- Dsp::Filter *f = NULL ;
109- float *filter_data[1 ];
110- filter_data[0 ] = data;
111-
112- switch (static_cast <FilterTypes> (filter_type))
113- {
114- case FilterTypes::BUTTERWORTH :
115- f = new Dsp::FilterDesign<Dsp::Butterworth::Design::BandPass<MAX_FILTER_ORDER >, 1 > ();
116- break ;
117- case FilterTypes::BUTTERWORTH_ZERO_PHASE :
118- f = new Dsp::FilterDesign<Dsp::Butterworth::Design::BandPass<MAX_FILTER_ORDER >, 1 > ();
119- break ;
120- case FilterTypes::CHEBYSHEV_TYPE_1 :
121- f = new Dsp::FilterDesign<Dsp::ChebyshevI::Design::BandPass<MAX_FILTER_ORDER >, 1 > ();
122- break ;
123- case FilterTypes::CHEBYSHEV_TYPE_1_ZERO_PHASE :
124- f = new Dsp::FilterDesign<Dsp::ChebyshevI::Design::BandPass<MAX_FILTER_ORDER >, 1 > ();
125- break ;
126- case FilterTypes::BESSEL :
127- f = new Dsp::FilterDesign<Dsp::Bessel::Design::BandPass<MAX_FILTER_ORDER >, 1 > ();
128- break ;
129- case FilterTypes::BESSEL_ZERO_PHASE :
130- f = new Dsp::FilterDesign<Dsp::Bessel::Design::BandPass<MAX_FILTER_ORDER >, 1 > ();
131- break ;
132- }
133-
134- Dsp::Params params;
135- params[0 ] = sampling_rate; // sample rate
136- params[1 ] = order; // order
137- params[2 ] = center_freq; // center freq
138- params[3 ] = band_width;
139- if ((filter_type == (int )FilterTypes::CHEBYSHEV_TYPE_1 ) ||
140- (filter_type == (int )FilterTypes::CHEBYSHEV_TYPE_1_ZERO_PHASE ))
141- {
142- params[3 ] = ripple; // ripple
143- }
144- f->setParams (params);
145- f->process (data_len, filter_data);
146- if ((filter_type == (int )FilterTypes::BUTTERWORTH_ZERO_PHASE ) ||
147- (filter_type == (int )FilterTypes::CHEBYSHEV_TYPE_1_ZERO_PHASE ) ||
148- (filter_type == (int )FilterTypes::BESSEL_ZERO_PHASE ))
149- {
150- reverse_array (data, data_len);
151- f->process (data_len, filter_data);
152- reverse_array (data, data_len);
153- }
154- delete f;
155118}
0 commit comments