Homec4science
Diffusion Event-Based Adaptive ECG Sampling (public) Browse

R8502/experiments.py
f5c0730d2d6apublic

Last Change
Raw File

experiments.py
View Options

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
# pylint: disable-msg= C0103, E0401
"""
Created on Thu Feb 01 13:42:46 2019
Script performing the dynamic adaptation of the approximation threshold to reduce the sampling
while keeping QRS detection abilities.
@author: T. Teijeiro
"""
from operator import attrgetter
import bisect
import wfdb
from wfdb.processing import gqrs_detect, compare_annotations
import numpy as np
def WD_polyapprox(arr, epsilon=0.25):
"""
Performs a piecewise linear approximation to *arr* by means of the Wall-Danielsson algorithm,
described in "A fast sequential method for polygonal approximation of digitized curves". It
includes a correction to be more sensitive to prominent peaks in the signal, that is also
described in the same paper.
Parameters
----------
arr:
Sequence object (list, numpy array...) containing the signal.
epsilon:
Error threshold accepted in the approximation. It corresponds to the area deviation between
the signal and the approximation.
Returns
-------
out:
Numpy array with the indices of the samples selected in the approximation.
"""
if len(arr) < 3:
return np.arange(len(arr))
dx = 1.0
result = [0]
f = x = y = p = L = 0
for i in range(1, len(arr)):
dy = arr[i]-arr[i-1]
x += dx
y += dy
f += x*dy - y*dx
#Peak-test correction
Ln = abs(y)+x
p = (Ln < L) * ((p == 0)*(i-1) + p)
L = Ln
if abs(f) > epsilon:
t = i-1 if p == 0 else p
result.append(t)
f, x, y = 0, (i-t)*dx, arr[i]-arr[t]
p, L = 0, abs(y) + x
result.append(len(arr)-1)
return np.array(result)
if __name__ == "__main__":
import sys
DB_DIR = '/home/teijeiro/Documents/databases/mitdb/'
rname = '100' if len(sys.argv) < 2 else sys.argv[1]
lead = 'MLII'
rec = wfdb.rdrecord(DB_DIR+rname, channel_names=[lead], physical=False)
from construe.utils.units_helper import set_sampling_freq
#Sampling frequency set has to be before all other imports, where constants are defined.
set_sampling_freq(rec.fs)
import construe.acquisition.signal_buffer as SIG
from construe.knowledge.observables import RDeflection
from construe.knowledge.abstraction_patterns.segmentation.QRS import QRS_PATTERN
from construe.model.abstraction_pattern import AbstractionPattern
from construe.model.constraint_network import InconsistencyError
from construe.utils.MIT import MITAnnotation, save_annotations
#Reusable variables for abstraction pattern execution
trans = QRS_PATTERN.abstractions[RDeflection][0]
r = RDeflection()
def delineate_qrs(ann):
"""
Performs the QRS delineation from an annotation time point, using the knowledge defined
in the *Construe* QRS abstraction pattern. It simulates an abduction operation on a beat
annotation.
"""
p = AbstractionPattern(QRS_PATTERN)
r.time.set(ann, ann)
p.istate = trans.istate
p.fstate = trans.fstate
p.trseq.append((trans, r))
p.evidence[type(r)].append(r)
trans.tconst(p, r)
p.check_temporal_consistency()
trans.gconst(p, r)
return p.hypothesis
#Empty the global buffer
SIG.reset()
#Original digital signal
s = rec.d_signal.reshape((1, -1))[0]
#Parameters to build the signal fragment to run the gqrs algorithm for the next beat
W_BEFORE = 5*rec.fs
W_AFTER = 2*rec.fs
W_DELIN = int(0.4*rec.fs)
#Adaptive sampling threshold. The minimum value is 0.5mm2 in common ECG scales.
thres = INI_THRES = 150
#RR thresholds obtained from the gqrs default configuration.
RR_MIN = 0.28*rec.fs
RR_MAX = 2.4*rec.fs
#Result measures
hbeats = []
thseq = []
nsamples = 0
#Upper and lower limits for the next R peak.
llim, ulim = 0, int(RR_MAX)
#Last detected R peak
b = -int(RR_MIN)
while len(SIG._SIGNAL[lead]) < len(s)-1:
sbuf = SIG._SIGNAL[lead]
#Obtain the next window of the signal, resample it and attach it to the global buffer.
w_search = slice(len(sbuf), ulim + W_DELIN + W_AFTER)
x = WD_polyapprox(s[w_search], thres) + w_search.start
#The signal between W_DELIN and W_AFTER is just for gqrs proper working, will be discarded.
attch_pt = min(bisect.bisect_right(x, ulim + W_DELIN), len(x)-1)
nsamples += attch_pt
#We need an interpolated version of the subsampled signal to run the gqrs algorithm
s2 = np.interp(np.arange(x[0], x[-1]+1), x, s[x])
SIG._SIGNAL[lead] = np.append(sbuf, s2[:x[attch_pt]-w_search.start])
#Some previous signal is needed also for gqrs.
s_bef = sbuf[b-W_BEFORE:]
beats = w_search.start - len(s_bef) + gqrs_detect(d_sig=np.concatenate((s_bef, s2)),
fs=rec.fs, adc_gain=rec.adc_gain[0],
adc_zero=rec.adc_zero[0])
try:
#We are only interested in the next beat
nxt_b = next(q for q in beats if b+RR_MIN <= q <= ulim)
qrs = delineate_qrs(nxt_b)
#Update the reference to the delineated peak
b = int(qrs.time.start)
hbeats.append(b)
#Get the largest wave of the QRS complex, and the new proposed threshold
wmax = max(qrs.shape[lead].waves, key=attrgetter('amp'))
nthres = abs(wmax.amp*wmax.dur/8)
#Threshold reductions are always applied, but increases are done only in a conservative
#way (90%-10%) and if the rhythm was kept.
thres = nthres if nthres < thres else 0.9*thres + 0.1*nthres if b >= llim else thres
#Static RR limits.
llim, ulim = int(b + RR_MIN), int(b + RR_MAX)
if len(hbeats) > 2:
#Check for a regular rr variation
rr1, rr = np.diff(hbeats[-3:])
#If the rhythm is "regular", we set limits on the dynamic evolution
if abs(rr-rr1) <= 0.2*rec.fs or 0.8 <= rr/rr1 <= 1.2:
llim = int(b + max(RR_MIN, min(0.8*rr, rr-0.2*rec.fs)))
ulim = int(b + min(RR_MAX, max(1.2*rr, rr+0.2*rec.fs)))
except (StopIteration, InconsistencyError):
#Detection or delineation fails. We lower the threshold and resample the same length
thres = INI_THRES
b, ulim = ulim, 2*ulim-b
thseq.append(thres)
# print('{}, {}'.format(b, int(thres)))
print('Mean sampling rate for record {}: {:.2f}'.format(rname, nsamples*rec.fs/len(s)))
#Comparison detections are performed on the full regular and resampled signals.
gbeats = gqrs_detect(d_sig=s, fs=rec.fs, adc_gain=rec.adc_gain[0], adc_zero=rec.adc_zero[0])
hbeats = gqrs_detect(d_sig=SIG._SIGNAL[lead], fs=rec.fs, adc_gain=rec.adc_gain[0],
adc_zero=rec.adc_zero[0])
#Comparison between the detections in the raw signal and in the subsampled signal
compare_annotations(gbeats, hbeats, 0.15*rec.fs).print_summary()
save_annotations([MITAnnotation.MITAnnotation(1, b) for b in hbeats], DB_DIR+rname+'.wds')
save_annotations([MITAnnotation.MITAnnotation(1, b) for b in gbeats], DB_DIR+rname+'.gqrs')
c4science · Help