From 77e270d356542e341ea4a69d16486735e62298c3 Mon Sep 17 00:00:00 2001 From: Erik Andresen Date: Sun, 27 Nov 2022 19:05:40 +0100 Subject: [PATCH] hrmmar: Add fft elimination algorithm to remove MA The fft with size=256 is calculated over a window of 8 seconds with step size 2 seconds. Sample rate of acceleration (12.5Hz) and PPG sensor (25Hz) is unchanged, so the algorithm runs with 12.5Hz. Generally the maximum peak in fft spectrum in an interval of -5..+10BPM of last one is used. The low fft resolution limits the accuracy to ~3BPM, but mean error is ~5-10BPM. If firmware reports confidence >= 90% its value is used. --- apps/hrmmar/README.md | 11 +++ apps/hrmmar/app.png | Bin 0 -> 1192 bytes apps/hrmmar/boot.js | 40 ++++++++ apps/hrmmar/fftelim.js | 190 ++++++++++++++++++++++++++++++++++++++ apps/hrmmar/metadata.json | 18 ++++ apps/hrmmar/settings.js | 26 ++++++ 6 files changed, 285 insertions(+) create mode 100644 apps/hrmmar/README.md create mode 100644 apps/hrmmar/app.png create mode 100644 apps/hrmmar/boot.js create mode 100644 apps/hrmmar/fftelim.js create mode 100644 apps/hrmmar/metadata.json create mode 100644 apps/hrmmar/settings.js diff --git a/apps/hrmmar/README.md b/apps/hrmmar/README.md new file mode 100644 index 000000000..ff90d9156 --- /dev/null +++ b/apps/hrmmar/README.md @@ -0,0 +1,11 @@ +# HRM Motion Artifacts removal + +Measurements from the build in PPG-Sensor (Photoplethysmograph) is sensitive to motion and can be corrupted with Motion Artifacts (MA). This module allows to remove these. + +## Settings + +* **MA removal** + +Select the algorithm to Remove Motion artifacts: + - None: (default) No Motion Artifact removal. + - fft elim: (*experimental*) Remove Motion Artifacts by cutting out the frequencies from the HRM frequency spectrum that are noisy in acceleration spectrum. Under motion this can report a heart rate that is closer to the real one but will fail if motion frequency and heart rate overlap. diff --git a/apps/hrmmar/app.png b/apps/hrmmar/app.png new file mode 100644 index 0000000000000000000000000000000000000000..9db19be376fcfbd57a41ed0a04eee292b34c828e GIT binary patch literal 1192 zcmV;Z1XufsP)7N8&wpCzj16QO=3?ye$>y9#7#&f2&F`l*fYj2&35nxJ-PlfTCr)jTJuJo^JBi0_eQgSo zztx>NbI-ZInRD)AV8n69>sIT^kJmm)b~9ScdZ`Vr^6Uo>pp4bk56z6G0yK#WB2lO0 z1Cj1W(agtUDnKOXC!J+xU)=7udpBe=gQ)-!LW?Opf&4IBpTIf{r=e&XS72`#u1_!& z0l|edlg=_o8fZOdYGDSRLQ`@X*;B-3SCQSNLAyV}zy*XBQbyYy9DW@PTtM(_+Gs2P z2LZCDc;vcMfG5;lz^Um&CLM>Op#ebw-bj|{xxbC;Y1aX=r%0o&yH!l4((v;sB1;L8 z&;~*RqAzdapUW86p|g7aZ`__k)4=I2VxOp;=oem0{x<~7uWac8Bn@b%&-R@eN_!3@ zsry$5Eu@)QN+4M@k???k0Pjo|*;8bEvV!YO0nq1vfK<#+E~)PW-Nyu+_7!j_CB~;J z_#(UghBxS2K=j+yW+^Qcp90m#P59c`4c#=&h`+cFUL_@V;u{N4@ zpx;Yx_yNeJnqMH0N0LC(z%te#6iuTwQ{Oj}W%;$g9y|D6mdY_qGp{lPV4 zZ-Md2D$xD{g`9(0#fE*NhQn2&n0MfdWU-&Bl1nKx>K43_EU&-$GXOWge1r8{(H_yd zXN1r!RI6Ia-*+KN`r=f`Iyrbar5iPgT-pVIyIBXV^x3%=B7P?VfYqk5ci+YOuk(G4 zL3kz3Z2OFMs`>U?wCcbR&@ABOT8Q}Fi1FMn-=FV$*D5UF!Isx(pAe8}*MZHha%ScM z4b?(=!%rseo>bpjNR_$T^TbMHXL)0g$*CGVIxe*x+G6=Yh!+@5fC>QkWaP$jU{6NW2-$V26CuRt_4tQ-Cbonay+>f>iUcb19=^0g9nUF9l%aU z0S}6y{8s6twCA9DXs1@S@i23G;Qog$K-2W~ppbPA#ea_!l8WZGrz4+INQk9R!OYM@ zfObpyCn%|B-;7bpAep0|u?SIms20)F`Vv$eX_ z`Vtb~h^Et{{fRz$NkufiwOF)QZ+!JdOeFq@5hF$nkADCxICO5y(}2|g0000 { + bpm_corrected = bpm; + }; + + Bangle.on('HRM', (hrm) => { + if (bpm_corrected > 0) { + // replace bpm data in event + hrm.bpm_orig = hrm.bpm; + hrm.confidence_orig = hrm.confidence; + hrm.bpm = bpm_corrected; + hrm.confidence = 0; + } + }); + + let run = () => { + const settings = Object.assign({ + mAremoval: 0 + }, require("Storage").readJSON("hrmmar.json", true) || {}); + + // select motion artifact removal algorithm + switch(settings.mAremoval) { + case 1: + require("hrmfftelim").run(settings, updateHrm); + break; + } + } + + // override setHRMPower so we can run our code on HRM enable + const oldSetHRMPower = Bangle.setHRMPower; + Bangle.setHRMPower = function(on, id) { + if (on && run !== undefined) { + run(); + run = undefined; // Make sure we run only once + } + return oldSetHRMPower(on, id); + }; +} diff --git a/apps/hrmmar/fftelim.js b/apps/hrmmar/fftelim.js new file mode 100644 index 000000000..98b7f33ad --- /dev/null +++ b/apps/hrmmar/fftelim.js @@ -0,0 +1,190 @@ +exports.run = (settings, updateHrm) => { + const SAMPLE_RATE = 12.5; + const NUM_POINTS = 256; // fft size + const ACC_PEAKS = 2; // remove this number of ACC peaks + + // ringbuffers + const hrmvalues = new Int16Array(8*SAMPLE_RATE); + const accvalues = new Int16Array(8*SAMPLE_RATE); + // fft buffers + const hrmfftbuf = new Int16Array(NUM_POINTS); + const accfftbuf = new Int16Array(NUM_POINTS); + let BPM_est_1 = 0; + let BPM_est_2 = 0; + + let hrmdata; + let idx=0, wraps=0; + + // init settings + Bangle.setOptions({hrmPollInterval: 40, powerSave: false}); // hrm=25Hz + Bangle.setPollInterval(80); // 12.5Hz + + calcfft = (values, idx, normalize, fftbuf) => { + fftbuf.fill(0); + let i_out=0; + let avg = 0; + if (normalize) { + const sum = values.reduce((a, b) => a + b, 0); + avg = sum/values.length; + } + // sort ringbuffer to fft buffer + for(let i_in=idx; i_in { + let maxVal = -Number.MAX_VALUE; + let maxIdx = 0; + + values.forEach((value,i) => { + if (value > maxVal) { + maxVal = value; + maxIdx = i; + } + }); + return {idx: maxIdx, val: maxVal}; + }; + + getSign = (value) => { + return value < 0 ? -1 : 1; + }; + + // idx in fft buffer to frequency + getFftFreq = (idx, rate, size) => { + return idx*rate/(size-1); + }; + + // frequency to idx in fft buffer + getFftIdx = (freq, rate, size) => { + return Math.round(freq*(size-1)/rate); + }; + + calc2ndDeriative = (values) => { + const result = new Int16Array(values.length-2); + for(let i=1; i { + // fft + const ppg_fft = calcfft(hrmvalues, idx, true, hrmfftbuf).subarray(minFreqIdx, maxFreqIdx+1); + const acc_fft = calcfft(accvalues, idx, false, accfftbuf).subarray(minFreqIdx, maxFreqIdx+1); + + // remove spectrum that have peaks in acc fft from ppg fft + const accGlobalMax = getMax(acc_fft); + const acc2nddiff = calc2ndDeriative(acc_fft); // calculate second derivative + for(let iClean=0; iClean < ACC_PEAKS; iClean++) { + // get max peak in ACC + const accMax = getMax(acc_fft); + + if (accMax.val >= 10 && accMax.val/accGlobalMax.val > 0.75) { + // set all values in PPG FFT to zero until second derivative of ACC has zero crossing + for (let k = accMax.idx-1; k>=0; k--) { + ppg_fft[k] = 0; + acc_fft[k] = -Math.abs(acc_fft[k]); // max(acc_fft) should no longer find this + if (k-2 > 0 && getSign(acc2nddiff[k-1-2]) != getSign(acc2nddiff[k-2]) && Math.abs(acc_fft[k]) < accMax.val*0.75) { + break; + } + } + // set all values in PPG FFT to zero until second derivative of ACC has zero crossing + for (let k = accMax.idx; k < acc_fft.length-1; k++) { + ppg_fft[k] = 0; + acc_fft[k] = -Math.abs(acc_fft[k]); // max(acc_fft) should no longer find this + if (k-2 >= 0 && getSign(acc2nddiff[k+1-2]) != getSign(acc2nddiff[k-2]) && Math.abs(acc_fft[k]) < accMax.val*0.75) { + break; + } + } + } + } + + // bpm result is maximum peak in PPG fft + const hrRangeMax = getMax(ppg_fft.subarray(rangeIdx[0], rangeIdx[1])); + const hrTotalMax = getMax(ppg_fft); + const maxDiff = hrTotalMax.val/hrRangeMax.val; + let idxMaxPPG = hrRangeMax.idx+rangeIdx[0]; // offset range limit + + if ((maxDiff > 3 && idxMaxPPG != hrTotalMax.idx) || hrRangeMax.val === 0) { // prevent tracking from loosing the real heart rate by checking the full spectrum + if (hrTotalMax.idx > idxMaxPPG) { + idxMaxPPG = idxMaxPPG+Math.ceil(6/freqStep); // step 6 BPM up into the direction of max peak + } else { + idxMaxPPG = idxMaxPPG-Math.ceil(2/freqStep); // step 2 BPM down into the direction of max peak + } + } + + idxMaxPPG = idxMaxPPG + minFreqIdx; + const BPM_est_0 = getFftFreq(idxMaxPPG, SAMPLE_RATE, NUM_POINTS)*60; + + // smooth with moving average + let BPM_est_res; + if (BPM_est_2 > 0) { + BPM_est_res = 0.9*BPM_est_0 + 0.05*BPM_est_1 + 0.05*BPM_est_2; + } else { + BPM_est_res = BPM_est_0; + } + + return BPM_est_res.toFixed(1); + }; + + Bangle.on('HRM-raw', (hrm) => { + hrmdata = hrm; + }); + + Bangle.on('accel', (acc) => { + if (hrmdata !== undefined) { + hrmvalues[idx] = hrmdata.filt; + accvalues[idx] = acc.x*1000 + acc.y*1000 + acc.z*1000; + idx++; + if (idx >= 8*SAMPLE_RATE) { + idx = 0; + wraps++; + } + + if (idx % (SAMPLE_RATE*2) == 0) { // every two seconds + if (wraps === 0) { // use rate of firmware until hrmvalues buffer is filled + updateHrm(undefined); + BPM_est_2 = BPM_est_1; + BPM_est_1 = hrmdata.bpm; + } else { + let bpm_result; + if (hrmdata.confidence >= 90) { // display firmware value if good + bpm_result = hrmdata.bpm; + updateHrm(undefined); + } else { + bpm_result = calculate(idx); + bpm_corrected = bpm_result; + updateHrm(bpm_result); + } + BPM_est_2 = BPM_est_1; + BPM_est_1 = bpm_result; + + // set search range of next BPM + const est_res_idx = getFftIdx(bpm_result/60, SAMPLE_RATE, NUM_POINTS)-minFreqIdx; + rangeIdx = [est_res_idx-maxBpmDiffIdxDown, est_res_idx+maxBpmDiffIdxUp]; + if (rangeIdx[0] < 0) { + rangeIdx[0] = 0; + } + if (rangeIdx[1] > maxFreqIdx-minFreqIdx) { + rangeIdx[1] = maxFreqIdx-minFreqIdx; + } + } + } + } + }); +}; diff --git a/apps/hrmmar/metadata.json b/apps/hrmmar/metadata.json new file mode 100644 index 000000000..232ff64a7 --- /dev/null +++ b/apps/hrmmar/metadata.json @@ -0,0 +1,18 @@ +{ + "id": "hrmmar", + "name": "HRM Motion Artifacts removal", + "shortName":"HRM MA removal", + "icon": "app.png", + "version":"0.01", + "description": "Removes Motion Artifacts in Bangle.js's heart rate sensor data.", + "type": "bootloader", + "tags": "health", + "supports": ["BANGLEJS","BANGLEJS2"], + "readme": "README.md", + "storage": [ + {"name":"hrmmar.boot.js","url":"boot.js"}, + {"name":"hrmfftelim","url":"fftelim.js"}, + {"name":"hrmmar.settings.js","url":"settings.js"} + ], + "data": [{"name":"hrmmar.json"}] +} diff --git a/apps/hrmmar/settings.js b/apps/hrmmar/settings.js new file mode 100644 index 000000000..3c6e62c91 --- /dev/null +++ b/apps/hrmmar/settings.js @@ -0,0 +1,26 @@ +(function(back) { + var FILE = "hrmmar.json"; + // Load settings + var settings = Object.assign({ + mAremoval: 0, + }, require('Storage').readJSON(FILE, true) || {}); + + function writeSettings() { + require('Storage').writeJSON(FILE, settings); + } + + // Show the menu + E.showMenu({ + "" : { "title" : "HRM MA removal" }, + "< Back" : () => back(), + 'MA removal': { + value: settings.mAremoval, + min: 0, max: 1, + format: v => ["None", "fft elim."][v], + onchange: v => { + settings.mAremoval = v; + writeSettings(); + } + }, + }); +})