Merre tovább? – Where to go?

Elérkeztem az ösztöndíjas időszakom utolsó hónapjába és további műszaki részletek helyett szeretnék megosztani néhány következtetést.

Az elmúl hónapokban alkalmam volt megismerni jó néhány jelfeldolgozási módszert, amelyek hozzásegítettek megismerni a nem-lineáris jeleket. Ugyanakkor elmélyültem az elektrokardiogramok tanulmányozásában is – megismertem úgy a fiziológiai hátterét, mint a feldolgozásának a lehetőségeit is. A befektetett munka eredményeként részt vettem egy nemzetközi konferencián, egy nemzetközi diákkonferencián, megszületett két cikk és nem utolsó sorban a mesteri képzésemhez tartozó disszertáció.

Természetesen a folyamat nem volt zökkenőmentes. Az időnkénti alkotói válságtól eltekintve egyéb akadályba is ütköztem az ösztöndíjas időszak alatt. A legnagyobb gondot az okozott, hogy nem minden esetben tudtam beszerezni a szükséges EKG jeleket. Bár rendelkezésemre állt a PhysioNet.org online adatbázis, a benne fellelhető minták nem minden esetben feleltek meg az elvárásaimnak. Legtöbbjüket utólagosan digitalizáltak, ezért a mintavételi frekvenciájuk igencsak szűk korlátok között mozog, ami probléma volt, főleg a Hilbert-Huang transzformáció alkalmazásakor.

Arra a felismerésre jutottam hát, hogy érdemes lenne időt és energiát fektetni egy olyan projektbe, aminek a célja egy egységes, nemzetközi kollaboráción alapuló adatbázis létrehozása. Egy ilyen adatbázis felállításának lehetséges lépései:

  • Az adatbázis felépítésének megtervezése. Figyelembe kell venni a jövőbeni felhasználók és felhasználások igényeit. Így az adatbázis az EKG jeleken kívül információt tartalmazna a jelekről: diagnózis, P, T, Q, R, S hullámok helye, esetleg intervallumok és szakaszok helye
  • Az adatbázissal való kommunikációt megengedő szoftver kifejlesztése. Ez a fejlesztés orvosok bevonásával kell megtörténjen, ugyanis az adatbegyűjtés fázisában ők lesznek az elsődleges felhasználói.
  • Az adatbázis feltöltése EKG jelekkel. Ebben a fázisban az adatbázissal való kommunikációt megvalósító szoftver kiosztásra kerül kardiológus szakorvosok között és megkezdődik az adatok begyűjtése. Ezt a lépést esetlegesen megelőzheti egy illesztőprogram-család megírása, amely közvetítene a különböző EKG begyűjtő eszközök és a saját fejlesztésű szoftver között. Ez a lépés igencsak megkönnyítené az orvosok munkáját, illetve kevésbé bonyolítaná azt.
  • Az EKG-k feldolgozásának megkezdése. Amikor megfelelő mennyiségű EKG összegyűlt (főleg egy adott betegségre koncentráltan), meg lehet kezdeni a kutatási tevékenységet.

A projekt végcélja, hogy létrejöjjön egy olyan adatbázis, amit nemzetközi kutatócsoportok is használhatnak, és bővítése folytonosan történik. Így mindig fellelhetőek lennének olyan felvételek, amelyek a legújabb adatbegyűjtő eszközökből származnak.

Nemrég jelent meg egy érdekes tájékoztatás a gyártástrend.hu oldalon arról, hogy Magyarországon kifejlesztésre kerül egy teleradiológia rendszer, amelynek célja összekötni különböző radiológiai intézeteket, elősegítve ezáltal az együttműködést és megkönnyítve a megfelelő kezelés meghatározását. Az EKG adatbázis létrehozását hasonló módon képzeltem el, egy regionális vagy országos kiterjedésű hálózatként.

A pályázat legfontosabb következtetése az, hogy nem lehet elektrofiziológiai jelek feldolgozását kutatni anélkül, hogy szoros együttműködés alakuljon ki orvosok és kutatók között, illetve, hogy ne álljon a kutatók rendelkezésére egy megfelelő adatbázis.

grad_red_right

This is the last month of my scholarship and I thought that instead of presenting further technical details of my research, I would rather draw some conclusions.

In the last months I had the chance to learn about quite a few signal processing methods, which helped me to understand non-stationary signals. At the same time I dug myself into the study of electrocardiograms too – I found out more about both its physiological background and the possible processing methods. As a result of my work I participated at an international conference, an international students’ conference, two research articles were born and last, but not least, I wrote my master’s dissertation.

Obviously, the process wasn’t as smooth as I expected it to be. Apart from the creative crisis that appeared from time to time I faced other difficulties as well. The biggest problem was that I couldn’t get the signals I wanted for my research. Although I had access to the PhysioNet online database, the samples it contains were mostly digitized after the signals were captured on paper. This process took place quite some years ago and, as a result, their sampling frequency is low in most of the cases. The low sampling frequency caused a lot of trouble, especially when making experiments with the Hilbert-Huang transform.

So I was thinking that it would be worth investing time and energy in a project which aims to construction a unified ECG database, based on international collaboration. The possible steps towards such a database are:

  • The design of the database. Both the users’ and applications’ needs should be taken into consideration. This way, the database would contain not only ECG signals but also additional information on the diagnose given based on the ECG signals, the position of the P, T, Q, R and S waves or even the position of the different segments and intervals.
  • The development of the software which enables interfacing with the database. The development should be performed in close collaboration with doctors, since they will be the main users of the software in the phase of gathering signals from patients.
  • Filling the database with ECG signals. In this phase the software is installed on the computers of the doctors participating in the study. In order to speed up the process of collecting ECG signals, driver software can be developed so that the different ECG recording devices could be integrated in the interfacing software.
  • Beginning to process the ECG signals. When there is enough ECG data (reflecting mostly one certain heart condition), the research to develop new processing methods can be started. Naturally, the process of ECG gathering would continue.

The final aim of the project is to create a database which can be used by international research groups and which is extended continuously. This way there would always be recordings in the database that were recorded with the latest recording devices.

