Gettin low pass filter into arduino sketch

[ESCHICK]

Hello!

Can someone help me to get the low pass filter into my arduino sketch? Maybe without libarys, simply c++ code?

DC Removal

There are two things you should notice in the graph (figure 3):

1. The graph is oscillating slightly
2. It has a DC offset of 50 000 units

To properly be able to read the heart rate and SaO2 we need to remove the DC signal and leave only the AC part.

It is actually very simple and can be done using these two equations:

w(t)=x(t)+∝∗w(t−1)

y(t)=w(t)−w(t−1)

y(t): is the output of the filter
x(t): current input/value
w(t): intermediate value, acts like the history of the DC value
α: is the response constant of the filter
If α = 1 then everything passes through
If α = 0 then nothing passes through
for DC removal you want the α as rather close to 1. I’ll be using α = 0.95

If you want to read more about DC removal, here is a good tutorial and much more detailed description of how it functions: Everyday DSP for Programmers: DC and Impulsive Noise Removal

Here is the filter implemented in a code:

struct fifo_t {
  uint16_t rawIR;
  uint16_t rawRed;
};
dcFilter_t MAX30100::dcRemoval(float x, float prev_w, float alpha)
{
  dcFilter_t filtered;
  filtered.w = x + alpha * prev_w;
  filtered.result = filtered.w - prev_w;

  return filtered;
}

Once, we pass the signal through the DC removal filter, we should get a signal similar to the one in figure 4:

MORF - Coding And Engineering

[ESCHICK]

The prev_w is the filtered value:

Hello,

It's been a while since I've looked at this code, but it looks like a simple low-pass type filter:
Where X = is your current value, and prev_w is the previous filtered value. Essentially DC filter works on principle of filtering against itself.

dcFilter = dcRemoval( (float)YOUR_NEW_DATA_POINT, dcFilter.w, ALPHA );

And your filtered result will be in this structure.
dcFilter.result

Conveniently I have mentioned the original article I used for DC filter:
Everyday DSP for Programmers: DC and Impulsive Noise Removal

So something like this?

float x;
float prev_w;
float alpha;
float filtered.result;
float filtered.w;

x = rawRed;
alpha = 0.89;

filtered.w = x + alpha * prev_w;
filtered.result = filtered.w - prev_w;
prev_w = filtered.result;

[PaulRB]

Can you please tell the forum what you are trying to achieve? Use non-technical English as far as possible.

The articles you quoted don't seem to match your request for a low pass filter, maybe the opposite.

[ESCHICK]

Hello Paul! Sure. I have a signal (AC and DC) and i want to temove the DC value from the signal, so i only get the AC value.

[PaulRB]

Ok, then maybe the articles you quoted are describing a suitable method. But this is not a low-pass filter. A low pass filter would pass only low frequency components of the signal, in other words remove the AC part.

Please post your full attempt at the sketch, even if it is not working or even compiling. Please read the forum guide first so you know how to post code correctly.

[PaulRB]

@ESCHICK

float filtered.result;

You cannot put "." in the name of a variable in C/C++.

A "." is used to separate the name of an object or structure from the name of an attribute or field inside the object/structure.

Objects or structures might be a useful technique to use if you want to use more than one or two low-pass or high-pass filters, but you seem to be a coding beginner so let's keep things simple as possible for now.

[Koepel]

You might be overthinking this.
The math notation is sometimes different from the source code.

Calculating the average with a low-pass filter is like this:

Average_(T) = 0.99 * Avarage_(T-1) + 0.01 * newValue

The newValue influences the average for 1%.
In the source code, the existing average can be taken, then add 1% of the new value and then store it in the same average variable.

average = 0.99 * average + 0.01 * newValue;

That average can be seen as the DC-offset, so the result is:

result = newValue - average;

Add a setup() and loop() and the sketch is finished:

float average;
const float strength = 0.99; 

void setup() 
{
}

void loop() 
{
  float newValue = (float) analogRead(A0);
  average = (strength * average) + ((1.0-strength) * newValue);
  float result = newValue - average;

  delay(10);
}

---

To make it visual, the Serial Plotter of the Arduino IDE can be used. I used the Wokwi simulator with this sketch:

// Forum: https://forum.arduino.cc/t/gettin-low-pass-filter-into-arduino-sketch/1105548
// This Wokwi project:
//
// The heart beat sensor is not released yet, it may change.
//
// Not using:
//    w(t)=x(t)+∝∗w(t−1)
//    y(t)=w(t)−w(t−1)
//    y(t): is the output of the filter
//    x(t): current input/value
//    w(t): intermediate value, acts like the history of the DC value
//    α: is the response constant of the filter
//    If α = 1 then everything passes through
//    If α = 0 then nothing passes through
//
// Using a standard low-pass filter instead:
//    average = strength * average + (1 - strength) * newVale

float average;
const float strength = 0.995; 

void setup() 
{
  Serial.begin(115200);

  // Print 1 second where all the signals are zero.
  // This is only to make it look better on the Serial Plotter.
  // Not for the real application.
  for(int i=0; i<100; i++)
  {
    Serial.println( "0.0,0.0,0.0,0.0");
    delay(10);
  }
}

void loop() 
{
  float newValue = (float) analogRead(A0);
  average = (strength * average) + ((1.0-strength) * newValue);
  float result = newValue - average;

  Serial.print( "0.0");        // Green, show where the X-axis is
  Serial.print( ",");
  Serial.print( newValue);     // Orange
  Serial.print( ",");
  Serial.print( average);      // Violet
  Serial.print( ",");
  Serial.print( result);       // Cyan
  Serial.println();

  // A quick and dirty way to sample at regular intervals
  delay(10);
}

