{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# GQRS Algorithm\n", "### An implementation for event-based sampled signal" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In this notebook, we're going to present a ri-formulation of the gqrs algorithm by George B. Moody (george@mit.edu), present in the [WFDB](https://www.physionet.org/physiotools/wag/gqrs-1.htm) library.\n", "\n", "The classical gQRS algorithm achieves above-the-average performances in the detection of the QRS complex.\n", "In this notebook, we will give an alternative formulation that can work with a sub-class of the event-based sampled signal while achieving good results (when compared to the original ones).\n", "\n", "In particular, the event-based sampled signals were generated from the ECG signals in the [MIT-BIH arrhythmia database](https://physionet.org/physiobank/database/mitdb/). For each signal, the sub-sampling is based on the idea to keep the error between the real(original) curve and the linearly interpolated signal below a given threshold. \n", "\n", "This gives us a strong prior: each pair of points in the sub-sampled signal can be linearly interpolated while keeping the error below the selected threshold. We will not execute a linear interpolation during this work but we will use this fact to obtain an efficient re-formulation.\n", "\n", "It's important to notice the total absence of the original smoothing filter present in the gQRS algorithm. This is because the sub-sampling can be thought of as a strong pice-wise smoothing.\n", "\n", "The original python implementation from which this work departed was originally posted on Github from: [MIT Laboratory for Computational Physiology](https://github.com/MIT-LCP/wfdb-python) " ] }, { "cell_type": "code", "execution_count": null, "metadata": { "ExecuteTime": { "end_time": "2020-01-27T15:20:10.682502Z", "start_time": "2020-01-27T15:20:10.455120Z" } }, "outputs": [], "source": [ "#Let's import some usefull standard library\n", "import numpy as np\n", "import matplotlib.pyplot as plt\n", "import pandas as pd\n", "import collections\n", "import copy\n", "import sys\n", "import os\n", "import time\n", "import pprint\n", "pp = pprint.PrettyPrinter(indent=4)\n", "\n", "#In order to write the results found to a standard ANNOTATION file that can be checked with the WFDB bxb tool\n", "scripts = '../data_parsing'\n", "if scripts not in sys.path:\n", " sys.path.insert(0,scripts)\n", "import MITAnnotation as MITA" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## File opening and showing" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Each event-based-sampled file is saved inside the data directory, inside the directory of the corresponding threshold.\n", "First let's select a threshold, after that we will be able to select the interested file(s).\n", "If we want to select all thresholds (and, after, all files) we can insert -1.\n", "\n", "The file format is a standard CSV file with 2 columns: time and value" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "ExecuteTime": { "end_time": "2020-01-27T15:20:10.686967Z", "start_time": "2020-01-27T15:20:10.683405Z" } }, "outputs": [], "source": [ "# We find the folder in 'dataOut' containing the data to use as input for the gQRS algorithm\n", "# change \"dataRaw\" with the output folder used for the ADC\n", "Folder = \"../data/dataOut\"\n", "dirOut = []\n", "\n", "for f in os.listdir(Folder):\n", " path = os.path.join(Folder,f)\n", " if 'lvlCrossingBits' in f and os.path.isdir(path):\n", " dirOut.append(path)\n", "\n", "print(\"Directories:\")\n", "print(\"Select using the corresponding number (-1 to use all the directories):\")\n", "for i,name in zip(range(len(dirOut)),dirOut):\n", " print(str(i)+\" --> \"+str(name))\n", "\n", "dataDict = dict()" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "ExecuteTime": { "end_time": "2020-01-27T15:20:10.693191Z", "start_time": "2020-01-27T15:20:10.688209Z" } }, "outputs": [], "source": [ "# Default folder variable, to test different folder, change here (-1 to run a complete test on all folders,\n", "# WARNING: a complete test will require a more time):\n", "defaultFolder = 0" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "ExecuteTime": { "end_time": "2020-01-27T15:20:10.700389Z", "start_time": "2020-01-27T15:20:10.694648Z" }, "scrolled": true }, "outputs": [], "source": [ "selectedFolder = defaultFolder\n", "if selectedFolder == -1:\n", " for i,name in zip(range(len(dirOut)),dirOut):\n", " #Find number of bits and histeresis lvl\n", " str1 = 'lvlCrossingBits'\n", " str2 = 'Hist'\n", " b = int(name[name.find(str1)+len(str1):name.find(str2)])\n", " h = int(name[name.find(str2)+len(str2):])\n", " #Create a dict\n", " dataDict[i]={'dir':name,'bits':b,'hist':h, 'data':dict()}\n", "elif selectedFolder >= len(dirOut):\n", " print(\"Invalid, used the firs as default\")\n", " name = dirOut[0]\n", " #Find number of bits and histeresis lvl\n", " str1 = 'lvlCrossingBits'\n", " str2 = 'Hist'\n", " b = int(name[name.find(str1)+len(str1):name.find(str2)])\n", " h = int(name[name.find(str2)+len(str2):])\n", " #Create a dict\n", " dataDict[0]={'dir':name,'bits':b,'hist':h, 'data' :dict()}\n", "else:\n", " name = dirOut[selectedFolder]\n", " #Find number of bits and histeresis lvl\n", " str1 = 'lvlCrossingBits'\n", " str2 = 'Hist'\n", " b = int(name[name.find(str1)+len(str1):name.find(str2)])\n", " h = int(name[name.find(str2)+len(str2):])\n", " #Create a dict\n", " dataDict[0]={'dir':name,'bits':b,'hist':h, 'data' :dict()}\n", "print(\"Using:\")\n", "print(dataDict)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "ExecuteTime": { "end_time": "2020-01-27T15:20:10.707933Z", "start_time": "2020-01-27T15:20:10.701293Z" } }, "outputs": [], "source": [ "fileList = [os.path.splitext(name)[0] for name in os.listdir(dataDict[0]['dir'])]\n", "print(\"Files:\")\n", "print(\"Select file using the file number (-1 to use all the directories):\")\n", "for name in fileList:\n", " print(str(name))" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "ExecuteTime": { "end_time": "2020-01-27T15:20:10.713288Z", "start_time": "2020-01-27T15:20:10.708790Z" } }, "outputs": [], "source": [ "# Default File variable, to test different files, change here (-1 to run a complete test on all files,\n", "# WARNING: a complete test will require a consistent ammount of time):\n", "defaultFile = 115" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "ExecuteTime": { "end_time": "2020-01-27T15:20:10.735680Z", "start_time": "2020-01-27T15:20:10.714314Z" }, "scrolled": true }, "outputs": [], "source": [ "select = defaultFile\n", "if select != -1 and str(select) not in fileList:\n", " select = fileList[0]\n", " print('not present, used: ',select)\n", "select = str(select)\n", "\n", "selectedFileNames = []\n", "if select != '-1':\n", " for i in range(len(dataDict)):\n", " fileThisThreshold = os.path.join(dataDict[i]['dir'],select+'.csv')\n", " data = pd.read_csv(fileThisThreshold,header = None)\n", " dataDict[i]['data'][str(select)] = {'t':data[0].values,'v':data[1].values}\n", " selectedFileNames.append(select)\n", "elif select == '-1':\n", " for i in range(len(dataDict)):\n", " for f in fileList:\n", " fileThisThreshold = os.path.join(dataDict[i]['dir'],f+'.csv')\n", " data = pd.read_csv(fileThisThreshold,header = None)\n", " dataDict[i]['data'][f] = {'t':data[0].values,'v':data[1].values}\n", " for f in fileList:\n", " selectedFileNames.append(f)\n", "\n", "pp.pprint(\"Slected file(s): \"+str(selectedFileNames))" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "ExecuteTime": { "end_time": "2020-01-27T15:20:10.739880Z", "start_time": "2020-01-27T15:20:10.736656Z" } }, "outputs": [], "source": [ "dataDict" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Let's observe one file:\n", "We ordered the \"data\" structure as it follows:\n", "\n", "dictionary[progressive number]:\n", " * ['dir'] = directory of the file \n", " * ['bits'] = number of bits \n", " * ['hist'] = hysteresis level \n", " * ['data']\n", " * ['$fileName']\n", " * ['t'] = list of times\n", " * ['v'] = list of values\n", "[Threshold][File][DataFrame]\n", "\n", "Select a starting and a stop time (in seconds)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "ExecuteTime": { "end_time": "2020-01-27T15:20:10.745628Z", "start_time": "2020-01-27T15:20:10.740880Z" } }, "outputs": [], "source": [ "print('Start: ',end = \" \")\n", "start = 10*360\n", "print(\"{} seconds\".format(start/360))\n", "print('Stop: ', end = \"\")\n", "stop = 15*360\n", "print(\"{} seconds\".format(stop/360))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now we define a function that returns us the indexes to use given the time vector:" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "ExecuteTime": { "end_time": "2020-01-27T15:20:10.751666Z", "start_time": "2020-01-27T15:20:10.746459Z" } }, "outputs": [], "source": [ "def getIdxs(time, start, stop):\n", " startIdx = 0\n", " stopIdx = len(t)-1\n", " for i in range(len(time)-1):\n", " if time[i]<=start and time[i+1]> start:\n", " startIdx = i\n", " if time[i]<=stop and time[i+1]> stop:\n", " stopIdx = i\n", " \n", " return startIdx,stopIdx" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "ExecuteTime": { "end_time": "2020-01-27T15:20:11.077342Z", "start_time": "2020-01-27T15:20:10.752966Z" }, "scrolled": false }, "outputs": [], "source": [ "fig, ax = plt.subplots(1, 1, sharey='row',figsize=(30,15))\n", "i = 0\n", "ax.set_title('ECG using ADC with bits = ' + str(dataDict[0]['bits']), fontsize=35)\n", "ax.set_xlabel(\"Time\", fontsize=35)\n", "ax.set_ylabel(\"ADC value\", fontsize=35)\n", "ax.tick_params(labelsize = 35)\n", "\n", "t = dataDict[0]['data'][selectedFileNames[0]]['t']\n", "v = dataDict[0]['data'][selectedFileNames[0]]['v']\n", "startI,stopI = getIdxs(t,start,stop)\n", "ax.plot(t[startI:stopI]/360,v[startI:stopI])\n", "ax.grid()\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "---\n", "---" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## gQRS Method\n", "Let's now see the gQRS method implementation:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Basics of the standard gQRS detector:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The standard gQRS is composed of 4 main phases. during these phases, a set of parameters is constantly re-adapted to find the (estimated and local) optimal values. The 4 phases can be summarized as:\n", "\n", "1. Smoothing and amplification\n", "2. \"Matched\" filter\n", "3. Peak detection and saving\n", "4. Peak lookback\n", "\n", "Each of these phases is important to understand in order to be able to re-formulate this algorithm as we like.\n", "\n", "Before going into details, is useful to discuss some meaningful parameter that the algorithm use to detect the QRS complex (Some learned, some fixed). All of them are described in \"number of samples\" instead of physical values, this is not a problem: the initial values, given by the user, are physical values that get transformed in sample numbers by the algorithm using the information about the ADC used :\n", "\n", "* rrmean = Average RR interval\n", "* rrdev = Standard deviation of the RR interval\n", "* rrmin = Minimum RR interval\n", "* rrmax = Maximum RR interval\n", "* rrinc = Increment step of the estimated RR interval (for learning)\n", "* dt = One-fourth of the standard duration of the QRS complex, other 3 variables like this are present (1, 1/2, 1/3 of the duration of the QRS complex)\n", "* pthr = Threshold for a cuspid to be considered peak\n", "* qthr = Threshold for a peak to be considered a candidate for QRS\n", "* pthmin = Lower possible value for pthr\n", "* qthmin = Lower possible value for qthr\n", "* tamean = Average amplitude of a T-Wave\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ " Smoothing and amplification\n", "During this phase the signal pass through a trapezoidal (R-C) filter and gets amplified by the fixed delta time that the algorithm suppose to be a standard QRS duration. This filter mainly try to remove the noise present" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### \"Matched\" filter\n", "Even if it's not immediate to get from the code, the \"Matched\" filter is a FIR filter that operates upon the integrated smoothed signal, it's shape is customly designed to obtain a strong output if a QRS complex is present and almost zero otherwise (except for strong T-Waves). We won't get any further in this section mainly because this filter is one of the keys to the good behavior of this algorithm and was one of the main blocks that needed to be re-formulated. We will see it in detail in the next section" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Peak detection and saving\n", "The peak detection is a standard check that controls the left and right (previous and next) point of the filtered signal and checks if it's higher than both and higher than pthr.\n", "\n", "Every point that passes this condition is saved in a circular structure that will be analyzed during the lookback phase: in this algorithm, the detection of the QRS complex is not immediate but happen after a variable (but bounded) time delay that we will see in the next phase" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Peak lookback\n", "Whenever a peak is detected the algorithm look inside the circular peak buffer analyzing all the peaks that fall between the last QRS complex found until the time of the present peak minus the standard time between QRS and T wave (this exclude the present detected peak from the analysis and the previous QRS complex if we detected a t-wave).\n", "\n", "If the peak analyzed is higher than qthr, is the dominant peak in his direct neighborhood and happened after rrmin from the last QRS complex then the peak is considered to be the next QRS complex, in this case, the algorithm also check the next peak in the buffer to search for a possible T-Wave. All the parameters (estimated intervals and thresholds) get updated in a classical way." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Customization of the gQRS\n", "Let's now discuss the problems of this algorithm and how they were solved in order to adapt it to our environment.\n", "\n", "As discussed in the introduction of this notebook, the idea of this study is to find a method to apply this state-of-the-art algorithm to event-based sampled signals. In particular, the prior that the event-based sampled signal (from now on, the signal) can be expressed as a sequence of piece-wise linear segments between each point without exceding a fixed error, will be extremely helpful." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Smoothing filter\n", "Because of the method used for the linear sampling, no smoothing was required. The only operation performed was just the amplification by $4\\delta$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### The problem with event-based sampling\n", "The original gQRS algorithm relay on a subsequent summation of the output of the matched filter:\n", "\n", "$y[t] = \\sum_{i=4\\delta}^{t} \\sum_{j=-4}^{4}h[j]x[i-j\\delta]$\n", "\n", "where $\\delta = \\frac{QRS duration}{4}$ is the delay of the FIR filter and $h[j]$ is the $j^{th}$ tap of the filter. The Impulsive response of the filter is the following:\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "ExecuteTime": { "end_time": "2020-01-27T15:20:11.201703Z", "start_time": "2020-01-27T15:20:11.078294Z" } }, "outputs": [], "source": [ "taps = [1,-2,-4,10,0,-10,4,2,-1]\n", "plt.plot(taps)\n", "f = np.fft.fftshift(np.fft.fft(taps)) \n", "fig = plt.gcf()\n", "fig.set_size_inches(10.5, 6.5)\n", "plt.grid()\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The filter composed of 9 taps: each of them delayed $\\delta$ from the previous one, covering a time-span of $8\\delta$ around the point under examination" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This requires to know all the points of the original signal: to know the value of the filter at $t$ we must know the output of it at $t-1$. It is possible to show that, if we simply apply the filter to each point of the sub-sampled signal (linearly interpolating the 8 points around the interesting one) the results, on average, diverge.\n", "\n", "One solution could be to simply interpolate linearly all the points between 2 samples and apply the algorithm to the signal obtained in this way. Instead, we want to have fewer possible operations. We will search, hence, to find a method that use only the needed points (so, 8 points for each new sample)\n", "\n", "In order to solve this, we can re-write our filter as:\n", "\n", "$y[t] = \\sum_{j=-4}^{4} h[j] \\sum_{i=4\\delta}^{t} x[i-j\\delta]$\n", "\n", "Now we define $s_{j}[t] = \\sum_{i=4\\delta}^{t} x[i-j\\delta]$\n", "\n", "If we assume to start from $-8\\delta$ with a signal filled with 0, then we can wrtie:\n", "\n", "$s_{j}[t] = \\sum_{i=0}^{t+j\\delta} x[i]$\n", "\n", "Wich is the numerical integration of the signal from 0 to $t+j\\delta$\n", "\n", "Hence: $y[t] = \\sum_{j=-4}^{4} h[j] s_{j}[t]$\n", "\n", "From here we can see that the result of the filtering is a FIR filter applied to the integrated signal.\n", "\n", "The idea from here is to accumulate and save the integral of the signal at each new point that arrives and apply the filter to this saved signal in order to avoid the need to accumulate the output of the filter. This will require to interpolate only 8 points from the integrated signal. Of course we will not be able to perform a linear interpolation of it and neither to find an interpolator formula as we will show shortly. Let's write the generic numerical integration of a digital signal:\n", "\n", "$i[t] = \\sum_{i=0}^{t-1} x[i]$\n", "\n", "(Note: The sum need to reach $t-1$ otherwise $i[t]+i[t:t+k] \\neq i[t+k]$)\n", "\n", "we can now use the prior that our signal is pice-wise linear and find a formula that give us the generic integral value between 2 given points on a linear signal:\n", "\n", "$i[a:b] = \\sum_{i=a}^{b-1} m*i+q = \\frac{m}{2}(b(b+1)-a(a+1))+q(b-a)$\n", "\n", "Now we have a vector of $i[t]$ where the $t_{s}$ are not equaly separeted. Now we only need to find a method to obtain the 8 points (spaced by $\\delta$) around any point of our signal\n", "\n", "Is it possible, given $i[t_{1}]$ and $i[t_{2}]$ to find $i[t_{k}]$ with $t_{1} \n", " # [-4dt,-3dt,-2dt,-1dt,0,+1dt,+2dt,+3dt,+4dt]\n", " #print(pts)\n", " \n", " #FILTER PART (qrst shape)\n", " dv1 = int(pts[-4] - pts[3])#dt1\n", " dv = dv1 << 1\n", " dv -= int(pts[-3]- pts[2])#dt2\n", " dv = dv << 1\n", " dv += dv1\n", " dv -= int(pts[-2]- pts[1])#dt3\n", " dv = dv << 1\n", " dv += int(pts[-1]- pts[0])#dt4\n", " \n", " \n", " self.v1 = dv\n", " \n", " v0 = int(self.v1 / self.c.v1norm)\n", " self.OUTFILTER.append(v0**2)\n", " self.TIMES.append(self.data.t[self.i])\n", " self.qfv_put(self.j, v0 * v0) \n", " \n", " def add_peak(self, peak_time, peak_amp, type):\n", " p = self.current_peak.next_peak\n", " p.time = peak_time\n", " p.amp = peak_amp\n", " p.type = type\n", " self.current_peak = p\n", " p.next_peak.amp = 0 \n", " def peaktype(self, p):\n", " # peaktype() returns 1 if p is the most prominent peak in its neighborhood, 2\n", " # otherwise. The neighborhood consists of all other peaks within rrmin.\n", " # Normally, \"most prominent\" is equivalent to \"largest in amplitude\", but this\n", " # is not always true. For example, consider three consecutive peaks a, b, c\n", " # such that a and b share a neighborhood, b and c share a neighborhood, but a\n", " # and c do not; and suppose that amp(a) > amp(b) > amp(c). In this case, if\n", " # there are no other peaks, a is the most prominent peak in the (a, b)\n", " # neighborhood. Since b is thus identified as a non-prominent peak, c becomes\n", " # the most prominent peak in the (b, c) neighborhood. This is necessary to\n", " # permit detection of low-amplitude beats that closely precede or follow beats\n", " # with large secondary peaks (as, for example, in R-on-T PVCs).\n", " if p.type:\n", " return p.type\n", " else:\n", " a = p.amp\n", " t0 = p.time - self.c.rrmin\n", " t1 = p.time + self.c.rrmin\n", "\n", " if t0 < 0:\n", " t0 = 0\n", "\n", " pp = p.prev_peak\n", " while t0 < pp.time and pp.time < pp.next_peak.time:\n", " if pp.amp == 0:\n", " break\n", " if a < pp.amp and self.peaktype(pp) == 1:\n", " p.type = 2\n", " return p.type\n", " # end:\n", " pp = pp.prev_peak\n", "\n", " pp = p.next_peak\n", " while pp.time < t1 and pp.time > pp.prev_peak.time:\n", " if pp.amp == 0:\n", " break\n", " if a < pp.amp and self.peaktype(pp) == 1:\n", " p.type = 2\n", " return p.type\n", " # end:\n", " pp = pp.next_peak\n", "\n", " p.type = 1\n", " return p.type\n", " def find_missing(self, r, p):\n", " if r is None or p is None:\n", " return None\n", "\n", " minrrerr = p.time - r.time\n", "\n", " s = None\n", " q = r.next_peak\n", " while q.time < p.time:\n", " if self.peaktype(q) == 1:\n", " rrtmp = q.time - r.time\n", " rrerr = rrtmp - self.c.rrmean\n", " if rrerr < 0:\n", " rrerr = -rrerr\n", " if rrerr < minrrerr:\n", " minrrerr = rrerr\n", " s = q\n", " # end:\n", " q = q.next_peak\n", "\n", " return s\n", " \n", " class Conf(object):\n", " \"\"\"\n", " Initial signal configuration object for this qrs detector\n", " \"\"\"\n", " #adc gain is the lvl/mV, can be found in the .hea file, default is 200\n", " #adc zero is in the header file, after the adc resolution (usualy a low number, in our case,11). default = 1024\n", " def __init__(self, fs = 360, adc_gain = 200, hr=75,\n", " RRdelta=0.2, RRmin=0.28, RRmax=2.4,\n", " QS=0.07, QT=0.35,\n", " RTmin=0.25, RTmax=0.33,\n", " QRSa=750, QRSamin=150,\n", " thresh=1.0):\n", " self.fs = fs\n", "\n", " self.sps = int(1 * fs + 0.5)\n", " self.spm = int(60* fs + 0.5)\n", "\n", " self.hr = hr\n", " self.RR = 60.0 / self.hr\n", " self.RRdelta = RRdelta\n", " self.RRmin = RRmin\n", " self.RRmax = RRmax\n", " self.QS = QS\n", " self.QT = QT\n", " self.RTmin = RTmin\n", " self.RTmax = RTmax\n", " self.QRSa = QRSa\n", " self.QRSamin = QRSamin\n", " self.thresh = thresh\n", "\n", " self._NORMAL = 1 # normal beat\n", " self._ARFCT = 16 # isolated QRS-like artifact\n", " self._NOTE = 22 # comment annotation\n", " self._TWAVE = 27 # T-wave peak\n", " self._NPEAKS = 10 # number of peaks buffered (per signal)\n", " self._BUFLN = 32768 # must be a power of 2, see qf()\n", "\n", " self.rrmean = int(self.RR * self.sps)\n", " self.rrdev = int(self.RRdelta * self.sps)\n", " self.rrmin = int(self.RRmin * self.sps)\n", " self.rrmax = int(self.RRmax * self.sps)\n", "\n", " self.rrinc = int(self.rrmean / 40)\n", " if self.rrinc < 1:\n", " self.rrinc = 1\n", "\n", " self.dt = int(self.QS * self.sps / 4)\n", " if self.dt < 1:\n", " raise Exception('Sampling rate is too low. Unable to use signal.')\n", "\n", " self.rtmin = int(self.RTmin * self.sps)\n", " self.rtmean = int(0.75 * self.QT * self.sps)\n", " self.rtmax = int(self.RTmax * self.sps)\n", " \n", " #dv: min peak to peak val for QRS complex in uVolts\n", " #QRSamin: min peak to peak val for QRS complex in samples\n", " dv = adc_gain * self.QRSamin * 0.001\n", " self.pthr = int((self.thresh * dv * dv) / 6)\n", " self.qthr = self.pthr << 1\n", " self.pthmin = self.pthr >> 2\n", " self.qthmin = int((self.pthmin << 2) / 3)\n", " self.tamean = self.qthr # initial value for mean T-wave amplitude\n", "\n", " # Filter constants and thresholds.\n", " self.dt2 = 2 * self.dt\n", " self.dt3 = 3 * self.dt\n", " self.dt4 = 4 * self.dt\n", "\n", " self.smdt = self.dt\n", " self.v1norm = self.smdt * self.dt * 64\n", "\n", " self.smt = 0\n", " self.smt0 = 0 + self.smdt \n", " class Peak(object):\n", " def __init__(self, peak_time, peak_amp, peak_type):\n", " self.time = peak_time\n", " self.amp = peak_amp\n", " self.type = peak_type\n", " self.next_peak = None\n", " self.prev_peak = None\n", " def __str__(self):\n", " return \"Peak object:\\\n", " \\nT: \"+str(self.time)+\\\n", " \"\\nAmp: \"+str(self.amp)+\\\n", " \"\\nTyppe: \"+str(self.type)\n", " class Annotation(object):\n", " def __init__(self, ann_time, ann_type, ann_subtype, ann_num):\n", " self.time = ann_time\n", " self.type = ann_type\n", " self.subtype = ann_subtype\n", " self.num = ann_num\n", " def __str__(self):\n", " return \"Annotation object:\\\n", " \\nT: \"+str(self.time)+\\\n", " \"\\nNum: \"+str(self.num)+\\\n", " \"\\nTyppe: \"+str(self.type) \n", " class DataInteg:\n", " def __init__(self,length):\n", " self.i = collections.deque( maxlen=length)\n", " self.m = collections.deque( maxlen=length)\n", " self.q = collections.deque( maxlen=length)\n", " self.t = collections.deque( maxlen=length)\n", " self.allInteg = []\n", " self.lastVal = None\n", " self.idx = 0\n", " self.len = length\n", " \n", " def put(self,t,v):\n", " presentVal = [t,v]\n", " if self.lastVal == None:\n", " self.lastVal = presentVal\n", " return 0\n", " #print(self.lastVal)\n", " #print(presentVal)\n", " integ,m,q = self.computeIMQ(self.lastVal,presentVal)\n", " self.lastVal = presentVal\n", " #print(integ)\n", " \n", " \n", " if len(self.i)>0:\n", " newInteg = self.i[len(self.i)-1]+integ\n", " else:\n", " newInteg = integ\n", " #print(t,' --> M: ',m,' Q: ',q, 'Integ: ',integ,' total integ: ',newInteg)\n", " self.allInteg.append(newInteg)\n", " self.i.append(newInteg)\n", " self.m.append(m)\n", " self.q.append(q)\n", " self.t.append(t)\n", " return len(self.i)\n", " \n", " \n", " def computeIMQ(self,pt1,pt2):\n", " m = (pt2[1]-pt1[1])/(pt2[0]-pt1[0])\n", " q = pt1[1]-pt1[0]*m\n", " t0 = pt1[0]\n", " t1 = pt2[0]\n", " #integral version not considering the fist point (for accumulation )\n", " integ = m/2*(t1*(t1+1)-t0*(t0+1)) + q*(t1-t0)\n", " #integ = q*(t1-t0+1)+m/2*(t1**2-t0**2+t1+t0)\n", " #print(t0,' --> M: ',m,' Q: ',q)\n", " return integ,m,q\n", " \n", " def reSampler(self, index, dt):\n", " t = self.t[index]\n", " newTimes = [t-4*dt]+[t-3*dt]+[t-2*dt]+[t-1*dt]+[t]+[t+1*dt]+[t+2*dt]+[t+3*dt]+[t+4*dt]\n", " \n", " newSamples = []\n", " firstIdx = None\n", " lastIdx = None\n", " \n", " for t in newTimes:\n", " idx = 0\n", " #can be optimized for t<0\n", " if t <= self.t[-2]:\n", " present = False\n", " for i in range(len(self.t)-1):\n", " if self.t[i] ==t:\n", " idx = i\n", " present = True\n", " break\n", " elif self.t[i+1] ==t:\n", " idx = i+1\n", " present = True\n", " break\n", " elif self.t[i] <=t and self.t[i+1]>t:\n", " idx = i\n", " break\n", " if present:\n", " new = self.i[idx]\n", " else:\n", " lastInteg = self.i[idx]\n", " t0 = self.t[idx]\n", " m = self.m[idx+1]\n", " q = self.q[idx+1]\n", " new = lastInteg+m/2*(t*(t+1)-t0*(t0+1)) + q*(t-t0)\n", " if firstIdx == None:\n", " firstIdx = idx\n", " if t == newTimes[-1]:\n", " lastIdx = idx\n", " else:\n", " if firstIdx == None:\n", " firstIdx = -1\n", " if t == newTimes[-1]:\n", " lastIdx = -1\n", " new = self.i[-1]\n", " t0 = self.t[idx]\n", " m = self.m[idx+1]\n", " q = self.q[idx+1]\n", " #print(t,' ---> ',t0,' -- ',self.t[idx+1],('M: ',m,'Q: ',q))\n", " newSamples.append(int(new))\n", "# print(firstIdx,lastIdx)\n", "# print(\"Saved times: \",list(self.t)[firstIdx:lastIdx+2])\n", "# print(\"Saved samples: \",list(self.i)[firstIdx:lastIdx+2])\n", "# print(\"Found samples: \",newSamples)\n", "# print(\"(for times t = \",newTimes,\")\")\n", " return newSamples\n", " \n", " \n", " #ENTRY POINT, NEED A CONF OBJ FOR INSIDE THE GQRS CLASS\n", " #adc zero is in the header file, after the adc resolution (usualy a low number, in our case,11). default = 1024\n", " #adc_zero in this case is 0 because we already removed the average\n", " \n", " def __init__(self, adc_zero = 1024):\n", " self.c = gQRS.Conf()\n", " \n", " self.annotations = []\n", " self.sample_valid = False\n", " self.adc_zero = adc_zero\n", "\n", " self.qfv = np.zeros((self.c._BUFLN), dtype=\"int64\")\n", " self.smv = np.zeros((self.c._BUFLN), dtype=\"int64\")\n", " self.OUTFILTER = []\n", " self.TIMES = []\n", " self.v1 = 0\n", " \n", " #self.data = self.Data(self.c.dt4*4+6)\n", " self.data = self.DataInteg(self.c.dt4*4+6)\n", " \n", "\n", " t0 = 0\n", " self.t = -self.c.dt4\n", " self.i = 0\n", " self.j = 0\n", "\n", " self.annot = gQRS.Annotation(-self.c.rrmin, \"NOTE\", 0, 0)\n", " # Cicular buffer of Peaks\n", " first_peak = gQRS.Peak(0, 0, 0)\n", " tmp = first_peak\n", " for _ in range(1, self.c._NPEAKS):\n", " tmp.next_peak = gQRS.Peak(0, 0, 0)\n", " tmp.next_peak.prev_peak = tmp\n", " tmp = tmp.next_peak\n", " tmp.next_peak = first_peak\n", " first_peak.prev_peak = tmp\n", " self.current_peak = first_peak\n", " \n", "\n", " self.state = \"RUNNING\"\n", "# self.gqrs(t0, self.tf)\n", "# return self.annotations\n", "\n", " #State variable\n", " self.q0 = None\n", " self.q1 = 0\n", " self.q2 = 0\n", " self.rr = None\n", " self.rrd = None\n", " self.rt = None\n", " self.rtd = None\n", " self.rtdmin = None\n", "\n", " self.p = None # (Peak)\n", " self.q = None # (Peak)\n", " self.r = None # (Peak)\n", " self.tw = None # (Peak)\n", "\n", " self.last_peak = 0\n", " self.last_qrs = 0\n", " self.lastdV = 0\n", " def isPeak(self,t,val,ver=False):\n", " #print('-'*40)\n", " #print(t)\n", " # VALUE TO RETURN\n", " # [integral, square, time, reDo, ok] is the classical return\n", " # for gqrs we wil return \n", " # [type, 0, time, False, ok]\n", " # type = dominant\n", " ok = False\n", " noneFlag = False\n", " rVal = 0\n", " noneVal = 0\n", " tVal = 0\n", " toReturn = []\n", " \n", " \n", " val = (val-1024)*25\n", " dataLen = self.data.put(t,val)\n", " if dataLen >= self.data.len//2:\n", " self.i = dataLen//2\n", " #self.qf()\n", " self.qfInteg()\n", " q0 = self.qfv_at(self.j)\n", " q1 = self.qfv_at(self.j - 1)\n", " q2 = self.qfv_at(self.j - 2)\n", " tPres = self.data.t[self.i-1]\n", " if ver:\n", " print(q0, q1, q2)\n", " pass\n", " \n", " if q1 > self.c.pthr and q2 < q1 and q1 >= q0 :\n", " #print(\"T: \", tPres, \"Filter Values: \", q0, q1, q2,\"\\tThreshold: \",self.c.qthr)\n", " tPeak = self.data.t[self.i-1]\n", " self.add_peak(tPeak, q1, 0)\n", " self.last_peak = tPeak\n", " #print(self.current_peak)\n", " p = self.current_peak.next_peak\n", " if ver:\n", " print(\"\\t\\tPEAK\")\n", " print(\"TIME: \",tPeak)\n", " print(\"THRESHOLD: \",self.c.pthr)\n", " print(\"VALUE: \",q1)\n", " print(\"TIME NEXT: \",p.time ,' Tpres-rtMax: ',tPres - self.c.rtmax)\n", " l = 0\n", " while p.time < tPres - self.c.rtmax:\n", " l+=1\n", " if ver:\n", " #print('*-*-'*10)\n", " #print(p.time,' < ' ,tPres, '-', self.c.rtmax)\n", " pass\n", " \n", " \n", "\n", " \n", " #peaks after the last annotation to the current one (dominant)\n", " if p.time >= self.annot.time + self.c.rrmin and self.peaktype(p) == 1:\n", " \n", " if p.amp > self.c.qthr:\n", " ok = True\n", " if ver:\n", " print(\"p.amp:\",p.amp,\" > self.c.qthr:\",self.c.qthr) \n", " self.rr = p.time - self.annot.time\n", " q = self.find_missing(self.r, p)\n", " if self.rr > self.c.rrmean + 2 * self.c.rrdev and \\\n", " self.rr > 2 * (self.c.rrmean - self.c.rrdev) and \\\n", " q is not None:\n", " p = q\n", " self.rr = p.time - self.annot.time\n", " self.annot.subtype = 1\n", " self.rrd = self.rr - self.c.rrmean\n", " if self.rrd < 0:\n", " self.rrd = -self.rrd\n", " self.c.rrdev += (self.rrd - self.c.rrdev) >> 3\n", " if self.rrd > self.c.rrinc:\n", " self.rrd = self.c.rrinc\n", " if self.rr > self.c.rrmean:\n", " self.c.rrmean += self.rrd\n", " else:\n", " self.c.rrmean -= self.rrd\n", " if p.amp > self.c.qthr * 4:\n", " self.c.qthr += 1\n", " #print(\"QTHR INCREASED: \",self.c.qthr)\n", " elif p.amp < self.c.qthr:\n", " self.c.qthr -= 1\n", " #print(\"QTHR DECREASED: \",self.c.qthr)\n", " if self.c.qthr > self.c.pthr * 20:\n", " self.c.qthr = self.c.pthr * 20\n", " #print(\"QTHR FIXED: \",self.c.qthr)\n", " self.last_qrs = p.time\n", "\n", "##################################################################################################################\n", " #ANNOTATION, HERE THE VALUE TO RETURN\n", " self.annot.time = p.time - self.c.dt2\n", " #self.annot.time = p.time - self.c.dt2\n", " self.annot.type = \"NORMAL\"\n", " qsize = int(p.amp * 10.0 / self.c.qthr)\n", " if qsize > 127:\n", " qsize = 127\n", " self.annot.num = qsize\n", " self.putann(self.annot)\n", " toReturn.append([\"NORMAL\", p.type, self.annot.time, False, ok])\n", " #self.annot.time += self.c.dt2\n", " self.annot.time += self.c.dt2\n", " # [type, 0, time, False, ok]\n", "##################################################################################################################\n", " \n", " #Look for this beat's T-wave\n", " self.tw = None\n", " self.rtdmin = self.c.rtmean\n", " q = p.next_peak\n", " while q.time > self.annot.time:\n", " self.rt = q.time - self.annot.time - self.c.dt2\n", " if self.rt < self.c.rtmin:\n", " # end:\n", " q = q.next_peak\n", " continue\n", " if self.rt > self.c.rtmax:\n", " break\n", " self.rtd = self.rt - self.c.rtmean\n", " if self.rtd < 0:\n", " self.rtd = -self.rtd\n", " if self.rtd < self.rtdmin:\n", " self.rtdmin = self.rtd\n", " self.tw = q\n", " # end:\n", " q = q.next_peak\n", " if self.tw is not None:\n", " tmp_time = self.tw.time - self.c.dt2\n", " tann = gQRS.Annotation(tmp_time, \"TWAVE\",\n", " 1 if tmp_time > self.annot.time + self.c.rtmean else 0,\n", " self.rtdmin)\n", " # if self.state == \"RUNNING\":\n", " # self.putann(tann)\n", " self.rt = tann.time - self.annot.time\n", " self.c.rtmean += (self.rt - self.c.rtmean) >> 4\n", " if self.c.rtmean > self.c.rtmax:\n", " self.c.rtmean = self.c.rtmax\n", " elif self.c.rtmean < self.c.rtmin:\n", " self.c.rtmean = self.c.rrmin\n", " self.tw.type = 2 # mark T-wave as secondary\n", " \n", " self.r = p\n", " self.q = None\n", " self.annot.subtype = 0\n", " elif tPres - self.last_qrs > self.c.rrmax and self.c.qthr > self.c.qthmin:\n", " self.c.qthr -= (self.c.qthr >> 4)\n", " #print(\"QTHR DECREASED: \",self.c.qthr)\n", " # end:\n", " p = p.next_peak\n", " \n", " elif tPres - self.last_peak > self.c.rrmax and self.c.pthr > self.c.pthmin:\n", " self.c.pthr -= (self.c.pthr >> 4)\n", " #print(\"PTHR DECREASED: \",self.c.pthr)\n", " self.j += 1 \n", " else:\n", " self.j = 0\n", " self.i = 0\n", " #print(t)\n", " \n", " if len(toReturn)==0:\n", " toReturn.append([0,noneVal,0,noneFlag,ok])\n", " \n", " return toReturn" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "ExecuteTime": { "end_time": "2020-01-27T15:20:21.025779Z", "start_time": "2020-01-27T15:20:11.273486Z" }, "code_folding": [], "scrolled": true }, "outputs": [], "source": [ "p = gQRS(adc_zero = 1024)\n", "#foundAllTh = []\n", "foundAllTh = dict()\n", "l = 1\n", "padT = 0\n", "padV = 0\n", "tailLen = 100\n", "tai=[]\n", "poi = []\n", "for bits in dataDict:\n", " print('---------',l,\"/\",len(dataDict))\n", " l+=1\n", " foundAll = []\n", " foundAll = dict()\n", " n=1\n", " verbose = False\n", " for name in dataDict[bits]['data']:\n", " print(n,\"/\",len(dataDict[bits]['data']), end = \" \")\n", " startT = time.time()\n", " n+=1\n", " p.__init__()\n", " times,points = dataDict[bits]['data'][name]['t'],dataDict[bits]['data'][name]['v']\n", " found = []\n", " #for tm,pt in tqdm(zip(times,points),total=len(times)):\n", " for tm,pt in zip(times,points):\n", " #print('-----------------------------------')\n", " \n", " a = p.isPeak(tm,pt,ver=verbose)\n", " if verbose:\n", " print(\"*\"*50)\n", " #print(a)\n", " if type(a[0]) != list:\n", " if a[3]:\n", " found[-1] = a\n", " elif a[4] == True:\n", " found.append(a)\n", " else:\n", " for elem in a:\n", " if elem[4] == True:\n", " found.append(elem)\n", " #TAIL\n", " for i in range(1,tailLen):\n", " tm = times[-1]+10*i\n", " pt = points[-1]\n", " tai.append(tm)\n", " poi.append(pt)\n", " a = p.isPeak(tm,pt,ver=verbose)\n", " if verbose:\n", " print(\"*\"*50)\n", " #print(a)\n", " if type(a[0]) != list:\n", " if a[3]:\n", " found[-1] = a\n", " elif a[4] == True:\n", " found.append(a)\n", " else:\n", " for elem in a:\n", " if elem[4] == True:\n", " found.append(elem)\n", " lastT = tai[-1]\n", " for i in range(1,60):\n", " tm = lastT+i*20\n", " pt = points[-1]+50*(-1)**i\n", " tai.append(tm)\n", " poi.append(pt)\n", " a = p.isPeak(tm,pt,ver=verbose)\n", " if verbose:\n", " print(\"*\"*50)\n", " #print(a)\n", " if type(a[0]) != list:\n", " if a[3]:\n", " found[-1] = a\n", " elif a[4] == True:\n", " found.append(a)\n", " else:\n", " for elem in a:\n", " if elem[4] == True:\n", " found.append(elem)\n", " \n", " mita = []\n", " for elem in found:\n", " mita.append(MITA.MITAnnotation(time = elem[2]))\n", " foundAll[name] = {'raw':found,'MITA':mita}\n", " end = time.time()\n", " print(\"---> Elapsed: \",end - startT)\n", " \n", " #foundAllTh.append(foundAll)\n", " foundAllTh[bits] = {'data': foundAll}\n", " foundAllTh[bits]['bits'] = dataDict[bits]['bits']\n", " foundAllTh[bits]['hist'] = dataDict[bits]['hist']\n", " \n", "foundSampleOneComplete = foundAllTh[0]['data'][selectedFileNames[0]]['raw']\n", "print(\"FINISHED\")" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "ExecuteTime": { "end_time": "2020-01-27T15:20:21.029317Z", "start_time": "2020-01-27T15:20:21.027091Z" } }, "outputs": [], "source": [ "len(foundSampleOneComplete)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "ExecuteTime": { "end_time": "2020-01-27T15:20:21.204544Z", "start_time": "2020-01-27T15:20:21.030255Z" } }, "outputs": [], "source": [ "t = dataDict[0]['data'][selectedFileNames[0]]['t'] \n", "v = dataDict[0]['data'][selectedFileNames[0]]['v'] \n", "\n", "tPeak = np.array(np.asarray(foundAllTh[0]['data'][selectedFileNames[0]]['raw'])[:,2], dtype=np.int64)\n", "startIdxMarks = (np.where(tPeak>= start) [0]).min()\n", "stopIdxMarks = (np.where(tPeak <= stop) [0]).max()\n", "foundSample = foundSampleOneComplete[startIdxMarks:stopIdxMarks+1]\n", "\n", "fig = plt.gcf()\n", "fig.set_size_inches(18.5, 10.5)\n", "\n", "plt.plot(t[startI:stopI]/360,v[startI:stopI])\n", "for i in range(len(foundSample)):\n", " plt.axvline(x=foundSample[i][2]/360,color = 'red',linewidth = .5)\n", " plt.axvline(x=foundSample[i][2]/360+0.075,color = 'green',linewidth = .5)\n", " plt.axvline(x=foundSample[i][2]/360-0.075,color = 'green',linewidth = .5)\n", "plt.grid()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Save and evaluate" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "ExecuteTime": { "end_time": "2020-01-27T15:20:21.209311Z", "start_time": "2020-01-27T15:20:21.207291Z" } }, "outputs": [], "source": [ "cwd = os.path.dirname(os.getcwd())\n", "sourceDir = cwd+'/data/dataRaw/'\n", "outDir = cwd+'/data/results/'" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "ExecuteTime": { "end_time": "2020-01-27T15:20:34.199029Z", "start_time": "2020-01-27T15:20:21.211829Z" } }, "outputs": [], "source": [ "# All the Files, independent number of thresholds\n", "l = 0\n", "for th in foundAllTh:\n", " print(\"----------------\",l+1,\"/\",len(foundAllTh))\n", " print(\"Bits: \"+str(foundAllTh[th]['bits'])+\" Hist: \"+str(foundAllTh[th]['hist']))\n", " ext = 'out'+str(foundAllTh[th]['bits'])+'Bits'+str(foundAllTh[th]['hist'])+'Hyst'\n", " evalFile = str(foundAllTh[th]['bits'])+'Bits'+str(foundAllTh[th]['hist'])+'HystEval.txt'\n", " for name in foundAllTh[th]['data']:\n", " pathOut = sourceDir+name+'.'+ext\n", " #print(\"saving \",pathOut)\n", " sample = foundAllTh[th]['data'][name]['MITA']\n", " print(pathOut)\n", " MITA.save_annotations(sample,pathOut)\n", " print(\"---------------- EVALUATING ----------------\")\n", " command = cwd+'/data_parsing/eval_mitdb.sh -f ' + evalFile +' -o ' + ext\n", " command += ' -s '+sourceDir+' -d '+outDir\n", " r = os.system(command)\n", " print('Executed: ',command,'\\nResult: ',r)\n", " l+=1\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "ExecuteTime": { "end_time": "2020-01-27T15:20:34.206160Z", "start_time": "2020-01-27T15:20:34.200364Z" } }, "outputs": [], "source": [ "outFullPath = os.path.join(outDir,evalFile)\n", "pos = 0\n", "fScore = 0\n", "print(\"File\\tSensistivity\\tPositive predictivity\\n\")\n", "with open(outFullPath) as f:\n", " for line in f.readlines():\n", " a = line.split()\n", " #print(str(a)+str(len(a)))\n", " #print(\"#############\")\n", " if pos != 0 and len(a) == 17:\n", " print(\"{}\\t{}\\t\\t{}\".format(a[0],a[12],a[13]))\n", " if len(a)>0:\n", " if a[0] == \"Average\":\n", " print(\"\\n-------------------- AVERAGE --------------------\\n\")\n", " print(\"AVG\\t{}\\t\\t{}\\n\".format(a[1],a[2]))\n", " se = float(a[1])\n", " pp = float(a[2])\n", " fScore = 2*se*pp/(se+pp)\n", " print(\"F1 Score: {}\".format(fScore))\n", " pos += 1" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Re-sample" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "ExecuteTime": { "end_time": "2020-01-27T15:20:34.213387Z", "start_time": "2020-01-27T15:20:34.208003Z" } }, "outputs": [], "source": [ "#RESAMPLER\n", "def getMidPoints(pt1,pt2):\n", " m = (pt2[1]-pt1[1])/(pt2[0]-pt1[0])\n", " q = pt1[1]-pt1[0]*m\n", " t1 = pt1[0]\n", " t2 = pt2[0]\n", " points = []\n", " for t in np.arange(t1,t2):\n", " points.append(m*t+q)\n", " return points" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "ExecuteTime": { "end_time": "2020-01-27T15:20:35.011440Z", "start_time": "2020-01-27T15:20:34.214382Z" } }, "outputs": [], "source": [ "for th in dataDict:\n", " print(\"th\")\n", " dataDict[th]['dataResamp'] = {}\n", " for file in dataDict[th]['data']:\n", " dataFileResamp = []\n", " p = np.array([dataDict[th]['data'][file]['t'],dataDict[th]['data'][file]['v']]).T\n", " for i in range(1,len(p)):\n", " p1 = p[i-1]\n", " p2 = p[i]\n", " res = getMidPoints(p1,p2)\n", " dataFileResamp.extend(res)\n", " if(len(dataFileResamp)<649999):\n", " p1 = p[-1]\n", " p2 = [649999,p1[1]]\n", " res = getMidPoints(p1,p2)\n", " dataFileResamp.extend(res)\n", " dataDict[th]['dataResamp'][file] = {'t':np.array(list(range(len(dataFileResamp)))), \n", " 'v':np.array(dataFileResamp)}\n", " \n", "# dataAllFileResamp.append(np.array([list(range(len(dataFileResamp))),dataFileResamp]))\n", "# dataResam.append(np.array(dataAllFileResamp))\n", "# dataResam = np.array(dataResam)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "ExecuteTime": { "end_time": "2020-01-27T15:20:35.116598Z", "start_time": "2020-01-27T15:20:35.012484Z" } }, "outputs": [], "source": [ "t = dataDict[0]['dataResamp'][selectedFileNames[0]]['t'][0*360:10*360] \n", "v = dataDict[0]['dataResamp'][selectedFileNames[0]]['v'][0*360:10*360]\n", "plt.plot(t,v)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "ExecuteTime": { "end_time": "2020-01-27T15:20:35.171172Z", "start_time": "2020-01-27T15:20:35.118098Z" }, "code_folding": [ 0 ] }, "outputs": [], "source": [ "# ORIGINAL METHOD\n", "def time_to_sample_number(seconds, frequency):\n", " return seconds * frequency + 0.5\n", "def relu(num):\n", " if num<0:\n", " return 0\n", " return num\n", "\n", "class GQRSorig(object):\n", " \"\"\"\n", " GQRS detection class\n", " \"\"\"\n", " class Conf(object):\n", " \"\"\"\n", " Initial signal configuration object for this qrs detector\n", " \"\"\"\n", " #adc gain is the lvl/mV, can be found in the .hea file, default is 200\n", " #adc zero is in the header file, after the adc resolution (usualy a low number, in our case,11). default = 1024\n", " def __init__(self, fs = 360, adc_gain = 200, hr=75,\n", " RRdelta=0.2, RRmin=0.28, RRmax=2.4,\n", " QS=0.07, QT=0.35,\n", " RTmin=0.25, RTmax=0.33,\n", " QRSa=750, QRSamin=150,\n", " thresh=1.0):\n", " self.fs = fs\n", "\n", " self.sps = int(time_to_sample_number(1, fs))\n", " self.spm = int(time_to_sample_number(60, fs))\n", "\n", " self.hr = hr\n", " self.RR = 60.0 / self.hr\n", " self.RRdelta = RRdelta\n", " self.RRmin = RRmin\n", " self.RRmax = RRmax\n", " self.QS = QS\n", " self.QT = QT\n", " self.RTmin = RTmin\n", " self.RTmax = RTmax\n", " self.QRSa = QRSa\n", " self.QRSamin = QRSamin\n", " self.thresh = thresh\n", "\n", " self._NORMAL = 1 # normal beat\n", " self._ARFCT = 16 # isolated QRS-like artifact\n", " self._NOTE = 22 # comment annotation\n", " self._TWAVE = 27 # T-wave peak\n", " self._NPEAKS = 10 # number of peaks buffered (per signal)\n", " self._BUFLN = 32768 # must be a power of 2, see qf()\n", "\n", " self.rrmean = int(self.RR * self.sps)\n", " self.rrdev = int(self.RRdelta * self.sps)\n", " self.rrmin = int(self.RRmin * self.sps)\n", " self.rrmax = int(self.RRmax * self.sps)\n", "\n", " self.rrinc = int(self.rrmean / 40)\n", " if self.rrinc < 1:\n", " self.rrinc = 1\n", "\n", " self.dt = int(self.QS * self.sps / 4)\n", " if self.dt < 1:\n", " raise Exception('Sampling rate is too low. Unable to use signal.')\n", "\n", " self.rtmin = int(self.RTmin * self.sps)\n", " self.rtmean = int(0.75 * self.QT * self.sps)\n", " self.rtmax = int(self.RTmax * self.sps)\n", " \n", " #dv: min peak to peak val for QRS complex in uVolts\n", " #QRSamin: min peak to peak val for QRS complex in samples\n", " dv = adc_gain * self.QRSamin * 0.001\n", " self.pthr = int((self.thresh * dv * dv) / 6)\n", " self.qthr = self.pthr << 1\n", " self.pthmin = self.pthr >> 2\n", " self.qthmin = int((self.pthmin << 2) / 3)\n", " self.tamean = self.qthr # initial value for mean T-wave amplitude\n", "\n", " # Filter constants and thresholds.\n", " self.dt2 = 2 * self.dt\n", " self.dt3 = 3 * self.dt\n", " self.dt4 = 4 * self.dt\n", "\n", " self.smdt = self.dt\n", " self.v1norm = self.smdt * self.dt * 64\n", "\n", " self.smt = 0\n", " self.smt0 = 0 + self.smdt\n", "\n", " class Peak(object):\n", " def __init__(self, peak_time, peak_amp, peak_type):\n", " self.time = peak_time\n", " self.amp = peak_amp\n", " self.type = peak_type\n", " self.next_peak = None\n", " self.prev_peak = None\n", " def __str__(self):\n", " return \"Peak object:\\\n", " \\nT: \"+str(self.time)+\\\n", " \"\\nAmp: \"+str(self.amp)+\\\n", " \"\\nTyppe: \"+str(self.type)\n", "\n", " class Annotation(object):\n", " def __init__(self, ann_time, ann_type, ann_subtype, ann_num):\n", " self.time = ann_time\n", " self.type = ann_type\n", " self.subtype = ann_subtype\n", " self.num = ann_num\n", " def __str__(self):\n", " return \"Annotation object:\\\n", " \\nT: \"+str(self.time)+\\\n", " \"\\nNum: \"+str(self.num)+\\\n", " \"\\nTyppe: \"+str(self.type)\n", "\n", "\n", " def putann(self, annotation):\n", " self.annotations.append(copy.deepcopy(annotation))\n", " \n", " def pad(self,lenFinal,signal):\n", " s = list(copy.deepcopy(signal))\n", " delta = lenFinal-len(s)\n", " if delta>0:\n", " pad = list(np.zeros(delta//2))\n", " if delta%2 == 0:\n", " s = pad+s+pad\n", " else:\n", " s = pad+s+pad+[0.0]\n", " else:\n", " delta = -delta\n", " rem = delta//2\n", " s = s[rem:rem+lenFinal]\n", " return s\n", " \n", "\n", " #ENTRY POINT, NEED A CONF OBJ FOR INSIDE THE GQRS CLASS\n", " #adc zero is in the header file, after the adc resolution (usualy a low number, in our case,11). default = 1024\n", " def detect(self, x, conf, adc_zero = 1024):\n", " \"\"\"\n", " Run detection. x is digital signal\n", " \"\"\"\n", " self.c = conf\n", " self.annotations = []\n", " self.sample_valid = False\n", " self.outTimes = [] # self.outTimes[i] = [0,0,T] for compatibility\n", "\n", " if len(x) < 1:\n", " return []\n", "\n", " self.x = x\n", " self.adc_zero = adc_zero\n", "\n", " self.qfv = np.zeros((self.c._BUFLN), dtype=\"int64\")\n", " self.smv = np.zeros((self.c._BUFLN), dtype=\"int64\")\n", " self.v1 = 0\n", " \n", "\n", " t0 = 0\n", " self.tf = len(x) - 1\n", " self.t = 0 - self.c.dt4\n", "\n", " self.annot = GQRSorig.Annotation(-self.c.rrmin, \"NOTE\", 0, 0)\n", " # Cicular buffer of Peaks\n", " first_peak = GQRSorig.Peak(0, 0, 0)\n", " tmp = first_peak\n", " for _ in range(1, self.c._NPEAKS):\n", " tmp.next_peak = GQRSorig.Peak(0, 0, 0)\n", " tmp.next_peak.prev_peak = tmp\n", " tmp = tmp.next_peak\n", " tmp.next_peak = first_peak\n", " first_peak.prev_peak = tmp\n", " self.current_peak = first_peak\n", "\n", " if self.c.spm > self.c._BUFLN:\n", " if self.tf - t0 > self.c._BUFLN:\n", " tf_learn = t0 + self.c._BUFLN - self.c.dt4\n", " else:\n", " tf_learn = self.tf - self.c.dt4\n", " else:\n", " if self.tf - t0 > self.c.spm:\n", " tf_learn = t0 + self.c.spm - self.c.dt4\n", " else:\n", " tf_learn = self.tf - self.c.dt4\n", "\n", "# self.countdown = -1\n", "# self.state = \"LEARNING\"\n", "# self.gqrs(t0, tf_learn)\n", "\n", " self.rewind_gqrs()\n", "\n", " self.state = \"RUNNING\"\n", " self.t = t0 - self.c.dt4\n", " self.gqrs(t0, self.tf)\n", " return self.outTimes\n", "\n", " def rewind_gqrs(self):\n", " self.countdown = -1\n", " self.at(self.t)\n", " self.annot.time = -self.c.rrmin\n", " self.annot.type = \"NORMAL\"\n", " self.annot.subtype = 0\n", " self.annot.num = 0\n", " p = self.current_peak\n", " for _ in range(self.c._NPEAKS):\n", " p.time = 0\n", " p.type = 0\n", " p.amp = 0\n", " p = p.next_peak\n", "\n", " def at(self, t):\n", " if t < 0:\n", " self.sample_valid = True\n", " return self.x[0]\n", " if t > len(self.x) - 1:\n", " self.sample_valid = False\n", " return self.x[-1]\n", " self.sample_valid = True\n", " return self.x[t]\n", "\n", " #Because _BUFLN is a power of 2 (2**N), _BUFLN-1 is (in binary)\n", " #all 1 from 0 to N-1, hence, self.smv[t & (self.c._BUFLN - 1)] \n", " #is just a fancy pants method for self.smv[t%_BUFLN] \n", " def smv_at(self, t):\n", " return self.smv[t & (self.c._BUFLN - 1)]\n", "\n", " def smv_put(self, t, v):\n", " self.smv[t & (self.c._BUFLN - 1)] = v\n", "\n", " def qfv_at(self, t):\n", " return self.qfv[t & (self.c._BUFLN - 1)]\n", "\n", " def qfv_put(self, t, v):\n", " self.qfv[t & (self.c._BUFLN - 1)] = v\n", " \n", "\n", "\n", " def sm(self, at_t):\n", " # implements a trapezoidal low pass (smoothing) filter\n", " # (with a gain of 4*smdt) applied to input signal sig\n", " # before the QRS matched filter qf().\n", " # Before attempting to 'rewind' by more than BUFLN-smdt\n", " # samples, reset smt and smt0.\n", "\n", " # Calculate samp values from self.smt to at_t.\n", " smt = self.c.smt\n", " smdt = int(self.c.smdt)\n", "\n", " v = 0\n", " while at_t > smt:\n", " smt += 1\n", " # from dt+1 onwards\n", " if smt > int(self.c.smt0):\n", " tmp = int(self.smv_at(smt - 1) + \\\n", " self.at(smt + smdt) + self.at(smt + smdt - 1) - \\\n", " self.at(smt - smdt) - self.at(smt - smdt - 1))\n", "# print('Pos: ',self.at(smt + smdt), self.at(smt + smdt - 1), \\\n", "# 'Neg: ',self.at(smt - smdt), self.at(smt - smdt - 1), \\\n", "# 'Pos-Neg: ',self.at(smt + smdt) + self.at(smt + smdt - 1) - \\\n", "# self.at(smt - smdt) - self.at(smt - smdt - 1))\n", " \n", " self.smv_put(smt, tmp)\n", " self.SIG_SMOOTH.append(tmp)\n", " # from 1 to dt. 0 is never calculated.\n", " else:\n", " v = int(self.at(smt))\n", " for j in range(1, smdt):\n", " smtpj = self.at(smt + j)\n", " smtlj = self.at(smt - j)\n", " v += int(smtpj + smtlj)\n", " self.smv_put(smt, (v << 1) + self.at(smt + j+1) + self.at(smt - j-1) - \\\n", " self.adc_zero * (smdt << 2))\n", "\n", " self.SIG_SMOOTH.append((v << 1) + self.at(smt + j+1) + self.at(smt - j-1) - \\\n", " self.adc_zero * (smdt << 2))\n", " self.c.smt = smt\n", "\n", " return self.smv_at(at_t)\n", "\n", " def qf(self):\n", " # evaluate the QRS detector filter for the next sample\n", "\n", " # do this first, to ensure that all of the other smoothed values needed below are in the buffer\n", " \n", " \n", " dv2 = self.sm(self.t + self.c.dt4)\n", " dv2 -= self.smv_at(self.t - self.c.dt4)\n", " dv1 = int(self.smv_at(self.t + self.c.dt) - self.smv_at(self.t - self.c.dt))\n", " dv = dv1 << 1\n", " dv -= int(self.smv_at(self.t + self.c.dt2) - self.smv_at(self.t - self.c.dt2))\n", " dv = dv << 1\n", " dv += dv1\n", " dv -= int(self.smv_at(self.t + self.c.dt3) - self.smv_at(self.t - self.c.dt3))\n", " dv = dv << 1\n", " dv += dv2\n", " self.v1 += dv\n", " v0 = int(self.v1 / self.c.v1norm)\n", " self.qfv_put(self.t, v0 * v0)\n", " \n", " self.SIG_QRS.append(v0 ** 2)\n", "# #DEBUG \n", "# temp = []\n", "# temp.append(self.smv_at(self.t - self.c.dt4))\n", "# temp.append(self.smv_at(self.t - self.c.dt3))\n", "# temp.append(self.smv_at(self.t - self.c.dt2))\n", "# temp.append(self.smv_at(self.t - self.c.dt))\n", "# temp.append(self.smv_at(self.t))\n", "# temp.append(self.smv_at(self.t + self.c.dt))\n", "# temp.append(self.smv_at(self.t + self.c.dt2))\n", "# temp.append(self.smv_at(self.t + self.c.dt3))\n", "# temp.append(self.smv_at(self.t + self.c.dt4))\n", "# print('T: ',self.t,' >>>>',temp,' ----> ', dv,' (', self.v1,')')\n", "# #DEBUG\n", "\n", "\n", " def gqrs(self, from_sample, to_sample):\n", " q0 = None\n", " q1 = 0\n", " q2 = 0\n", " rr = None\n", " rrd = None\n", " rt = None\n", " rtd = None\n", " rtdmin = None\n", "\n", " p = None # (Peak)\n", " q = None # (Peak)\n", " r = None # (Peak)\n", " tw = None # (Peak)\n", "\n", " last_peak = from_sample\n", " last_qrs = from_sample\n", "\n", " self.SIG_SMOOTH = []\n", " self.SIG_QRS = []\n", "\n", " def add_peak(peak_time, peak_amp, type):\n", " p = self.current_peak.next_peak\n", " p.time = peak_time\n", " p.amp = peak_amp\n", " p.type = type\n", " self.current_peak = p\n", " p.next_peak.amp = 0\n", "\n", " def peaktype(p):\n", " # peaktype() returns 1 if p is the most prominent peak in its neighborhood, 2\n", " # otherwise. The neighborhood consists of all other peaks within rrmin.\n", " # Normally, \"most prominent\" is equivalent to \"largest in amplitude\", but this\n", " # is not always true. For example, consider three consecutive peaks a, b, c\n", " # such that a and b share a neighborhood, b and c share a neighborhood, but a\n", " # and c do not; and suppose that amp(a) > amp(b) > amp(c). In this case, if\n", " # there are no other peaks, a is the most prominent peak in the (a, b)\n", " # neighborhood. Since b is thus identified as a non-prominent peak, c becomes\n", " # the most prominent peak in the (b, c) neighborhood. This is necessary to\n", " # permit detection of low-amplitude beats that closely precede or follow beats\n", " # with large secondary peaks (as, for example, in R-on-T PVCs).\n", " if p.type:\n", " return p.type\n", " else:\n", " a = p.amp\n", " t0 = p.time - self.c.rrmin\n", " t1 = p.time + self.c.rrmin\n", "\n", " if t0 < 0:\n", " t0 = 0\n", "\n", " pp = p.prev_peak\n", " while t0 < pp.time and pp.time < pp.next_peak.time:\n", " if pp.amp == 0:\n", " break\n", " if a < pp.amp and peaktype(pp) == 1:\n", " p.type = 2\n", " return p.type\n", " # end:\n", " pp = pp.prev_peak\n", "\n", " pp = p.next_peak\n", " while pp.time < t1 and pp.time > pp.prev_peak.time:\n", " if pp.amp == 0:\n", " break\n", " if a < pp.amp and peaktype(pp) == 1:\n", " p.type = 2\n", " return p.type\n", " # end:\n", " pp = pp.next_peak\n", "\n", " p.type = 1\n", " return p.type\n", "\n", " def find_missing(r, p):\n", " if r is None or p is None:\n", " return None\n", "\n", " minrrerr = p.time - r.time\n", "\n", " s = None\n", " q = r.next_peak\n", " while q.time < p.time:\n", " if peaktype(q) == 1:\n", " rrtmp = q.time - r.time\n", " rrerr = rrtmp - self.c.rrmean\n", " if rrerr < 0:\n", " rrerr = -rrerr\n", " if rrerr < minrrerr:\n", " minrrerr = rrerr\n", " s = q\n", " # end:\n", " q = q.next_peak\n", "\n", " return s\n", "\n", " r = None\n", " next_minute = 0\n", " minutes = 0\n", " while self.t <= to_sample + self.c.sps:\n", "# print('-'*40)\n", "# print(self.t)\n", " if self.countdown < 0:\n", " if self.sample_valid:\n", " self.qf()\n", " else:\n", " self.countdown = int(time_to_sample_number(1, self.c.fs))\n", " self.state = \"CLEANUP\"\n", " else:\n", " self.countdown -= 1\n", " if self.countdown < 0:\n", " break\n", " #print(self.t,self.countdown,self.sample_valid)\n", "\n", " q0 = self.qfv_at(self.t)\n", " q1 = self.qfv_at(self.t - 1)\n", " q2 = self.qfv_at(self.t - 2)\n", " #print(q0,q1,q2)\n", " # state == RUNNING only\n", " \n", " if q1 > self.c.pthr and q2 < q1 and q1 >= q0 and self.t > self.c.dt4:\n", " \n", " \n", " add_peak(self.t - 1, q1, 0)\n", " last_peak = self.t - 1\n", " p = self.current_peak.next_peak\n", "# if debug:\n", "# print(\"\\t\\tPEAK\")\n", "# print(\"TIME: \",self.t-1)\n", "# print(\"THRESHOLD: \",self.c.pthr)\n", "# print(\"VALUE: \",q1)\n", "# print(\"TIME NEXT: \",p.time ,' LIMIT TIME: ',self.t - self.c.rtmax)\n", " \n", " while p.time < self.t - self.c.rtmax:\n", "# if debug:\n", "# print(p.time,'\\t' ,self.annot.time , self.c.rrmin,'\\t',p.amp,' ---> ',self.c.qthr)\n", " if p.time >= self.annot.time + self.c.rrmin and peaktype(p) == 1:\n", " if p.amp > self.c.qthr:\n", "# if debug:\n", "# print(\"p.amp:\",p.amp,\" > self.c.qthr:\",self.c.qthr) \n", " rr = p.time - relu(self.annot.time)\n", " q = find_missing(r, p)\n", " if rr > self.c.rrmean + 2 * self.c.rrdev and \\\n", " rr > 2 * (self.c.rrmean - self.c.rrdev) and \\\n", " q is not None:\n", " p = q\n", " rr = p.time - self.annot.time\n", " self.annot.subtype = 1\n", " rrd = rr - self.c.rrmean\n", " if rrd < 0:\n", " rrd = -rrd\n", " self.c.rrdev += (rrd - self.c.rrdev) >> 3\n", " if rrd > self.c.rrinc:\n", " rrd = self.c.rrinc\n", " if rr > self.c.rrmean:\n", " self.c.rrmean += rrd\n", " else:\n", " self.c.rrmean -= rrd\n", " if p.amp > self.c.qthr * 4:\n", " self.c.qthr += 1\n", " elif p.amp < self.c.qthr:\n", " self.c.qthr -= 1\n", " if self.c.qthr > self.c.pthr * 20:\n", " self.c.qthr = self.c.pthr * 20\n", " last_qrs = p.time\n", "\n", " if self.state == \"RUNNING\":\n", " self.annot.time = p.time - self.c.dt2\n", " self.annot.type = \"NORMAL\"\n", " qsize = int(p.amp * 10.0 / self.c.qthr)\n", " if qsize > 127:\n", " qsize = 127\n", " self.annot.num = qsize\n", " self.putann(self.annot)\n", "###############################################################################################################\n", " ##OUTPUT HERE##\n", " self.outTimes.append([0,0,self.annot.time])\n", "###############################################################################################################\n", " self.annot.time += self.c.dt2\n", "\n", " #Look for this beat's T-wave\n", " tw = None\n", " rtdmin = self.c.rtmean\n", " q = p.next_peak\n", " while q.time > self.annot.time:\n", " rt = q.time - self.annot.time - self.c.dt2\n", " if rt < self.c.rtmin:\n", " # end:\n", " q = q.next_peak\n", " continue\n", " if rt > self.c.rtmax:\n", " break\n", " rtd = rt - self.c.rtmean\n", " if rtd < 0:\n", " rtd = -rtd\n", " if rtd < rtdmin:\n", " rtdmin = rtd\n", " tw = q\n", " # end:\n", " q = q.next_peak\n", " if tw is not None:\n", " tmp_time = tw.time - self.c.dt2\n", " tann = GQRSorig.Annotation(tmp_time, \"TWAVE\",\n", " 1 if tmp_time > self.annot.time + self.c.rtmean else 0,\n", " rtdmin)\n", " # if self.state == \"RUNNING\":\n", " # self.putann(tann)\n", " rt = tann.time - self.annot.time\n", " self.c.rtmean += (rt - self.c.rtmean) >> 4\n", " if self.c.rtmean > self.c.rtmax:\n", " self.c.rtmean = self.c.rtmax\n", " elif self.c.rtmean < self.c.rtmin:\n", " self.c.rtmean = self.c.rrmin\n", " tw.type = 2 # mark T-wave as secondary\n", " r = p\n", " q = None\n", " self.annot.subtype = 0\n", " elif self.t - last_qrs > self.c.rrmax and self.c.qthr > self.c.qthmin:\n", " self.c.qthr -= (self.c.qthr >> 4)\n", " # end:\n", " p = p.next_peak\n", " \n", " elif self.t - last_peak > self.c.rrmax and self.c.pthr > self.c.pthmin:\n", " self.c.pthr -= (self.c.pthr >> 4)\n", "\n", " self.t += 1\n", " if self.t >= next_minute:\n", " next_minute += self.c.spm\n", " minutes += 1\n", " if minutes >= 60:\n", " minutes = 0\n", "\n", " if self.state == \"LEARNING\":\n", " return\n", "\n", " # Mark the last beat or two.\n", " p = self.current_peak.next_peak\n", " i = 0\n", " while p.time <= p.next_peak.time:\n", " if p.time >= self.annot.time + self.c.rrmin and p.time < self.tf and peaktype(p) == 1:\n", " self.annot.type = \"NORMAL\"\n", " self.annot.time = p.time\n", " self.putann(self.annot)\n", "###############################################################################################################\n", " ##OUTPUT HERE##\n", " self.outTimes.append([0,0,self.annot.time])\n", "###############################################################################################################\n", " # end:\n", " p = p.next_peak\n", " if i>= self.c._NPEAKS:\n", " break\n", " i+=1" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "ExecuteTime": { "end_time": "2020-01-27T15:20:35.180231Z", "start_time": "2020-01-27T15:20:35.172177Z" } }, "outputs": [], "source": [ "foundAllTh = dict()\n", "for bit in dataDict:\n", " foundAllTh[bit] = {'dataResamp':dict()}\n", " foundAllTh[bit]['bits'] = dataDict[bit]['bits']\n", " foundAllTh[bit]['hist'] = dataDict[bit]['hist']" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "ExecuteTime": { "end_time": "2020-01-27T15:20:42.192067Z", "start_time": "2020-01-27T15:20:35.181143Z" }, "scrolled": true }, "outputs": [], "source": [ "l = 1\n", "for bit in dataDict:\n", " print('---------',l,\"/\",len(dataDict))\n", " l+=1\n", " foundAll = []\n", " n=1\n", " verbose = False\n", "\n", " for name in dataDict[bit]['dataResamp']:\n", " startT = time.time()\n", " c = GQRSorig.Conf()\n", " detector = GQRSorig()\n", " print(n,\"/\",len(dataDict[bit]['dataResamp']),end = \" \")\n", " n+=1\n", " times,points = dataDict[bit]['dataResamp'][name]['t'],dataDict[bit]['dataResamp'][name]['v']\n", " found = detector.detect(points,c,1024)\n", " stopT = time.time()\n", " print(\" ---> \",stopT-startT)\n", " \n", " mita = []\n", " for elem in found:\n", " mita.append(MITA.MITAnnotation(time = elem[2]))\n", " foundAllTh[bit]['dataResamp'][name] = {'raw':found,'MITA':mita}" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "ExecuteTime": { "end_time": "2020-01-27T15:20:55.262921Z", "start_time": "2020-01-27T15:20:42.193614Z" } }, "outputs": [], "source": [ "# All the Files, independent number of thresholds\n", "cwd = os.path.dirname(os.getcwd())\n", "sourceDir = cwd+'/data/dataRaw/'\n", "outDir = cwd+'/data/results/'\n", "\n", "l = 0\n", "for th in foundAllTh:\n", " print(\"----------------\",l+1,\"/\",len(foundAllTh))\n", " print(\"Bits: \"+str(foundAllTh[th]['bits'])+\" Hist: \"+str(foundAllTh[th]['hist']))\n", " ext = 'out'+str(foundAllTh[th]['bits'])+'Bits'+str(foundAllTh[th]['hist'])+'Hyst_Orig'\n", " evalFile = str(foundAllTh[th]['bits'])+'Bits'+str(foundAllTh[th]['hist'])+'HystEval_Orig.txt'\n", " for name in foundAllTh[th]['dataResamp']:\n", " pathOut = sourceDir+name+'.'+ext\n", " #print(\"saving \",pathOut)\n", " sample = foundAllTh[th]['dataResamp'][name]['MITA']\n", " print(pathOut)\n", " MITA.save_annotations(sample,pathOut)\n", " print(\"---------------- EVALUATING ----------------\")\n", " command = cwd+'/data_parsing/eval_mitdb.sh -f ' + evalFile +' -o ' + ext\n", " command += ' -s '+sourceDir+' -d '+outDir\n", " r = os.system(command)\n", " print('Executed: ',command,'\\nResult: ',r)\n", " l+=1" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "ExecuteTime": { "end_time": "2020-01-27T15:20:55.268781Z", "start_time": "2020-01-27T15:20:55.264121Z" } }, "outputs": [], "source": [ "outFullPath = os.path.join(outDir,evalFile)\n", "pos = 0\n", "fScore = 0\n", "print(\"File\\tSensistivity\\tPositive predictivity\\n\")\n", "with open(outFullPath) as f:\n", " for line in f.readlines():\n", " a = line.split()\n", " #print(str(a)+str(len(a)))\n", " #print(\"#############\")\n", " if pos != 0 and len(a) == 17:\n", " print(\"{}\\t{}\\t\\t{}\".format(a[0],a[12],a[13]))\n", " if len(a)>0:\n", " if a[0] == \"Average\":\n", " print(\"\\n-------------------- AVERAGE --------------------\\n\")\n", " print(\"AVG\\t{}\\t\\t{}\\n\".format(a[1],a[2]))\n", " se = float(a[1])\n", " pp = float(a[2])\n", " fScore = 2*se*pp/(se+pp)\n", " print(\"F1 Score: {}\".format(fScore))\n", " pos += 1" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.7.1" }, "varInspector": { "cols": { "lenName": 16, "lenType": 16, "lenVar": 40 }, "kernels_config": { "python": { "delete_cmd_postfix": "", "delete_cmd_prefix": "del ", "library": "var_list.py", "varRefreshCmd": "print(var_dic_list())" }, "r": { "delete_cmd_postfix": ") ", "delete_cmd_prefix": "rm(", "library": "var_list.r", "varRefreshCmd": "cat(var_dic_list()) " } }, "types_to_exclude": [ "module", "function", "builtin_function_or_method", "instance", "_Feature" ], "window_display": false } }, "nbformat": 4, "nbformat_minor": 2 }