Recently an interesting post appeared on the site gyartastrend.hu about how a teleradiology system is getting developed in Hungary. Its goal is to link different radiology institutions, thus helping collaboration and diagnostic decision making easier. I thought about the ECG database in somewhat similar manner.

However, the most important and relevant conclusion of my scholarship period is that the research of electrophysiological signals cannot happen without the help of highly skilled medical professionals. It is also vital for the researchers to have a reliable signal database for daily testing and validation.

Reklámok

Aritmia-analízis Poincaré-ábrával – Arrhythmia analysis with Poincaré maps

Az eddigi bejegyzéseimben próbáltam rávilágítani arra, hogy miként lehetséges az EKG egyes paramétereit azonosítani és megjelölni. Egy EKG-feldolgozó rendszer azonban többet nyújt, mint csupán azonosítást (előfeldolgozást), azaz a felismert jel-részleteket diagnózis felállítására kell felhasználni. Ennek értelmében megvizsgáltam a P, R és T hullámok hasznosíthatóságát. A következőkben az R hullámok helyzetének ismeretére támaszkodó módszert mutatok be, amely alkalmas szívritmus-zavarok (aritmiák) detektálására.

Bármilyen, a normális szívritmustól eltérő szívaktivitást aritmiának nevezünk. Fő kiváltó oka az aritmiáknak a szív elektromos ingerületvezető pályáinak az elváltozása. A különböző típusú aritmiákat jól lehet jellemezni az R csúcsok közötti távolsággal: RR_{i}=R_{i+1}-R_{i}. Az így kiszámolt távolságok sorozatát szívfrekvencia változékonyságnak nevezik (heart rate variability). Ideális esetben a távolságok egy igen szűk intervallumban mozognak, ami normális, hiszen az egészséges szív nyugalmi állapotban megközelítően konstans szívfrekvenciával működik. Abban a pillanatban, amikor az R hullámok közötti távolság erősen változni kezd, aritmiáról van szó.

Az R hullámok közötti távolságokat és azok változását jól lehet ábrázolni egy Poincaré-ábra segítségével. Ez egy olyan két dimenziós ábra, amely páronként ábrázolja az RR távolságokat: (RR_{i}, RR_{i+1}). Az így ábrázolt pontok különböző formákba csoportosulnak. A fenti megfigyelések a következőképpen fordíthatók le a Poincaré-térkép nyelvére:

  • A normál szinuszritmus Poincaré-ábráján az első szögfelező mentén elhelyezkedő pontokat fogunk látni; ezek kör- vagy ellipszisformába csoportosulnak
  • Bármely más szívritmusnak megfelelő Poincaré-ábra egy vagy több formát is tartalmazhat az első szögfelezőn elhelyezkedő mellett.

A két esetet az 1-es és a 2-es ábra szemlélteti. A Poincare-ábra nyújtotta lehetőségeket T hullám alternáns felismerésére és értékelésére is lehet használni. Ebben az esetben az ábrát a T hullám-párok amplitúdói adják és az elemzés során a normálistól eltérő alakzatokat keresünk. A T hullám alternáns detektálása hasznos lehet a hirtelen szívhalál elkerülését segítő defibrillátorok beültetéséhez.

grad_red_right

Up until now I tried to shed some light on why and how the parameters of an ECG can be identified. However, an ECG processing system offers more than just identifying and localizing ECG features, i.e. the extracted features have to be used for diagnosing the patient. For this I have looked for the usefulness of the P, R and T waves. In the following section I will present a very interesting and useful method that makes use of the position of the R waves in order to detect arrhythmia.

Basically, any kind of heart activity that differs from the normal sinus rhythm, is called arrhythmia. In the background of most arrhythmia there is a problem of the heart’s electric circuitry. Different types of arrhythmia can be characterized by the time elapsed between two R peaks: RR_{i}=R_{i+1}-R_{i}. The series obtained from these values is called the heart rate variability (HRV). Ideally, in the case of sinus rhythm, these values vary in a relatively narrow interval, which is logical, because in resting position the human heart beats in an approximately uniform manner. When the interval of RR values enters a larger interval, we can surely talk about some kind of arrhythmia.

The distances between R waves and their variability can be displayed on Poincaré-plots. Such a plot is a two dimensional representation of RR distance pairs: (RR_{i}, RR_{i+1}). There is no temporal information in these plots. The points on a Poincaré-plot are organized into different shapes. The observations stated in the previous paragraph translate to the Poincare-plot world as follows:

  • The normal sinus rhythm appears on the Poincare-plot on the first bisector as a symmetrical blob;
  • Any other blobs that appear on the plot are a sign of arrhythmia. For example, a principal blob and a side blob may suggest an alternating rhythm pattern.

The two cases (i.e. normal and abnormal) are shown on fig. 1 and 2.

Poincare-plot based analysis can also be applied in T-wave alternans detection. In this case pairs corresponding to amplitudes of T-waves are shown on the plot and malformed shapes are being searched after. Detecting T-wave alternans can be really useful to avoid sudden cardiac death by implanting a defibrillator.

Ajánlott olvasmány – Recommended reading

  1. Li Fei, Zhao Jie, Jia Hui-Lin, Zhang Chun-Yun and Zhu Xiao-Lei, Poincare Mapping: A Potential Method for Detection of T-wave Alternans, Procedia Environmental Sciences, vol. 8, pp. 575-581, 2011;
  2. Agnieszka Kitlas Golińska, Poincaré Plots in Analysis of Selected Biomedical Signals, Studies in Logic, Grammar and Rhetoric, vol. 35 (48), pp. 117-127, 2013.

Poincare_NSR

1. ábra – figure 1. – Normális szinuszritmusnak megfelelő Poincaré-ábra; Poincaré-plot corresponding to normal sinus rhythm

Poincare_Afibfig

2. ábra – figure 2. – Pitvari fibrilláció egy rövid intervallumának megfelelő Poincaré-ábra; Poincaré-plot corresponding to a short interval of atrial fibrillation

Magzati EKG és anyai EKG szétválasztása – Separating fetal ECG and maternal ECG

