Spracovanie signálu ECG za účelom rozpoznávania fázy cyklu srdca

18. December, 2009, Autor článku: Války Gabriel, Elektrotechnika, Medicína, Študentské práce
Ročník 2, číslo 12 This page as PDF Pridať príspevok

Picture 2Práca sa zaoberá spracovaním signálu ECG v reálnom čase. Vychádza z hardvérového a softwérového riešenia prezentovaného v minulom ročníku ŠVOČ. Za cieľ si kladie vytvorenie fázového cyklu srdca ktoré sa použije ako zdroj časovacích signálov pre elektronický synezátor hudby. Za predpokladu prirodzenej synchronizácie ľudského tepu na vonkajší zdroj by modifikácia fázového cyklu srdca mohla byť využitá na korekcie srdcových rytmických porúch.

1. Úvod

Zariadenie [1] (obrázok 1, 2) na interaktívne sledovanie a spracovanie signálov EKG predstavuje komplexnú platformu, ktorá formou cenovej dostupnosti a otvorenej architektúry umožňuje v širokom rozsahu experimentovať a ľubovoľne spracovávať ľudské biosignály. Toto zariadenie pozostáva z aktívnych elektród, [5] ktoré sú pripevnené na tele. Elektródy sú prepojené s modulom, ktorý vykonáva zosilnenie signálu z elektród na požadovanú úroveň, vzorkovanie a prevod vzoriek do digitálnej podoby. Takýto dátový tok je vysielaný prostredníctvom technológie bluetooth do nadradeného počítača alebo mobilného telefónu. Viaceré vedecké štúdie poukazujú na vplyv hudby ako súboru harmonickej, melodickej a rytmickej štruktúry na respiráciu a srdcovú aktivitu [2, 3].

V práci sa zaoberám otázkou, aký bude subjektívny efekt posluchu hudby, ktorá bude s poslucháčom zosynchronizovaná. A taktiež aký je potenciál takto syntetizovanej hudby na cielenú zmenu aktivity srdca. Hlavný predpoklad pre takéto experimentovanie je zostavenie fázového cyklu, ktoré v každom okamihu prijatia vzorky reprezentujúcej elektrický potenciál bude vedieť určiť, v akej fáze periodického cyklu sa srdce nachádza. V súvislosti s riešením tohto problému vznikol nástroj, ktorý metódou autokorelácie vie s dostatočnou presnosťou určiť priebežnú fázu periodického signálu bez znalosti jeho podoby.

2. Fázový cyklus srdca

Predpokladajúc, že priebežne vzorkovaný ECG signál je periodická funkcia s predpisom y = card (\varphi) . Snažíme sa navrhnúť rýchlu numerickú metódu, ktorá bude transformovať vzorky y diskrétne v čase na uhol \varphi = card-1 (y) . Meranie ukázalo, že tvar ECG signálu závisí od zvolenej polohy elektród a v ohraničenom čase merania sa nemení. Napriek tomu sa však mení perióda signálu, pričom tvar jednotlivých charakteristických segmentov ECG signálu ostáva nemenný.

Na rozdiel od bežných inverzných funkcií periodických funkcií (\arccos, \arcsin ), my budeme od našej požadovať výstup v celom intervale 0..2 \pi . Keďže vstupný signál predstavuje neinjektívne zobrazenie, jednému vstupnému argumentu funkcie môže prislúchať viac výsledkov v závislosti od predošlého výsledku operácie.

valky 1
Obr. 1. Modul na spracovanie signálov ECG

valky 2
Obr. 2. Kompletné zariadenie na spracovanie signálov ECG s aktívnymi elektródami

valky 3
Obr. 3. Získaný a spracovávaný záznam

valky 4
Obr. 4. Záznam získaný iným rozložením elektród

Nebude sa teda jednať o konvenčnú matematickú funkciu, ale o algoritmus s pamäťou transformujúci absolútnu hodnotu na uhol. Keďže pracujeme s reálnymi signálmi, algoritmus bude musieť byť odolný voči nasuperponovanému šumu, pomaly meniacemu sa DC offsetu a amplitúdou ktorá sa tiež môže s časom meniť. Dôraz bude taktiež kladený na minimálnu matematickú náročnosť výpočtu, keďže bude tento vykonávaný približne 300 krát za sekundu, čo je zvolená rýchlosť vzorkovania.

