Latest Monbaby demo – full stack: device to iphone to cloud
Monbaby is a baby health monitor with a form factor that it could be attached anywhere on baby clothing.
It measures movements, detects motions, transmits signal over Bluetooth 4.0 (Low Energy) to an IPhone app, which analyzes and displays the data and re-transmits signals further into the cloud, where more sophisticated analysis could be done on a backend. Here is the video of the full technology pipeline:
http://youtu.be/J1diQt1irEQ
Read MoreFourth Iteration of Monbaby Prototype
Prototypes are meant to be molded and refined until it embodies viable set of features, given the present constraints. Prototyping a device is one thing, but refining it to fit the numerous criteria of technical feasibility, manufacturing, and final product design can be challenging.
PCB image with CR2032 battery holder to show relative size of the device
That said, we are making very good progress. And these latest images show the device at its 4th iteration in its prototype stage. We are getting close, and all we need to do is to pull all the pieces together to make Monbaby into truly integrated device!
3D image of Monbaby prototype PCB
Monbaby movement detection via IPhone app
Demo showing Monbaby movement and orientation detection on an IPhone app
Read MoreMonbaby Android BLE demo
Here is a link to a quite crude, but functional demo of a monbaby sensor transmitting accelerometer measurements wirelessly over BLE to a smartphone that routes it to the web/cloud.
Read MoreHeart Rate Measurements via Pulsesensor
The Sensor
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.
Read More
Recent Comments