Skip to content

Commit 45b8188

Browse files
cx1111cx1111
committed
Tidy and clarify several aspects. Rename gqrs variables and other processing variables, move and change heart rate and gqrs demo. Update function documentation. Correct all cases of dac and adc. Update readme to capture new changes.
1 parent a01250b commit 45b8188

File tree

9 files changed

+443
-264
lines changed

9 files changed

+443
-264
lines changed

‎README.rst

Lines changed: 139 additions & 47 deletions
Large diffs are not rendered by default.

‎demo.ipynb

Lines changed: 119 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,73 @@
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+
"outputs": [],
450+
"source": [
451+
"def peaks_hr(x, peak_indices, fs, title, figsize=(20, 10), saveto=None):\n",
452+
" \n",
453+
" # Calculate heart rate\n",
454+
" hrs = wfdb.processing.compute_hr(siglen=x.shape[0], peak_indices=peak_indices, fs=fs)\n",
455+
" \n",
456+
" N = x.shape[0]\n",
457+
" \n",
458+
" fig, ax_left = plt.subplots(figsize=figsize)\n",
459+
" ax_right = ax_left.twinx()\n",
460+
" \n",
461+
" ax_left.plot(x, color='#3979f0', label='Signal')\n",
462+
" ax_left.plot(peak_indices, x[peak_indices], 'rx', marker='x', color='#8b0000', label='Peak', markersize=12)\n",
463+
" ax_right.plot(np.arange(N), hrs, label='Heart rate', color='m', linewidth=2)\n",
464+
"\n",
465+
" ax_left.set_title(title)\n",
466+
"\n",
467+
" ax_left.set_xlabel('Time (ms)')\n",
468+
" ax_left.set_ylabel('ECG (mV)', color='#3979f0')\n",
469+
" ax_right.set_ylabel('Heart rate (bpm)', color='m')\n",
470+
" # Make the y-axis label, ticks and tick labels match the line color.\n",
471+
" ax_left.tick_params('y', colors='#3979f0')\n",
472+
" ax_right.tick_params('y', colors='m')\n",
473+
" if saveto is not None:\n",
474+
" plt.savefig(saveto, dpi=600)\n",
475+
" plt.show()\n",
476+
"\n",
477+
"\n",
478+
"recordname = 'sampledata/100'\n",
479+
"\n",
480+
"def gqrs_plot(recordname, t0=0, tf=10000):\n",
481+
" # Load the wfdb record and the physical samples\n",
482+
" record = wfdb.rdsamp(recordname, sampfrom=t0, sampto=tf, channels=[0])\n",
483+
" \n",
484+
" # Use the gqrs algorithm to find peaks in the first channel\n",
485+
" # The gqrs_detect argument expects a digital signal for the first argument.\n",
486+
" d_signal = record.adc()[:,0]\n",
487+
" peak_indices = wfdb.processing.gqrs_detect(d_signal, fs=record.fs, adcgain=record.adcgain[0], adczero=record.adczero[0], threshold=1.0)\n",
488+
" print('gqrs detected peak indices:', peak_indices)\n",
489+
" peaks_hr(x=record.p_signals, peak_indices=peak_indices, fs=record.fs, title=\"GQRS peak detection on sampledata/100\")\n",
490+
" \n",
491+
" # Do peak detection with constraints\n",
492+
" min_bpm = 20\n",
493+
" max_bpm = 230\n",
494+
" min_gap = record.fs*60/min_bpm\n",
495+
" max_gap = record.fs*60/max_bpm\n",
496+
" peak_indices = wfdb.processing.correct_peaks(d_signal, peak_indices=peak_indices, min_gap=min_gap, max_gap=max_gap, smooth_window=150)\n",
497+
" print('corrected gqrs detected peak indices:', sorted(peak_indices))\n",
498+
" peaks_hr(x=record.p_signals, peak_indices=sorted(peak_indices), fs=record.fs, title=\"Corrected GQRS peak detection on sampledata/100\")\n",
499+
"\n",
500+
"gqrs_plot(recordname)"
501+
]
502+
},
401503
{
402504
"cell_type": "code",
403505
"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)