Kutatásom következő lépését az anyai és magzati EKG jelek szétválsztása jelentette. Sok esetben a terhesség során fel kell mérni a magzat szívműküdésének az állapotát. Többféle eljárás is létezik erre a célra, a legismertebb talán az ultrahangos vizsgálat során mintavételezett EKG. Azonban ezen a módon csak a vizsgálat pillanatában lehetséges betekintést nyerni a szív működésébe, szélesebb körű információval nem szolgál a módszer. Ennek érdekében hosszú távú, hordozható EKG monitort csatolnak az anya hasára (több elektróda kerül bizonyos pontokra – az elektódák száma változó) és ez 24 órás megfigyelést biztosít, akár otthon is. A monitorozás végén az EKG jeleket kiértékelik és ilyenkor lehet megállapítani, hogy van-e a magzatnak epizodikus jellegű vezetési (vagy másmilyen) zavara.

Az anya és a magzat EKG-jának szétválasztására való törekvésemet az igazolja, hogy a hasi EKG tartalmazza úgy a magzat mint az anya szivének elektromos aktivitását, azaz a két jel egymásra tevődik a felvett jeleken. Ahhoz, hogy szétválaszthatók legyenek a jelek, legalább két mért jelre van szükség, ha abból a feltételezésből indulunk ki, hogy az anyai és a magzati szív elektormos aktivitása egymástól statisztikailag független. Legalább két mért jelre van szükség a két EKG szétválasztásához – ha feltételezzük, hogy a mért EKG jeleket nem zavarta meg semmilyen zaj a mintavételezéskor. Valós körülmények között ez a feltételezés nem helytálló, ezért a kísérletekben 4 kanálisos hasi EKG felvételeket használtam. Így a négy jel mindegyikében a következő négy eredeti jel lineáris kombinációja található meg: magazti EKG, anyai EKG, a gyerek mozgásából adódó elektromiogram és egyéb zaj.

Az algoritmus, amely megoldja a feladatot, a következő:

  1. Az alapvonal-vándorlás eltávolítása
  2. Szűrés hiánya vagy nagyon óvatos használata
  3. Független-komponens analízis használata a  négy módosított kanálisra
  4. A magzat EKG-jának megjelölése és szívfrekvenciájának meghatározása

Ezzel az algoritmussal 98%-os érzékenységet és 99%-os pozitív prediktivitást sikerült elérni. A következőkben sorra bemutatom az egyes  lépéseket:

1. Az alapvonal-vándorlás eltávolítása volt az első lépés az eredeti jelek uniformizálása fele. Enélkül a lépés nélkül a jelek nem minden esetben összehasonlíthatóak a független komponens felbontás szemszögéből. Az eltávolítást medián szűrővel végeztem el. Annak függvényében, hogy mekora ablakot használtam a medián szűrőhöz, a független komponensek között hol a magzati EKG, hol az anyai EKG hangsúlyozódott ki. Kisebb ablakméret mellett sikerült a legjobb eredményeket elérni.

2. Sem az alapvonal-vándorlás előtt, sem utána nem ajánlott a jel konvencionális értelemben vett szűrése. Alul- és felüláteresztő szűrők egyaránt erősen módosítják a jelek valószínűség-sűrűség függvényét, amely hamis független komponensek kivonását eredményezné.

3. A független komponens analízis különböző rendű statisztikák függetlenségét próbálja garantálni a szétválasztott jelek között. Ily módon a statisztikai függetlenséget nemcsak a korreláció hiánya határozza meg, hanem magasabb rendű statisztikák (szórás, ferdeség, lapultság) függetlensége is. Ezt tükrözi az alábbi egyenlet is: E[x^p]E[y^q]=E[x^py^q]. Egyik módszer az ilyen értelemben vett függetlenség elnyerésére a lapultság maximalizálása. Ebben az esetben felrajzolják a jelek által alkotott egyesített valószínűség-sűrűség függvényt és rajta megkeresik a legnagyobb lapultságot. Az ebbe az irányba mutató vektor segítségével ki lehet vonni a jejelk alkotta mátrixból az egyik független komponenst. A műveletet megismétlik amíg az utolsó kívánt független komponenst is megkapjuk. Jelen esetben a föggetlen komponens analízis sikeresen elkülöníti a magzati EKG-t az anyaitól. A másik két független komponens elektromiogramot és zajt tartalmaz.

4. Mivel a független komponensek sorrendje nem követ semmilyen szabályt, szükség volt egy mechanizmusra, amely képes megkülönböztetni az anya és a magzat EKG-ját. Ennek érdekében a kapott független komponensekben megkerestem az R hullámokat. Amelyik jelben több R hullám volt, az felelt meg a magzati EKG-nak.

A módszernek nyilvánvaló hiányossága is van. Abban az esetben, amikor a magzat bradycardiában szenved, az R hullámokat felhasználó diszkrimináló mechanizmus nem hatékony. Minden más esetben használható, ugyanis a magzat szívfrekvenciája nagyobb az anyáénál.

A legnagyobb problémát az okozza, hogy a hasi EKG-t felhasználó vizsgálati módszer még nem túl elterjedt és nem is szabványosított, ezért nagyon sok féle jellel találkoztam, amelyek más-más pontokon voltak mintavételezve és sokszor nemcsak négy független komponenst tartalmaznak. A független komponens analízis annyi független komponenst szolgáltat, ahány bemenő jel volt. Ha a bemenő jelekben több eredeti jel keveredik, mint ahány bemenő jel van, ezek nem lesznek helyesen szétválasztva.

grad_red_right

The separation of the fetal and maternal ECG constitutes the next step of my research project. In many cases the heart activity of the heart has to be assessed during pregnancy. Naturally, there are already some procedures for this, the most known being the ECG sampled during an ultrasound examination. This method offers only a momentary picture about the fetal heart. For a long term monitoring a holter ECG device can be used which has electrodes placed on the mother’s abdomen. This way episodic appearances of diseases of the fetal heart’s electrical system can be detected – even in the mother’s home environment.

The need for the separation of the maternal and fetal ECG is justified by the very nature of the abdominal ECGs: these contain both the mother’s and the fetus’ ECG. In order for the ECG signals to be separable we need at least two measured abdominal signals, presuming that the ECG of the mother and the one of the fetus are statistically independent.In real world cases two measured signals won’t be enough to extract the fetal and the maternal ECG signals since there are other type of signals that add up to the two ECGs (e.g. power system noise, electromyogram). During the tests I used 4 channeled abdominal ECGs so that every channel contained some amount of the mother’s ECG, the fetal ECG, the mother’s electromyogram and some noise.

