/* Custom CPAP/Oximetry Calculations Header
*
* Copyright (c) 2011-2014 Mark Watkins <jedimark@users.sourceforge.net>
*
* This file is subject to the terms and conditions of the GNU General Public
* License. See the file COPYING in the main directory of the Linux
* distribution for more details. */
#include <QMutex>
#include <QFile>
#include <QDataStream>
#include <QTextStream>
#include <cmath>
#include <algorithm>
#include "calcs.h"
#include "profiles.h"
bool SearchEvent(Session * session, ChannelID code, qint64 time, int dur, bool update=true)
{
qint64 t, start;
QHash<ChannelID, QVector<EventList *> >::iterator it;
it = session->eventlist.find(code);
quint32 *tptr;
EventStoreType *dptr;
int cnt;
//qint64 rate;
// bool fixdurations = (session->machine()->loaderName() != STR_MACH_ResMed);
if (!p_profile->cpap->resyncFromUserFlagging()) {
update=false;
}
QHash<ChannelID, QVector<EventList *> >::iterator evend = session->eventlist.end();
if (it != evend) {
int el_size=it.value().size();
for (int i = 0; i < el_size; i++) {
EventList *el = it.value()[i];
// rate=el->rate();
cnt = el->count();
// why would this be necessary???
if (el->type() == EVL_Waveform) {
qDebug() << "Called SearchEvent on a waveform object!";
return false;
} else {
start = el->first();
tptr = el->rawTime();
dptr = el->rawData();
for (int j = 0; j < cnt; j++) {
t = start + *tptr;
// Move the position and set the duration
qint64 end1 = time + 5000L;
qint64 start1 = end1 - quint64(dur+10)*1000L;
qint64 end2 = t + 5000L;
qint64 start2 = end2 - quint64(*dptr+10)*1000L;
bool overlap = (start1 <= end2) && (start2 <= end1);
if (overlap) {
if (update) {
qint32 delta = time-start;
if (delta >= 0) {
*tptr = delta;
*dptr = (EventStoreType)dur;
}
}
return true;
}
tptr++;
dptr++;
}
}
}
}
return false;
}
bool SearchApnea(Session *session, qint64 time, double dur)
{
if (SearchEvent(session, CPAP_Obstructive, time, dur))
return true;
if (SearchEvent(session, CPAP_Apnea, time, dur))
return true;
if (SearchEvent(session, CPAP_ClearAirway, time, dur))
return true;
if (SearchEvent(session, CPAP_Hypopnea, time, dur))
return true;
if (SearchEvent(session, CPAP_UserFlag1, time, dur, false))
return true;
if (SearchEvent(session, CPAP_UserFlag2, time, dur, false))
return true;
return false;
}
// Sort BreathPeak by peak index
bool operator<(const BreathPeak &p1, const BreathPeak &p2)
{
return p1.start < p2.start;
}
//! \brief Filters input to output with a percentile filter with supplied width.
//! \param samples Number of samples
//! \param width number of surrounding samples to consider
//! \param percentile fractional percentage, between 0 and 1
void percentileFilter(EventDataType *input, EventDataType *output, int samples, int width,
EventDataType percentile)
{
if (samples <= 0) {
return;
}
if (percentile > 1) {
percentile = 1;
}
QVector<EventDataType> buf(width);
int s, e;
int z1 = width / 2;
int z2 = z1 + (width % 2);
int nm1 = samples - 1;
//int j;
// Scan through all of input
for (int k = 0; k < samples; k++) {
s = k - z1;
e = k + z2;
// Cap bounds
if (s < 0) { s = 0; }
if (e > nm1) { e = nm1; }
//
int j = 0;
for (int i = s; i < e; i++) {
buf[j++] = input[i];
}
j--;
EventDataType val = j * percentile;
EventDataType fl = floor(val);
// If even percentile, or already max value..
if ((val == fl) || (j >= width - 1)) {
nth_element(buf.begin(), buf.begin() + j, buf.begin() + width - 1);
val = buf[j];
} else {
// Percentile lies between two points, interpolate.
double v1, v2;
nth_element(buf.begin(), buf.begin() + j, buf.begin() + width - 1);
v1 = buf[j];
nth_element(buf.begin(), buf.begin() + j + 1, buf.begin() + width - 1);
v2 = buf[j + 1];
val = v1 + (v2 - v1) * (val - fl);
}
output[k] = val;
}
}
void xpassFilter(EventDataType *input, EventDataType *output, int samples, EventDataType weight)
{
// prime the first value
output[0] = input[0];
for (int i = 1; i < samples; i++) {
output[i] = weight * input[i] + (1.0 - weight) * output[i - 1];
}
//output[samples-1]=input[samples-1];
}
FlowParser::FlowParser()
{
m_session = nullptr;
m_flow = nullptr;
m_filtered = nullptr;
m_gain = 1;
m_samples = 0;
m_startsUpper = true;
// Allocate filter chain buffers..
m_filtered = (EventDataType *) malloc(max_filter_buf_size);
}
FlowParser::~FlowParser()
{
free(m_filtered);
// for (int i=0;i<num_filter_buffers;i++) {
// free(m_buffers[i]);
// }
}
void FlowParser::clearFilters()
{
m_filters.clear();
}
EventDataType *FlowParser::applyFilters(EventDataType *data, int samples)
{
EventDataType *in = nullptr, *out = nullptr;
if (m_filters.size() == 0) {
//qDebug() << "Trying to apply empty filter list in FlowParser..";
return nullptr;
}
for (int i = 0; i < num_filter_buffers; i++) {
m_buffers[i] = (EventDataType *) malloc(max_filter_buf_size);
}
int numfilt = m_filters.size();
for (int i = 0; i < numfilt; i++) {
if (i == 0) {
in = data;
out = m_buffers[0];
if (in == out) {
//qDebug() << "Error: If you need to use internal m_buffers as initial input, use the second one. No filters were applied";
return nullptr;
}
} else {
in = m_buffers[(i + 1) % num_filter_buffers];
out = m_buffers[i % num_filter_buffers];
}
// If final link in chain, pass it back out to input data
if (i == numfilt - 1) {
out = data;
}
Filter &filter = m_filters[i];
if (filter.type == FilterNone) {
// Just copy it..
memcpy(in, out, samples * sizeof(EventDataType));
} else if (filter.type == FilterPercentile) {
percentileFilter(in, out, samples, filter.param1, filter.param2);
} else if (filter.type == FilterXPass) {
xpassFilter(in, out, samples, filter.param1);
}
}
for (int i = 0; i < num_filter_buffers; i++) {
free(m_buffers[i]);
}
return out;
}
// Opens the flow rate EventList, applies the input filter chain, and calculates peaks
void FlowParser::openFlow(Session *session, EventList *flow)
{
if (!flow) {
qDebug() << "called FlowParser::processFlow() with a null EventList!";
return;
}
m_session = session;
m_flow = flow;
m_gain = flow->gain();
m_rate = flow->rate();
m_samples = flow->count();
EventStoreType *inraw = flow->rawData();
// Make sure we won't overflow internal buffers
if (m_samples > max_filter_buf_size) {
qDebug() << "Error: Sample size exceeds max_filter_buf_size in FlowParser::openFlow().. Capping!!!";
m_samples = max_filter_buf_size;
}
// Begin with the second internal buffer
EventDataType *buf = m_filtered;
// Apply gain to waveform
EventStoreType *eptr = inraw + m_samples;
// Convert from store type to floats..
for (; inraw < eptr; ++inraw) {
*buf++ = EventDataType(*inraw) * m_gain;
}
// Apply the rest of the filters chain
buf = applyFilters(m_filtered, m_samples);
// Scan for and create an index of each breath
calcPeaks(m_filtered, m_samples);
}
// Calculates breath upper & lower peaks for a chunk of EventList data
void FlowParser::calcPeaks(EventDataType *input, int samples)
{
if (samples <= 0) {
return;
}
EventDataType min = 0, max = 0, c, lastc = 0;
EventDataType zeroline = 0;
double rate = m_flow->rate();
double flowstart = m_flow->first();
double lasttime, time;
double peakmax = flowstart, peakmin = flowstart;
time = lasttime = flowstart;
breaths.clear();
// Estimate storage space needed using typical average breaths per minute.
m_minutes = double(m_flow->last() - m_flow->first()) / 60000.0;
const double avgbpm = 20; // average breaths per minute of a standard human
int guestimate = m_minutes * avgbpm;
// reserve some memory
breaths.reserve(guestimate);
// Prime min & max, and see which side of the zero line we are starting from.
c = input[0];
min = max = c;
lastc = c;
m_startsUpper = (c >= zeroline);
qint32 start = 0, middle = 0;
int sps = 1000 / m_rate;
int len = 0, lastk = 0;
qint64 sttime = time;
// For each samples, find the peak upper and lower breath components
for (int k = 0; k < samples; k++) {
c = input[k];
if (c >= zeroline) {
// Did we just cross the zero line going up?
if (lastc < zeroline) {
// This helps filter out dirty breaths..
len = k - start;
if ((max > 3) && ((max - min) > 8) && (len > sps) && (middle > start)) {
// peak detection may not be needed..
breaths.push_back(BreathPeak(min, max, start, middle, k));
// Set max for start of the upper breath cycle
max = c;
peakmax = time;
// Starting point of next breath cycle
start = k;
sttime = time;
}
} else if (c > max) {
// Update upper breath peak
max = c;
peakmax = time;
}
}
if (c < zeroline) {
// Did we just cross the zero line going down?
if (lastc >= zeroline) {
// Set min for start of the lower breath cycle
min = c;
peakmin = time;
middle = k;
} else if (c < min) {
// Update lower breath peak
min = c;
peakmin = time;
}
}
lasttime = time;
time += rate;
lastc = c;
lastk = k;
}
}
//! \brief Calculate Respiratory Rate, TidalVolume, Minute Ventilation, Ti & Te..
// These are grouped together because, a) it's faster, and b) some of these calculations rely on others.
void FlowParser::calc(bool calcResp, bool calcTv, bool calcTi, bool calcTe, bool calcMv)
{
if (!m_session) {
return;
}
// Don't even bother if only a few breaths in this chunk
const int lowthresh = 4;
int nm = breaths.size();
if (nm < lowthresh) {
return;
}
const qint64 minute = 60000;
double start = m_flow->first();
double time = start;
int bs, be, bm;
double st, et, mt;
/////////////////////////////////////////////////////////////////////////////////
// Respiratory Rate setup
/////////////////////////////////////////////////////////////////////////////////
EventDataType minrr = 0, maxrr = 0;
EventList *RR = nullptr;
quint32 *rr_tptr = nullptr;
EventStoreType *rr_dptr = nullptr;
if (calcResp) {
RR = m_session->AddEventList(CPAP_RespRate, EVL_Event);
minrr = RR->Min(), maxrr = RR->Max();
RR->setGain(0.2F);
RR->setFirst(time + minute);
RR->getData().resize(nm);
RR->getTime().resize(nm);
rr_tptr = RR->rawTime();
rr_dptr = RR->rawData();
}
int rr_count = 0;
double len, st2, et2, adj, stmin, b, rr = 0;
double len2;
/////////////////////////////////////////////////////////////////////////////////
// Inspiratory / Expiratory Time setup
/////////////////////////////////////////////////////////////////////////////////
double lastte2 = 0, lastti2 = 0, lastte = 0, lastti = 0, te, ti, ti1, te1, c;
EventList *Te = nullptr, * Ti = nullptr;
if (calcTi) {
Ti = m_session->AddEventList(CPAP_Ti, EVL_Event);
Ti->setGain(0.02F);
}
if (calcTe) {
Te = m_session->AddEventList(CPAP_Te, EVL_Event);
Te->setGain(0.02F);
}
/////////////////////////////////////////////////////////////////////////////////
// Tidal Volume setup
/////////////////////////////////////////////////////////////////////////////////
EventList *TV = nullptr;
EventDataType mintv = 0, maxtv = 0, tv = 0;
double val1, val2;
quint32 *tv_tptr = nullptr;
EventStoreType *tv_dptr = nullptr;
int tv_count = 0;
if (calcTv) {
TV = m_session->AddEventList(CPAP_TidalVolume, EVL_Event);
mintv = TV->Min(), maxtv = TV->Max();
TV->setGain(20);
TV->setFirst(start);
TV->getData().resize(nm);
TV->getTime().resize(nm);
tv_tptr = TV->rawTime();
tv_dptr = TV->rawData();
}
/////////////////////////////////////////////////////////////////////////////////
// Minute Ventilation setup
/////////////////////////////////////////////////////////////////////////////////
EventList *MV = nullptr;
EventDataType mv;
if (calcMv) {
MV = m_session->AddEventList(CPAP_MinuteVent, EVL_Event);
MV->setGain(0.125); // gain set to 1/8th
}
EventDataType sps = (1000.0 / m_rate); // Samples Per Second
qint32 timeval = 0; // Time relative to start
BreathPeak * bpstr = breaths.data();
BreathPeak * bpend = bpstr + nm;
// For each breath...
for (BreathPeak * bp = bpstr; bp != bpend; ++bp) {
bs = bp->start;
bm = bp->middle;
be = bp->end;
// Calculate start, middle and end time of this breath
st = start + bs * m_rate;
mt = start + bm * m_rate;
timeval = be * m_rate;
et = start + timeval;
/////////////////////////////////////////////////////////////////////
// Calculate Inspiratory Time (Ti) for this breath
/////////////////////////////////////////////////////////////////////
if (calcTi) {
// Ti is simply the time between the start of a breath and it's midpoint
// Note the 50.0 is the gain value, to give it better resolution
// (and mt and st are in milliseconds)
ti = ((mt - st) / 1000.0) * 50.0;
// Average for a little smoothing
ti1 = (lastti2 + lastti + ti * 2) / 4.0;
Ti->AddEvent(mt, ti1); // Add the event
// Track the last two values to use for averaging
lastti2 = lastti;
lastti = ti;
}
/////////////////////////////////////////////////////////////////////
// Calculate Expiratory Time (Te) for this breath
/////////////////////////////////////////////////////////////////////
if (calcTe) {
// Te is simply the time between the second half of the breath
te = ((et - mt) / 1000.0) * 50.0;
// Average last three values..
te1 = (lastte2 + lastte + te * 2 ) / 4.0;
Te->AddEvent(mt, te1);
lastte2 = lastte;
lastte = te;
}
/////////////////////////////////////////////////////////////////////
// Calculate TidalVolume
/////////////////////////////////////////////////////////////////////
if (calcTv) {
val1 = 0, val2 = 0;
// Scan the upper breath
for (int j = bs; j < bm; j++) {
// convert flow to ml/s to L/min and divide by samples per second
c = double(qAbs(m_filtered[j])) * 1000.0 / 60.0 / sps;
val2 += c;
//val2+=c*c; // for RMS
}
tv = val2;
bool usebothhalves = false;
if (usebothhalves) {
for (int j = bm; j < be; j++) {
// convert flow to ml/s to L/min and divide by samples per second
c = double(qAbs(m_filtered[j])) * 1000.0 / 60.0 / sps;
val1 += c;
//val1 += c*c; // for RMS
}
tv = (qAbs(val2) + qAbs(val1)) / 2;
}
// Add the other half here and average might make it more accurate
// But last time I tried it was pretty messy
// Perhaps needs a bit of history averaging..
// calculate root mean square
//double n=bm-bs;
//double q=(1/n)*val2;
//double x=sqrt(q)*2;
//val2=x;
if (tv < mintv) { mintv = tv; }
if (tv > maxtv) { maxtv = tv; }
*tv_tptr++ = timeval;
*tv_dptr++ = tv / 20.0;
tv_count++;
}
/////////////////////////////////////////////////////////////////////
// Respiratory Rate Calculations
/////////////////////////////////////////////////////////////////////
if (calcResp) {
stmin = et - minute;
if (stmin < start) {
stmin = start;
}
len = et - stmin;
rr = 0;
len2 = 0;
// Step back through last minute and count breaths
BreathPeak *bpstr1 = bpstr-1;
for (BreathPeak *p = bp; p != bpstr1; --p) {
st2 = start + double(p->start) * m_rate;
et2 = start + double(p->end) * m_rate;
if (et2 < stmin) {
break;
}
len = et2 - st2;
if (st2 < stmin) {
// Partial breath
st2 = stmin;
adj = et2 - st2;
b = (1.0 / len) * adj;
len2 += adj;
} else {
b = 1;
len2 += len;
}
rr += b;
}
if (len2 < minute) {
rr *= minute / len2;
}
// Calculate min & max
if (rr < minrr) {
minrr = rr;
}
if (rr > maxrr) {
maxrr = rr;
}
// Add manually.. (much quicker)
*rr_tptr++ = timeval;
// Use the same gains as ResMed..
*rr_dptr++ = rr * 5.0;
rr_count++;
}
if (calcMv && calcResp && calcTv) {
// Minute Ventilation is tidal volume times respiratory rate
mv = (tv / 1000.0) * rr;
// The 8.0 is the gain of the MV EventList to boost resolution
MV->AddEvent(et, mv * 8.0);
}
}
/////////////////////////////////////////////////////////////////////
// Respiratory Rate post filtering
/////////////////////////////////////////////////////////////////////
if (calcResp) {
RR->setMin(minrr);
RR->setMax(maxrr);
RR->setFirst(start);
RR->setLast(et);
RR->setCount(rr_count);
}
/////////////////////////////////////////////////////////////////////
// Tidal Volume post filtering
/////////////////////////////////////////////////////////////////////
if (calcTv) {
TV->setMin(mintv);
TV->setMax(maxtv);
TV->setFirst(start);
TV->setLast(et);
TV->setCount(tv_count);
}
}
void FlowParser::flagUserEvents(ChannelID code, EventDataType restriction, EventDataType duration)
{
int numbreaths = breaths.size();
EventDataType val, mx, mn;
QVector<EventDataType> br;
QVector<qint32> bstart;
QVector<qint32> bend;
// Allocate some memory beforehand so it doesn't have to slow it down mid calculations
bstart.reserve(numbreaths * 2);
bend.reserve(numbreaths * 2);
br.reserve(numbreaths * 2);
double start = m_flow->first();
double st, et, dur;
qint64 len;
bool allowDuplicates = p_profile->cpap->userEventDuplicates();
// Get the Breath list, which is calculated by the previously run breath marker algorithm.
BreathPeak *bpstr = breaths.data();
BreathPeak *bpend = bpstr + numbreaths;
// Create a list containing the abs of min and max waveform flow for each breath
for (BreathPeak *p = bpstr; p != bpend; ++p) {
br.push_back(qAbs(p->max));
br.push_back(qAbs(p->min));
}
// The following ignores the outliers to get a cleaner cutoff value
// 60th percentile was chosen because it's a little more than median, and, well, reasons I can't remember specifically.. :-}
// Look for the 60th percentile of the abs'ed min/max values
const EventDataType perc = 0.6F;
int idx = float(br.size()) * perc;
nth_element(br.begin(), br.begin() + idx, br.end() - 1);
// Take this value as the peak
EventDataType peak = br[idx]; ;
// Scale the restriction percentage to the peak to find the cutoff value
EventDataType cutoffval = peak * (restriction / 100.0F);
int bs, bm, be, bs1, bm1, be1;
// For each Breath, search for flow restrictions
for (BreathPeak *p = bpstr; p != bpend; ++p) {
// Todo: Define these markers in the comments better
bs = p->start; // breath start
bm = p->middle; // breath middle
be = p->end; // breath end
mx = p->max;
mn = p->min;
val = mx - mn; // the total height from top to bottom of this breath
// Scan the breath in the flow data and stop at the first location more than the cutoff value
// (Only really needs to scan to the middle.. I'm not sure why I made it go all the way to the end.)
for (bs1 = bs; bs1 < be; bs1++) {
if (qAbs(m_filtered[bs1]) > cutoffval) {
break;
}
}
// if bs1 has reached the end, this means the entire marked breath is within the restriction
// Scan backwards from the middle to the start, stopping at the first value past the cutoff value
for (bm1 = bm; bm1 > bs; bm1--) {
if (qAbs(m_filtered[bm1]) > cutoffval) {
break;
}
}
// Now check if a value was found in the first half of the breath above the cutoff value
if (bm1 >= bs1) {
// Good breath... And add it as a beginning/end marker for the next stage
bstart.push_back(bs1);
bend.push_back(bm1);
} // else we crossed over and the first half of the breath is under the threshold, therefore has a flow restricted
// Now do the other half of the breath....
// Scan from middle to end of breath, stopping at first cutoff value
for (bm1 = bm; bm1 < be; bm1++) {
if (qAbs(m_filtered[bm1]) > cutoffval) {
break;
}
}
// Scan backwards from the end to the middle of the breath, stopping at the first cutoff value
for (be1 = be; be1 > bm; be1--) {
if (qAbs(m_filtered[be1]) > cutoffval) {
break;
}
}
// Check for crossover again
if (be1 >= bm1) {
// Good strong healthy breath.. add the beginning and end to the breath markers.
bstart.push_back(bm1);
bend.push_back(be1);
} // else crossed over again.. breathe damn you!
// okay, this looks like cruft..
// st = start + bs1 * m_rate;
// et = start + be1 * m_rate;
}
int bsize = bstart.size(); // Number of breath components over cutoff threshold
EventList *uf = nullptr;
// For each good breath marker, look at the gaps in between
for (int i = 0; i < bsize - 1; i++) {
bs = bend[i]; // start at the end of the healthy breath
be = bstart[i + 1]; // look ahead to the beginning of the next one
// Calculate the start and end timestamps
st = start + bs * m_rate;
et = start + be * m_rate;
// Calculate the length of the flow restriction
len = et - st;
dur = len / 1000.0; // (scale to seconds, not milliseconds)
if (dur >= duration) { // is the event greater than the duration threshold?
// Unless duplicates have been specifically allowed, scan for any apnea's already detected by the machine
if (allowDuplicates || !SearchApnea(m_session, et, dur)) {
if (!uf) {
// Create event list if not already done
uf = m_session->AddEventList(code, EVL_Event);
}
// Add the user flag at the end
uf->AddEvent(et, dur);
}
}
}
}
void FlowParser::flagEvents()
{
if (!p_profile->cpap->userEventFlagging()) { return; }
int numbreaths = breaths.size();
if (numbreaths < 5) { return; }
flagUserEvents(CPAP_UserFlag1, p_profile->cpap->userFlowRestriction(), p_profile->cpap->userEventDuration());
flagUserEvents(CPAP_UserFlag2, p_profile->cpap->userFlowRestriction2(), p_profile->cpap->userEventDuration2());
}
void calcRespRate(Session *session, FlowParser *flowparser)
{
if (session->machine()->type() != MT_CPAP) { return; }
// if (session->machine()->loaderName() != STR_MACH_PRS1) return;
if (!session->eventlist.contains(CPAP_FlowRate)) {
//qDebug() << "calcRespRate called without FlowRate waveform available";
return; //need flow waveform
}
bool trashfp;
if (!flowparser) {
flowparser = new FlowParser();
trashfp = true;
//qDebug() << "calcRespRate called without valid FlowParser object.. using a slow throw-away!";
//return;
} else {
trashfp = false;
}
bool calcResp = !session->eventlist.contains(CPAP_RespRate);
bool calcTv = !session->eventlist.contains(CPAP_TidalVolume);
bool calcTi = !session->eventlist.contains(CPAP_Ti);
bool calcTe = !session->eventlist.contains(CPAP_Te);
bool calcMv = !session->eventlist.contains(CPAP_MinuteVent);
int z = (calcResp ? 1 : 0) + (calcTv ? 1 : 0) + (calcMv ? 1 : 0);
// If any of these three missing, remove all, and switch all on
if (z > 0 && z < 3) {
if (!calcResp && !calcTv && !calcMv) {
calcTv = calcMv = calcResp = true;
}
QVector<EventList *> &list = session->eventlist[CPAP_RespRate];
int size = list.size();
for (int i = 0; i < size; ++i) {
delete list[i];
}
session->eventlist[CPAP_RespRate].clear();
QVector<EventList *> &list2 = session->eventlist[CPAP_TidalVolume];
size = list2.size();
for (int i = 0; i < size; ++i) {
delete list2[i];
}
session->eventlist[CPAP_TidalVolume].clear();
QVector<EventList *> &list3 = session->eventlist[CPAP_MinuteVent];
size = list3.size();
for (int i = 0; i < size; ++i) {
delete list3[i];
}
session->eventlist[CPAP_MinuteVent].clear();
}
flowparser->clearFilters();
// No filters works rather well with the new peak detection algorithm..
// Although the output could use filtering.
//flowparser->addFilter(FilterPercentile,7,0.5);
//flowparser->addFilter(FilterPercentile,5,0.5);
//flowparser->addFilter(FilterXPass,0.5);
EventList *flow;
QVector<EventList *> &EVL = session->eventlist[CPAP_FlowRate];
int size = EVL.size();
for (int ws = 0; ws < size; ++ws) {
flow = EVL[ws];
if (flow->count() > 20) {
flowparser->openFlow(session, flow);
flowparser->calc(calcResp, calcTv, calcTi , calcTe, calcMv);
flowparser->flagEvents();
}
}
if (trashfp) {
delete flowparser;
}
}
EventDataType calcAHI(Session *session, qint64 start, qint64 end)
{
bool rdi = p_profile->general->calculateRDI();
double hours, ahi, cnt;
if (start < 0) {
// much faster..
hours = session->hours();
cnt = session->count(CPAP_Obstructive)
+ session->count(CPAP_Hypopnea)
+ session->count(CPAP_ClearAirway)
+ session->count(CPAP_Apnea);
if (rdi) {
cnt += session->count(CPAP_RERA);
}
ahi = cnt / hours;
} else {
hours = double(end - start) / 3600000L;
if (hours == 0) { return 0; }
cnt = session->rangeCount(CPAP_Obstructive, start, end)
+ session->rangeCount(CPAP_Hypopnea, start, end)
+ session->rangeCount(CPAP_ClearAirway, start, end)
+ session->rangeCount(CPAP_Apnea, start, end);
if (rdi) {
cnt += session->rangeCount(CPAP_RERA, start, end);
}
ahi = cnt / hours;
}
return ahi;
}
int calcAHIGraph(Session *session)
{
bool calcrdi = session->machine()->loaderName() == "PRS1";
const qint64 window_step = 30000; // 30 second windows
double window_size = p_profile->cpap->AHIWindow();
qint64 window_size_ms = window_size * 60000L;
bool zeroreset = p_profile->cpap->AHIReset();
if (session->machine()->type() != MT_CPAP) { return 0; }
bool hasahi = session->eventlist.contains(CPAP_AHI);
bool hasrdi = session->eventlist.contains(CPAP_RDI);
if (hasahi && hasrdi) {
return 0; // abort if already there
}
if (!(!hasahi && !hasrdi)) {
session->destroyEvent(CPAP_AHI);
session->destroyEvent(CPAP_RDI);
}
if (!session->channelExists(CPAP_Obstructive) &&
!session->channelExists(CPAP_Hypopnea) &&
!session->channelExists(CPAP_Apnea) &&
!session->channelExists(CPAP_ClearAirway) &&
!session->channelExists(CPAP_RERA)
) { return 0; }
qint64 first = session->first(),
last = session->last(),
f;
EventList *AHI = new EventList(EVL_Event);
AHI->setGain(0.02F);
session->eventlist[CPAP_AHI].push_back(AHI);
EventList *RDI = nullptr;
if (calcrdi) {
RDI = new EventList(EVL_Event);
RDI->setGain(0.02F);
session->eventlist[CPAP_RDI].push_back(RDI);
}
EventDataType ahi, rdi;
qint64 ti = first, lastti = first;
double avgahi = 0;
double avgrdi = 0;
int cnt = 0;
double events;
double hours = (window_size / 60.0F);
if (zeroreset) {
// I personally don't see the point of resetting each hour.
do {
// For each window, in 30 second increments
for (qint64 t = ti; t < ti + window_size_ms; t += window_step) {
if (t > last) {
break;
}
events = session->rangeCount(CPAP_Obstructive, ti, t)
+ session->rangeCount(CPAP_Hypopnea, ti, t)
+ session->rangeCount(CPAP_ClearAirway, ti, t)
+ session->rangeCount(CPAP_Apnea, ti, t);
ahi = events / hours;
AHI->AddEvent(t, ahi * 50);
avgahi += ahi;
if (calcrdi) {
events += session->rangeCount(CPAP_RERA, ti, t);
rdi = events / hours;
RDI->AddEvent(t, rdi * 50);
avgrdi += rdi;
}
cnt++;
}
lastti = ti;
ti += window_size_ms;
} while (ti < last);
} else {
for (ti = first; ti < last; ti += window_step) {
f = ti - window_size_ms;
//hours=window_size; //double(ti-f)/3600000L;
events = session->rangeCount(CPAP_Obstructive, f, ti)
+ session->rangeCount(CPAP_Hypopnea, f, ti)
+ session->rangeCount(CPAP_ClearAirway, f, ti)
+ session->rangeCount(CPAP_Apnea, f, ti);
ahi = events / hours;
avgahi += ahi;
AHI->AddEvent(ti, ahi * 50);
if (calcrdi) {
events += session->rangeCount(CPAP_RERA, f, ti);
rdi = events / hours;
RDI->AddEvent(ti, rdi * 50);
avgrdi += rdi;
}
cnt++;
lastti = ti;
ti += window_step;
}
}
AHI->AddEvent(last, 0);
if (calcrdi) {
RDI->AddEvent(last, 0);
}
if (!cnt) {
avgahi = 0;
avgrdi = 0;
} else {
avgahi /= double(cnt);
avgrdi /= double(cnt);
}
cnt++;
session->setAvg(CPAP_AHI, avgahi);
if (calcrdi) {
session->setAvg(CPAP_RDI, avgrdi);
}
return cnt;
}
struct TimeValue {
TimeValue() {
time = 0;
value = 0;
}
TimeValue(qint64 t, EventStoreType v) {
time = t;
value = v;
}
TimeValue(const TimeValue ©) {
time = copy.time;
value = copy.value;
}
qint64 time;
EventStoreType value;
};
struct zMaskProfile {
public:
zMaskProfile(MaskType type, QString name);
~zMaskProfile();
void reset() { pressureleaks.clear(); }
void scanLeaks(Session *session);
void scanPressure(Session *session);
void updatePressureMin();
void updateProfile(Session *session);
void load(Profile *profile);
void save();
QMap<EventStoreType, EventDataType> pressuremax;
QMap<EventStoreType, EventDataType> pressuremin;
QMap<EventStoreType, EventDataType> pressurecount;
QMap<EventStoreType, EventDataType> pressuretotal;
QMap<EventStoreType, EventDataType> pressuremean;
QMap<EventStoreType, EventDataType> pressurestddev;
QVector<TimeValue> Pressure;
EventDataType calcLeak(EventStoreType pressure);
protected:
static const quint32 version = 0;
void scanLeakList(EventList *leak);
void scanPressureList(EventList *el);
MaskType m_type;
QString m_name;
Profile *m_profile;
QString m_filename;
QMap<EventStoreType, QMap<EventStoreType, quint32> > pressureleaks;
EventDataType maxP, minP, maxL, minL, m_factor;
};
bool operator<(const TimeValue &p1, const TimeValue &p2)
{
return p1.time < p2.time;
}
zMaskProfile::zMaskProfile(MaskType type, QString name)
: m_type(type), m_name(name)
{
}
zMaskProfile::~zMaskProfile()
{
save();
}
void zMaskProfile::load(Profile *profile)
{
m_profile = profile;
m_filename = m_profile->Get("{" + STR_GEN_DataFolder + "}/MaskProfile.mp");
QFile f(m_filename);
if (!f.open(QFile::ReadOnly)) {
return;
}
QDataStream in(&f);
in.setVersion(QDataStream::Qt_4_6);
in.setByteOrder(QDataStream::LittleEndian);
quint32 m, v;
in >> m;
in >> v;
if (m != magic) {
qDebug() << "Magic wrong in zMaskProfile::load";
}
in >> pressureleaks;
f.close();
}
void zMaskProfile::save()
{
if (m_filename.isEmpty()) {
return;
}
QFile f(m_filename);
if (!f.open(QFile::WriteOnly)) {
return;
}
QDataStream out(&f);
out.setVersion(QDataStream::Qt_4_6);
out.setByteOrder(QDataStream::LittleEndian);
out << (quint32)magic;
out << (quint32)version;
out << pressureleaks;
f.close();
}
void zMaskProfile::scanPressureList(EventList *el)
{
qint64 start = el->first();
int count = el->count();
EventStoreType *dptr = el->rawData();
EventStoreType *eptr = dptr + count;
quint32 *tptr = el->rawTime();
qint64 time;
EventStoreType pressure;
for (; dptr < eptr; dptr++) {
time = start + *tptr++;
pressure = *dptr;
Pressure.push_back(TimeValue(time, pressure));
}
}
void zMaskProfile::scanPressure(Session *session)
{
Pressure.clear();
int prescnt = 0;
//EventStoreType pressure;
if (session->eventlist.contains(CPAP_Pressure)) {
prescnt = session->count(CPAP_Pressure);
Pressure.reserve(prescnt);
// WTF IS THIS DOING??? WHY THE HECK DID I PUT AN INNER LOOP HERE??
QVector<EventList *> &EVL=session->eventlist[CPAP_Pressure];
int size = EVL.size();
for (int j = 0; j < size; ++j) {
scanPressureList(EVL[j]);
// QVector<EventList *> &el = session->eventlist[CPAP_Pressure];
// for (int e = 0; e < el.size(); e++) {
// scanPressureList(el[e]);
// }
}
} else if (session->eventlist.contains(CPAP_IPAP)) {
prescnt = session->count(CPAP_IPAP);
Pressure.reserve(prescnt);
QVector<EventList *> &EVL=session->eventlist[CPAP_IPAP];
int size = EVL.size();
for (int j = 0; j < size; ++j) {
scanPressureList(EVL[j]);
}
}
qSort(Pressure);
}
void zMaskProfile::scanLeakList(EventList *el)
{
qint64 start = el->first();
int count = el->count();
EventStoreType *dptr = el->rawData();
EventStoreType *eptr = dptr + count;
quint32 *tptr = el->rawTime();
//EventDataType gain=el->gain();
EventStoreType pressure, leak;
//EventDataType fleak;
QMap<EventStoreType, EventDataType>::iterator pmin;
qint64 ti;
bool found;
int psize = Pressure.size();
if (psize == 0) return;
TimeValue *tvstr = Pressure.data();
TimeValue *tvend = tvstr + (psize - 1);
TimeValue *p1, *p2;
for (; dptr < eptr; dptr++) {
leak = *dptr;
ti = start + *tptr++;
found = false;
pressure = Pressure[0].value;
if (psize > 1) {
for (p1 = tvstr; p1 != tvend; ++p1) {
p2 = p1+1;
if ((p2->time > ti) && (p1->time <= ti)) {
pressure = p1->value;
found = true;
break;
} else if (p2->time == ti) {
pressure = p2->value;
found = true;
break;
}
// for (int i = 0; i < Pressure.size() - 1; i++) {
// const TimeValue &p1 = Pressure[i];
// const TimeValue &p2 = Pressure[i + 1];
// if ((p2.time > ti) && (p1.time <= ti)) {
// pressure = p1.value;
// found = true;
// break;
// } else if (p2.time == ti) {
// pressure = p2.value;
// found = true;
// break;
// }
}
} else {
found = true;
}
if (found) {
pressureleaks[pressure][leak]++;
// pmin=pressuremin.find(pressure);
// fleak=EventDataType(leak) * gain;
// if (pmin==pressuremin.end()) {
// pressuremin[pressure]=fleak;
// pressurecount[pressure]=pressureleaks[pressure][leak];
// } else {
// if (pmin.value() > fleak) {
// pmin.value()=fleak;
// pressurecount[pressure]=pressureleaks[pressure][leak];
// }
// }
} else {
//int i=5;
}
}
}
void zMaskProfile::scanLeaks(Session *session)
{
QVector<EventList *> &elv = session->eventlist[CPAP_LeakTotal];
int size=elv.size();
if (!size)
return;
QVector<EventList *>::iterator end = elv.end();
QVector<EventList *>::iterator it;
for (it = elv.begin(); it != end; ++it) {
scanLeakList(*it);
}
}
void zMaskProfile::updatePressureMin()
{
QMap<EventStoreType, QMap<EventStoreType, quint32> >::iterator it;
EventStoreType pressure;
//EventDataType perc=0.1;
//EventDataType tmp,l1,l2;
//int idx1, idx2,
int lks;
double SN;
double percentile = 0.40;
double p = 100.0 * percentile;
double nth, nthi;
int sum1, sum2, w1, w2, N, k;
QMap<EventStoreType, QMap<EventStoreType, quint32> >::iterator plend = pressureleaks.end();
QMap<EventStoreType, quint32>::iterator lmend;
for (it = pressureleaks.begin(); it != plend; it++) {
pressure = it.key();
QMap<EventStoreType, quint32> &leakmap = it.value();
lks = leakmap.size();
SN = 0;
// First sum total counts of all leaks
lmend = leakmap.end();
for (QMap<EventStoreType, quint32>::iterator lit = leakmap.begin(); lit != lmend; ++lit) {
SN += lit.value();
}
nth = double(SN) * percentile; // index of the position in the unweighted set would be
nthi = floor(nth);
sum1 = 0, sum2 = 0;
w1 = 0, w2 = 0;
EventDataType v1 = 0, v2;
N = leakmap.size();
k = 0;
bool found = false;
double total = 0;
// why do this effectively twice? and k = size
for (QMap<EventStoreType, quint32>::iterator lit = leakmap.begin(); lit != lmend; ++lit, ++k) {
total += lit.value();
}
pressuretotal[pressure] = total;
for (QMap<EventStoreType, quint32>::iterator lit = leakmap.begin(); lit != lmend; ++lit, ++k) {
//for (k=0;k < N;k++) {
v1 = lit.key();
w1 = lit.value();
sum1 += w1;
if (sum1 > nthi) {
pressuremin[pressure] = v1;
pressurecount[pressure] = w1;
found = true;
break;
}
if (sum1 == nthi) {
break; // boundary condition
}
}
if (found) { continue; }
if (k >= N) {
pressuremin[pressure] = v1;
pressurecount[pressure] = w1;
continue;
}
v2 = (leakmap.begin() + (k + 1)).key();
w2 = (leakmap.begin() + (k + 1)).value();
sum2 = sum1 + w2;
// value lies between v1 and v2
double px = 100.0 / double(SN); // Percentile represented by one full value
// calculate percentile ranks
double p1 = px * (double(sum1) - (double(w1) / 2.0));
double p2 = px * (double(sum2) - (double(w2) / 2.0));
// calculate linear interpolation
double v = v1 + ((p - p1) / (p2 - p1)) * (v2 - v1);
pressuremin[pressure] = v;
pressurecount[pressure] = w1;
}
}
EventDataType zMaskProfile::calcLeak(EventStoreType pressure)
{
if (maxP == minP) {
return pressuremin[pressure];
}
EventDataType leak = (pressure - minP) * (m_factor) + minL;
return leak;
}
void zMaskProfile::updateProfile(Session *session)
{
scanPressure(session);
scanLeaks(session);
updatePressureMin();
if (pressuremin.size() <= 1) {
maxP = minP = 0;
maxL = minL = 0;
m_factor = 0;
return;
}
EventDataType p, l, tmp, mean, sum;
minP = 250, maxP = 0;
minL = 1000, maxL = 0;
long cnt = 0, vcnt;
int n;
EventDataType maxcnt, maxval, lastval, lastcnt;
QMap<EventStoreType, QMap<EventStoreType, quint32> >::iterator plend = pressureleaks.end();
QMap<EventStoreType, QMap<EventStoreType, quint32> >::iterator it = pressureleaks.begin();
QMap<EventStoreType, quint32>::iterator lit;
QMap<EventStoreType, quint32>::iterator lvend;
for (; it != plend; ++it) {
p = it.key();
l = pressuremin[p];
QMap<EventStoreType, quint32> &leakval = it.value();
cnt = 0;
vcnt = -1;
n = leakval.size();
sum = 0;
maxcnt = 0, maxval = 0, lastval = 0, lastcnt = 0;
lvend = leakval.end();
for (lit = leakval.begin(); lit != lvend; ++lit) {
cnt += lit.value();
if (lit.value() > maxcnt) {
lastcnt = maxcnt;
maxcnt = lit.value();
lastval = maxval;
maxval = lit.key();
}
sum += lit.key() * lit.value();
if (lit.key() == (EventStoreType)l) {
vcnt = lit.value();
}
}
pressuremean[p] = mean = sum / EventDataType(cnt);
if (lastval > 0) {
maxval = qMax(maxval, lastval); //((maxval*maxcnt)+(lastval*lastcnt)) / (maxcnt+lastcnt);
}
pressuremax[p] = lastval;
sum = 0;
for (lit = leakval.begin(); lit != lvend; lit++) {
tmp = lit.key() - mean;
sum += tmp * tmp;
}
pressurestddev[p] = tmp = sqrt(sum / EventDataType(n));
if (vcnt >= 0) {
tmp = 1.0 / EventDataType(cnt) * EventDataType(vcnt);
}
}
QMap<EventStoreType, EventDataType> pressureval;
QMap<EventStoreType, EventDataType> pressureval2;
EventDataType max = 0, tmp2, tmp3;
QMap<EventStoreType, EventDataType>::iterator ptit;
QMap<EventStoreType, EventDataType>::iterator ptend = pressuretotal.end();
for (ptit = pressuretotal.begin(); ptit != ptend; ++ptit) {
if (max < ptit.value()) { max = ptit.value(); }
}
for (ptit = pressuretotal.begin(); ptit != pressuretotal.end(); ptit++) {
p = ptit.key();
tmp = pressurecount[p];
tmp2 = ptit.value();
tmp3 = (tmp / tmp2) * (tmp2 / max);
if (tmp3 > 0.05) {
tmp = pressuremean[p];
if (tmp > 0) {
pressureval[p] = tmp; // / tmp2;
if (p < minP) { minP = p; }
if (p > maxP) { maxP = p; }
if (tmp < minL) { minL = tmp; }
if (tmp > maxL) { maxL = tmp; }
}
}
}
if ((maxP == minP) || (minL == maxL) || (minP == 250) || (minL == 1000)) {
// crappy data set
maxP = minP = 0;
maxL = minL = 0;
m_factor = 0;
return;
}
m_factor = (maxL - minL) / (maxP - minP);
// for (QMap<EventStoreType, EventDataType>::iterator it=pressuremin.begin();it!=pressuremin.end()-1;it++) {
// p1=it.key();
// p2=(it+1).key();
// l1=it.value();
// l2=(it+1).value();
// if (l2 > l1) {
// factor=(l2 - l1) / (p2 - p1);
// sum+=factor;
// cnt++;
// }
// }
// m_factor=sum/double(cnt);
// int i=0;
// if (i==1) {
// QFile f("/home/mark/leaks.csv");
// f.open(QFile::WriteOnly);
// QTextStream out(&f);
// EventDataType p,l,c;
// QString fmt;
// for (QMap<EventStoreType, EventDataType>::iterator it=pressuremin.begin();it!=pressuremin.end();it++) {
// p=EventDataType(it.key()/10.0);
// l=it.value();
// fmt=QString("%1,%2\n").arg(p,0,'f',1).arg(l);
// out << fmt;
// }
// cruft
// for (QMap<EventStoreType, QMap<EventStoreType,quint32> >::iterator it=pressureleaks.begin();it!=pressureleaks.end();it++) {
// QMap<EventStoreType,quint32> & leakval=it.value();
// for (QMap<EventStoreType,quint32>::iterator lit=leakval.begin();lit!=leakval.end();lit++) {
// l=lit.key();
// c=lit.value();
// fmt=QString("%1,%2,%3\n").arg(p,0,'f',2).arg(l).arg(c);
// out << fmt;
// }
// }
// f.close();
// }
}
void flagLargeLeaks(Session *session)
{
// Already contains?
if (session->eventlist.contains(CPAP_LargeLeak))
return;
if (!session->eventlist.contains(CPAP_Leak))
return;
EventDataType threshold = p_profile->cpap->leakRedline();
if (threshold <= 0) {
return;
}
QVector<EventList *> & EVL = session->eventlist[CPAP_Leak];
int evlsize = EVL.size();
if (evlsize == 0)
return;
EventList * LL = nullptr;
qint64 time = 0;
EventDataType lastvalue;
qint64 leaktime;
for (int ec = 0; ec < evlsize; ++ec) {
EventList &el = *EVL[ec];
int count = el.count();
if (!count) continue;
leaktime = 0;
EventDataType leakvalue = 0;
lastvalue = -1;
for (int i=0; i < count; ++i) {
time = el.time(i);
EventDataType value = el.data(i);
if (value >= threshold) {
if (lastvalue < threshold) {
leaktime = time;
leakvalue = value;
}
} else if (lastvalue > threshold) {
if (!LL) {
LL=session->AddEventList(CPAP_LargeLeak, EVL_Event);
}
int duration = (time - leaktime) / 1000L;
LL->AddEvent(time, duration);
}
lastvalue = value;
}
}
if (lastvalue > threshold) {
if (!LL) {
LL=session->AddEventList(CPAP_LargeLeak, EVL_Event);
}
int duration = (time - leaktime) / 1000L;
LL->AddEvent(time, duration);
}
}
QMutex zMaskmutex;
zMaskProfile mmaskProfile(Mask_NasalPillows, "ResMed Swift FX");
bool mmaskFirst = true;
int calcLeaks(Session *session)
{
if (session->machine()->type() != MT_CPAP) { return 0; }
if (session->eventlist.contains(CPAP_Leak)) { return 0; } // abort if already there
if (!session->eventlist.contains(CPAP_LeakTotal)) { return 0; } // can't calculate without this..
zMaskmutex.lock();
zMaskProfile *maskProfile = &mmaskProfile;
if (mmaskFirst) {
mmaskFirst = false;
maskProfile->load(p_profile);
}
// if (!maskProfile) {
// maskProfile=new zMaskProfile(Mask_NasalPillows,"ResMed Swift FX");
// }
maskProfile->reset();
maskProfile->updateProfile(session);
EventList *leak = session->AddEventList(CPAP_Leak, EVL_Event, 1);
QVector<EventList *> & EVL = session->eventlist[CPAP_LeakTotal];
int evlsize = EVL.size();
TimeValue *p2, *pstr, *pend;
// can this go out of the loop?
int mppressize = maskProfile->Pressure.size();
pstr = maskProfile->Pressure.data();
pend = maskProfile->Pressure.data()+(mppressize-1);
for (int i = 0; i < evlsize; ++i) {
EventList &el = *EVL[i];
EventDataType gain = el.gain(), tmp, val;
int count = el.count();
EventStoreType *dptr = el.rawData();
EventStoreType *eptr = dptr + count;
quint32 *tptr = el.rawTime();
qint64 start = el.first(), ti;
EventStoreType pressure;
tptr = el.rawTime();
start = el.first();
bool found;
for (; dptr < eptr; ++dptr) {
tmp = EventDataType(*dptr) * gain;
ti = start + *tptr++;
found = false;
pressure = pstr->value;
for (TimeValue *p1 = pstr; p1 != pend; ++p1) {
p2 = p1+1;
if ((p2->time > ti) && (p1->time <= ti)) {
pressure = p1->value;
found = true;
break;
} else if (p2->time == ti) {
pressure = p2->value;
found = true;
break;
}
}
// for (int i = 0; i < mppressize - 1; ++i) {
// const TimeValue &p1 = &maskProfile->Pressure[i];
// const TimeValue &p2 = maskProfile->Pressure[i + 1];
// if ((p2.time > ti) && (p1.time <= ti)) {
// pressure = p1.value;
// found = true;
// break;
// } else if (p2.time == ti) {
// pressure = p2.value;
// found = true;
// break;
// }
// }
if (found) {
val = tmp - maskProfile->calcLeak(pressure);
if (val < 0) {
val = 0;
}
leak->AddEvent(ti, val);
}
}
}
zMaskmutex.unlock();
return leak->count();
}
int calcPulseChange(Session *session)
{
if (session->eventlist.contains(OXI_PulseChange)) { return 0; }
QHash<ChannelID, QVector<EventList *> >::iterator it = session->eventlist.find(OXI_Pulse);
if (it == session->eventlist.end()) { return 0; }
EventDataType val, val2, change, tmp;
qint64 time, time2;
qint64 window = p_profile->oxi->pulseChangeDuration();
window *= 1000;
change = p_profile->oxi->pulseChangeBPM();
EventList *pc = new EventList(EVL_Event, 1, 0, 0, 0, 0, true);
pc->setFirst(session->first(OXI_Pulse));
qint64 lastt;
EventDataType lv = 0;
int li = 0;
int max;
int size = it.value().size();
for (int e = 0; e < size; ++e) {
EventList &el = *(it.value()[e]);
int elcount=el.count();
for (int i = 0; i < elcount; ++i) {
val = el.data(i);
time = el.time(i);
lastt = 0;
lv = change;
max = 0;
for (int j = i + 1; j < elcount; ++j) { // scan ahead in the window
time2 = el.time(j);
if (time2 > time + window) { break; }
val2 = el.data(j);
tmp = qAbs(val2 - val);
if (tmp > lv) {
lastt = time2;
if (tmp > max) { max = tmp; }
//lv=tmp;
li = j;
}
}
if (lastt > 0) {
qint64 len = (lastt - time) / 1000.0;
pc->AddEvent(lastt, len, tmp);
i = li;
}
}
}
if (pc->count() == 0) {
delete pc;
return 0;
}
session->eventlist[OXI_PulseChange].push_back(pc);
session->setMin(OXI_PulseChange, pc->Min());
session->setMax(OXI_PulseChange, pc->Max());
session->setCount(OXI_PulseChange, pc->count());
session->setFirst(OXI_PulseChange, pc->first());
session->setLast(OXI_PulseChange, pc->last());
return pc->count();
}
int calcSPO2Drop(Session *session)
{
if (session->eventlist.contains(OXI_SPO2Drop)) { return 0; }
QHash<ChannelID, QVector<EventList *> >::iterator it = session->eventlist.find(OXI_SPO2);
if (it == session->eventlist.end()) { return 0; }
EventDataType val, val2, change, tmp;
qint64 time, time2;
qint64 window = p_profile->oxi->spO2DropDuration();
window *= 1000;
change = p_profile->oxi->spO2DropPercentage();
EventList *pc = new EventList(EVL_Event, 1, 0, 0, 0, 0, true);
qint64 lastt;
EventDataType lv = 0;
int li = 0;
// Fix me.. Time scale varies.
//const unsigned ringsize=30;
//EventDataType ring[ringsize]={0};
//qint64 rtime[ringsize]={0};
//int rp=0;
int min;
int cnt = 0;
tmp = 0;
qint64 start = 0;
// Calculate median baseline
QList<EventDataType> med;
int evsize = it.value().size();
for (int e = 0; e < evsize; ++e) {
EventList &el = *(it.value()[e]);
int elcount = el.count();
for (int i = 0; i < elcount; i++) {
val = el.data(i);
time = el.time(i);
if (val > 0) { med.push_back(val); }
if (!start) { start = time; }
if (time > start + 3600000) { break; } // just look at the first hour
tmp += val;
cnt++;
}
}
EventDataType baseline = 0;
if (med.size() > 0) {
qSort(med);
int midx = float(med.size()) * 0.90;
if (midx > med.size() - 1) { midx = med.size() - 1; }
if (midx < 0) { midx = 0; }
baseline = med[midx];
}
session->settings[OXI_SPO2Drop] = baseline;
//EventDataType baseline=round(tmp/EventDataType(cnt));
EventDataType current;
qDebug() << "Calculated baseline" << baseline;
for (int e = 0; e < evsize; ++e) {
EventList &el = *(it.value()[e]);
int elcount = el.count();
for (int i = 0; i < elcount; ++i) {
current = el.data(i);
if (!current) { continue; }
time = el.time(i);
/*ring[rp]=val;
rtime[rp]=time;
rp++;
rp=rp % ringsize;
if (i<ringsize) {
for (unsigned j=i;j<ringsize;j++) {
ring[j]=val;
rtime[j]=0;
}
}
tmp=0;
cnt=0;
for (unsigned j=0;j<ringsize;j++) {
if (rtime[j] > time-300000) { // only look at recent entries..
tmp+=ring[j];
cnt++;
}
}
if (!cnt) {
unsigned j=abs((rp-1) % ringsize);
tmp=(ring[j]+val)/2;
} else tmp/=EventDataType(cnt); */
val = baseline;
lastt = 0;
lv = val;
min = val;
for (int j = i; j < elcount; ++j) { // scan ahead in the window
time2 = el.time(j);
//if (time2 > time+window) break;
val2 = el.data(j);
if (val2 > baseline - change) { break; }
lastt = time2;
li = j + 1;
}
if (lastt > 0) {
qint64 len = (lastt - time);
if (len >= window) {
pc->AddEvent(lastt, len / 1000, val - min);
i = li;
}
}
}
}
if (pc->count() == 0) {
delete pc;
return 0;
}
session->eventlist[OXI_SPO2Drop].push_back(pc);
session->setMin(OXI_SPO2Drop, pc->Min());
session->setMax(OXI_SPO2Drop, pc->Max());
session->setCount(OXI_SPO2Drop, pc->count());
session->setFirst(OXI_SPO2Drop, pc->first());
session->setLast(OXI_SPO2Drop, pc->last());
return pc->count();
}