Heart Rate Measurements via Pulsesensor

Posted on Dec 23, 2011 | 6 comments

The Sensor

pulsesensor Pulsesensor is an open-source heart rate monitoring sensor built to work with the Arduino. The underlying mechanism uses periodic sampling of ADC value from the light intensity sensor, and signal filtering to detect the peak-to-peak length in the data. The current Arduino code requires relatively high rate of sample polling, and thus prone to heavier power consumption compared to low-powered designs.

While the method used to detect the pulse is not new, the pulsesensor is the first open source product we’ve seen for heart rate sensors. As such, we felt the product was a good fit for our first prototype design for Monbaby. However, the platform we will be using isn’t Arduino. Hence, we will need to examine the algorithm used, and understand its pros and cons. There should be other signal filtering algorithms that are appropriate in finding out peak-to-peak length in the noisy raw data from the light intensity sensor. Once the signal processing algorithm is fully understood, we can then reimplement the logic within the CC2540 software stack.

The Algorithm

The latest release of pulsesensor code is posted in its entirety below. This code is meant to compile under Arduino IDE. Hence, it’s missing headers and normal format of typical AVR firmware C source. It uses the timer interrupt every 1ms to read the analog value from the sensor. It averages out the readings at every 300 run, and de-means the raw signal by subtracting it by the calculated mean.

This normalized (de-meaned) data is then refined further. A signal processing method called “Butterworth filter” is used to reduce the noise in the normalized raw data. I have yet to read more about the algorithm, but it seems akin to ARIMA time series model, where the predicted value (Y_hat) is a function of past values and the forecast error term (residual) between the real and prediction values. The way Butterworth filter is used, however, is through “Band-Pass” filtering process. It defines low-pass and high-pass bands that eliminate frequencies that are outside of those bounds. We will be reading more about these signal processing techniques in the coming days.

/*
>> Pulse Sensor Digital Filter <<
This code is the library prototype for Pulse Sensor www.pulsesensor.com. 
    >>> Pulse Sensor purple wire goes to Analog Pin 0 <<<
Pulse Sensor sample aquisition and processing happens in the background via Timer 1 interrupt. 1mS sample rate.
The following variables are automatically updated:
Pulse :     boolean that is true when a heartbeat is sensed then false in time with pin13 LED going out.
Signal :    int that holds the analog signal data straight from the sensor. updated every 1mS.
HRV  :      int that holds the time between the last two beats. 1mS resolution.
B  :        boolean that is made true whenever HRV is updated. User must reset.
BPM  :      int that holds the heart rate value. derived from averaging HRV every 10 pulses.
QS  :       boolean that is made true whenever BPM is updated. User must reset.
Scale  :    int that works abit like gain. use to change the amplitude of the digital filter output. useful range 12<>20 : high<>low default = 12
FSignal  :  int that holds the output of the digital filter/amplifier. updated every 1mS.

See the README for detailed information and known issues.
Joel Murphy  December 2011  Happy New Year! 
*/

long Hxv[4]; // these arrays are used in the digital filter
long Hyv[4]; // H for highpass, L for lowpass
long Lxv[4];
long Lyv[4];

unsigned long readings; // used to help normalize the signal
unsigned long peakTime; // used to time the start of the heart pulse
unsigned long lastPeakTime = 0;// used to find the time between beats
volatile int Peak;     // used to locate the highest point in positive phase of heart beat waveform
int rate;              // used to help determine pulse rate
volatile int BPM;      // used to hold the pulse rate
int offset = 0;        // used to normalize the raw data
int sampleCounter;     // used to determine pulse timing
int beatCounter = 1;   // used to keep track of pulses
volatile int Signal;   // holds the incoming raw data
int NSignal;           // holds the normalized signal 
volatile int FSignal;  // holds result of the bandpass filter
volatile int HRV;      // holds the time between beats
volatile int Scale = 15;  // used to scale the result of the digital filter. range 12<>20 : high<>low amplification
volatile int Fade = 0;