The following are the steps that I followed to separate the signals:

  1. Remove baseline drift
  2. Do not or very carefully use filtering methods
  3. Do the independent component analysis of the four channels
  4. Identify fetal ECG and determine its heart rate

I have reached a sensibility of 98% and a positive predictability of 99%. In the next sections each of the above steps are presented in more detail.

1. The baseline drift removal is the first step toward the uniformization of the measured signals. Without this step the signals wouldn’t be comparable in all cases, which would be disadvantageous for the independent component analysis algorithm. I removed the baseline drift using median filters. Depending on the size of the filtering window, either the mother’s or the fetus’ ECG became accentuated after running independent component analysis. In the end, smaller window widths proved to offer the best results.

2. Neither before removing the baseline drift, nor after it the signal should not be filtered. Low pass and high pass filters both cause a modification of the joint probability density function of the measured signals. This modification can cause an alteration of the results of the independent component analysis compared to the unfiltered case.

3. Independent component analysis tries to guarantee the statistical independence of extracted components. This way statistical independence is not only defined by no correlation, but also by the independence of higher order statistics (deviance, skewness, kurtosis). This is described by the equation: E[x^p]E[y^q]=E[x^py^q], where E[g] represents the expected value of the signal g. A method to achieve statistical independence between the extracted components is done by maximizing kurtosis. In this case the joint probability density function is drawn and a vector is rotated until it reaches the maximum of the kurtosis. Then, using this vector, the original component can be extracted from the mixtures. Then these steps are repeated until all components are extracted.

4. The extracted independent components are not ordered in any way. This is the reason for which there is a need of discriminating the fetal ECG from the other signals. I have used an R wave detector based discriminator. First, R waves were identified in the extracted signals. The signal with the most R waves was chosen as the fetal ECG.

The method has an obvious disadvantage. When the fetus suffers from bradycardia (i.e. low heart frequency), the R wave based identification may not be suitable. In all other cases it proved to be efficient, as normally the heart rate of the fetus is higher than the mother’s heart rate.

The biggest problem I’ve faced during the development of the separation algorithm is that there are no standard points where electrodes are to be placed during an abdominal ECG examination  – this may be due to the limited use of this technique. I’ve met groups of channels, that contained more or less than 4 channels and they were sampled at different points of the abdomen. This “diversity” of signals means that the algorithm necessitates considerable changes and cannot be described as a generic one.

4 ponton mintavételezett hasi EKG. Jól látszanak az anya QRS komplexusain kívül megjelenő kisebb amplitúdójú magzati QRS komplexusok. - 4 channel abdomincal ECG recording. Lower amplitude fetal QRS complexes can be observed between the mother's QRS complexes.

4 ponton mintavételezett hasi EKG. Jól látszanak az anya QRS komplexusain kívül megjelenő kisebb amplitúdójú magzati QRS komplexusok. – 4 channel abdomincal ECG recording. Lower amplitude fetal QRS complexes can be observed between the mother’s QRS complexes.

A független komponens analízis (rész)eredményei: az anyai (lent) és a részleges magzati (fent) EKG. Bár a zaj jelentős a magzati jelben, az R hullámok megmaradtak a helyükön, lehetővé téve a szívfrekvencia meghatározását. - The (partial) results of the independent component analysis: fetal (top) and maternal (bottom) ECG. The fetal ECG is heavily corrupted by noise but the R waves remain visible making possible the determination of the fetal heart rate.

A független komponens analízis (rész)eredményei: az anyai (lent) és a részleges magzati (fent) EKG. Bár a zaj jelentős a magzati jelben, az R hullámok megmaradtak a helyükön, lehetővé téve a szívfrekvencia meghatározását. – The (partial) results of the independent component analysis: fetal (top) and maternal (bottom) ECG. The fetal ECG is heavily corrupted by noise but the R waves remain visible making possible the determination of the fetal heart rate.

R hullámok lokalizálása – Identifying R waves

Előző bejegyzésemhez hasonló módon most az R hullámok azonosításáról szeretnék néhány szót ejteni. Az R hullámok a QRS komplexusok szerves és többnyire állandó részei, olyan hullámok, amelyek először “felfele” tartanak és utána “lefele” – megjelenésük a kamrai depolarizációhoz köthető. Mivel a kamrák depolarizációja a legsúlyosabb pitvari ritmuszavarok esetén is kialakul, az R hullámok közötti távolságot szokás használni a szívfrekvencia mérésére [1].

Ami az R hullám lokalizációját illeti, a Matlab honlapján is található olyan bemutató, amely ezt a célt szolgálja [2]. Nevezetesen a findpeaks függvénnyel lehetséges az EKG méréseken megjelölni az R hullámok pozícióját. Ez a megközelítés “több sebből is vérzik”:

  • a findpeaks függvény általános felhasználásra készült, semmilyen EKG-ra vonatkozó speciális információval nem rendelkezik;
  • mivel az általános felhsználhatóság volt a cél, a függvény paraméterezhető; a MinPeakHeight paraméterrel a detektált csúcsok minimális amplitúdója (“magassága”) állítható – nyilvánvaló, hogy ez a paraméter a rögzítés módjától és az EKG-géptől is függ, az elemzést megelőzően ismerni kell;
  • a MinPeakDistance paraméter két csúcs közötti minimális távolságot hivatott rögzíteni, ezt az értéket is szükséges ismerni még az elemzés előtti fázisban.

