[Matlab] Pan Tompkins Algoritması Kullanılarak EKG’de QRS Complex Tespiti ve Ritm Analizi

Biomedical Signal Analysis dersinde Timur Düzenli ve Salih Aslan ile beraber yaptığımız bir proje bu. EKG verileri üzerinden QRS Complex kısımlarının tespitinden bahsedeceğim. Öncelikle nedir ve neden önemlidir buna değinelim.

Çoğumuzun bildiği üzere EKG (ECG - Electrocardiogram) kalp aktivitelerinin elektriksel olarak kaydıdır. Electrocardiograph ile ölçülerek kaydedilir.

EKG analiz edilerek kalp hakkında oldukça detaylı bilgiler elde edilebilir. Tıp okumadığım için çok detaylı bilemem tabii ki. Fakat EKG sinyalinin birazdan bahsedeceğim kısımlarının zamanlaması hastalıklar hakkında oldukça yararlı bilgiler sağlamaktadır.

Tipik bir EKG sinyali aşağıdaki şekilde gösterilebilir.

Burada görüleceği gibi EKG sinyali belirli tipik kısımlar içerir. Bunlar P,Q,R,S ve T kısımlarıdır. Hepsi kalp hakkında önemli bilgiler taşır. Bu yazıda, sinyalde QRS Comlex’in yerinin bulunmasından bahsedeceğim. Bu bilgi önemlidir çünkü QRS Complex’in kaç milisaniye sürdüğü çıplak gözle anlaşılamaz.

Konuya girmeden önce şundan da bahsetmek isterim. Doktorların, ellerine gelen EKG kaydına bakarak karar vermesinde bahsedeceğim konu oldukça önem taşır. EKG sinyalleri elektriksel olarak zayıf sinyaller olduğu için kaydedilmesi ve çevre gürültülerden arındırılması işi oldukça zahmetlidir. Demek istediğim, ham bir EKG kaydından hiç bir doktor hiç bir teşhis koyamaz. O şeritteki net ve anlaşılır EKG’nin nasıl o hale geldiğini ele alacağız. Fakat tüm aşamaları koca bir kitap yapan bu analizin sadece küçük bir kısmı bu yazıca mevcut.

Analiz için kullanacağım önceden kaydedilmiş EKG verileri elimde mevcut. 4000 örnekten oluşan veriyi şuradan alabilirsiniz.

Analize başlayalım,

close all
clear all
clc

%% Degisken ve Sabitlerin Tanimlanmasi
fs=200; % sample rate

%% Sinyalin Matlab Ortamina Alinmasi
hamsinyal=load('ECG3');
plot(hamsinyal);
title('Ham Sinyal');
figure

Görüldüğü gibi sinyal oldukça fazla gürültü içeriyor. Hem yüksek frekanslı hem de şebekeden kaynaklanan ek bileşenler sinyale karışmış durumda. Öncelikle sinyali DC sıfır seviyesine çekelim

%% DC Bilesenlerin Atilmasi
dcsizsinyal=(hamsinyal-mean(hamsinyal));
plot(dcsizsinyal);
title('DC Bilesenleri Atilan Sinyal');
figure

Şimdi de sinyalden yüksek frekanslı gürültüleri temizleyelim. Bunun için bir alçak geçiren (low pass) filtre tasarlayıp sinyali bu filtreden geçirmemiz gerekli. Filtre tipi olarak 10 point avarage filtre kullanalım.

%% Filtre
% 10 point avarage filter
B=(1/10)*ones(1,10);
A=1;
freqz(B,A);
title('10 Point Moving Avarage Filtre');
figure
avaragefiltrelisinyal=filter(B,A,dcsizsinyal);
plot(avaragefiltrelisinyal)
title('Moving Avarage (Low Pass) Filtreden Gecmis Sinyal');
figure

Yüksek frekanslı gürültülerin azaldığını görebiliriz. Fakat hala kusursuz değil. Şimdi de düşük frekanslı gürültüleri yok etmek için bir yüksek geçiren (high pass) filtre tasarlayıp sinyali bu filtreden geçirelim. Bunun için de Derivative Based filtre oluşturacağız.

%% Derivative Based Filter
B=(1/1.0025)*[1 -1];% 1.0025 normalizasyon degeri
A=[1 -0.995];
freqz(B,A);
title(’Derivative Based Filter’);
figure
derivativefiltrelisinyal=filter(B,A,avaragefiltrelisinyal);
plot(derivativefiltrelisinyal)
title(’Derivative Based (High Pass) Filtreden Gecmis Sinyal’);
figure

Sinyalden atmamız gereken bir de 60Hz’lik şebeke gürültüsü var. Kayıt esnasında cihazın şehir şebekesinden kaptığı gütültü de sinyalle beraber kaydedilmiş durumda. Zaten bakıldığında çıplak gözle de görülebiliyor. Bu arada kullandığımız örnek yabancı bir kaynaktan olduğu için şebeke gürültüsü 60Hz. Bizim şebekemiz 50Hz. Onu da hatırlatayım.

60Hz’lik gürültüyü atmak için bir filtre tasarlanacak ve önceki filtreden çıkan sinyal bu filtreden geçirilecek. Önemli bir sorun da, şebekeden kapılan gürültünün yalnızca 60Hz değil aynı zamanda bunun harmonikleri olarak ortaya çıkması. Bu yüzden yalnızca 60Hz’i bastıran bir filtre işimiz görmeyecek. Bunun için 60Hz ve tüm harmoniklerini bastırmamızı sağlayacak bir Comb filtre tasarlamamız gerekir. Aşağıda yazdığım filtre katsayılarının nasıl bulunduğuna değinmeyeceğim, zira çok uzun mesele Zaten ben de hazır kullandım. Bulmaya çalışmak apayrı bir mesele.

