Skip to content

Commit 356b873

Browse files
cx1111cx1111
committed
fix inifite xqrs backsearch
1 parent 1c9069b commit 356b873

File tree

7 files changed

+78
-106
lines changed

7 files changed

+78
-106
lines changed

‎demo.ipynb

Lines changed: 16 additions & 13 deletions
Original file line numberDiff line numberDiff line change
@@ -470,7 +470,10 @@
470470
},
471471
"outputs": [],
472472
"source": [
473-
"from wfdb import processing"
473+
"import wfdb\n",
474+
"from wfdb import processing\n",
475+
"import matplotlib.pyplot as plt\n",
476+
"import numpy as np"
474477
]
475478
},
476479
{
@@ -481,18 +484,18 @@
481484
"source": [
482485
"# Demo 19 - Use the gqrs detection algorithm.\n",
483486
"\n",
484-
"def peaks_hr(x, peak_indices, fs, title, figsize=(20, 10), saveto=None):\n",
487+
"def peaks_hr(sig, peak_indices, fs, title, figsize=(20, 10), saveto=None):\n",
485488
" \"Plot a signal with its peaks and heart rate\"\n",
486489
" # Calculate heart rate\n",
487-
" hrs = processing.compute_hr(sig_len=x.shape[0], qrs_inds=peak_indices, fs=fs)\n",
490+
" hrs = processing.compute_hr(sig_len=sig.shape[0], qrs_inds=peak_indices, fs=fs)\n",
488491
" \n",
489-
" N = x.shape[0]\n",
492+
" N = sig.shape[0]\n",
490493
" \n",
491494
" fig, ax_left = plt.subplots(figsize=figsize)\n",
492495
" ax_right = ax_left.twinx()\n",
493496
" \n",
494-
" ax_left.plot(x, color='#3979f0', label='Signal')\n",
495-
" ax_left.plot(peak_indices, x[peak_indices], 'rx', marker='x', color='#8b0000', label='Peak', markersize=12)\n",
497+
" ax_left.plot(sig, color='#3979f0', label='Signal')\n",
498+
" ax_left.plot(peak_indices, sig[peak_indices], 'rx', marker='x', color='#8b0000', label='Peak', markersize=12)\n",
496499
" ax_right.plot(np.arange(N), hrs, label='Heart rate', color='m', linewidth=2)\n",
497500
"\n",
498501
" ax_left.set_title(title)\n",
@@ -514,20 +517,20 @@
514517
"qrs_locs = processing.gqrs_detect(sig=record.p_signal[:,0], fs=record.fs)\n",
515518
"\n",
516519
"# Plot results\n",
517-
"peaks_hr(x=record.p_signal, peak_indices=qrs_locs, fs=record.fs,\n",
518-
" title=\"GQRS peak detection on record 100\")\n",
520+
"# peaks_hr(sig=record.p_signal, peak_indices=qrs_locs, fs=record.fs,\n",
521+
"# title=\"GQRS peak detection on record 100\")\n",
519522
" \n",
520523
"# Correct the peaks by applying constraints\n",
521524
"min_bpm = 20\n",
522525
"max_bpm = 230\n",
523-
"min_gap = record.fs*60/min_bpm\n",
524-
"max_gap = record.fs*60/max_bpm\n",
525-
"peak_indices = processing.correct_peaks(record.p_signal[:,0], peak_indices=qrs_locs, min_gap=min_gap,\n",
526-
" max_gap=max_gap, smooth_window=150)\n",
526+
"#in_gap = record.fs * 60 / min_bpm\n",
527+
"max_gap = record.fs * 60 / max_bpm\n",
528+
"peak_indices = processing.correct_peaks(record.p_signal[:,0], peak_inds=qrs_locs,\n",
529+
" max_gap=max_gap, smooth_window_size=150)\n",
527530
"\n",
528531
"# Display results\n",
529532
"print('Corrected gqrs detected peak indices:', sorted(peak_indices))\n",
530-
"peaks_hr(x=record.p_signal, peak_indices=sorted(peak_indices), fs=record.fs,\n",
533+
"peaks_hr(sig=record.p_signal, peak_indices=sorted(peak_indices), fs=record.fs,\n",
531534
" title=\"Corrected GQRS peak detection on sampledata/100\")\n",
532535
" "
533536
]

‎tests/test_processing.py

Lines changed: 8 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -5,7 +5,9 @@
55

66

77
class test_processing():
8-
8+
"""
9+
Test processing functions
10+
"""
911
def test_resample_single(self):
1012
sig, fields = wfdb.rdsamp('sample-data/100')
1113
ann = wfdb.rdann('sample-data/100', 'atr')
@@ -33,7 +35,7 @@ def test_resample_multi(self):
3335
assert new_sig.shape[0] == expected_length
3436
assert new_sig.shape[1] == sig.shape[1]
3537

36-
def test_normalize(self):
38+
def test_normalize_bound(self):
3739
sig, _ = wfdb.rdsamp('sample-data/100')
3840
lb = -5
3941
ub = 15
@@ -80,10 +82,11 @@ def test_correct_peaks(self):
8082
min_bpm = 10
8183
max_bpm = 350
8284
min_gap = fs*60/min_bpm
83-
max_gap = fs*60/max_bpm
85+
max_gap = fs * 60 / max_bpm
8486

85-
y_idxs = processing.correct_peaks(sig[:,0], ann.sample, min_gap,
86-
max_gap, smooth_window=150)
87+
y_idxs = processing.correct_peaks(sig=sig[:,0], peak_inds=ann.sample,
88+
max_gap=max_gap,
89+
smooth_window_size=150)
8790

8891
yz = numpy.zeros(sig.shape[0])
8992
yz[y_idxs] = 1

‎wfdb/io/download.py

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -107,7 +107,7 @@ def get_dbs():
107107

108108
# ---- Helper functions for downloading physiobank files ------- #
109109

110-
def get_record_list(dburl, records):
110+
def get_record_list(dburl, records='all'):
111111
# Check for a RECORDS file
112112
if records == 'all':
113113
r = requests.get(posixpath.join(dburl, 'RECORDS'))

‎wfdb/processing/__init__.py

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -2,7 +2,7 @@
22
The processing subpackage contains signal-processing tools.
33
"""
44
from .basic import (resample_ann, resample_sig, resample_singlechan,
5-
resample_multichan, normalize_bound)
5+
resample_multichan, get_filter_gain)
66
from .evaluate import Comparitor, compare_annotations
77
from .gqrs import gqrs_detect
88
from .hr import compute_hr

‎wfdb/processing/evaluate.py

Lines changed: 1 addition & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -141,8 +141,7 @@ def compare(self):
141141
# Get the next closest sample for this reference sample.
142142
# Can this be empty? Need to catch case where nothing left?
143143
closest_samp_num, smallest_samp_diff = (
144-
self._get_closest_samp_num(ref_samp_num, test_samp_num,
145-
closest_samp_num))
144+
self._get_closest_samp_num(ref_samp_num, test_samp_num))
146145

147146
# If no clash, it is straightforward.
148147

‎wfdb/processing/peaks.py

Lines changed: 49 additions & 81 deletions
Original file line numberDiff line numberDiff line change
@@ -3,52 +3,55 @@
33

44
from .basic import smooth
55

6-
7-
def find_peaks(x):
6+
import pdb
7+
def find_peaks(sig):
88
"""
99
Find hard peaks and soft peaks in a signal, defined as follows:
1010
- Hard peak: a peak that is either /\ or \/
1111
- Soft peak: a peak that is either /-*\ or \-*/ (In that case we define the
1212
middle of it as the peak)
13-
13+
1414
Parameters
1515
----------
16-
x : np array
17-
The signal array
16+
sig : np array
17+
The 1d signal array
1818
1919
Returns
2020
-------
2121
hard_peaks : np array
22-
Array containing the indices of the hard peaks:
22+
Array containing the indices of the hard peaks:
2323
soft_peaks : np array
2424
Array containing the indices of the soft peaks
25-
25+
2626
"""
27-
if len(x) == 0:
27+
if len(sig) == 0:
2828
return np.empty([0]), np.empty([0])
2929

30-
tmp = x[1:]
31-
tmp = np.append(tmp, [x[-1]])
32-
tmp = x-tmp
33-
tmp[np.where(tmp>0)] = +1
30+
tmp = sig[1:]
31+
tmp = np.append(tmp, [sig[-1]])
32+
tmp = sig - tmp
33+
tmp[np.where(tmp>0)] = 1
3434
tmp[np.where(tmp==0)] = 0
3535
tmp[np.where(tmp<0)] = -1
3636
tmp2 = tmp[1:]
3737
tmp2 = np.append(tmp2, [0])
3838
tmp = tmp-tmp2
39-
hard_peaks = np.where(np.logical_or(tmp==-2,tmp==+2))[0]+1
39+
40+
hard_peaks = np.where(np.logical_or(tmp==-2, tmp==+2))[0] + 1
4041
soft_peaks = []
42+
4143
for iv in np.where(np.logical_or(tmp==-1,tmp==+1))[0]:
4244
t = tmp[iv]
4345
i = iv+1
4446
while True:
4547
if i==len(tmp) or tmp[i] == -t or tmp[i] == -2 or tmp[i] == 2:
4648
break
4749
if tmp[i] == t:
48-
soft_peaks.append(int(iv+(i-iv)/2))
50+
soft_peaks.append(int(iv + (i - iv)/2))
4951
break
5052
i += 1
51-
soft_peaks = np.asarray(soft_peaks)+1
53+
soft_peaks = np.array(soft_peaks, dtype='int') + 1
54+
5255
return hard_peaks, soft_peaks
5356

5457

@@ -89,74 +92,39 @@ def find_local_peaks(sig, radius):
8992
return (np.array(peak_inds))
9093

9194

92-
def correct_peaks(x, peak_indices, min_gap, max_gap, smooth_window):
95+
def correct_peaks(sig, old_peak_inds, search_radius, min_gap,
96+
smooth_window_size=None):
9397
"""
98+
Adjust a set of detected peaks to coincide with local signal maxima,
99+
and
100+
101+
Parameters
102+
----------
103+
sig : numpy array
104+
The 1d signal array
105+
peak_inds : np array
106+
Array of the original peak indices
107+
max_gap : int
108+
The radius within which the original peaks may be shifted.
109+
smooth_window_size : int
110+
The window size of the moving average filter to apply on the
111+
signal. The smoothed signal
112+
113+
Returns
114+
-------
115+
corrected_peak_inds : numpy array
116+
Array of the corrected peak indices
117+
94118
"""
95-
N = x.shape[0]
119+
sig_len = sig.shape[0]
120+
n_peaks = len(peak_inds)
96121

97-
rpeaks = np.zeros(N)
98-
rpeaks[peak_indices] = 1.0
122+
# Peak ranges. What for?
123+
peak_ranges = [[peak_inds[i], peak_inds[i+1]] for i in range(n_peaks - 1)]
124+
sig_smoothed = smooth(sig=sig, window_size=smooth_window_size)
99125

100-
rpeaks = rpeaks.astype('int32')
101126

102-
# 1- Extract ranges where we have one or many ones side by side
103-
rpeaks_ranges = []
104-
tmp_idx = 0
105-
for i in range(1, len(rpeaks)):
106-
if rpeaks[i-1] == 1:
107-
if rpeaks[i] == 0:
108-
rpeaks_ranges.append((tmp_idx, i-1))
109-
else:
110-
if rpeaks[i] == 1:
111-
tmp_idx = i
112-
113-
smoothed = smooth(x, smooth_window)
114-
115-
# Compute signal's peaks
116-
hard_peaks, soft_peaks = find_peaks(x=x)
117-
all_peak_idxs = np.concatenate((hard_peaks, soft_peaks)).astype('int64')
118-
119-
# Replace each range of ones by the index of the best value in it
120-
tmp = set()
121-
for rp_range in rpeaks_ranges:
122-
r = np.arange(rp_range[0], rp_range[1]+1, dtype='int64')
123-
vals = x[r]
124-
smoothed_vals = smoothed[r]
125-
p = r[np.argmax(np.absolute(np.asarray(vals)-smoothed_vals))]
126-
tmp.add(p)
127-
128-
# Replace all peaks by the peak within x-max_gap < x < x+max_gap which have
129-
# the bigget distance from smooth curve
130-
dist = np.absolute(x-smoothed) # Peak distance from the smoothed mean
131-
rpeak_indices = set()
132-
for p in tmp:
133-
a = max(0, p-max_gap)
134-
b = min(N, p+max_gap)
135-
r = np.arange(a, b, dtype='int64')
136-
idx_best = r[np.argmax(dist[r])]
137-
rpeak_indices.add(idx_best)
138-
139-
rpeak_indices = list(rpeak_indices)
140-
141-
# Prevent multiple peaks to appear in the max bpm range (max_gap)
142-
# If we found more than one peak in this interval, then we choose the peak
143-
# with the maximum amplitude compared to the mean of the signal
144-
tmp = np.asarray(rpeak_indices)
145-
to_remove = {}
146-
for idx in rpeak_indices:
147-
if idx in to_remove:
148-
continue
149-
r = tmp[np.where(np.absolute(tmp-idx)<=max_gap)[0]]
150-
if len(r) == 1:
151-
continue
152-
rr = r.astype('int64')
153-
vals = x[rr]
154-
smoo = smoothed[rr]
155-
the_one = r[np.argmax(np.absolute(vals-smoo))]
156-
for i in r:
157-
if i != the_one:
158-
to_remove[i] = True
159-
for v, _ in to_remove.items():
160-
rpeak_indices.remove(v)
161-
162-
return sorted(rpeak_indices)
127+
128+
129+
return corrected_peak_inds
130+

‎wfdb/processing/qrs.py

Lines changed: 2 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -473,10 +473,9 @@ def _backsearch(self):
473473
"""
474474
for peak_num in range(self.last_qrs_peak_num + 1, self.peak_num + 1):
475475
if self._is_qrs(peak_num=peak_num, backsearch=True):
476-
self._update_qrs_params(peak_num=peak_num, backsearch=True)
476+
self._update_qrs(peak_num=peak_num, backsearch=True)
477477
# No need to update noise parameters if it was classified as
478478
# noise. It would have already been updated.
479-
return
480479

481480
def _run_detection(self):
482481
"""
@@ -501,7 +500,7 @@ def _run_detection(self):
501500

502501
# Before continuing to the next peak, do backsearch if
503502
# necessary
504-
while self._require_backsearch():
503+
if self._require_backsearch():
505504
self._backsearch()
506505

507506
# Detected indices are relative to starting sample

0 commit comments

Comments
 (0)