boolean first = true; // reminds us to seed the filter on the first go
volatile boolean Pulse = false;  // becomes true when there is a heart pulse
volatile boolean B = false;     // becomes true when there is a heart pulse
volatile boolean QS = false;      // becomes true when pulse rate is determined. every 20 pulses

int pulsePin = 0;  // pulse sensor purple wire connected to analog pin 0

void setup(){
pinMode(13,OUTPUT);    // pin 13 will blink to your heartbeat!
Serial.begin(115200); // we agree to talk fast!
// this next bit will wind up in the library. it initializes Timer1 to throw an interrupt every 1mS.
TCCR1A = 0x00; // DISABLE OUTPUTS AND BREAK PWM ON DIGITAL PINS 9 & 10
TCCR1B = 0x11; // GO INTO 'PHASE AND FREQUENCY CORRECT' MODE, NO PRESCALER
TCCR1C = 0x00; // DON'T FORCE COMPARE
TIMSK1 = 0x01; // ENABLE OVERFLOW INTERRUPT (TOIE1)
ICR1 = 8000;   // TRIGGER TIMER INTERRUPT EVERY 1mS  
sei();         // MAKE SURE GLOBAL INTERRUPTS ARE ENABLED

}

void loop(){
  Serial.print("S");          // S tells processing that the following string is sensor data
  Serial.println(Signal);    //  filterdSignal holds the latest filtered Pulse Sensor signal
  if (B == true){             //  B is true when arduino finds the heart beat
    Serial.print("B");        // 'B' tells Processing the following string is HRV data (time between beats in mS)
    Serial.println(HRV);      //  HRV holds the time between this pulse and the last pulse in mS
    B = false;                // reseting the QS for next time
  }
  if (QS == true){            //  QS is true when arduino derives the heart rate by averaging HRV over 20 beats
    Serial.print("Q");        //  'QS' tells Processing that the following string is heart rate data
    Serial.println(BPM);      //  BPM holds the heart rate in beats per minute
    QS = false;               //  reset the B for next time
  }
  Fade -= 15;
  Fade = constrain(Fade,0,255);
  analogWrite(11,Fade);
  
delay(20);                    //  take a break

}

