Notes

# Gerekli kütüphaneler
import csv # veri seti için
import math #matematik hesaplamalar
import pandas as pd # csv dosyasını okumak için
import numpy as np #lineer cebir , dizi işlemleri için
import matplotlib.pyplot as plt # Grafik çizdirme
import scipy.fftpack # Sinyalin FFT'sini almak için
import scipy.signal as signal # sinyal işlemleri için

# veri setinin okunması
dataset = pd.read_csv("noise.csv")
'''verisetinin içerisinde olan "hart" verilerini for ile y dizisine
alır
'''
y = [e for e in dataset.hart]
''' ECG bir sinyal olduğu için sürekli zamandan örnekleme yapılarak
sinyal görüntülenir.
'''
N = len(y)
# örnek aralığı
Fs = 1000
T = 1.0 / Fs # periyodu
#x eksenini hesaplanması. Sinyalin çizdirilmesi için (periyot ve sinyal)
x = np.linspace(0.0, N*T, N)

'''zaman ekseninden fouirer eksenine geçilir
fft alınır(fast fouirer transform) '''
yf = scipy.fftpack.fft(y)
#FFTsi alınan sinyalin x ekseni oluşturulur.
xf = np.linspace(0.0, 1.0//(2.0*T), N//2) # Bu işlem sinyal çizdirme ile alaklıdır
# x'in frekansa donusumu işlemidir standart fft- sinyal işlemidir.

# Zaman ve frekans ekseni için çizimler yapılır.
fig_td = plt.figure()
fig_td.canvas.set_window_title('Time domain signals')
fig_fd = plt.figure()
fig_fd.canvas.set_window_title('Frequency domain signals')
ax1 = fig_td.add_subplot(211)
ax1.set_title('Before filtering')
ax2 = fig_td.add_subplot(212)
ax2.set_title('After filtering')
ax3 = fig_fd.add_subplot(211)
ax3.set_title('Before filtering')
ax4 = fig_fd.add_subplot(212)
ax4.set_title('After filtering')

#Filtrelenmemiş girişleri çizilir
ax1.plot(x,y, color='r', linewidth=0.7)
ax3.plot(xf, 2.0/N * np.abs(yf[:N//2]), color='r', linewidth=0.7, label='raw')
ax3.set_ylim([0 , 0.2])

###50 hzlil gürültüsüyü ortadan kaldırmak için filtreleme katsayılarını hesaplayın###

#band_filt = np.array([45, 55])
#b, a = signal.butter(2, band_filt/(Fs/2), 'bandstop', analog=False)

''' bandpass filtre uygulanı 4 filtre sırası
50/(Fs/2) : kritik frekanslar
low: alçak geçiren filtre
'''
b, a = signal.butter(4, 50/(Fs/2), 'low')

###ax3.plot(w, 20 * np.log10(abs(h)))
#Compute filtered signal

''' Sinyale bir ileri bir geri dijital filtre uygulanır
b: Filtrenin pay katsayısı vektörü.
a : Filtrenin payda katsayısı vektörü
y:Filtrelenecek veri dizisi.
'''
tempf = signal.filtfilt(b,a, y)
#b, a = signal.butter(1, band_filt/(Fs/2), 'bandstop')
tempf = signal.filtfilt(b,a, y)
#Dijital filteleme yapılmış sinyalin fftsi alınır
yff = scipy.fftpack.fft(tempf)

### Temel kayma gürültüsünü ortadan kaldırmak için Kaiser pencere katsayıları hesaplanır ###
nyq_rate = Fs/ 2.0
# filtrede arzı edilen geçirme büyüklüğü
width = 5.0/nyq_rate
# Durdurma bandında istenen zayıflama, dB cinsinden.
ripple_db = 60.0
# FIR filtresi için sıra ve Kaiser parametresini hesaplanması
O, beta = signal.kaiserord(ripple_db, width)
#Filtrenin kesme frekansı
cutoff_hz = 4.0
###Alçak geçiren bir FIR filtresi oluşturmak için bir Kaiser penceresiyle firwin kullanılır.###

taps = signal.firwin(O, cutoff_hz/nyq_rate, window=('kaiser', beta), pass_zero=False)
# FIR filtresiyle x'i filtrelemek için lfilter kullanılır.
y_filt = signal.lfilter(taps, 1.0, tempf)
yff = scipy.fftpack.fft(y_filt)
# filtrelenmiş çıktının çizdirilnmesesi
ax4.plot(xf, 2.0/N * np.abs(yff[:N//2]), color='g', linewidth=0.7)
ax4.set_ylim([0 , 0.2])
ax2.plot(x,y_filt, color='g', linewidth=0.7);



# kalt atımının-vuruşlarının hesaplanması
dataset['filt']=y_filt
'''Her iki yönde 0,75'lik hareketli ortalamayı hesaplanır
ardından do veri kümesini eklenir
'''
hrw = 1 #Örnekleme frekansının oranı olarak tek taraflı pencere boyutu
fs = 333 #Örnek veri seti 300Hz'de kaydedildi
#mov_avg = pd.rolling_mean(dataset.filt, window=(int(hrw*fs))) #Calculate moving average
mov_avg = dataset.filt.rolling(int(hrw * fs)).mean()

'''x hrw'nin bulunduğu sinyalin başlangıcı olan NaN değerini döndürülmesi '''
avg_hr = (np.mean(dataset.filt))
mov_avg = [avg_hr if math.isnan(x) else x for x in mov_avg]
mov_avg = [(0.5+x) for x in mov_avg]

#ikincil kalp kasılmasının araya girmesini önlemek için ortalam %20 artırılmış

mov_avg = [x*1.2 for x in mov_avg]
dataset['filt_rollingmean'] = mov_avg # hareketli ortalama dataframe eklenmiş

''' #ilgili bölgelerin işaretlenmesi için işlemler yapılır.
uygulanmış filterler for ile datapoinete atılır ve bunlar append ile
windowun içine atılır. Bu işlem ile r komples adı verilen etkinlik kontorl edilecektir.
bu işlemlerin tamamı ECG sinyalinin analizi ile alakalıdır.
'''
window = []
peaklist = [] # tepe noktaları eklenecek -kalp atışının maks noktaları
listpos = 0 #Farklı veri sütunları üzerinde hareket edebilmek için bir sayaç kullanılmış
for datapoint in dataset.filt:
rollingmean = dataset.filt_rollingmean[listpos] #yerel ortalamayı al

#Tespit edilebilir bir R-kompleks etkinliği yoksa -> hiçbir şey yapmayın
if (datapoint < rollingmean) and (len(window) < 1):
listpos += 1
#Sinyal yerel ortalamanın üzerine çıkarsa, ROI(region of interset -ilgi alanı) işaretlenir.
elif (datapoint > rollingmean):

window.append(datapoint)
listpos += 1

else: #Sinyal yerel ortalamanın altına düşerse -> en yüksek noktayı belirle
maximum = max(window)
#X eksenindeki noktanın konumunu not eder
beatposition = listpos - len(window) + (window.index(max(window)))
peaklist.append(beatposition) #Algılanan zirveyi listeye ekle
window = [] #İşaretli ROI'yi temizle
listpos += 1
#Çizim amacıyla tüm tepe noktalarının y değeri alınır
ybeat = [dataset.filt[x] for x in peaklist]

# yılarıda elde edilen sonuçlar çizdirilir
fig_hr= plt.figure()
fig_hr.canvas.set_window_title('Peak detector')
ax5 = fig_hr.add_subplot(111)
ax5.set_title("Detected peaks in signal")
#ax6.set_xlim(0,2500)
ax5.plot(dataset.filt, alpha=0.5, color='blue')
ax5.plot(mov_avg, color ='green') #Plot hareketli ortalama
ax5.scatter(peaklist, ybeat, color='red') #Plot tespit edilen tepe noktaları

#kalp atış hızını hesapla
RR_list = [] #kalp atış hızının eklenmesi için boş
cnt = 0
while (cnt < (len(peaklist)-1)):
#Örnek sayısı olarak vuruşlar arasındaki mesafeyi hesaplayın
RR_interval = (peaklist[cnt+1] - peaklist[cnt]) #Örnek sayısı olarak vuruşlar arasındaki mesafeyi hesaplayın
ms_dist = ((RR_interval / fs) * 1000.0) #Örnek mesafeleri ms mesafelerine dönüştürün
RR_list.append(ms_dist) #listeye eklenir
cnt += 1

bpm = 60000 / np.mean(RR_list) # 1 dakikalık sürede ortamala R sinyal sıklığı (60000 ms= 1dakika)
print("nnnAverage Heart Beat is: %.01fn" %(bpm)) # çıkan sonuç ondalık basamağa yuvarlanır
print("No of peaks in sample are {0}".format(len(peaklist)))

plt.show()