Do úvahy spadajú nasledovné metódy:

  • Pamäťová inverzná funkcia, ktorá bude mať predlohu očakávaného signálu a podľa posledného známeho uhlu „napasuje“ súčasnú nameranú vzorku na definovaný graf signálu a určí nový uhol
    • \surd Jednoduchá a veľmi rýchla metóda
    \boxtimes Slabá odolnosť voči šumu
    \boxtimes Nepoužiteľná ak vopred nepoznáme sledovaný priebeh funkcie
  • Využitím Fourierovej transformácie, pri zvolení vhodného tvaru okna a jeho dĺžky podobnej perióde meraného signálu, argument druhého koeficientu transformácie určuje fázu dominantnej zložky.
    • \surd Odolnosť voči šumu
    \surd Vylúčená nutnosť poznať predlohu signálu
    \boxtimes Pomalá – matematicky náročná metóda
  • Korelácia – korelovaním krátkeho úseku posledných nameraných hodnôt s predlohou získanou počas predošlej periódy
    • \surd Metóda odolná voči šumu
    \surd Vylúčená nutnosť poznať predlohu signálu
    \surd Možnosť optimalizácie algoritmu na rýchlosť

3. Autokorelátor so zotrvačnosťou

valky 5
Obr. 5. Princíp autokorelátora

Algoritmus prijíma vzorky a ukladá si ich do pamäte. Jeho prvou úlohou je nájsť význačné body v nameranom priebehu. Ako význačný bod je definovaný ten, ktorý presiahol 90% hodnoty z plávajúceho maxima. Takto zachytíme čelo impulzu charakteristickej QRS vlny. Priestor medzi dvomi poslednými charakteristickými bodmi označíme za aktívny priestor, tento bude pre nás predstavovať matricu – referenčný priebeh signálu. V rámci aktívneho priestoru vyhradíme časový región dĺžky m vzoriek, v ktorom budeme počítať koreláciu so súborom posledných n-nameraných vzoriek (na obr. 5 vpravo). Tento región je vypočítaný ako okolie bodu T1-T0. Pri pravidelnom dokonale periodickom signáli by sa vzorka v tomto bode mala zhodovať s poslednou nameranou vzorkou.

valky 6
Obr. 6. Znázornenie korelácie

Teraz nasleduje samotný numerický výpočet korelácie. Máme dva vektory – jeden dlhší ktorý reprezentuje vybratý región v aktívnom priestore, a jeden kratší ktorý obsahuje posledné namerané vzorky. Obidva vektory vzájomne posúvame v x-ovej aj y-ovej osi a hľadáme kedy sa čo najdokonalejšie prekryjú. Ako hodnotiace kritérium prekrytia som použil osvedčenú metódu najmenších štvorcov.

4. Optimalizácia korelačnej funkcie

Priama implementácia myšlienky posúvania jedného vektora vzhľadom na druhý v rámci možných rozsahov je nepraktická. Keď ešte uvážime, že v každom bode vykonáme výpočet korelácie ako sumu druhých mocnín rozdielov hodnôt, je táto metóda pomalá. Je tu však priestor na optimalizáciu, je zrejmé že pre konštantnú hodnotu posunu v x-ovej osi je možné ideálnu hodnotu y-ového posunu vypočítať.

Predstavme si, že máme zvolenú hodnotu posunu na x-ovej osi a ideme počítať zhodu vektorov metódou najmenších štvorcov. Označme si kratší vektor písmenom A a zodpovedajúci interval dlhšieho vektora písmenom B. Obidva budú rovnakej dĺžky n. Výsledná suma hovoriaca a zhode vektorov bude mať tvar:

\frac{\displaystyle \sum ^n _{i = 1} (A_i - B_i)^2}{\displaystyle n} (1)

Ak by sme chceli vzájomne posúvať tieto vektory v y-ovej osi, matematicky by to znamenalo zaradenie ďalšieho súčtového člena v zátvorke. Sumu by sme mohli testovať pre rôzne hodnoty spomenutého súčtového člena a čakať, kedy bude výsledok minimálny. Čím menšia hodnota sumy nám výjde, tým sú vektory preložené lepšie. V prípade, že by boli identické A_i = B_i, \quad \forall i = 1..n výsledok by bol rovný nule. Hľadáme teda ako sú hodnoty vektorov vzájomne posunuté. Definujme si, že vzorka A_i = a_i + q . Týmto sme ju vlastne posunuli o q, a vzorka B_i = b_i . Dosadením do vzorca č.1 nám výde:

\frac{\displaystyle \sum ^n _{i = 1} (A_i - B_i)^2}{\displaystyle n} = \frac{\displaystyle \sum ^n _{i = 1} ((a_i + q) - b_i )^2}{\displaystyle n} =
= \frac{\displaystyle \sum ^n _{i = 1} (a_i - b_i)^2}{\displaystyle n} + \frac{\displaystyle 2 q \sum ^n _{i = 1} (A_i + B_i )}{\displaystyle n} - q^2 (2)