// THIS IS THE TIMER 1 INTERRUPT SERVICE ROUTINE. IT WILL BE PUT INTO THE LIBRARY
ISR(TIMER1_OVF_vect){ // triggered every time Timer 1 overflows
// Timer 1 makes sure that we take a reading every milisecond
Signal = analogRead(pulsePin);

// First normailize the waveform around 0
readings += Signal; // take a running total
sampleCounter++;     // we do this every milisecond. this timer is used as a clock
if ((sampleCounter %300) == 0){   // adjust as needed
  offset = readings / 300;        // average the running total
  readings = 0;                   // reset running total
}
NSignal = Signal - offset;        // normalizing here

// IF IT'S THE FIRST TIME THROUGH THE SKETCH, SEED THE FILTER WITH CURRENT DATA
if(first = true){
  for (int i=0; i<4; i++){
    Lxv[i] = Lyv[i] = NSignal <<10;  // seed the lowpass filter
    Hxv[i] = Hyv[i] = NSignal <<10;  // seed the highpass filter
  }
first = false;      // only seed once please
}
// THIS IS THE BANDPAS FILTER. GENERATED AT www-users.cs.york.ac.uk/~fisher/mkfilter/trad.html
//  BUTTERWORTH LOWPASS ORDER = 3; SAMPLERATE = 1mS; CORNER = 5Hz
    Lxv[0] = Lxv[1]; Lxv[1] = Lxv[2]; Lxv[2] = Lxv[3];
    Lxv[3] = NSignal<<10;    // insert the normalized data into the lowpass filter
    Lyv[0] = Lyv[1]; Lyv[1] = Lyv[2]; Lyv[2] = Lyv[3];
    Lyv[3] = (Lxv[0] + Lxv[3]) + 3 * (Lxv[1] + Lxv[2])
          + (3846 * Lyv[0]) + (-11781 * Lyv[1]) + (12031 * Lyv[2]);
//  Butterworth; Highpass; Order = 3; Sample Rate = 1mS; Corner = .8Hz
    Hxv[0] = Hxv[1]; Hxv[1] = Hxv[2]; Hxv[2] = Hxv[3];
    Hxv[3] = Lyv[3] / 4116; // insert lowpass result into highpass filter
    Hyv[0] = Hyv[1]; Hyv[1] = Hyv[2]; Hyv[2] = Hyv[3];
    Hyv[3] = (Hxv[3]-Hxv[0]) + 3 * (Hxv[1] - Hxv[2])
          + (8110 * Hyv[0]) + (-12206 * Hyv[1]) + (12031 * Hyv[2]);
FSignal = Hyv[3] >> Scale;  // result of highpass shift-scaled
//PLAY AROUND WITH THE SHIFT VALUE TO SCALE THE OUTPUT ~12 <> ~20 = High <> Low Amplification.

if (FSignal >= Peak && Pulse == false){  // heart beat causes ADC readings to surge down in value.  
  Peak = FSignal;                        // finding the moment when the downward pulse starts
  peakTime = sampleCounter;              // recodrd the time to derive HRV. 
}
//  NOW IT'S TIME TO LOOK FOR THE HEART BEAT
if ((sampleCounter %20) == 0){// only look for the beat every 20mS. This clears out alot of high frequency noise.
  if (FSignal < 0 && Pulse == false){  // signal surges down in value every time there is a pulse
     Pulse = true;                     // Pulse will stay true as long as pulse signal < 0
     digitalWrite(13,HIGH);            // pin 13 will stay high as long as pulse signal < 0  
     Fade = 255;                       // set the fade value to highest for fading LED on pin 11 (optional)   
     HRV = peakTime - lastPeakTime;    // measure time between beats
     lastPeakTime = peakTime;          // keep track of time for next pulse
     B = true;                         // set the Quantified Self flag when HRV gets updated. NOT cleared inside this ISR     
     rate += HRV;                      // add to the running total of HRV used to determine heart rate
     beatCounter++;                     // beatCounter times when to calculate bpm by averaging the beat time values
     if (beatCounter == 10){            // derive heart rate every 10 beats. adjust as needed
       rate /= beatCounter;             // averaging time between beats
       BPM = 60000/rate;                // how many beats can fit into a minute?
       beatCounter = 1;                 // reset counter
       rate = 0;                        // reset running total
       QS = true;                       // set Beat flag when BPM gets updated. NOT cleared inside this ISR
     }
  }
  if (FSignal > 0 && Pulse == true){    // when the values are going up, it's the time between beats
    digitalWrite(13,LOW);               // so turn off the pin 13 LED
    Pulse = false;                      // reset these variables so we can do it again!
    Peak = 0;                           // 
  }
}

}// end isr

A Quick Test

There are few packages in R that implements the Butterworth filtering on timeseries data. I’ve used mFilter to test out the theory. First I generated a cosine wave with random noise:


> x = cos(pi*seq(100)/5)+rnorm(100, sd=0.5)
> plot(x, type='l')

And then I run it through the bwfilter function to see how the algorithm smoothes out the noisy data. It should do a reasonable job of getting back our original cosine wave with the noise removed.


> f3 = bwfilter(x, freq=3)
> f4 = bwfilter(x, freq=4)
> f5 = bwfilter(x, freq=5)
> f6 = bwfilter(x, freq=6)
> plot(f3$x, col=1, type='l')
> lines(f3$trend, col=2)
> lines(f4$trend, col=3)
> lines(f5$trend, col=4)
> lines(f6$trend, col=5)
> legend("topleft", legend=c("rawdata","freq=3","freq=4","freq=5","freq=6"), col=1:5, lty=rep(1,5))

The result shows the effect of “smoothing” that Butterworth algorithm carries out on the raw data as you increase the number of sample count per cycle. In our case, the number of samples per sinusoidal wave was 10 ( pi * 100 / 5 = 20 pi = 10 waves ). With 100 total samples, that gives you 10 samples per wave. As you throttle the freq argument from 3 to 6, you can see how the filtering is taking out the random noise I added in.