Egy valamennyire is pontos R hullám detektor az EKG alaptulajdonságait felhsználva kell képes legyen a kamrai depolarizációk időpillanatát azonosítani. Az általam javasolt megoldás lépései a következők:

  1. Az alapvonal-vándorlás eltávolítása. Ez a lépés az EKG uniformizálásához szükséges. Míg előtte az alacsony frekvenciás komponensek az R hullámok amplitúdójának lényeges változását eredményezték, az izoelektromos vonalra való visszahelyezés nagyban segíti az R hullámok amplitújának szűkebb intervallumba való helyezését;
  2. Ahhoz, hogy az R hullámok felismerése során ne zavarjanak más, magas frekvenciás komponensek, egy nem-lineáris zajszűrést végeztem el: kiszámoltam a jelnek megfelelő hisztogramot, azonosítottam a legnépesebb halmazt a 10 elemű hisztogramon, majd a jelben minden olyan értéket, amely a legnépesebb halmaz felső határánál kisebb volt, nullával helyettesítettem; így nem csak a nagyfrekvenciás komponenseket sikerült kiszűrni, hanem minden más értéket is, amely a jel átlaga körül volt – ezáltal nagyrészt az R hullámok maradtak meg;
  3. A megmaradt csúcsok között lelhetők fel az R hullámok. Mivel a kamrai szívizom képtelen 200-300 ms-nál rövidebb idő alatt repolarizálódni és újra depolarizálódni [3], elég, ha ennek az időintervallumnak megfelelő számú mintánként vizsgálom a jelet R hullámokat keresve.

Az algoritmus a legtöbb esetben működőképesnek bizonyult, ám nyilvánvalóan eredményezett néhány hamis pozitív pozíciót is. Az R hullámok felismerésének egyik haszna a már említettk szívfrekvencia meghatározhatósága, de segíthet különbüző arritmiák felismerésében is.

grad_red_right

In my previous post I presented a method for detecting T waves in ECG recordings. This time, Iád like to discuss some aspects on the detection of R waves.  R waves are part of the QRS complexes and they mostly appear in every ECG. They are waves that start with a rising front followed by a falling one. Their appearance is related to ventricular depolarization. Because the ventricles are getting depolarized even in the most severe cases of atrial arrhythmia, the distance between the R waves is used to determine the heart rate [1].

Matlab offers a demonstration on their website regarding the detection of R waves [2] by using the findpeaks function. Using this function has a series of disadvantages:

  • the findpeaks function was designed for general use, it does not take into consideration any particularities of the ECG signal;
  •  because of its general scope, the function takes quite a lot of parameters. The MinPeakHeight parameter, for instance can be used to set a lower limit for the amplitude of the peaks to be detected (i.e. the R waves). It is clear that the value of this parameter has to be known beforehand;
  • the MinPeakDistance parameter sets the minimum threshold for the distance between two consecutive peaks (i.e. R waves); similarly to the MinPeakHeight parameter, this values has to be known also before the whole detection process can be started.

A somewhat useful R wave detector should be able to mark the position of the R waves based on the characteristics of the ECG signal. Considering this criterion, I constructed an algorithm that goes through the following steps:

  1. Remove the baseline-wandering. This step is needed for the uniformization of the ECG. While initially the low frequency components caused the significant variation of the amplitude of the R waves, taking the ECG back to the isoelectric line helps to reduce this variation to a narrow interval;
  2. In order to recognize R waves, high frequency components have to be filtered. I used a non-linear filtering approach. I computed the histogram of the signal in 10 bins. Then I set to zero all the values in the signal that were less than the higher limit of the most populated bin. This way the high frequency components were filtered out as well as any other values that were close to the mean of the signal;
  3. The remaining values contain the R waves. Because the myocardium of the venctricles cannot depolarize and repolarize in less than 200-300 ms [3], I counted every peak considering a number of samples corresponding to this time interval.

The algorithm offered satisfactory results most of the time but naturally produced some false positive results too. The use of knowing the position of R waves is not only limited to the determination of the heart rate but it is also used to identify some types of arrhythmia.

Könyvészet – References

[1] A. R. Houghton and D. Gray, Making sense of the ECG, 2nd edition. CRC Press, 2003.

[2] Peak Analysis – http://www.mathworks.com/help/signal/examples/peak-analysis.html

