Skip to content

Commit d58286c

Browse files
cx1111cx1111
authored
Merge pull request MIT-LCP#77 from MIT-LCP/clarify
Clarify
2 parents a01250b + d7a7cbb commit d58286c

File tree

9 files changed

+459
-286
lines changed

9 files changed

+459
-286
lines changed

‎README.rst

Lines changed: 153 additions & 69 deletions
Large diffs are not rendered by default.

‎demo.ipynb

Lines changed: 121 additions & 17 deletions
Original file line numberDiff line numberDiff line change
@@ -12,12 +12,15 @@
1212
{
1313
"cell_type": "code",
1414
"execution_count": null,
15-
"metadata": {},
15+
"metadata": {
16+
"collapsed": true
17+
},
1618
"outputs": [],
1719
"source": [
1820
"import wfdb\n",
1921
"import numpy as np\n",
2022
"import os\n",
23+
"import matplotlib.pyplot as plt\n",
2124
"from IPython.display import display"
2225
]
2326
},
@@ -46,7 +49,9 @@
4649
{
4750
"cell_type": "code",
4851
"execution_count": null,
49-
"metadata": {},
52+
"metadata": {
53+
"collapsed": true
54+
},
5055
"outputs": [],
5156
"source": [
5257
"# Demo 1 - Read a wfdb record using the 'rdsamp' function into a wfdb.Record object.\n",
@@ -64,7 +69,9 @@
6469
{
6570
"cell_type": "code",
6671
"execution_count": null,
67-
"metadata": {},
72+
"metadata": {
73+
"collapsed": true
74+
},
6875
"outputs": [],
6976
"source": [
7077
"# Demo 2 - Read certain channels and sections of the WFDB record using the simplified 'srdsamp' function\n",
@@ -80,7 +87,9 @@
8087
{
8188
"cell_type": "code",
8289
"execution_count": null,
83-
"metadata": {},
90+
"metadata": {
91+
"collapsed": true
92+
},
8493
"outputs": [],
8594
"source": [
8695
"# Demo 3 - Read a WFDB header file only (without the signals)\n",
@@ -94,7 +103,9 @@
94103
{
95104
"cell_type": "code",
96105
"execution_count": null,
97-
"metadata": {},
106+
"metadata": {
107+
"collapsed": true
108+
},
98109
"outputs": [],
99110
"source": [
100111
"# Demo 4 - Read part of a WFDB annotation file into a wfdb.Annotation object, and plot the samples\n",
@@ -109,7 +120,9 @@
109120
{
110121
"cell_type": "code",
111122
"execution_count": null,
112-
"metadata": {},
123+
"metadata": {
124+
"collapsed": true
125+
},
113126
"outputs": [],
114127
"source": [
115128
"# Demo 5 - Read a WFDB record and annotation. Plot all channels, and the annotation on top of channel 0.\n",
@@ -134,7 +147,9 @@
134147
{
135148
"cell_type": "code",
136149
"execution_count": null,
137-
"metadata": {},
150+
"metadata": {
151+
"collapsed": true
152+
},
138153
"outputs": [],
139154
"source": [
140155
"# Demo 6 - Read the multi-segment record and plot waveforms from the MIMIC matched waveform database. \n",
@@ -151,7 +166,9 @@
151166
{
152167
"cell_type": "code",
153168
"execution_count": null,
154-
"metadata": {},
169+
"metadata": {
170+
"collapsed": true
171+
},
155172
"outputs": [],
156173
"source": [
157174
"# Demo 7 - Read the multi-segment record and plot waveforms from the MIMIC matched waveform database.\n",
@@ -186,7 +203,9 @@
186203
{
187204
"cell_type": "code",
188205
"execution_count": null,
189-
"metadata": {},
206+
"metadata": {
207+
"collapsed": true
208+
},
190209
"outputs": [],
191210
"source": [
192211
"# Demo 8 - Read a wfdb record in which one channel has multiple samples/frame. Return a smoothed uniform array.\n",
@@ -197,7 +216,9 @@
197216
{
198217
"cell_type": "code",
199218
"execution_count": null,
200-
"metadata": {},
219+
"metadata": {
220+
"collapsed": true
221+
},
201222
"outputs": [],
202223
"source": [
203224
"# Demo 9 - Read a wfdb record in which one channel has multiple samples/frame. Return a list of all the expanded samples.\n",
@@ -220,7 +241,9 @@
220241
{
221242
"cell_type": "code",
222243
"execution_count": null,
223-
"metadata": {},
244+
"metadata": {
245+
"collapsed": true
246+
},
224247
"outputs": [],
225248
"source": [
226249
"# Demo 10 - Read a WFDB record's digital samples and create a copy via the wrsamp() instance method \n",
@@ -261,7 +284,9 @@
261284
{
262285
"cell_type": "code",
263286
"execution_count": null,
264-
"metadata": {},
287+
"metadata": {
288+
"collapsed": true
289+
},
265290
"outputs": [],
266291
"source": [
267292
"# Demo 12 - Write a WFDB record with multiple samples/frame in a channel\n",
@@ -280,7 +305,9 @@
280305
{
281306
"cell_type": "code",
282307
"execution_count": null,
283-
"metadata": {},
308+
"metadata": {
309+
"collapsed": true
310+
},
284311
"outputs": [],
285312
"source": [
286313
"# Demo 13 - Read a WFDB annotation file and create a copy via the wrann() instance method\n",
@@ -320,7 +347,9 @@
320347
{
321348
"cell_type": "code",
322349
"execution_count": null,
323-
"metadata": {},
350+
"metadata": {
351+
"collapsed": true
352+
},
324353
"outputs": [],
325354
"source": [
326355
"# Demo 15 - View what the 'anntype' symbols mean in the standard WFDB library\n",
@@ -341,7 +370,9 @@
341370
{
342371
"cell_type": "code",
343372
"execution_count": null,
344-
"metadata": {},
373+
"metadata": {
374+
"collapsed": true
375+
},
345376
"outputs": [],
346377
"source": [
347378
"# Demo 16 - List the Physiobank Databases\n",
@@ -353,7 +384,9 @@
353384
{
354385
"cell_type": "code",
355386
"execution_count": null,
356-
"metadata": {},
387+
"metadata": {
388+
"collapsed": true
389+
},
357390
"outputs": [],
358391
"source": [
359392
"# Demo 17 - Download all the WFDB records and annotations from a small Physiobank Database\n",
@@ -375,7 +408,9 @@
375408
{
376409
"cell_type": "code",
377410
"execution_count": null,
378-
"metadata": {},
411+
"metadata": {
412+
"collapsed": true
413+
},
379414
"outputs": [],
380415
"source": [
381416
"# Demo 18 - Download specified files from a Physiobank database\n",
@@ -398,6 +433,75 @@
398433
"display(os.listdir(os.path.join(dldir, 'data')))"
399434
]
400435
},
436+
{
437+
"cell_type": "markdown",
438+
"metadata": {
439+
"collapsed": true
440+
},
441+
"source": [
442+
"## ECG Peak Detection"
443+
]
444+
},
445+
{
446+
"cell_type": "code",
447+
"execution_count": null,
448+
"metadata": {
449+
"collapsed": true
450+
},
451+
"outputs": [],
452+
"source": [
453+
"def peaks_hr(x, peak_indices, fs, title, figsize=(20, 10), saveto=None):\n",
454+
" \n",
455+
" # Calculate heart rate\n",
456+
" hrs = wfdb.processing.compute_hr(siglen=x.shape[0], peak_indices=peak_indices, fs=fs)\n",
457+
" \n",
458+
" N = x.shape[0]\n",
459+
" \n",
460+
" fig, ax_left = plt.subplots(figsize=figsize)\n",
461+
" ax_right = ax_left.twinx()\n",
462+
" \n",
463+
" ax_left.plot(x, color='#3979f0', label='Signal')\n",
464+
" ax_left.plot(peak_indices, x[peak_indices], 'rx', marker='x', color='#8b0000', label='Peak', markersize=12)\n",
465+
" ax_right.plot(np.arange(N), hrs, label='Heart rate', color='m', linewidth=2)\n",
466+
"\n",
467+
" ax_left.set_title(title)\n",
468+
"\n",
469+
" ax_left.set_xlabel('Time (ms)')\n",
470+
" ax_left.set_ylabel('ECG (mV)', color='#3979f0')\n",
471+
" ax_right.set_ylabel('Heart rate (bpm)', color='m')\n",
472+
" # Make the y-axis label, ticks and tick labels match the line color.\n",
473+
" ax_left.tick_params('y', colors='#3979f0')\n",
474+
" ax_right.tick_params('y', colors='m')\n",
475+
" if saveto is not None:\n",
476+
" plt.savefig(saveto, dpi=600)\n",
477+
" plt.show()\n",
478+
"\n",
479+
"\n",
480+
"recordname = 'sampledata/100'\n",
481+
"\n",
482+
"def gqrs_plot(recordname, t0=0, tf=10000):\n",
483+
" # Load the wfdb record and the physical samples\n",
484+
" record = wfdb.rdsamp(recordname, sampfrom=t0, sampto=tf, channels=[0])\n",
485+
" \n",
486+
" # Use the gqrs algorithm to find peaks in the first channel\n",
487+
" # The gqrs_detect argument expects a digital signal for the first argument.\n",
488+
" d_signal = record.adc()[:,0]\n",
489+
" peak_indices = wfdb.processing.gqrs_detect(d_signal, fs=record.fs, adcgain=record.adcgain[0], adczero=record.adczero[0], threshold=1.0)\n",
490+
" print('gqrs detected peak indices:', peak_indices)\n",
491+
" peaks_hr(x=record.p_signals, peak_indices=peak_indices, fs=record.fs, title=\"GQRS peak detection on sampledata/100\")\n",
492+
" \n",
493+
" # Correct the peaks by applying constraints\n",
494+
" min_bpm = 20\n",
495+
" max_bpm = 230\n",
496+
" min_gap = record.fs*60/min_bpm\n",
497+
" max_gap = record.fs*60/max_bpm\n",
498+
" peak_indices = wfdb.processing.correct_peaks(d_signal, peak_indices=peak_indices, min_gap=min_gap, max_gap=max_gap, smooth_window=150)\n",
499+
" print('corrected gqrs detected peak indices:', sorted(peak_indices))\n",
500+
" peaks_hr(x=record.p_signals, peak_indices=sorted(peak_indices), fs=record.fs, title=\"Corrected GQRS peak detection on sampledata/100\")\n",
501+
"\n",
502+
"gqrs_plot(recordname)"
503+
]
504+
},
401505
{
402506
"cell_type": "code",
403507
"execution_count": null,

‎examples/Heart Rate.ipynb

Lines changed: 0 additions & 120 deletions
This file was deleted.

‎wfdb/processing/basic.py

Lines changed: 3 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -6,9 +6,9 @@
66

77
def resample_ann(tt, sample):
88
# tt: numpy.array as returned by signal.resample
9-
# sample: numpy.array containing indexes of annotations (Annotation.sample)
9+
# sample: numpy.array containing indices of annotations (Annotation.sample)
1010

11-
# Compute the new annotation indexes
11+
# Compute the new annotation indices
1212

1313
tmp = numpy.zeros(len(tt), dtype='int16')
1414
j = 0
@@ -86,7 +86,7 @@ def resample_multichan(xs, ann, fs, fs_target, resamp_ann_chan=0):
8686
# ann: an Annotation object
8787
# fs: the current frequency
8888
# fs_target: the target frequency
89-
# resample_ann_channel: the signal channel that is used to compute new annotation indexes
89+
# resample_ann_channel: the signal channel that is used to compute new annotation indices
9090

9191
# Resample multiple channels with their annotations
9292

‎wfdb/processing/gqrs.py

Lines changed: 26 additions & 20 deletions
Original file line numberDiff line numberDiff line change
@@ -453,31 +453,37 @@ def find_missing(r, p):
453453
p = p.next_peak
454454

455455

456-
def gqrs_detect(x, frequency, adcgain, adczero, threshold=1.0,
456+
def gqrs_detect(x, fs, adcgain, adczero, threshold=1.0,
457457
hr=75, RRdelta=0.2, RRmin=0.28, RRmax=2.4,
458458
QS=0.07, QT=0.35, RTmin=0.25, RTmax=0.33,
459459
QRSa=750, QRSamin=130):
460460
"""
461-
A signal with a frequency of only 50Hz is not handled by the original algorithm,
462-
thus it is not recommended to use this algorithm for this case.
463-
464-
* x: The signal as an array
465-
* frequency: The signal frequency
466-
* adcgain: The gain of the signal (the number of adus (q.v.) per physical unit)
467-
* adczero: The value produced by the ADC given a 0 volt input.
468-
* threshold: The threshold for detection
469-
* hr: Typical heart rate, in beats per minute
470-
* RRdelta: Typical difference between successive RR intervals in seconds
471-
* RRmin: Minimum RR interval ("refractory period"), in seconds
472-
* RRmax: Maximum RR interval, in seconds; thresholds will be adjusted if no peaks are detected within this interval
473-
* QS: Typical QRS duration, in seconds
474-
* QT: Typical QT interval, in seconds
475-
* RTmin: Minimum interval between R and T peaks, in seconds
476-
* RTmax: Maximum interval between R and T peaks, in seconds
477-
* QRSa: Typical QRS peak-to-peak amplitude, in microvolts
478-
* QRSamin: Minimum QRS peak-to-peak amplitude, in microvolts
461+
Detect qrs locations in a single channel ecg.
462+
463+
Functionally, a direct port of the gqrs algorithm from the original
464+
wfdb package. Therefore written to accept wfdb record fields.
465+
466+
Input arguments:
467+
- x (required): The digital signal as a numpy array
468+
- fs (required): The sampling frequency of the signal
469+
- adcgain: The gain of the signal (the number of adus (q.v.) per physical unit)
470+
- adczero (required): The value produced by the ADC given a 0 volt input.
471+
- threshold (default=1.0): The threshold for detection
472+
- hr (default=75): Typical heart rate, in beats per minute
473+
- RRdelta (default=0.2): Typical difference between successive RR intervals in seconds
474+
- RRmin (default=0.28): Minimum RR interval ("refractory period"), in seconds
475+
- RRmax (default=2.4): Maximum RR interval, in seconds; thresholds will be adjusted
476+
if no peaks are detected within this interval
477+
- QS (default=0.07): Typical QRS duration, in seconds
478+
- QT (default=0.35): Typical QT interval, in seconds
479+
- RTmin (default=0.25): Minimum interval between R and T peaks, in seconds
480+
- RTmax (default=0.33): Maximum interval between R and T peaks, in seconds
481+
- QRSa (default=750): Typical QRS peak-to-peak amplitude, in microvolts
482+
- QRSamin (default=130): Minimum QRS peak-to-peak amplitude, in microvolts
483+
484+
Note: This function should not be used for signals with fs <= 50Hz
479485
"""
480-
conf = Conf(freq=frequency, gain=adcgain, hr=hr,
486+
conf = Conf(freq=fs, gain=adcgain, hr=hr,
481487
RRdelta=RRdelta, RRmin=RRmin, RRmax=RRmax,
482488
QS=QS, QT=QT,
483489
RTmin=RTmin, RTmax=RTmax,

‎wfdb/processing/hr.py

Lines changed: 29 additions & 13 deletions
Original file line numberDiff line numberDiff line change
@@ -1,21 +1,37 @@
1-
import numpy
1+
import numpy as np
22

33

4-
def compute_hr(length, peaks_indexes, fs):
5-
result = numpy.full(length, numpy.nan, dtype='float32')
4+
def compute_hr(siglen, peak_indices, fs):
5+
"""
6+
Compute instantaneous heart rate from peak indices.
67
7-
if len(peaks_indexes) < 2:
8-
return result
8+
Usage:
9+
heart_rate = compute_hr(siglen, peak_indices, fs)
910
10-
current_hr = numpy.nan
11+
Input argumnets:
12+
- siglen (required): The length of the corresponding signal
13+
- peak_indices (required): The peak indices.
14+
- fs (required): The corresponding signal's sampling frequency.
1115
12-
for i in range(0, len(peaks_indexes)-2):
13-
a = peaks_indexes[i]
14-
b = peaks_indexes[i+1]
15-
c = peaks_indexes[i+2]
16+
Output arguments:
17+
- heart_rate: A numpy array of the instantaneous heart rate, with the length
18+
of the corresponding signal. Contains numpy.nan where heart rate could not be computed.
19+
"""
20+
heart_rate = np.full(siglen, np.nan, dtype='float32')
21+
22+
if len(peak_indices) < 2:
23+
return heart_rate
24+
25+
current_hr = np.nan
26+
27+
for i in range(0, len(peak_indices)-2):
28+
a = peak_indices[i]
29+
b = peak_indices[i+1]
30+
c = peak_indices[i+2]
1631
RR = (b-a) * (1.0 / fs) * 1000
1732
hr = 60000.0 / RR
18-
result[b+1:c+1] = hr
19-
result[peaks_indexes[-1]:] = result[peaks_indexes[-1]]
33+
heart_rate[b+1:c+1] = hr
34+
35+
heart_rate[peak_indices[-1]:] = heart_rate[peak_indices[-1]]
2036

21-
return result
37+
return heart_rate

0 commit comments

Comments
 (0)