By default, the bwfilter function uses order 2. The order basically determines how “complex” your model can get. Higher the order, more accurate the filtering will be, but at a computational cost. Below is the printout of bwfilter R function.

function (x, freq = NULL, nfix = NULL, drift = FALSE) 
{
    if (is.null(drift)) 
        drift <- TRUE
    xname = deparse(substitute(x))
    if (is.ts(x)) 
        frx = frequency(x)
    else frx = 1
    if (is.null(freq)) {
        if (frx > 1) 
            freq = trunc(frx * 2.5)
        else freq = 2
    }
    if (is.null(nfix)) 
        nfix = 2
    xo = x
    x = as.matrix(x)
    if (drift) 
        x = undrift(x)
    n = length(x)
    cut.off = 2 * pi/freq
    mu = (1/tan(cut.off/2))^(2 * nfix)
    imat = diag(n)
    Ln = rbind(matrix(0, 1, n), diag(1, n - 1, n))
    Ln = imat - Ln
    if (nfix > 1) {
        for (i in 1:(nfix - 1)) Ln = (imat - Ln) %*% Ln
    }
    Q = t(Ln[3:n, ])
    SIGMA.R = t(Q) %*% Q
    SIGMA.n = abs(SIGMA.R)
    g = t(Q) %*% as.matrix(x)
    b = solve(SIGMA.n + mu * SIGMA.R, g)
    x.cycle = c(mu * Q %*% b)
    x.trend = x - x.cycle
    if (is.ts(xo)) {
        tsp.x = tsp(xo)
        x.cycle = ts(x.cycle, star = tsp.x[1], frequency = tsp.x[3])
        x.trend = ts(x.trend, star = tsp.x[1], frequency = tsp.x[3])
        x = ts(x, star = tsp.x[1], frequency = tsp.x[3])
    }
    A = mu * Q %*% solve(SIGMA.n + mu * SIGMA.R) %*% t(Q)
    if (is.ts(xo)) {
        tsp.x = tsp(xo)
        x.cycle = ts(x.cycle, star = tsp.x[1], frequency = tsp.x[3])
        x.trend = ts(x.trend, star = tsp.x[1], frequency = tsp.x[3])
        x = ts(x, star = tsp.x[1], frequency = tsp.x[3])
    }
    res <- list(cycle = x.cycle, trend = x.trend, fmatrix = A, 
        title = "Butterworth Filter", xname = xname, call = as.call(match.call()), 
        type = "asymmetric", lambda = mu, nfix = nfix, freq = freq, 
        method = "bwfilter", x = x)
    return(structure(res, class = "mFilter"))
}

Hmm.. the guts of the function looks somewhat familiar to me. It contains matrix manipulation that seems to construct the coefficient matrix that models the relationship among variables on a "moving window". The matrix multiplication is carried out between this coefficient matrix and the actual data, which gives the values to the right-hand side of the equation : Ax = b. The problem is then treated as an inverse problem. The function "solve" is used when you have the model and the output, but you do not know which inputs you need to make the output hold true given the relationships defined in your linear model. The coefficients, or the loadings of the variables are calculated this way, and it is a common technique used in numerical computing when the model parameters are fit to the existing dataset.

Further investigation is necessary to understand the system of equations used and what objective function this filter is trying to optimize.

6 Responses to “Heart Rate Measurements via Pulsesensor”

  1. Michael Sison says:

    Good day!

    I have a question about the code. In normalizing the waveform, why is it that there is a constant 300? What is its relation to time? Maybe, a mathematical conversion could help. I’m confused about the value, 300. How do we know what value to input if it is adjustable?

    • 300 constant is meant to calculate an offset every 300 ms. You can replace this in the code with an adjustable constant

  2. shah riza says:

    how to calculate low frequency and high frequency ratio from the value that we get from pulse sensor

  3. Mohd Izzad says:

    Hi, I need help. Why does the BPM value did not appear at the serial monitor? Besides that, how can I implement Kalman filter for the pulse sensor?