Heart Rate Data Analysis
I needed a way to collect real heart rate sample data for more in-depth analysis. I wanted to not only validate the Butterworth filtering algorithm, but also understand the sensitivity of its parameters and its robustness against other exogenous noise. Furthermore, having a good dataset will enable me to examine other possible algorithms in detecting the pulse accurately. The CC2540 implementation has a bug. Hopefully, I can use the dataset to identify the bug and modify the algorithm to make it work under my new implementation.
I’ve written some small python scripts to extract the raw data from both the Arduino and my implementation on CC2540. Fortunately, the SmartRF05EB board comes with access to UART so I was able to spit out the same data out to PC. From my linux box, I listen to the serial port over USB and just write the data out to raw text. Here’s the script that dumps data from the serial port.
#!/usr/bin/env python
import sys
import time
sys.path.insert(0,'/nas/dev/c/uC/8051/projects/cc2540/scripts/pyserial-2.6')
import serial
def readlines(ser):
starttime = time.time()
while True:
text = ser.readline();
difftime = time.time() - starttime
sys.stdout.write("%f %s" % (difftime,text))
ser = serial.Serial(port='/dev/ttyUSB0', baudrate=115200, timeout=None)
readlines(ser)
ser.flush()
ser.close()
Once the raw text data is collected, another python script parses this data and writes it into tabular csv data. The csv data is much more ammenable for experimentation, and importing it is a snap on any statistical or computing environment of your choice. (for me it’s R) This is the script that reads from the serial port.
#!/usr/bin/env python
import sys
import re
import time
import os
path = sys.argv[1]
f = open(path, "r");
gotS = False
gotB = False
gotQ = False
data = []
row = [ None, None, None, None]
ts = None
print ",".join(['ts','signal','filtered','hrv','bpm'])
while f:
line = f.readline()
tok = line.split()
n = len(line)
if n == 0:
break
if len(tok) < 2:
continue
s = tok[1]
# when we detect the first S, the rest is part of the same row
if s[0] == 'S' and gotS == True:
gotS = False
print ",".join(str(i) for i in [ts]+ row)
# data.append(row[:]) # copies
if s[0] == 'S' and gotS == False:
ts = float(tok[0])
row[0] = int(s[1:])
gotS = True
if s[0] == 'F' and gotS == True:
sign = -1 if s[1] == 'N' else 1
row[1] = sign*int(s[2:])
if s[0] == 'B' and gotS == True:
row[2] = int(s[1:])
if s[0] == 'Q' and gotS == True:
row[3] = int(s[1:])
The pulse is detected whenever the filtered data line crosses zero. The algorithm refines the raw data by filtering it so that pulse detection is easier.
The data from CC2540 gives almost the same result.
Another difference to note is the scale of the raw data. The ADC values (signal) from Arduino hovers around ~900, whereas the CC2540 values hover around ~300. The differences are due to the reference voltage, and the operating voltage among other things. Arduino uses 5V by default for the reference voltage, whereas CC2540 uses 1.25V by default. Arduino runs at 5V, CC2540 runs at 3V.
Read MoreArticle on DIY galvanic response, skin conductance
DIY galvanic response, skin conductance
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 MoreKickstarter is loosing it’s appeal
Kickstarter is loosing it’s appleal,
but according to comments, so is gizmodo
Read MoreY-Combinator – Frighteningly Ambitious Startup Ideas
Very interesting Paul Graham article on Frighteningly Ambitious Startup Ideas, please note idea #7.
Read MoreA new X-prize competition
A new X-prize <a title=”x-prize” href=”http://www.engadget.com/2012/05/25/nokia-and-x-prize-put-medical-sensors-on-the-spot-for-next-chall/”>competition </a>have been announced with collaboration with Nokia, prize is $2.25 million.
Read MoreMetawatch with BLE support
A new version of <a title=”metawatch” href=”http://www.metawatch.org/forums/thread/455/using-the-metawatch-as-bluetooth-host-device-ble”>metawatch </a>came out that supports BLE now.
Read Morescanadu.com developing real life medical tricoder
Remeber Star Trek <a title=”tricoder” href=”http://en.wikipedia.org/wiki/Tricorder”>tricoder</a>? Scanadu developing a real life <a title=”scanadu tricoder” href=”http://techcrunch.com/2011/11/08/scanadu-raises-2m-check-your-body-as-often-as-your-email/”>one</a>.
Read MoreTI CC2541 launch
TI’s second Bluetooth low energy single-mode SoC IC, <a title=”CC2541″ href=”http://e2e.ti.com/blogs_/b/bluetooth_low_energy_blog/archive/2012/02/15/cc2541-now-released-including-support-for-bluetooth-low-energy-and-proprietary.aspx”>CC2541</a>, have been launched.
Read MoreVenture lab final presentation
Monbaby is a baby health monitor, wireless wearable monitor specifically designed for babies and parents. It keeps track of the baby’s vital signals and alerts you upon critical events. Our product aim is to reduce parents anxiety about their babies, and to help parents to monitor their babies’ health with highest reliance but least amount of effort.
We are building a wireless monitor based on Bluetooth LE wearable as an ankle bracelet. It measures baby orientation and other data, then transmits measured signal to a Bluetooth capable receiver, such as IPhone, or Android smartphone or a Bluetooth Smart baby monitor. Unlike existing product on the market that send alert when baby stops moving, at which point it maybe too late to react, Monbaby warns parents if baby turned around and sleeping face down, giving parents plenty of time to go and turn baby face up.
We are former colleagues and NJ neighbors, and we met with an idea to pursue projects outside of work that combine wireless hardware, software and firmware programming, as we share love and passion about those technologies. By day, we are engineers and quantitative program traders, who are experienced in diversified software and hardware development, quantitative analysis and problem solving with combined 20 years of experience in Wall Street finance. We love to come up with new ideas, and convert these ideas into real life applications.
The idea for Monbaby product came as an improvement on existing video/audio baby monitors. We remember waking in the middle of the night just to check on the baby. We have submitted several ideas to bluetooth competition but Monbaby was noticed the most.
We got a pretty good feedback on this idea so far, as we were finalists in last year bluetooth competition, http://remnart.com/2012/02/bluetooth-iwc-competition/ and we were noticed and invited to CES 2012, http://remnart.com/2012/01/monbabyinces_2012/ to Freescale booth.
Now we want to take this product further. We were excited to participate in this class to get feedback from unique blend of entrepreneurial community and get a taste of silicon valley spirit.
Now we want to see Monbaby to fruition and to take it from a prototype into a product.
Developing for CC2540: Part II
A Bit About GAP and GATT
The entire BLE stack uses HAL to communicate with the hardware and the associated peripherals such as LED’s, LCD, and push buttons. The OSAL layer is the glue that lets different tasks and layers of software to interact with each other. The innards of the Bluetooth LE protocol is largely divided into two software layers:
- GAP (Generic Access Profile)
- GATT (Generic Attribute Profile)
The combination of those two layers implement the connectivity and the inter-device communication standards of the Bluetooth LE protocol. These two layers define what it is to be a bluetooth device, and by doing so, it ensures interoperability across various bluetooth-enabled devices that are built using the same specifications, this to carry a proper security device or any other app for monitoring ssolutions. Want to know the best Windows 10 Keylogger? Check out this review about the xnspy Keylogger for 2020
The CC2540 masks most of the underlying details of the link layer. The specifics of the wireless characteristics of bluetooth such as frequency hopping, secure pairing, and subsequent encryption process, are all implemented within a closed-source compiled library. This library is provided as a linkable binary only for IAR Workbench. And for this reason, the development for CC2540 is only possible with IAR Workbench. With open source code, it would have been possible port the software stack for other compilers such as SDCC.
Mind the GAP
The GAP layer defines specific set of roles each BLE device can assume. In the wireless device-to-device connection stage, each node can take one of the following roles (or profiles):
- Broadcaster
- Peripheral
- Observer
- Central
Without going too much into detail, it’s helpful to remember that the Observer/Central are the ones that will look for devices to connect to – they act as the master nodes in the pairing process. On the other hand, the Broadcaster/Peripheral node is the one that will sit there and advertise its existence, acting as the slave. The only difference to note among these profiles is that both Broadcaster and Observer are non-connectable.
What’s GATT?
The GATT layer sets up standards for data and inter-node communication flow between the nodes. It’s this layer that enables various bluetooth devices to share information with one another in a common “language”. Without this common language, devices manufactured by different vendors will speak their own langauges and will not be interoperable. The GATT specifications give definition and structure to both the data and functions each node in the connection can perform. Within this paradigm, any Bluetooth LE device can discover the identity and the functionality of any other Bluetooth LE devices.
BTTool
Here are the commands necessary to start receiving data from the heartrate monitor GATT profile.
* discover characteristic by UUID
Characteristic UUID = 37:2A (from HEARTRATE_MEAS_UUID define below)
* heart rate sensor comes back with “10 11 00 37 2A”
– 10 : means notify only
– 11 00 (0x0011) : is the handle for characteristic value
– 37 2A (0x2A37) : is the UUID
* in order to enable notifications, client need to write 0x0001 to
server characteristic config descriptor. The handle for this
config descriptor immediately follows the characteristic value’s
handle. So in this case, it is: 0x0012
– chracteristic value handle : 0x0012
– enter into write text box : 01:00
– once the write succeeds, the heart rate is transmitted
Mission for longer life
Life expectancy in US have reached 80 years and some experts claim that babies born in 2011 are 8 times more likely to reach a lifespan of 100 years. The world is getting older and more healthy and we at, at Remnart technologies, want to make sure that every newborn baby would have a chance for a long and healthy life.
We’ve been putting our time and effort into Monbaby product, which was envisioned and designed to help parents to monitor their newborn babies better than existing audio/video monitors.
It works like this: an ankle bracelet or a wristband with a clasp containing our sensor is worn around baby’s arm or an ankle, similar to existing bracelets worn by newborns to prevent baby theft. Clasp contains a sensor that monitors baby health and measures vital health signals. The signal is processed by a chip in a bracelet and send over Bluetooth wireless signal to a wireless reader in the vicinity. Wireless reader then records or transmits the data over WiFi and may trigger an alert if anything goes wrong.
This website and this blog covers our thoughts and our path to developing this product. We are trying to stay as open as possible in this endeavour and we are proponent of open source hardware movement. We will try to establish an acceptable solution and invite people to contribute to it.
As for the market potential of our product, recession or not, the field of medical health monitoring, and the medical electronic devices in general, is rapidly growing. According to DataBeans 2010 Medical Semiconductors report the market of medical electronic devices, is expected to reach US$191.1 billion by 2015 with 8% annual growth rate.
The reasons behind the growth are increased average life span, aging population and growing demand from emerging markets.
We are hoping that our product will be both useful, popular and profitable.
Our mission statement is a quest for longer life.
Read MoreBattery life and low power wireless tech
A very good reference about the battery life among different low power wireless technologies applications. Another reason why we use Bluetooth V4.0 (Bluetooth Low Energy) for our Monbaby product.
http://www.embedded.com/ContentEETimes/Documents/Schweber/C0924/C0924post.pdf
Read MoreComparison between ZigBee and Bluetooth Low Energy
Both ZigBee and Bluetooth Low Energy are designed for ultra-low power wireless communicating devices. Below are some of their major differences:
| ZigBee | Bluetooth Low Energy | |
| Range | 10-100m | 50m |
| Networking topology |
|
|
| Modulation |
|
|
| Bit rate | ~100kbps | ~200kbps |
| Setup time | 20ms | 3ms |
| Frequency |
|
|
| Hardware support | In some consumer electronics and remote controls | Most of the handheld devices support Bluetooth. BLE dual-mode is going to be included in all new devices. |
| Power consumption* | 185.9uW/bit | 0.153uW/bit |
* Reference: http://www.electronicscomponentsworld.com/articleView~idArticle~71867_2737432862011.html
Read MoreRadiation concern
One of the concern about the monbaby device is radiation safeness, especially it is going to be used on the little babies who are more sensitive to environmental stimulus.
Radiation of radio frequency (RF) devices are measured by Special absorption rate (SAR). It is a measure of the rate at which energy is absorbed by the body when exposed to a radio frequency electromagnetic field. It is defined as the power absorbed per mass of tissue and has units of watts per kilogram (W/kg). A general rule of thumb is, higher emission power would increase the SAR measurement.
Here are some sample facts and measurement of SAR for different RF devices:
- Cellphones in US are limited to 1.6W/kg by regulation
- Apple iPhone 4 : 1.17W/kg
- Baby monitors range between 0.01-0.08W/kg
- Bluetooth devices about 0.001W/kg
- There is no public measurement data found for BLE yet, but given BLE has even lower power emission than Bluetooth, the SAR number should be very insignificant to be compared.
Bluetooth technology has such an insignificant radiation that in fact, this is one of the reason why a lot of bluetooth earpieces came out in the market after cellphones become our daily routine, to reduce the radiation we are exposing to when making calls.
Bluetooth and Bluetooth low energy technology are both using 2.4-2.5GHz band, which is the same as the common FM radio.
Read More
The BLE Prototyping Hardware: TI CC2540 – DEV KIT
Texas Instrument manufactures an integrated BLE chip called CC2540. This chip contains not only the bluetooth LE radio, but also programmable 8051 core. They provide complete development tools, and they come in two forms:
The larger of the devkit CC2540DK contains general-purpose development / evaluation board called SmartRF05EB. These boards are also used for their line of wireless modules created by ChipCon (subsidiary of TI). As such, it’s compatible with other wireless evaluation modules CCxxxx.
The CC2540DK kit comes with two CC2540EM evaluation modules and two SmartRF05EB boards. The SmartRF05EB contains lots of peripherals such as LCD, Joystick, LEDs, Potentiometer, switches, etc.. In addition, it has necessary components to flash and debug the target chip. The target chip is any compatible CC device as previously mentioned, which is surface mounted onto the evaluation module PCBs.
A noteworthy aspects about this kit
- Only BLE hardware implementation that has coupled host microcontroller in a single chip (SoC: System On Chip)
- Significant sofware stack for interfacing with the hardware
- Development project workspaces available for testing on IAR Workbench
- Bluetooth v4.0 compliant profile implementations and examples
In order to test out different examples provided with TI’s ble stack download, several other tools need to be setup:
- SmartRF Flasher – flashes the chip with designated firmware (.hex file)
- SmartRF Protocol Packet Sniffer – sniffs the packets and provides helpful display
- SmartRF Studio – application that can control CCxxxx devices
- BTTool – comes with ble download and can be used to send BLE specific commands
- IAR Workbench for 8051– there’s 30 day trial version you can use to test
As for this writing, this is the only dev kit that provides both the hardware and the necessary software stack to test out BLE functionality in a single package. The only downside is that the software stack is written only for IAR Workbench, which requires license to create compiled hex beyond 4Kb (free version limit). We would love to see an opensource version of the software stack written for compilers such as SDCC.
The Chip
There are 2 variants of CC2540 : CC2540F128, CC2540F256. The only different between the two is the flash space as the name indicates. (128kB vs 256kB). This chip has all major peripherals and hardware capability of modern microcontrollers, as you can see from the functional diagram below. Some highlights:
- USB controller
- 2 hardware UART
- Hardware encryption engine
- bunch of timers
- ADC’s go upto 14bit resolution
Next, the CC2540′s pin labels in the current package:
The evaluation module that houses the actual CC2540 chip has all of the IO pins brought out to accessible debug headers. Most of the IO ports are routed to the first SMD socket (evaluation module socket), and is accessible from one of the debug headers located on P18 on SmartRF05EB.
Connectors for Peripherals
The pins that are brought out via evaluation module connectors are mapped to various peripherals on SmartRF05EB board. The mapping is as follows. On my rev 1.8 board, these pin outs correspond to board connectors P5 and P6, respectively. More details can be found on SmartRF05EB board’s user guide.
The probe connectors give access to these pins, but the mapping from the EM -> Probe pinouts are not straight forward. The probe pin mapping is as follows:
To access a specific pin on from CC2540, we need to do a little mapping of the ports specified from the tables shown above. For example, the IO pin P0.1 from CC2540 corresponds to EM_BUTTON1/EM_LED4_SOC at EM connector pin 5. The same pin corresponds to the pin 7 on the prob connector. I wish the mapping of the pins from the debug headers to the actual chip was more straight forward. I had to trudge through 3 different technical documentations to find this.
I’ll be connecting the pulsesensor to pin P0.7, which corresponds to EM_POT_R (potentiometer) on pin 17 @ P5 on the EM connector and pin 17 @ P18 on the probe connector. In order to disable the onboard potentiometer, I needed to take the jumper off on connector P10 on SmartRF05EB, at position 33-34. The connector P10 is where you can sever or link each individual peripheral that connects to the evaluation board housing CC2540.
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