PPGv2

[mark9064]

• New heart rate algorithm
  • Adaptive filtering to remove noise caused by motion
  • Harmonic analysis to match measured waveform with a real PPG waveform to reject noise
  • Frequency tracking and refinement over time
• Heart rate detection can now function in more challenging environments e.g. while running
• Full automatic gain and drive implementation for HRS3300 resulting in power saving in most scenarios while being able to adapt to full sunlight conditions on light skin, also adapts to darker skin tones better
  • Including no touch detection
• Much better noise tolerance: this algorithm reports failure instead of an incorrect heart rate value when the signal is poor. The heart rate values reported can be trusted a lot more
  • Algorithm accuracy validated against an ECG chest strap (Polar H10)
  • Also good performance on PPG-DaLiA (public PPG dataset): Mean absolute error (MAE) ~8BPM

Just to set expectations, while this is a solid improvement the HRS3300 sensor is really not that great so this implementation is nowhere the performance of flagship smartwatches today in terms of accuracy and recall

Big thanks to everyone over at wasp-os who investigated the sensor performance, this work wouldn't have been possible without

I'm not sure it will be possible to merge this due to the memory usage of the adaptive filter. It might conflict with the G7710 watchface (I haven't tested as I use a custom G7710 version with seconds, so the fonts are smaller and use less memory)

This changes the way heart rate values are presented. With this implementation, you will only see a number when the algorithm is actively tracking the heart rate, or if the last background measurement was successful and the algorithm does not yet have enough data to begin searching for a heart rate (8 seconds of data)

[mark9064]

Oh, this needs GCC 14 as well

[github-actions]

Build checks have not completed. Possible reasons for this are:

1. The checks need to be approved by a maintainer
2. The branch has conflicts
3. The firmware build has failed

[mark9064]

GCC 14 here #2372, not going to rebase on to it as I don't think that will fix CI

[tituscmd]

Ooooh, I like this a lot.

If this is ready to test on-device, I would be more than willing to daily drive it for a while and report its results!

[mark9064]