Try the sketch in Wokwi:

Remove the DC offset of a heart beat with the average - Wokwi ESP32, STM32,...

Green line: X-axis
Orange line: Heart Beat
Pink line: Average
Cyan line: Result: Heart Beat without DC-offset.

[ESCHICK]

Many thanks, this version works! But maybe we get the other solution also to run? (as far as i know it does not need so long the get the DC) I am at the moment not 100% sure how the first version works.

[camsysca]

For removing the DC, in simple terms, you can pre-seed the average with, for example, the average of the first n readings to 'mostly' short-circuit the slow convergence of a filter as described. It takes some fiddling because it will depend on the nature of the waveform being sampled, but it can work wonders in the right circumstances. If you know in advance, for example, that your primary 'disturbance' waveform is about 7 samples in duration, preseeding with an average of 7 samples taken at the start can 'snap' the waveform close to it's final average nicely. (Obviously, this has more relevance if you're using longer averages, for example averaging 1000 samples it's more useful than if averaging 10.) Sometimes, it's even sufficient to simply pre-load the 'average' with a first reading, especially if the DC component is large.
There's probably a proper term for the approach, but I've long since forgotten it.

[Koepel]

@camsysca

pre-seed the average

Of course, that is normally added for a 'average' filter. Add this at the end of setup() for a flying start:

  // Initialize the average value, using 10 samples with 10ms in between.
  average = 0.0;
  for( int i=0; i<10; i++)
  {
    average += (float) analogRead(A0);
    delay(10);
  }
  average /= 10.0;

[ESCHICK]

Really helpfull, thx! Its not the topic but maybe someone can help me here on "the fly". I need also the AC peak Value, what was a good solution? Maybe something like this?

pulse =  (float) particleSensor.getRed() - average;

red = pulse; 


 if (increasing)
        {
            if (red < redMax)
            {
                redMin = redMax;
                increasing = false;
                peakMax = redMax;
            }
            else
            {
                redMax = red;
            }
        }
        else
        {
            if (red > redMin)
            {
                redMax = redMin;
                increasing = true;
            }
            else
            {
                redMin = red;
            }
        }
 
AC = peakMax;

[gcjr]

@Koepel

average = 0.99 * average + 0.01 * newValue;

somewhat simpler
average += (newValue - average) * K; // K < 1

for N = 1/K, average will approach the value of newValue after 3N samples (similar to RC time constant)

@ESCHICK

I need also the AC peak Value,

if (max < abs(samp))
    max = abs(samp);
if (min > abs(samp))
    min = abs(samp);

(what is the purpose of increasing?

[ESCHICK]

Its a bool variable to get the peaks: Compiler Explorer (godbolt.org)

[gcjr]

ok. hopefully you see there no need for it in capturing the max/min?

[Koepel]

ESCHICK, are you trying to find random code online ? How do you know if that code is any good ?

The link to the code on godbolt has questionable code. If all the values in the table are between 100 and 200, then it does not work, and the 'bool' variable is not needed. There is also no explanation why such weird choice of coding is done. It is bad code, I give three

Good code starts with a good description. If you can describe something from the real world in a logical way, then half of the code is already written.

Getting the AC peak value is not possible because that has not enough information.
Do you want the minimum and maximum sample value from the start ?
Do you want the minimum and maximum value over the last few seconds ? Then you should specify how many seconds.
Do you want the minimum and maximum value from the previous peak ? Then you have to remove the noise and determine when a peak is high or deep enough to be accepted as a peak.

[alto777]

@gcjr

somewhat simpler
average += (newValue - average) * K; // K < 1

This

   average = 0.99 * average + 0.01 * newValue;

reads even for a non mathematician.

"the new average is mostly the old average with a bit of the new reading mixed in"

Without a bit of work, it isn't at all obvious that your version even computes the same thing.

You did save a multiplication, so it may be faster and better. Not simpler.

a7

[ESCHICK]

Hello! The Code was from stackflow so i tought it will okay, but of course, better code are always better! But i don't finde a better one.

After i cleare the DC with the code from here i get an okay wave form, i will searche for a filter to get a cleaner wave form and then i want to get the peak value from the last wave, maybe to get an better AC Peak value its a good idea to save 3 peaks and to divide it.

PS: I need some other values to, but its better to make it step by step. I also need the peak numbers from 10 seconds so get the puls (peaks per minute), but first i need the peak value (like in the image)

[Koepel]

The values from the last peaks are not easy, because the signal can be almost anything.

According to wikipedia the maximum is 480 pulses per minute. That is 8Hz. A low-pass filter of 8Hz would clean up the signal.
Detecting when the signal crosses the zero (x-axis) is a common way, together with the low-pass filter it might be a start.

What is the sample rate ?

Note: There are blood pressure measurement devices that easily fail to properly detect the pulse rate. I guess those have just simple code with a filter as well.

[ESCHICK]

@Koepel

...

Thanks but i would need the AC peak value explicitly. I could try to create a group that starts when the first 0 value is exceeded and ends when it falls below it. Then I determine the peak values in each case, I save the value that is just being read in via the sensor and then compare it with the new one, if the new one is smaller than the previous one, then that is a peak value. These are also all saved and compared with each other again. The highest value is then my peak value. Would that be possible?

The sample rate:

" However, the MAX30102 and MAX30101 have a maximum sampling rate of 3200sps ; when reading data under burst mode, the effective data write/read speed is 9600 per second, which is 9.6kHz, so 400kHz is sufficient for the component."

[Koepel]

I'm not sure yet about the heart rate calculation. As a human being, we see immediately what the curve is and where the peaks are.

About 100Hz sample rate over I2C is more realistic.