[3] J. Malmivuo and R. Plonsey, Bioelectromagnetism – Principles and Applications of Bioelectric and Biomagnetic Fields. Oxford University Press, New York, 1995 (http://www.bem.fi/book).

T hullámok lokalizálása – Identifying T waves

A májusi bejegyzésemben röviden bemutattam az EKG-t és azt, hogy miként lehet szűrni azt. Jelen bejegyzésben egy saját módszert szeretnék leírni, amely képes a T hullámok felismerésére, a jelben való lokalizálásukra.

A módszer az empirikus módusfelbontást (EMD) használja fel a T hullámok felismeréséhez. Mivel az EMD által szolgáltatott benső módusfüggvények mindegyikében fellelhető az QRS komplexus valamilyen mértékben, a módszer első fontos lépése ezen komplexusok kiküszöbölése a jelből. A keletkező “hézagot” egyszerűen lineáris interpolációval pótoljuk. Az 1. ábra szemlélteti ezt a műveletet.

Miután a jelet szabaddá tettük a QRS komplexusoktól, újra elvégezzük az empirikus módusfelbontást az módosított jelen és immár az 5., 6. és 7. benső módusföggvényekben észrevehetővé válik egy-egy olyan rezgés, amely alakra és pozícióra is hasonlít a jelben levő T hullámhoz. Az említett benső módusfüggvényeket összeadjuk, majd a QRS komplexusok helye után először előforduló rezgést jelöljük meg T hullámként. Ezen lépéseket a 2. ábra szemlélteti.

A módszer hatékonyságát bizonyítja, hogy 2555 T hullámból 2450-et ismert fel.Érzékenysége 99.96% és pozitív előrejelző képessége 95.99%.

A módszernek természetesen megvannak a maga hátrányai is:

  1. A helyes T hullámot jelképező rezgések körül kisebb amplitúdójú rezgések is megjelenhetnek a benső módusfüggvények összegében, amelyek megzavarhatják a T hullámok pozíciójának a felismerését. Kiküszöbölésükhöz elegendőnek bizonyult egy 4.-5. fokú átlagoló szűrő;
  2. Időnként megjelenik a módusfüggvények összegében egy, a T hullám “hasonmásával” azonos amplitúdójú rezgés, amelyet a módszer T hullámként ismer fel. Ebben az esetben az algoritmus hibás eredményt nyújt.

A bemutatott módszer részletesebb leírása megtalálható egy nemrég megjelent cikkben itt: http://scientificbulletin.upm.ro/papers/2014-1/10_Mozes_Szala.pdf

grad_red_right

In my last blog post I made a brief presentation of the ECG and how it can be preprocessed and filtered.  In this post I would like to describe a method which I contributed to and that is capable of recognizing and localizing T waves from the ECG.

The method makes use of the empirical mode decomposition (EMD) to detect T waves. Because all the intrinsic mode functions provided by the EMD contain a trace of the QRS complex in some form or another, the proposed method’s first aim is to eliminate QRS complexes from the ECG signal. The gaps that appear in the signal are filled by using linear interpolation between the starting point and the endpoint of the intervals from where the QRS complexes were cut out. Figure 1 illustrates the process of QRS complex elimination.

After the signal is freed of QRS complexes it undergoes EMD again. Then the 5th, 6th and sometimes 7th intrinsic mode functions will contain oscillations that resemble T waves both in their morphology and position. These intrinsic mode functions are summed then the oscillations that appear after the intervals which contained once the QRS complexes are marked. These steps can be seen on fig. 2.

The efficiency of the method is proven by the fact that 2450 T waves were identified out of 2555 during a test. The sensibility of the algorithm is 99.96% while its positive predictability is 95.99%.

However, the method has its disadvantages:

  1. The oscillations representing the T wave in the intrinsic mode functions may be preceded by smaller amplitude oscillations that can influence the outcome of the method. These can be filtered using a mean filter of order 4-5.;
  2. From time to time, oscillations with similar amplitude to the one we’re interested in appear right after the QRS interval and before the oscillation that looks like a T wave.  These cannot be filtered out and the method will prove to be wrong in these cases.

The presented method is described more in detail in an article that has only appeared recently: http://scientificbulletin.upm.ro/papers/2014-1/10_Mozes_Szala.pdf

shannon

1. ábra – figure 1. A QRS komplexusok azonosítása az EKG-ban. Fent: az eredeti EKG jel, középen: az EKG első három benső módusfüggvényének az összege, lent: az összeg Shannon energiaburkolója. A függőleges vonalak jelzik a meghatározott QRS komplexusok határait. – The identification of the QRS complexes in an ECG signal. Top: the original ECG signal, middle: the sum of the first three intrinsic mode functions of the ECG signal, bottom: the Shannon energy envelope of the sum. The vertical lines represent the boundaries of the identified QRS complexes.

waves

2. ábra – figure. 2. A T hullámok lokalizálása. Fent: az eredeti EKG jel, középen: az EKG jel a QRS komplexusok nélkül, lent: az 5. és 6. benső módusfüggvény összege. Jól látható, hogy a QRS intervallumok utáni első csúcs a T hullámoknak felel meg. – The localization of the T waves. Top: the original ECG signal, middle: the ECG signal without the QRS complexes, bottom: the sum of the 5th and 6th intrinsic mode functions. It can be seen that the first peak after the interval of a QRS complex corresponds to a T wave.

Az EKG és annak előfeldolgozása – The ECG and its preprocessing

Az elektrokardiogram (EKG) az emberi szív elektromos aktivitását hivatott megjeleníteni – papíron vagy digitális médiumon. Ennek elérése érdekében a bőr felszínére különböző pontokon elektródákat helyeznek és a bőr felszínén megjelenő potenciált mérik. A jól ismert EKG-hullámforma a különböző elektródák közötti potenciálkülönbségek megjelenítése [1]. Ezeket a különbségeket elvezetéseknek nevezzük. A leggyakoribb mérési módszer a 12 elvezetéses EKG, amelynek során tíz elektródát helyeznek a beteg testére: hatot a mellkasra és négyet a végtagokra (egyet-egyet mindegyikre).
A különöző elvezetésekre azért van szükség, mert ezek különböző szögből “néznek” a szív fele, így különböző részek aktivitását képesek kiemelni. Ebből következően egy a szívről akkor alakulhat ki teljes kép, ha minél több elvezetést vizsgálunk. Az 1. ábra a szív anatómiai szerkezetét szemlélteti.
A szívben történő ingerületvezetés módjából adódóan az EKG jelnek specifikus formája van, amely az 2. ábrán látható. A különböző hullámok jelentése:
P hullám: a pitvarok depolarizációjának felel meg
PR szakasz: a pitvari depolarizáció vége és a kamrai depolarizáció kezdete között eltelt idő; ekkor az elektromos impulzus eléri az AV csomót;
QRS komplexus: a kamrák depolarizációjának a megfelelője
ST szakasz: az S hullám végétől a T hullám elejéig tart, ez alatt az idő alatt a myocardium nem képes elektromos aktivitásra
T hullám: a kamrák repolarizációjának felel meg
QT szakasz: a kamrák teljes elektromos aktivitásának az időtartama
U hullám: igen ritkán jelenik meg az EKG-n, eredete még nem ismert

Természetesen a rögzített EKG jelek a legritkább esetben olyan “tiszták”, mint a 2. ábrán látható minta.

Általában gaussi zaj adódik a jelhez. Régebben nem volt ritka a hálózatból eredő 50 Hz-es komponens sem, ma már a legtöbb EKG készülék képes kiszűrni ezt. Egy másik nagyon gyakori zaj a páciens mozgása, légzése által keltett, úgynevezett alapvonalvándorlás [2]. Ez utóbbi műtermék nemcsak az orvost zavarja a felvétel vizuális tanulmányozása közben, hanem a további számítógépes feldolgozást is akadályozza.
A továbbiakban néhány előfeldolgozási módszerre szeretnék rávilágítani. Ezek elvégzése elengedhetetlen az EKG komplex feldolgozásához és paramétereinek kinyeréséhez.

Alapvonalvándorlás szűrése
Több módszer is létezik ennek a zajnak a kiszűrésére. Igen gyakran az alapvonalat alacsony frekvenciás komponensként fogják fel, amit lineáris, véges impulzusválaszú felüláteresztő szűrőkkel távolítanak el. Ezen szűrők vágási frekvenciája a 0.5 és 2 Hz közötti tartományban mozog. Egy másik lehetőséget a medián szűrő használata adja; ez egy nemlineáris szűrő, alapja a medián statisztikai mértéke. Egy számsor mediánját úgy kapjuk meg, hog a számokat rendezzük és kiválasztjuk a középső értéket. Attól függően, hogy mekkora ablakméretet választunk a mediánszűrőnek, különböző eredményekhez juthatunk. A legjobb, ha egy periódusnak megfelelő hosszúságú ablakot választunk (azaz akkorát, hogy egy P, QRS, T hullámegyüttes beleférjen). A 3. ábra bemutat egy eredeti minát és az alapvonalvándorlás eltörlése után kapott jelet.

Zajszűrés
A második nagyon fontos lépés az EKG előfeldolgozásában a nagyfrekvenciás komponensek kiszűrése. Ezek a komponensek általában zajnak tekinthetők, bár egyes esetekben óvatosan kell kezelnünk őket, ugyanis lehetnek a magzati EKG részei is. A legelterjedtebb szűrési eljárás a lineáris, aluláteresztő FIR szűrők használatán alapszik, amelyek vágási frekvenciája a 30 – 40 Hz-es tartományban van. Egy másik módszer szűrésre (úgy zajszűrésre, mint az alapvonalvádnorlás kiküszöbölésére) a wavelet transzformáció használata [3]. Bizonyos számű felbontás után kapott együtthatókat valamilyen küszöbérték szerint ki lehet hagyni a jel visszaépítésénél (finom vagy durva küszöbölés, illetve a Johnston-Donoho féle univerzális küszöb egyaránt használható).

A bemutatott módszerek alkalmazása után olyan jelhez jutunk, amelyet fel lehet használni összetettebb, fejlettebb elemző módszereknél, tipikusan paraméterek kinyerésének céljából.

grad_red_right

The electrocardiogram (ECG) is used to visualize the electrical activity of the human heart. This visualization can be either on paper or in a digital format. To achieve this, electrodes are mounted on the surface of a patient’s skin. These electrodes will measure the potential on the skin. The well-known ECG waveform is actually the result of the potential difference between different electrodes – also known as leads. The most widely used measuring method consists of 12 leads. During a measurement there are ten electrodes on the skin: six in the precordial part of the chest and four on the limbs – one on each of them. Multiple leads are needed in order to obtain different “views” of the heart, which may offer different information about the heart. Thus, a physician can have a whole view of the heart if there are enough leads to be examined. Figure 1. presents the anatomy of the heart. The ECG has a quite specific form, explained by the modality of the propagation of an electric impulse inside the heart (as illustrated in fig. 2.). The meaning of the different parts of the ECG is summarized below:
P wave: corresponds to the atrial depolarization;
PR interval: the time between the end of the atrial depolarization and the start of the ventricular depolarization; this is the time when the impulse arrives at the AV node;
QRS complex: the equivalent of the depolarization of the ventricles;
ST segment: the time between the start of wave S and the beginnig of wave T; during this time, the myocardium is free of any electrical activities;
T wave: corresponds to the repolarization of the ventricles;
QT interval: represents the length of the time while the ventricles are electrically active;
U wave: it only appears rarely; its origins are not yet known.

The recorded ECG signals are rarely as “clean” as the one presented in fig. 2. Generally gaussian noise, power line interference or baseline drift [2] also appear on the ECG. ECG recording equipment nowadays are able to eliminate power line interference. However, the other two artifacts can be disturbing not only for the examining physician but for computers as well, making it difficult to run advanced parameter extraction algorithms on the signal.
In the next sections two important preprocessing steps will be presented, that help transform the signal into a more useable form.

Baseline drift removal
There are several methods for removing the baseline drift (also known as baseline wandering). One of them consists of applying a linear, finite impulse response (FIR) high-pass filter of a cutting frequency between 0.5-2 Hz. Another possibility is offered by the median filter (which is a non-linear filter) The median filter makes use of the median, that is the middle value in an ordered statistical series of data. Depending on the width of the examining (filtering) window, different results can be obtained. An ideal value for the median filter would be the length of a set of P,QRS and T waves. The baseline obtained by the median filter is then subtracted from the signal to obtain a more cleaner signal. Figure 3. shows an original signal and the result of the baseline drift removal process.

Filtering
The second very important step in preprocessing the ECG is the removal of high frequency components from the signal. This happens by applying a linear, low-pass FIR filter with a cutting frequency from the interval [30, 40] Hz. A more general method for threshold filtering (both of the high frequency noise or the low frequency baseline drift removal; soft, hard and the Johnston-Donoho universal threshold can be used) is making use of the wavelet transform [3]. Failing to reconstruct the signal from all of its details will lead to a successful filtering.

If applied, the presented methods will offer a signal that can be passed on to more complex and advanced analysis methods, tipically with the final purpose of parameter extraction.

 

anatomy_of_the_heart

1. A szív anatómiai felépítése – The anatomycal structure of the heart (source: http://howmed.net/anatomy/gross-anatomy/aorta/)

ecg_form

2. Az ideális EKG egy periódusa a főbb paraméterekkel – One period of the ideal ECG signal with the main parameters (source: http://www.todayifoundout.com/index.php/2011/10/how-to-read-an-ekg-electrocardiograph/)

baseline

3. Eredeti (fent) és alapvonalvándorlás nélküli (lent) hasi EKG – Abdominal ECG with (top) and without (bottom) baseline wandering

Irodalomjegyzék – References

[1] Houghton, A., Gray, D. (2008), Making Sense of the ECG: A Hands-on Guide, Third Edition, CRC Press

[2] Kaur, M., Singh, B., Longowal, S. (2011), Comparisons of Different Approaches for Removal of Baseline Wander from ECG Signal, Proceedings of the 2nd International Conference and Workshop on Emerging Trends in Technology, pp. 30-36

[3] Borries, R. F., Pierluissi, J. H., Nazeran, H. (2005), Wavelet Transform-Based ECG Baseline Drift Removal for Body Surface Potential Mapping, Proceedings of the 2005 IEEE, pp. 3891-3894.

Egy lehetséges mód a pillanatnyi frekvencia meghatározására – A possible way of computing the instantaneous frequency

Folytatva a legutóbbi poszt gondolatmenetét, szeretném bemutatni a második utat amelyet követve abban reménykedtem, hogy sikerül egy egyszerűbb, pontosabb módon meghatározni a pillanatnyi frekvenciát.

Az alapötlet az empirikus módusfelbontásból származó benső módusfüggvények modellezése volt. Ezeket a függvényeket meg lehet feleltetni egy változó körfrekvenciájú és amplitúdójú körmozgás egyik vetületének. Ha ismernénk a másik vetületet, akkor lehetővé válna a pillanatnyi frekvencia kiszámítása. Hasonló módon működik a Hilbert transzformáció és az AM-FM felbontás is.

A körmozgás két vetülete között létezik egy összefüggés, amely egyenletes körmozgás esetén a jól ismert   \sin^2(t)+\cos^2(t)=1 egyenletre redukálódik. Általános esetben fennáll az   \displaystyle \dot{x}^2+\dot{y}^2=v^2 összefüggés, ahol x és y a körmozgás két vetülete és v a pillanatnyi sebességvektor modulusza. A két vetület deriváltja nem más, mint a sebességvektor komponensei a koordinátatengelyek szerint. Ahhoz, hogy kiküszöböljük v-t az egyenletből, felhasználhatjuk a két vetület általános alakját:   \displaystyle x(t)=A(t)\sin(\omega(t)), y(t)=A(t)\cos(\omega(t)). Innen   \displaystyle \dot{x}(t)=\dot{A}(t)\sin(\omega(t))+A(t)\dot{\omega}(t)\cos(\omega(t)) és \displaystyle \dot{y}(t)=\dot{A}(t)\cos(\omega(t))-A(t)\dot{\omega}(t)\sin(\omega(t)), ami az \displaystyle \dot{x}^2+\dot{y}^2=\dot{A}^2+\dot{\omega}^2A^2 eredményhez vezet. Ez nyilván egyenlő a sebességvektor abszolút értékének négyzetével.

A következő lépés az, hogy kifejezzük y-t x segítségével. Azonban ez a kifejezés nem enged közelebb a kvadratúra meghatározásához, mivel a kapott egyenletben két ismeretlen jelenik meg: az amplitúdó és a körfrekvencia. Ezért a következő próba a gyorsulás bevezetésére irányult, annak a reményében, hogy sikerül még egy egyenletet felírni és így meghatározni a kvadratúrát. A gyorsulás két komponensének (érintőleges és sugárirányú) kifejezése: \displaystyle a_t(t)=A(t)\dot{\omega}(t)+2\dot{A}(t)\omega(t) és \displaystyle a_r(t)=\ddot{A}(t)-A(t)\omega(t). Ha ezt a két komponenst felhasználjuk a gyorsulásvektor moduluszának kiszámításához, ahogy a két vetület (x és y) dupla deriváltjának négyzetösszegét is, egy többedrendű differenciálegyenlet lenne az eredmény, amibe x-et be lehet helyettesíteni. A legnagyobb hátránya ennek az egyenletnek az, hogy megoldásához szükség van néhány (az egyenlet rendjétől függően) kezdőértékre, ami a teljesen ismeretlen kvadratúra esetén nehézséget okoz.

A Wolfram Mathematica programcsomaggal végzett kísérletek után és tekintettel az egyenletrendszer a megszokottnál magasabb fokára, az ezen az úton való haladást elvetettem, a sokkal igéretesebb empirikus módusfelbontás javára.

grad_red_right

As a continuation of my last post, I’d like to present a second way I tried in my hope to compute the instantaneous frequency more efficiently and in a simpler manner.

The basic idea lies in the modeling of intrinsic mode functions that are the output of the empirical mode decomposition. These mode functions correspond to a projection of a non-uniform circular motion, i.e. a circular motion with variable frequency and amplitude. If the other projection of the circular motion would be known, the computation of the instantaneous frequency would become possible. The Hilbert transform or the AM-FM decomposition based methods work in a very similar manner.

A relation exist between the projections of a circular motion, which, in case of a uniform circular motion reduces to the well known equation \sin^2(t)+\cos^2(t)=1. In the general case the \displaystyle \dot{x}^2+\dot{y}^2=v^2 relation holds, where x and y are the two projections of the circular motion and v is the absolute value of the velocity vector. The derivative of the two projections are the components of the velocity vector with respect to the axes of coordinates. In order to eliminate v from the equation, the general form of the two projections have been used: \displaystyle x(t)=A(t)\sin(\omega(t)), y(t)=A(t)\cos(\omega(t)). From here, \displaystyle \dot{x}(t)=\dot{A}(t)\sin(\omega(t))+A(t)\dot{\omega}(t)\cos(\omega(t)) and \displaystyle \dot{y}(t)=\dot{A}(t)\cos(\omega(t))-A(t)\dot{\omega}(t)\sin(\omega(t)), which yields to \displaystyle \dot{x}^2+\dot{y}^2=\dot{A}^2+\dot{\omega}^2A^2. Naturally, this is equal to the square of the absolute value of the velocity vector. The next step is to find an expression for y by using x. However, this expression does not let us closer to determine the quadrature of the intrinsic mode function, because the equation contains two variables: the amplitude and the angular frequency.

Thus, the next step was to try to make use of the acceleration, hoping that this track will lead to another equation for y. The two components of the acceleration, i.e. the tangential and the radial, are: \displaystyle a_t(t)=A(t)\dot{\omega}(t)+2\dot{A}(t)\omega(t) and \displaystyle a_r(t)=\ddot{A}(t)-A(t)\omega(t). If we use these two components to compute the absolute value of the acceleration vector, as well as the second order derivatives of the two projections we would obtain a higher order, non-linear differential equation, from where y could be expressed. The major problem of this approach is that approximating the derivatives, some initial values are needed (depending on the order of the equation) which is an obvious disadvantage in the case of an unknown quadrature.

After experiments made in Wolfram Mathematica and being aware of the complexity of the equations, I gave up continuing to follow this path and turned towards the more promising empirical mode decomposition.

A körmozgás sebessége és gyorsulása. The velocity and acceleration of the circular motion. (Forrás/Source: http://upload.wikimedia.org/wikipedia/commons/2/22/Nonuniform_circular_motion.svg)