Yeah this is ready, I've been working on this for >2 years so I'm pretty happy with everything here (of course there's gonna be review feedback etc. but the design of the detector is done in my eyes)

[tituscmd]

I'm gonna merge both GCC 14 and this PR onto my main device and run it for a while. I'll report back here when something noteworthy happens!

[tituscmd]

I’m particularly excited about my gym session and the bike ride there tomorrow. Interested to see how it performs there!

[tituscmd]

Frequency tracking and refinement over time

What does this mean? I've had a few measurements spaced apart more than my background measurement interval setting which confused me.

[minacode]

Wow, good work!

Including no touch detection

Can we now reliably detect if the watch is not worn?

[mark9064]

Frequency tracking and refinement over time

What does this mean? I've had a few measurements spaced apart more than my background measurement interval setting which confused me.

The heart rate is tracked internally by the algorithm rather than smoothing over a buffer of the last values so accuracy is a bit better. There's no outward behaviour change

If measurements are spaced farther than the background interval, it's probably because a measurement failed

[mark9064]

Wow, good work!

Including no touch detection

Can we now reliably detect if the watch is not worn?

I'm not sure I'd call it super reliable, any surface reflects light in the same way as skin. But it does catch the sensor not being near any surface

[Andersama]

I should've known to check the PR's for interesting changes. I'd love to test this.

[minacode]

I would like to review this from a more theoretical perspective. Can you do a little write-up of the logic behind this? I tried to get it from the code, but got a little overwhelmed 😄

[mark9064]

Yeah it's quite a lot, and understanding each step exactly requires some DSP knowledge that can be quite tricky to learn (at least it was for me!)

The pipeline has two overall stages, the first of which runs for every ingested sample and the second of which is evaluated on the first stage output every half of a second.

Stage 1:

• Pass the incoming tuples of (reflectance, accx, accy, accz) through individual high-pass filters (high-pass filters transmit higher frequencies and attenuate lower ones; the cut off frequency and steepness of the cut off is a choice made when designing the filter)
  • Why? Removes DC baseline from the signal and other long term trends (low frequency components)
• Use the filtered acceleration channels acc to estimate noise in the reflectance signal and remove it (this is performed by the adaptive filter)
  • Why? Removes (some) noise caused by motion artefacts, which is the biggest contributor to the reflectance signal

We now have the filtered reflectance signal which we stream into an 8s long buffer. Twice per second:
(Stage 2)

• Energy normalise a copy of the buffer data
  • Why? Provides input scale invariance, i.e. tolerance against different skin reflectivities, different sensor drive/gain modes, different ambient lighting conditions
• Transform it to the frequency domain using the Fourier transform
  • Why? Heart rate BVP signal is periodic so easier to identify in the frequency domain
• Detect peaks in the signal with a prominence filter (bandpass the spectrum)
  • Why? Suppresses wideband spectral noise caused by discontinuties in the signal (usually transient motion artefacts)
• Filter the peaks based on whether the signal matches the well known shape of a heart rate signal at that frequency
  • Why? Attenuates peaks caused by things other than the BVP
• Threshold resulting filtered spectrum to find strong peaks
  • Why? Want to avoid tracking peaks with little energy (weak peaks) as these are more likely to be noise
• Accept strong peaks and track peaks over time
  • Why? Strength of the true HR peak will vary over time as signal quality changes, so once it has been detected tracking it over time provides much better recall (as opposed to only reporting a value when a strong enough peak is present)

Aside: The behaviour around the stopped state is a bit broken currently, will fix when I have some time

[mark9064]

(Will squash fixups down later, I know a few people have merged this in locally so trying to avoid rebasing for now)

[Andersama]

@mark9064 attempted a build merging my commit as well, it looks like there's a lot of compiler errors. I'll patch them up and give this a shot.

Ok from an initial quick overview: I'd replace the lines:

std::array<std::complex<float>, IntegerLog2(N)> result;

with

std::array<std::complex<float>, IntegerLog2(N)> result{};

I don't know about GCC, but I've had plenty of experience with constexpr functions failing because variables aren't initialized*. It's a bit silly but you should {} even though the values will be overwritten in the for loop.

and then I've never seen this as an expression, might be a gcc extension? I can't compile this with GCC on godbolt.

-2.i * std::numbers::pi

my best guess is

-2. * i * std::numbers::pi

I haven't touched FFTs in ages so I'd have to work out what's going on here.

I'd drop bringing in <numbers> just for the pi constant. That seems to resolve the issues in the FFT section.

Update: Assuming the guess above was ok, I made edits that got it to compile and build. I'm not sure that it was a correct assumption because the samples I'm getting back are anywhere between 150 to 200 bpm, which is a bit high considering I was maybe at 70-90 bpm.

To clarify this is buildable with a few small edits. The constexpr errors can be resolved by doing the calculations needed in a constructor for the std::complex<float>s.

[Andersama]

@mark9064 I'm going to guess I messed up the FFT in some way trying to patch up the errors around std::complex, it doesn't appear to me that the code would need a compiler version bump, for reference here were my edits off of this PR: https://github.com/Andersama/InfiniTime/tree/ppg_upgrade

[mark9064]

Does it compile for you with GCC14? It does on my machine

Upgrading the compiler is not a big deal, so I wouldn't worry too much about trying to workaround unimplemented functionality. We try to keep the compiler up to date anyway (as well as other dependencies) but haven't recently due to some linker warnings (which we now understand)

I appreciate you looking into this though :)

[mark9064]

Also anyone got any feedback? Criticism is welcome too, this approach is no good if it only works well for me :)

[Andersama]