%% Comb Filter
% 60Hz sebeke gurultusunu ve harmoniklerini bastiran filtre
B=conv([1 1],[0.6310 -0.2149 0.1512 -0.1288 0.1227 -0.1288 0.1512 -0.2149 0.6310]);
A=1;
freqz(B,A);
title(’Comb Filter’);
figure
comb=filter(B,A,derivativefiltrelisinyal);
plot(comb)
title(’Comb (60Hz ve Harmoniklerini Bastiran) Filtreden Gecmis Sinyal’);

Çok kusursuz olmasa da sinyal işlenebilecek hale geldi. Şimdi ise fark alarak sinyaldeki en keskin tepeleri bulalım. Bunlar bize R noktalarını verecek. Ardından da bu R noktasının sağında ve solunda en düşük seviyedeki noktaları bulalım. Bunlar da Q ve S noktalarını verecek.

%% Differentiator
differentiatorcikisi=diff(comb);

%% Squaring Operation
kare=differentiatorcikisi.*differentiatorcikisi;

%% Moving Window Integrator
window=ones(1,30); % N=30 moving window
integral= medfilt1(filter(window,1,kare),10);
delay = ceil(length(window)/2);
integral = integral(delay:length(integral));

%% QRS Arama
max_h = max(integral);
thresh = 0.2;
poss_reg = integral>(thresh*max_h);

sol = find(diff([0 poss_reg'])==1);
sag = find(diff([poss_reg' 0])==-1);

 for i=1:length(sol)
    [maxdeger(i) maxloc(i)] = max( comb(sol(i):sag(i)) );
    maxloc(i) = maxloc(i)-1+sol(i); % offset ekle
    [mindeger(i) minloc(i)] = min( comb(sol(i):sag(i)) );
    %mindeger Q noktasini verir
    minloc(i) = minloc(i)-1+sol(i); % offset ekle
    % Q noktasi icin offset
 end

% minpozisyon=ones(1,4000)*-max(comb);
% minpozisyon(1,minloc)=max(comb);
maxpozisyon=ones(1,4000)*(-max(comb));
maxpozisyon(1,maxloc)=max(comb);
% maxpozisyon(1,sol)=max(comb);
% maxpozisyon(1,sag)=max(comb);
% maxpozisyon(1,minloc)=max(comb);

figure
% plot(comb(1:(length(comb)/4)),’b-’)
plot(comb)
title(’R noktalari belirlenmis sinyal’);
hold on
% plot(minloc,mindeger,’r–’)
plot(maxloc,maxdeger,’g-’)

figure
plot(comb)
title(’R noktalari isaretlenmi� sinyal’);
hold on
%plot(minpozisyon,’r') % Q baslangici
plot(maxpozisyon,’r')
legend(’Comb Filtreden Gecen Sinyal’,'R’);

Görüldüğü gibi QRS Complex’in yeri artık tam olarak bilinebiliyor. İstenen ölçümler buradan çıkarılabilir.

Bu arada sondan üçüncü Complex’de bir hata var. Algoritma burada hatalı sonuç veriyor. Gürültülerin kusursuz bir şekilde atılamamasından kaynaklanan bir hata söz konusu. P noktası olması gerekenden daha büyük olduğu için R olarak algılanmış ve bir hataya sebebiyet vermiş. Fakat sinyalin geri kalanı istediğimiz veriyi elde etmemiz için yeterli veriyi sağlıyor.

Şimdi elde ettiğimiz verilerden bazı sonuçlara ulaşmaya çalışalım. Girilen veride kaç adet beat (vuruş) bulunduğunu tespit edelim. R noktalarını tespit ettiğimize göre ve R noktası kadar vuruş bulunduğunu bildiğime göre kolaylıkla sonuca ulaşabiliriz.

%% Beat Sayisi
beat=length(maxloc)
heartrate=beat/(length(hamsinyal)/fs)

Cevap olarak şu sonuçları elde ediyoruz.

beat =

25
heartrate =

1.2500

Analiz sırasında P noktalarından birisi de R olarak algılanmıştı. Aslında beat sayısı 24 olmalıydı.

Şimdi de RR aralığının ortalamasını bulalım.

%% RR Araligi Ortalamasi
ortalama=0;
for i=1:beat-1
    ortalama=ortalama+(maxloc(i+1)-maxloc(i)-1)/fs;
    u(i)=(maxloc(i+1)-maxloc(i)-1)/fs;
end
ortalama=1000*ortalama/beat
% ms cinsinden bulmak icin 1000 ile carpildi.

Sonuç olarak elde edilen çıktı,

ortalama =

745.8000

%% Standart Sapma
variance=var(u)*1000 %ms cinsinden

variance = 22.4692

QRS genişliklerinin tespiti için ise zaten bulmuş olduğumuz noktaların farkları alınabilir.

% QRS Genisligi
for i=1:beat
    qrs(i)=sag(i)-minloc(i)
end

Çıktısında da QRS genişliklerinin birbirine yakın değerde olduğu gözlenebilir.

Görüldüğü gibi hastalık teşhisi için canla başla çalışanlar yalnızca doktorlar değil Bizleri de unutmayın.

---

Possibly related posts: (automatically generated)

• Hz. Muhammed karikatürleri hakkında:
• Serbest Çizgi’nin D-Evrim sayısı
• Mekanik Obsesyon ve György Ligeti