Prvý člen výsledku je suma štvorcov vektorov bez pridaného konštantného offsetu q. Vo zvyšných dvoch členoch je skryté navýšenie výsledku kvôli tomuto offsetu. Pri bežnom počítaní je q skryté a jeho hodnotu nepoznáme, jednoducho vieme, že vektory sú nejako vzájomne posunuté ale nevieme ako. Poznáme vektory A a B, ale pre nás je zaujímavý iba prvý súčtový člen výsledku, ktorý najlepšie charakterizuje podobnosť vektorov so zanedbaním ich vzájomného posunutia. Zaujíma nás teda, ako dosiahnuť čo najviac zápornú hodnotu prídavných dvoch členov. Zderivovaním podľa q dostaneme:

\frac{\partial}{\partial q} \displaystyle \Bigg ( \displaystyle \frac{ 2 q \sum ^n _{i = 1} (A_i + B_i )}{\displaystyle n} - q^2 \displaystyle \Bigg ) =
= \displaystyle 2 \frac{ \sum ^n _{i = 1} (A_i - B_i)}{\displaystyle n} + \displaystyle \frac{\partial}{\partial q} (q) - \displaystyle \frac{\partial}{\partial q} (q^2) (3)

Hľadáme extrém:

2 \displaystyle \frac{ \sum ^n _{i = 1} (A_i - B_i)}{\displaystyle n} - 2 q = 0
q = + \displaystyle \frac{ \sum ^n _{i = 1} (A_i - B_i)}{\displaystyle n} (4)

Komplikovaným spôsobom sme zistili jednoduchý záver. Tento nám hovorí, že ak chceme dosiahnuť najväčšiu zhodu, treba vektory vzájomne posunúť opačne ako je ich priemerný rozdiel hodnôt. Môžeme to však využiť a dosadiť do rovnice (2).

\displaystyle \frac{ \sum ^n _{i = 1} (a_i - b_i)^2}{\displaystyle n} = \displaystyle \frac{ \sum ^n _{i = 1} (A_i - B_i)^2}{\displaystyle n} - \displaystyle \frac{\Big ( \sum ^n _{i = 1}\Big )^2 (A_i - B_i)}{\displaystyle n} (5)

Týmto sme získali nástroj na zistenie podobnosti vektorov so zanedbaním ich vzájomného posunutia. Rovnica predstavuje známy ukazovateľ z teórie matematickej štatistiky s názvom variancia [6]:

var (X) = E (X^2) - E^2 (X) (6)

Takúto sumu je možné vypočítať v jednej iterácii vektorov, čo môžeme považovať za veľmi rýchle riešenie a ideálne pre naše účely:

fSum = 0, fSumSq = 0
For i = 1 to n
{

fDelta = A[i]-B[i]
fSum += fDelta
fSumSq += fDelta^2

}
fSum /= n, fSumSq /= n
Return fSumSq - fSum*fSum

Autokorelátor teda počíta výsledok korelácie m – n + 1 krát, pričom zvýhodňuje také posunutie v x-ovej osi, ktoré bolo jeho výsledkom v predošlom volaní. Toto sa prejavuje ako zotrvačnosť, ktorá eliminuje vplyvy šumu. Výsledkom autokorelátora je práve posunutie v x-ovej osi, ktoré určuje samotnú fázu vzorky. Táto fáza je určená ako zodpovedajúca poloha nameranej vzorky v rámci aktívneho priestoru. Ak vzorka spadá na začiatok aktívneho priestoru, prislúcha jej fáza 0°, ak vzorka spadá do presného stredu priestoru, prislúcha jej fáza 180°.

5. Aplikácia fázového cyklu

Vytvorili sme fázový cyklus ktorý teraz môžeme prepojiť so syntezátorom hudby. Protokol MIDI, ktorý je zaužívaný komunikačný prostriedok hudobných nástrojov, vyžaduje pre každý takt 24 časovacích impulzov. Celý fázový cyklus rozdelíme podľa tejto požiadavky na 24 intervalov a vždy keď jedným takýmto intervalom prejde fázor fázového cyklu, vyšleme do MIDI zariadenia časovací signál.

valky 7
Obr. 7. Diagram zapojenia

Na tento experiment boli podľa obrázku 7 použité nasledovné zariadenia:

  • MIDI syntezátor Yamaha DJX-1
  • USB-MIDI prevodník M-Audio
  • Laptop s OS Windows a ECG softvérom
  • ECG modul s aktívnymi elektródami

