# set cell width to 100%
from IPython.core.display import display, HTML
display(HTML("<style>.container { width:100% !important; }</style>"))
Pobranie plików¶
Wybieramy bazę MIT-svbd Pobieram z niej pliki oraz je przetwarzam.
Pobieranie plików za pomocą biblioteki urllib¶
import urllib.request
# 850
url = 'https://physionet.org/physiobank/database/svdb/850.hea'
urllib.request.urlretrieve(url, '850.hea')
url = 'https://physionet.org/physiobank/database/svdb/850.dat'
urllib.request.urlretrieve(url, '850.dat')
url = 'https://physionet.org/physiobank/database/svdb/850.atr'
urllib.request.urlretrieve(url, '850.atr')
# 868
url = 'https://physionet.org/physiobank/database/svdb/868.hea'
urllib.request.urlretrieve(url, '868.hea')
url = 'https://physionet.org/physiobank/database/svdb/868.dat'
urllib.request.urlretrieve(url, '868.dat')
url = 'https://physionet.org/physiobank/database/svdb/868.atr'
urllib.request.urlretrieve(url, '868.atr')
Pliki zapisów sygnałów EKG zostały pobrane.
Przetwarzanie¶
import wfdb
import pywt
import matplotlib.pyplot as plt
import numpy as np
Wczytanie rekordu¶
%matplotlib notebook
record = wfdb.rdrecord('868')
# record = wfdb.rdrecord('850', sampfrom=800, channels=[1, 1])
display(record.__dict__)
Zapisanie rekordów do zmiennych.
signals_850, fields_850 = wfdb.rdsamp('850', sampfrom=0)
signals_868, fields_868 = wfdb.rdsamp('868', sampfrom=0)
Wyświetlenie części sygnału w celu przedstawienia poprawnego załadowania sygnału.
cutoff = 1000
plt.figure(1, figsize=(10, 5))
plt.plot(signals_850[:cutoff, 1])
plt.show()
Opcjonalne: Skrócenie rekordów na potrzeby szybszego przetwarzania¶
sig0 = signals_850[:, 0]
sig1 = signals_868[:, 0]
# sig0 = signals_850[500:5000, 0]
# sig1 = signals_868[500:5000, 0]
Ciągła transformata falkowa (CWT)¶
Używam falki mexican hat.
Pierwszy sygnał przetworzony odwrócono, aby ułatwić dalsze przetwarzanie. Drugi sygnał po transformacie (868) podniesiono do kwadratu, aby ułatwić progowanie.
coef0, freqs0 = pywt.cwt(sig0, np.arange(1, 50), 'mexh')
coef1, freqs1 = pywt.cwt(sig1, np.arange(1, 50), 'mexh')
my_coef_0 = coef0[3]
my_coef_1 = coef1[1]
my_coef_0 = -my_coef_0
my_coef_1 = np.square(my_coef_1)
plt.figure(2, figsize=(24, 8))
plt.subplot(2, 2, 1)
plt.plot(sig0[:cutoff])
plt.title('Original signal 850')
plt.grid(True)
plt.subplot(2, 2, 2)
plt.plot(sig1[:cutoff])
plt.title('Original signal 868')
plt.grid(True)
plt.subplot(2, 2, 3)
plt.plot(my_coef_0[:cutoff])
plt.title('CWT signal #0')
plt.grid(True)
plt.subplot(2, 2, 4)
plt.plot(my_coef_1[:cutoff])
plt.title('CWT signal #1')
plt.grid(True)
Progowanie
my_coef_0[my_coef_0 < 1.8] = 0
my_coef_1[my_coef_1 < 3.5] = 0
Przedstawienie sygnałów oraz ich przetworzonych wersji już po progowaniu.
plt.figure(3, figsize=(24, 8))
plt.subplot(2, 2, 1)
plt.plot(sig0[:cutoff])
plt.title('Original signal 850')
plt.grid(True)
plt.subplot(2, 2, 2)
plt.plot(sig1[:cutoff])
plt.title('Original signal 868')
plt.grid(True)
plt.subplot(2, 2, 3)
plt.plot(my_coef_0[:cutoff])
plt.title('850 CWT after threshold')
plt.grid(True)
plt.subplot(2, 2, 4)
plt.plot(my_coef_1[:cutoff])
plt.title('868 CWT after threshold')
plt.grid(True)
Wykrycie peaków, wyświetlenie ich oraz ich liczby.
from scipy.signal import find_peaks
peaks0 = find_peaks(my_coef_0, height=0)
peaks1 = find_peaks(my_coef_1, height=0)
print('QRSy dla sygnału 850:')
print(len(peaks0[0]))
print(peaks0)
print('\nQRSy dla sygnału 868:')
print(len(peaks1[0]))
print(peaks1)
Wyświetlenie sygnałów wraz z nałożonymi liniami pionowymi wskazującymi na wykryte załamki QRS (czerwone linie).
%matplotlib inline
plt.figure(4, figsize=(24, 8))
plt.subplot(2, 1, 1)
plt.plot(sig0[:cutoff])
plt.title('Original signal 850')
plt.grid(True)
for peak_pos in peaks0[0][:9]:
plt.axvline(x=peak_pos, c="r")
plt.subplot(2, 1, 2)
plt.plot(sig1[:cutoff])
plt.title('Original signal 868')
plt.grid(True)
for peak_pos in peaks1[0][:14]:
plt.axvline(x=peak_pos, c="r")
# Jakich falek można użyć
for i in pywt.wavelist(family=None, kind='continuous'): print(i, end=' ')
Wczytanie adnotacji dla obu sygnałów.
from multiprocessing import freeze_support
import wfdb
import os
ann_850 = wfdb.rdann('850', 'atr')
sampleAnnotation_850 = [(sample, ann_850.symbol[i]) for i, sample in enumerate(ann_850.sample)]
ann_868 = wfdb.rdann('868', 'atr')
sampleAnnotation_868 = [(sample, ann_868.symbol[i]) for i, sample in enumerate(ann_868.sample)]
Funkcja wyświetlająca czułość, precyzję oraz TP, FP, FN.
Sprawdzanie skuteczności detektora polega na przyjęciu przedziałów referencyjnych o szerokości 20 próbek, co przy częstotliwości próbkowania 128Hz odpowiada przedziale około 150ms. Środkami przedziałów referencyjnych są pobrane adnotacje (poprawne miejsca załamków). Następnie wyznaczono liczności zbiorów TP, FP, FN.
- TP: gdy wyliczony załamek trafia do przedziału
- FP: gdy wyliczony załamek
- FN:
def get_TPR_PPV(annotation, peaks):
TP, FP, FN = 0, 0, 0
for ann in annotation:
middle = ann[0]
# range: <start; end>
start = middle-10
end = middle+10
a = peaks[0]
hits = ((start <= a) & (a <= end)).sum()
if hits > 0:
TP += 1
if hits > 1:
FP += hits-1
if hits == 0:
FN += 1
print("|TP|={}, |FP|={}, |FN|={}\n".format(TP, FP, FN))
# Sensitivity / True positive value / czułość
TPR = TP / (TP + FN) * 100
# Precision / Positive predictive value / precyzja
PPV = TP / (TP + FP) * 100
print("Czułość: {}%".format(np.around(TPR, 2)))
print("Precyzja: {}%".format(np.around(PPV, 2)))
return TPR, PPV
Czułość to stosunek wyników prawdziwie dodatnich do sumy prawdziwie dodatnich i fałszywie ujemnych. Czułość 100% w przypadku testu medycznego oznaczałaby, że wszystkie osoby chore lub ogólnie z konkretnymi poszukiwanymi zaburzeniami zostaną rozpoznane. $${\displaystyle \mathrm {TPR} ={\frac {\mathrm {TP} }{P}}={\frac {\mathrm {TP} }{\mathrm {TP} +\mathrm {FN} }}}$$
Precyzja to stosunek prawdziwie pozytywnych wyników do wszystkich pozytywnych wyników. $${\displaystyle \mathrm {PPV} ={\frac {\mathrm {TP} }{\mathrm {TP} +\mathrm {FP} }}}$$
Wysoka czułość oznacza dużą ilość wykryć QRSów.
Wysoka precyzja oznacza dużą poprawną ilość wykryć QRSów.
get_TPR_PPV(sampleAnnotation_850, peaks0)
W przypadku badania 850 czułość wynosi 99.02%, a precyzja 100%.
Oznacza to, że algorytm znajduje prawie wszystkie załamki QRS, natomiast precyzja równa 100% oznacza że nie było błędnych wykryć QRSów.
get_TPR_PPV(sampleAnnotation_868, peaks1)
W przypadku badania 868 czułość wynosi 97.79%, a precyzja 93.9%.
W algorytmie precyzja jest mniejsza od czułości, gdyż niekiedy algorytm dodatkowo wykrywał QRSy. Algorytm wykrył prawie wszystkie załamki QRS. Niższa precyzja wskazuje, że nie wszystkie załamki jakie wykrył były poprawne. Możliwe, że kilka razy wykrywał ten sam załamek, co jest błędem.