Skip to content

Clarify #77

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
merged 3 commits into from
Oct 3, 2017
Merged

Clarify #77

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
45b8188
Tidy and clarify several aspects. Rename gqrs variables and other pro…
cx1111 Oct 3, 2017
33bcfc7
fix readme
cx1111 Oct 3, 2017
d7a7cbb
fix examples
cx1111 Oct 3, 2017
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
222 changes: 153 additions & 69 deletions README.rst
View file
Open in desktop

Large diffs are not rendered by default.

138 changes: 121 additions & 17 deletions demo.ipynb
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -12,12 +12,15 @@
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"import wfdb\n",
"import numpy as np\n",
"import os\n",
"import matplotlib.pyplot as plt\n",
"from IPython.display import display"
]
},
Expand Down Expand Up @@ -46,7 +49,9 @@
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"# Demo 1 - Read a wfdb record using the 'rdsamp' function into a wfdb.Record object.\n",
Expand All @@ -64,7 +69,9 @@
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"# Demo 2 - Read certain channels and sections of the WFDB record using the simplified 'srdsamp' function\n",
Expand All @@ -80,7 +87,9 @@
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"# Demo 3 - Read a WFDB header file only (without the signals)\n",
Expand All @@ -94,7 +103,9 @@
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"# Demo 4 - Read part of a WFDB annotation file into a wfdb.Annotation object, and plot the samples\n",
Expand All @@ -109,7 +120,9 @@
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"# Demo 5 - Read a WFDB record and annotation. Plot all channels, and the annotation on top of channel 0.\n",
Expand All @@ -134,7 +147,9 @@
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"# Demo 6 - Read the multi-segment record and plot waveforms from the MIMIC matched waveform database. \n",
Expand All @@ -151,7 +166,9 @@
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"# Demo 7 - Read the multi-segment record and plot waveforms from the MIMIC matched waveform database.\n",
Expand Down Expand Up @@ -186,7 +203,9 @@
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"# Demo 8 - Read a wfdb record in which one channel has multiple samples/frame. Return a smoothed uniform array.\n",
Expand All @@ -197,7 +216,9 @@
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"# Demo 9 - Read a wfdb record in which one channel has multiple samples/frame. Return a list of all the expanded samples.\n",
Expand All @@ -220,7 +241,9 @@
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"# Demo 10 - Read a WFDB record's digital samples and create a copy via the wrsamp() instance method \n",
Expand Down Expand Up @@ -261,7 +284,9 @@
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"# Demo 12 - Write a WFDB record with multiple samples/frame in a channel\n",
Expand All @@ -280,7 +305,9 @@
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"# Demo 13 - Read a WFDB annotation file and create a copy via the wrann() instance method\n",
Expand Down Expand Up @@ -320,7 +347,9 @@
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"# Demo 15 - View what the 'anntype' symbols mean in the standard WFDB library\n",
Expand All @@ -341,7 +370,9 @@
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"# Demo 16 - List the Physiobank Databases\n",
Expand All @@ -353,7 +384,9 @@
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"# Demo 17 - Download all the WFDB records and annotations from a small Physiobank Database\n",
Expand All @@ -375,7 +408,9 @@
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"# Demo 18 - Download specified files from a Physiobank database\n",
Expand All @@ -398,6 +433,75 @@
"display(os.listdir(os.path.join(dldir, 'data')))"
]
},
{
"cell_type": "markdown",
"metadata": {
"collapsed": true
},
"source": [
"## ECG Peak Detection"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"def peaks_hr(x, peak_indices, fs, title, figsize=(20, 10), saveto=None):\n",
" \n",
" # Calculate heart rate\n",
" hrs = wfdb.processing.compute_hr(siglen=x.shape[0], peak_indices=peak_indices, fs=fs)\n",
" \n",
" N = x.shape[0]\n",
" \n",
" fig, ax_left = plt.subplots(figsize=figsize)\n",
" ax_right = ax_left.twinx()\n",
" \n",
" ax_left.plot(x, color='#3979f0', label='Signal')\n",
" ax_left.plot(peak_indices, x[peak_indices], 'rx', marker='x', color='#8b0000', label='Peak', markersize=12)\n",
" ax_right.plot(np.arange(N), hrs, label='Heart rate', color='m', linewidth=2)\n",
"\n",
" ax_left.set_title(title)\n",
"\n",
" ax_left.set_xlabel('Time (ms)')\n",
" ax_left.set_ylabel('ECG (mV)', color='#3979f0')\n",
" ax_right.set_ylabel('Heart rate (bpm)', color='m')\n",
" # Make the y-axis label, ticks and tick labels match the line color.\n",
" ax_left.tick_params('y', colors='#3979f0')\n",
" ax_right.tick_params('y', colors='m')\n",
" if saveto is not None:\n",
" plt.savefig(saveto, dpi=600)\n",
" plt.show()\n",
"\n",
"\n",
"recordname = 'sampledata/100'\n",
"\n",
"def gqrs_plot(recordname, t0=0, tf=10000):\n",
" # Load the wfdb record and the physical samples\n",
" record = wfdb.rdsamp(recordname, sampfrom=t0, sampto=tf, channels=[0])\n",
" \n",
" # Use the gqrs algorithm to find peaks in the first channel\n",
" # The gqrs_detect argument expects a digital signal for the first argument.\n",
" d_signal = record.adc()[:,0]\n",
" peak_indices = wfdb.processing.gqrs_detect(d_signal, fs=record.fs, adcgain=record.adcgain[0], adczero=record.adczero[0], threshold=1.0)\n",
" print('gqrs detected peak indices:', peak_indices)\n",
" peaks_hr(x=record.p_signals, peak_indices=peak_indices, fs=record.fs, title=\"GQRS peak detection on sampledata/100\")\n",
" \n",
" # Correct the peaks by applying constraints\n",
" min_bpm = 20\n",
" max_bpm = 230\n",
" min_gap = record.fs*60/min_bpm\n",
" max_gap = record.fs*60/max_bpm\n",
" peak_indices = wfdb.processing.correct_peaks(d_signal, peak_indices=peak_indices, min_gap=min_gap, max_gap=max_gap, smooth_window=150)\n",
" print('corrected gqrs detected peak indices:', sorted(peak_indices))\n",
" peaks_hr(x=record.p_signals, peak_indices=sorted(peak_indices), fs=record.fs, title=\"Corrected GQRS peak detection on sampledata/100\")\n",
"\n",
"gqrs_plot(recordname)"
]
},
{
"cell_type": "code",
"execution_count": null,
Expand Down
120 changes: 0 additions & 120 deletions examples/Heart Rate.ipynb
View file
Open in desktop

This file was deleted.

6 changes: 3 additions & 3 deletions wfdb/processing/basic.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -6,9 +6,9 @@

def resample_ann(tt, sample):
# tt: numpy.array as returned by signal.resample
# sample: numpy.array containing indexes of annotations (Annotation.sample)
# sample: numpy.array containing indices of annotations (Annotation.sample)

# Compute the new annotation indexes
# Compute the new annotation indices

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

# Resample multiple channels with their annotations

Expand Down
46 changes: 26 additions & 20 deletions wfdb/processing/gqrs.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -453,31 +453,37 @@ def find_missing(r, p):
p = p.next_peak


def gqrs_detect(x, frequency, adcgain, adczero, threshold=1.0,
def gqrs_detect(x, fs, adcgain, adczero, threshold=1.0,
hr=75, RRdelta=0.2, RRmin=0.28, RRmax=2.4,
QS=0.07, QT=0.35, RTmin=0.25, RTmax=0.33,
QRSa=750, QRSamin=130):
"""
A signal with a frequency of only 50Hz is not handled by the original algorithm,
thus it is not recommended to use this algorithm for this case.

* x: The signal as an array
* frequency: The signal frequency
* adcgain: The gain of the signal (the number of adus (q.v.) per physical unit)
* adczero: The value produced by the ADC given a 0 volt input.
* threshold: The threshold for detection
* hr: Typical heart rate, in beats per minute
* RRdelta: Typical difference between successive RR intervals in seconds
* RRmin: Minimum RR interval ("refractory period"), in seconds
* RRmax: Maximum RR interval, in seconds; thresholds will be adjusted if no peaks are detected within this interval
* QS: Typical QRS duration, in seconds
* QT: Typical QT interval, in seconds
* RTmin: Minimum interval between R and T peaks, in seconds
* RTmax: Maximum interval between R and T peaks, in seconds
* QRSa: Typical QRS peak-to-peak amplitude, in microvolts
* QRSamin: Minimum QRS peak-to-peak amplitude, in microvolts
Detect qrs locations in a single channel ecg.

Functionally, a direct port of the gqrs algorithm from the original
wfdb package. Therefore written to accept wfdb record fields.

Input arguments:
- x (required): The digital signal as a numpy array
- fs (required): The sampling frequency of the signal
- adcgain: The gain of the signal (the number of adus (q.v.) per physical unit)
- adczero (required): The value produced by the ADC given a 0 volt input.
- threshold (default=1.0): The threshold for detection
- hr (default=75): Typical heart rate, in beats per minute
- RRdelta (default=0.2): Typical difference between successive RR intervals in seconds
- RRmin (default=0.28): Minimum RR interval ("refractory period"), in seconds
- RRmax (default=2.4): Maximum RR interval, in seconds; thresholds will be adjusted
if no peaks are detected within this interval
- QS (default=0.07): Typical QRS duration, in seconds
- QT (default=0.35): Typical QT interval, in seconds
- RTmin (default=0.25): Minimum interval between R and T peaks, in seconds
- RTmax (default=0.33): Maximum interval between R and T peaks, in seconds
- QRSa (default=750): Typical QRS peak-to-peak amplitude, in microvolts
- QRSamin (default=130): Minimum QRS peak-to-peak amplitude, in microvolts

Note: This function should not be used for signals with fs <= 50Hz
"""
conf = Conf(freq=frequency, gain=adcgain, hr=hr,
conf = Conf(freq=fs, gain=adcgain, hr=hr,
RRdelta=RRdelta, RRmin=RRmin, RRmax=RRmax,
QS=QS, QT=QT,
RTmin=RTmin, RTmax=RTmax,
Expand Down
42 changes: 29 additions & 13 deletions wfdb/processing/hr.py
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,21 +1,37 @@
import numpy
import numpy as np


def compute_hr(length, peaks_indexes, fs):
result = numpy.full(length, numpy.nan, dtype='float32')
def compute_hr(siglen, peak_indices, fs):
"""
Compute instantaneous heart rate from peak indices.

if len(peaks_indexes) < 2:
return result
Usage:
heart_rate = compute_hr(siglen, peak_indices, fs)

current_hr = numpy.nan
Input argumnets:
- siglen (required): The length of the corresponding signal
- peak_indices (required): The peak indices.
- fs (required): The corresponding signal's sampling frequency.

for i in range(0, len(peaks_indexes)-2):
a = peaks_indexes[i]
b = peaks_indexes[i+1]
c = peaks_indexes[i+2]
Output arguments:
- heart_rate: A numpy array of the instantaneous heart rate, with the length
of the corresponding signal. Contains numpy.nan where heart rate could not be computed.
"""
heart_rate = np.full(siglen, np.nan, dtype='float32')

if len(peak_indices) < 2:
return heart_rate

current_hr = np.nan

for i in range(0, len(peak_indices)-2):
a = peak_indices[i]
b = peak_indices[i+1]
c = peak_indices[i+2]
RR = (b-a) * (1.0 / fs) * 1000
hr = 60000.0 / RR
result[b+1:c+1] = hr
result[peaks_indexes[-1]:] = result[peaks_indexes[-1]]
heart_rate[b+1:c+1] = hr

heart_rate[peak_indices[-1]:] = heart_rate[peak_indices[-1]]

return result
return heart_rate
Loading