Na obrázku 8 je zobrazený bežný prevádzkový stav programu počas spustenej MIDI synchronizácie. Prvé fázové koleso znázorňuje fázu srdca, druhé fázu generovania MIDI synchronizačných impulzov, tieto sú vzájomne posunuté o hodnotu, ktorá koriguje offset korelátora a latenciu syntezátora.

valky 8
Obr. 8. Ukážka programu

Valky_QT
Ukážka programu – video so zvukom

6. Korekcia srdcového rytmu

Podľa zdroja [2] bol pozorovaný jav, keď uvoľnení kandidáti počúvali rôzne žánre hudby. Pri pomalých skladbách sa kandidátom spomalila respirácia ako aj srdcový rytmus. Podobne pri rýchlych skladbách sa im respirácia zrýchlila spolu aj s tepom. Z tohto usudzujem, že človek má prirodzenú schopnosť a tendenciu k vlastnej synchronizácii k vonkajším zdrojom.

Srdce môžeme chápať ako PLL (Phase Locked Loop) systém vo voľnobežnom stave, ktorý reaguje na aktuálnu požiadavku organizmu dodávať kyslík do orgánov v závislosti od jeho vyťaženia. Pri nájdení vhodného vonkajšieho zdroja sa tento systém snaží s ním v istej miere synchronizovať. Ak by sme zrelaxovanému kandidátovi ponúkli rytmus v rýchlosti blízkej vlastnej aktivite jeho srdca, je pravdepodobné, že by ho PLL systém srdca nasledoval.

Vychádzajúc z tejto domnienky som navrhol metódu na synchronizáciu srdca s elektronickým syntezátorom. Vkladaním korekčného offsetu do fázového cyklu podľa vzťahu 7, môžeme zvoliť cieľovú frekvenciu f_{Dest} .

\displaystyle \Delta = k \Big ( \frac{1}{f_{Dest}} - T_{Meas} \Big ) (7)

Výstupom syntezátora už nie je iba rytmus srdca kandidáta, ale rytmus s vneseným ľahkým driftom, ktorý je proporcionálny ku rozdielu súčasnej frekvencie srdca (definovaný poslednou periódou T_{Meas} ) a nami zvolenej frekvencie f_{Dest} . Ak srdce bije rýchlejšie ako zvolenou frekvenciou, takt generovaného 4/4 rytmu je vyvolaný s oneskorením. Podobne ak srdce bije pomalšie, takt je generovaný v predstihu. Takto reprodukovaná hudba má frekvenciu zhodnú resp. blízku súčasnej frekvencii srdca a jediná zmena je v jej fáze.

7. Záver

Popísaná metóda a algoritmus predstavuje zaujímavý prístup ku korekcii a ovplyvňovaniu srdcového rytmu. Tento efekt by mohol byť prínosný pri štúdiách a liečení srdcových arytmií.

Doterajšie štúdie sledovali vplyvy hudobných nahrávok na fyziologické prejavy organizmu, no žiadna z nich nevyužívala efekt spätnej väzby. Táto v sebe skrýva ďalší priestor na experimentovanie a sledovanie psychologických efektov posluchu hudby na organizmus.

Práca popisuje nekonvenčnú metódu hľadania inverznej funkcie autokoreláciou signálu. Metóda je matematicky aj časovo efektívna a vhodná na spracovanie signálov v reálnom čase.

8. Použitá literatúra

  1. Války, G.: Zariadenie pre dlhodobé monitorovanie aktivity ľudského srdca. In: ŠVOČ 2008 : Zborník víťazných prác. Bratislava, Slovak Republic, 23.4.2008. – Bratislava : STU v Bratislave FEI, 2008. – ISBN 978-80-227-2865-2. – CD-Rom, Dostupné na internete, [20.4.2009]: http://gabonator.yweb.sk/sim/?ecg
  2. Bernardi, L. et. al, „Cardiovascular, cerebrovascular, and respiratory changes induced by different types od music“, Heart, april 2006
  3. Pearson, R. M. „Up tempo, down tempo“, Heart, december 1990
  4. Wikipedia, Electrocardiogram, [online] Publikované 27.03.2008, [citované 12.04.2008], Dostupné z http://en.wikipedia.org/wiki/Electrocardiogram
  5. OpenEEG community, Active Electrodes for the modularEEG, [online] [citované 12.04.2008], Dostupné z http://openeeg.sourceforge.net/doc/hw/ae.html
  6. Volauf, P. „Numerické a štatistické výpočty v Matlabe“, FEI STU BA, 2005
Práca vznikla na Katedre rádioelektroniky, Slovenská technická univerzita, Fakulta elektrotechniky a informatiky STU, Ilkovičova 3, 812 19 Bratislava.

Napísať príspevok