@mark9064 I'm building off of the included docker image, so whichever GCC's included there, it wouldn't compile, but also I tested the same expression on godbolt (with GCC14 and trunk) https://godbolt.org/z/5hnEo5ej6 and it also wouldn't compile so I'm not clear on what -2.i * std::numbers::pi was supposed to do, because to me, that looks like an error, and it appears to be treated as an error. If it were me I'd just write whatever mathematically works out to calculate the twiddle values into the std::complex<float> constructor for the real and imaginary components, because float operations have been constexpr for quite a while and the constructor is simple enough it didn't appear GCC complained about that approach when I tried it. Although I probably messed up the math.

I did with some fiddling and edits and managed to get it to compile with whichever GCC's in use now...overall I'm just not 100% certain I did quite the right thing, because the heart rate values I was getting back were wildly off.

Also: It's the same errors that are appearing in the build checks for the PR here: https://github.com/InfiniTimeOrg/InfiniTime/actions/runs/19586678888/job/56096912190?pr=2371
I wanted to test this out so, I at least got past the compiler errors, but I'm not sure if I'm dealing with a semantic error.

[mark9064]

Use #2372 and rebuild the docker container from there

The compiler error in Godbolt tells you exactly what you need to add :) and you'll see that using namespace std::complex_literals; is in the original code too (this provides i so inline imaginary constants are easy to write)

[Andersama]

Ah, well serves me right for not paying attention to the godbolt logs, weird that it's not the same error. Ah, you know what I bet it's buried behind the constexpr failure. The simplified expression in godbolt I dropped constexpr in the surrounding function so that's probably the actual issue, no contexpr i operator.

Oh, and that probably explains my bad results, because I thought i was the loop variable...so my twiddle values in my code are probably garbage.

[Andersama]

@mark9064 slight update, fixed up the apparent math error, compiles and appears to work as intended. I'll have to grab my other heart rate sensor to verify, but the numbers make much more sense now, currently recording ~80 bpm. I rewrote the offending blocks like so:

      template <std::size_t N>
      static consteval std::array<std::complex<float>, IntegerLog2(N)> GenComplexTwiddle() {
        using namespace std::complex_literals;

        std::array<std::complex<float>, IntegerLog2(N)> result {};
        for (std::size_t i = 0; i < IntegerLog2(N); i++) {
          std::complex<double> tmp = std::complex<double> {0.0, -2.0} * std::numbers::pi / static_cast<double>(1 << (i + 1));
          std::complex<float> value = exp_consteval(std::complex<float> {(float) tmp.real(), (float) tmp.imag()});
          result[i] = value;
        }
        return result;
      }

      template <std::size_t N>
      static consteval std::array<std::complex<float>, (N / 4) - 1> GenRealTwiddle() {
        using namespace std::complex_literals;

        std::array<std::complex<float>, (N / 4) - 1> result {};
        for (std::size_t i = 0; i < ((N / 4) - 1); i++) {
          std::complex<double> tmp =
            std::complex<double> {0.0, -2.0} * std::numbers::pi * static_cast<double>(i + 1) / static_cast<double>(N);
          std::complex<float> value = exp_consteval(std::complex<float> {(float) tmp.real(), (float) tmp.imag()});
          result[i] = value;
        }
        return result;
      }

There were more errors that cropped up since the compiler didn't appreciate the conversion from std::complex<double> to std::complex<float> so I had to break the expression apart into steps. Apparently there's an option to use if as a literal operator to generate std::complex<float> from a literal instead of a double which probably would simplify this mess...

I suppose the only thing I would wonder is if a windowing function was used, it seemed like the old arduino fft based code just performed an fft with no windowing function at all, and I'm wondering if that would've helped.

Update, after sitting with a finger ppg sensor, this seems very accurate, I'd say it almost stays in complete sync with the finger sensor, except for the random blips where it drops to 0, will be testing it on a workout later.

Further update: Two different workouts later, did very poorly on an outdoor bike ride most of the readings were 0's virtually unusable for feedback for how I was doing. Second one indoors did much much better, only a few samples of 0s across a whole hour, and was responsive enough I actually managed to complete a workout with set target heartrates, so very nice...in good conditions. My best guess, lighting and potentially an amount of motion is a significant issue.