Analiza Conectivitatii Functionale Si Efective In Prelucrarea Informatiei Cerebrale

CUPRINS

Lista acronimelor

EEG – Electroencefalografie

EcoG – Electrocorticografie

ECG – Electrocardiografie

EMG – Electromiografie

MEG – Magnetoencefalografie

EGG – Electrogastrografie

CT – Computer Tomografie

MRI – Magnetic Resonance Imaging (Imagisticǎ prin Rezonanṭǎ Magneticǎ)

fMRI – functional MRI (Imagisticǎ Funcṭionalǎ prin Rezonanṭǎ Magneticǎ)

RMN – Rezonanṭǎ Magneticǎ Nuclearǎ

EVP – Evoked Potential (Potenṭial Evocat)

ERP – Event-Related Potential (Potenṭial generat de evenimente)

DTI – Diffusion-weighted Imaging Techniques (Tehnici de imagistică prin difuzie)

AR – Autoregresiv

MVAR – Multivariat autoregresiv

tiMVAR – time-invariant MVAR (MVAR invariant ȋn timp)

tvMVAR – time-variant MVAR (MVAR variant ȋn timp)

TFR – Transformata Fourier Rapidǎ

FFT – Fast Fourier Transform (Transformata Fourier Rapidǎ)

TMS – Transcranial Magnetic Stimulation (Stimulare Magneticǎ Transcranianǎ)

REM – Rapid Eye Movement (Miṣcare rapidǎ a ochilor)

PDC – Partial Directed Coherence (Coerenṭǎ parṭialǎ direcṭionatǎ)

PDC-I – Partial Directed Coherence – Invariant (PDC invariant)

GCI – Granger Causality Index (Indexul de cauzalitate Granger)

CPU – Central Processing Unit (Unitate Centrala de Procesare)

GPGPU – General Purpose Graphics Processing Unit (Unitate de procesare graficǎ pentru uz general)

GPU – Graphics Processing Unit (Unitate de procesare graficǎ)

CUDA – Compute Unified Device Architecture (arhitecturǎ CUDA)

Introducere

Întrebări ca „Ce reprezintă un gând?” „Ce reprezintă telepatia sau posibilitatea de a comunica prin mijloace ascunse?” „Este posibilă citirea gândurilor?” au preocupat timp de secole imaginația oamenilor de rând, în timp ce diverse științe apăreau și apuneau totodată, încercând să le răspundă. Timp de generații întregi, subiecte de discuție în cadrul cercurilor științifice au studiat posibilitatea interacțiunii cu creierul uman și a comunicării extra-senzoriale.

Recent, un alt subiect, considerat până nu demult o fantezie tinde să se îndepărteze de domeniul științelor oculte și să se integreze tot mai mult în abordarea tehnologică recentă: controlul cu puterea minții. Deși considerată o practică ocultă până recent, avansul studiilor și a tehnologiei a început să dezrădăcineze concepția populară, încercând să facă posibilă interacțiunea om-mașină, folosindu-se de concepte fizice și metode matematice complexe.

Definirea termenului de electroencefalografie (EEG) de către Hans Berger în 1924 a dus la crearea unei noi direcții în tagma cercetătorilor (neuro-știință cognitivă) și lărgirea tehnicilor de imagistică cerebrală. Descoperirea undelor cerebrale de către Berger și avansul tehnologic recent din domeniul semiconductoarelor și a materialelor electronice au dus la o adevărată revoluție în cadrul studiului creierului uman, deoarece precizia măsurătorilor cu care poate fi efectuată astăzi o electroencefalografie oferă tot mai multe informații.

Astăzi, studiul electroencefalogramei are o multitudine de aplicații, începând cu studiul clinic al bolilor psihice (e.g. epilepsia) și ajungând până la stimularea capacității de relaxare prin dispozitive comerciale. O întrebuințare nouă, aflată în faza de cercetare, este posibilitatea de redare a mobilității pentru persoanele cu dizabilități sau paralizii. Deși aplicații experimentale recente au reușit să permită anumitor pacienți controlul unor membre paralizate cu puterea minții, echipamentele prin care acestea funcționau aveau performanțe limitate sau dimensiuni considerabile, existând încă o distanță considerabilă între acestea și aplicabilitatea comercială (de uz larg). În plus, obținerea de performanțe sporite se bazează, în general, pe utilizarea ECoG (electrocorticogramei), tehnică invazivă ce prezintă riscuri și costuri ridicate, spre deosebire de EEG. Astfel, un nivel tot mai mare de efort de cercetare este concentrat în găsirea unor alternative viabile care să înlăture constrângerile menționate anterior.

De-a lungul ultimelor decade, în paralel cu studiul algoritmilor de prelucrare a informației cerebrale, se pune tot mai mult accentul pe eficientizarea acestora, pentru a se crea posibilitatea prelucrării seturilor mari de date în timp real. Această necesitate a apărut din dorința cercetătorilor de a lărgi orizontul studiului și al aplicațiilor conectivității cerebrale, eliminând constrângerile legate de costuri și putere de procesare. Astfel, deși noțiunile de statistică matematică complexă existau de ceva vreme, studiile recente au dus la apariția unor metode statistice de analiză mult mai eficiente.

Deși apărută în timpul anilor 1960, noțiunea de cauzalitate Granger și-a găsit o gama largă de aplicații treizeci de ani mai târziu în rândul științelor neuro-cognitive și în evaluarea conectivității cerebrale. Formalismul matematic dezvoltat de Clive Granger, un economist britanic, laureat al Premiului Nobel pentru Economie, a dus la simplificarea și îmbunătățirea algoritmilor statistici de evaluare a conectivității cerebrale, contribuind semnificativ la avansul recent al științelor cerebrale.

Deși printre cei mai folosiți algoritmi astăzi pentru determinarea conectivității sunt cei dezvoltați de Granger (pentru studiul în domeniul timp) sau Baccala (pentru studiul coerenței în domeniul frecvența), aplicabilitatea lor este puternic limitată de numărul proceselor analizate. Volumul de calcul necesar evaluării unui set complex de EEG nu permite aplicarea algoritmilor complecși, mai ales în cazul sistemelor de calcul uzuale.

Astfel, în prima parte a lucrării de față vor fi studiate performanțele unor algoritmi de complexitate medie și ridicată pentru studiul conectivității și prelucrării informației cerebrale.

Datorită imposibilității aplicării algoritmilor avansați pe sisteme EEG cu un număr mare de canale, a două parte a lucrării va trata posibilitatea reducerii timpului de calcul printr-o metodă hibridă de analiză a seturilor mari de date.

Evoluția recentă a tehnicii de calcul și migrarea către soluțiile de procesare paralelă și heterogene permit reducerea considerabilă a timpului de procesare în cazul calcului matriceal. Astfel, ultima parte a lucrării va trata posibilitatea reducerii considerabile a efortului computațional prin implementarea unui exemplu simplu de prelucrare masiv-paralelă a unui semnal, în vederea utilizării rezultatului în analiză complexă efectuată cu algoritmii de specialitate.

Fiziologie

În cadrul medicinei moderne se folosește o varietate de tehnici de imagistică a corpului uman. Din grupul de măsurători biologice fac parte electrocardiografia (ECG), electromiografia (EMG), electroencefalografia (EEG), magnetoencefalografia (MEG), electrogastrografia (EGG) și electrooptigrafia (EOG). Tehnicile de imagistică medicală se bazează pe diverse principii fizice și includ Computer Tomografia (CT), imagistică prin rezonanță magnetică (MRI), functional MRI (fMRI), etc.

Introducere in EEG

Electroencefalografia este o tehnică de imagistică medicală care se ocupă cu cuantizarea activității electrice generate de structurile cerebrale. Electroencefalograma (EEG) este definită ca o activitate electrică variabilă în timp, înregistrată la nivelul scalpului și măsurată cu ajutorul electrozilor de metal și a mediilor conductive. Această metodă este o procedura de investigație medicală neinvazivă, ce poate fi aplicată în mod repetat unui pacient, indiferent de vârstă sau de alte limitări.

La activarea celulelor nervoase (neuroni), se produc diferențe de potențial spontane ce duc la apariția unui curent de natură electrică între diverse structuri cerebrale. EEG măsoară, astfel, diferențele de potențial apărute în timpul excitațiilor sinaptice ale dendritelor multor neuroni piramidali din cortexul cerebral. Aceste diferențe de potențial sunt cauzate de însumarea potențialelor postsinaptice ale celulelor cerebrale piramidale care constituie dipoli între soma (corpul neuronului) și dendritele apicale (extremitățile neurale). Curentul electric cerebral are o natură de tip ionic și este constituit din mișcarea particulelor de Na+, K+, Ca++ și Cl- prin canalele membranelor neuronale în direcția indicată de diferența de potențial. Este de remarcat că activarea unui singur neuron este nesemnificativă pentru realizarea sau evidențierea într-un EEG. Astfel, numai populații largi de neuroni activi pot genera activitate electrică perceptibilă la nivelul scalpului, deoarece câmpul electric de la nivelul electrodului este atenuat considerabil de piele, craniu și alte țesuturi (e.g. dura mater) [Teplan, 2002].

Semnalele electrice de la suprafața scalpului trebuie, astfel, amplificate considerabil, pentru ca apoi să poată fi cuantizate și procesate.

Datorită capabilității de a evidenția atât activitatea normală cât și pe cea anormală a creierului, EEG-ul s-a dovedit a fi o tehnică foarte puternică în domeniul neurologiei și a neurofiziologiei clinice.

Un avantaj considerabil al EEG-ului este rezoluția temporală ridicată, echipamentele moderne având rate de eșantionare uzuale de 500 Hz, în timp ce anumite sisteme pentru cercetare pot ajunge și la 20 KHz. În schimb, spre deosebire de metodele de imagistică medicală de tip MRI, CT etc., electroencefalograma prezintă o rezoluție spațială redusă, fiind, ca și metodă, incapabilă să evidențieze punctual structuri sau grupuri restrânse de neuroni care se activează la un moment dat. Practic, rezoluția spațială a EEG-ului crește odată cu numărul de electrozi plasați pe cap. În practică, rar se întâmplă să se depășească un număr de 64 de electrozi plasați pe suprafață scalpului, ducând la o rezoluție de cel mult 3 cm. (mult mai mare față de alte metode cu rezoluții < 1 mm – RMN).

Electrocorticograma este un caz aparte de EEG din cauza plasării electrozilor pe creier. Un avantaj considerabil al electrocorticogramei (ECoG) este reprezentat de o mai bună precizie în măsurarea zonelor de interes. Beneficiile acestui tip de investigație medicală sunt, însă, compensate de caracterul său invaziv și riscant, deoarece ECoG presupune o incizie prin scalp și amplasarea electrozilor pe stratul cerebral. Tocmai din acest motiv, tendințele recente încearcă înlocuirea electrocorticogramei cu EEG, pentru reducerea riscurilor și costurilor unei intervenții. Deși progresul recent în domeniul electroencefalografiei permite investigarea mai multor afecțiuni sau situații, electrocorticografia rămâne o metodă de baza pentru determinarea precisă a zonelor de interes în cazul anumitor proceduri (e.g. excizii pentru ameliorarea crizelor epileptice grave sau, mai nou, operații cu implanturi electronice care să permită redarea mobilității în cazul unor persoane paralizate), în timp ce EEG-ul păstrează încă un rol cu scop consultativ, dar foarte important în determinarea următorilor pași de urmat pentru tratamentul pacientului.

Evoluția tehnologiei și a preciziei electrozilor și echipamentelor de măsură a dus la o răspândire largă a electroencefalografiei pentru studiul și evaluarea activității funcționale și efective ale creierului, permițând electroencefalografiei să înlocuiască tot mai mult metodele invazive convenționale.

Activitatea electrică a cortexului cerebral apare în intervalul săptămânilor 17 – 23 ale dezvoltării prenatale. Se presupune că la naștere numărul total de neuroni ai creierului atinge o valoare aproximativă de 100 de miliarde de celule, cu o densitate medie de 10.000 de neuroni per milimetru cub. Aceste celule nervoase sunt interconectate în rețele neuronale prin aproximativ 500 de trilioane de sinapse, dezvoltarea cerebrală postnatală bazându-se pe multiplicarea numărului de sinapse al fiecărui neuron și, astfel, pe dezvoltarea rețelei de interconectivitate neuronală.

EEG măsurat în lipsa unui stimul extern poartă denumirea de EEG spontan. Când se vorbește despre aplicarea unui stimul extern, atunci măsurătoarea preia numele de evoked potential (EVP) sau event-related posituații, electrocorticografia rămâne o metodă de baza pentru determinarea precisă a zonelor de interes în cazul anumitor proceduri (e.g. excizii pentru ameliorarea crizelor epileptice grave sau, mai nou, operații cu implanturi electronice care să permită redarea mobilității în cazul unor persoane paralizate), în timp ce EEG-ul păstrează încă un rol cu scop consultativ, dar foarte important în determinarea următorilor pași de urmat pentru tratamentul pacientului.

Evoluția tehnologiei și a preciziei electrozilor și echipamentelor de măsură a dus la o răspândire largă a electroencefalografiei pentru studiul și evaluarea activității funcționale și efective ale creierului, permițând electroencefalografiei să înlocuiască tot mai mult metodele invazive convenționale.

Activitatea electrică a cortexului cerebral apare în intervalul săptămânilor 17 – 23 ale dezvoltării prenatale. Se presupune că la naștere numărul total de neuroni ai creierului atinge o valoare aproximativă de 100 de miliarde de celule, cu o densitate medie de 10.000 de neuroni per milimetru cub. Aceste celule nervoase sunt interconectate în rețele neuronale prin aproximativ 500 de trilioane de sinapse, dezvoltarea cerebrală postnatală bazându-se pe multiplicarea numărului de sinapse al fiecărui neuron și, astfel, pe dezvoltarea rețelei de interconectivitate neuronală.

EEG măsurat în lipsa unui stimul extern poartă denumirea de EEG spontan. Când se vorbește despre aplicarea unui stimul extern, atunci măsurătoarea preia numele de evoked potential (EVP) sau event-related potential (ERP). Amplitudinea unui semnal EEG măsurat pentru un adult în stare conștientă variază între 10uV-100uV. În cazul unor episoade epileptice, amplitudinea EEG poate să crească cu un ordin de mărime. În interiorul cortexului cerebral, amplitudinea se situează în gama 500uV-1500 uV.

Ritmurile EEG

Semnalele înregistrate cu ajutorul electroencefalogramei au un caracter oscilant. De-a lungul timpului, studiile au relevat diverse tipuri de forme de unde caracteristice anumitor stări de concentrare, categorii de vârstă sau chiar rase de animale. Deși cele mai multe studii au fost realizate pe subiecți umani, multe descoperiri și caracterizări fiziologice au fost realizate și pe alte mamifere.

În urmă însumării și particularizării semnalelor EEG pentru situațiile analizate, spectrul undelor cerebrale a fost împărțit în câteva intervale de bază, denumite ritmuri EEG.

Următoarele tipuri de ritmuri au fost clasificate în cazul electroencefalografiei: delta (0.5-4 Hz), theta (4 – 8 Hz), alpha (8-13 Hz), beta (13 – 30 Hz) și gamma (>30 Hz). Componentele gamma ale EEG sunt dificil de înregistrat de electrozii plasați pe scalp, iar frecvența acestora nu depășește 45 Hz. Cu toate acestea, în cazul electrocorticogramei (electrozi intra-cranieni), pot fi înregistrate semnale cu până la 100 Hz. Contribuția diferitelor ritmuri la compunerea electroencefalografiei depinde de vârstă și starea comportamentală a pacientului, foarte mult contând starea de atenție și concentrare (alterness). Avansurile medicinei neuro-cerebrale și a teoriilor și formalismului matematic privind pattern-urile de conectivitate funcțională și efectivă au permis clasificarea stărilor omului pe baza EEG-ului. Cu toate acestea, există diferențe semnificative în compunerea EEG-urilor diferiților subiecți, deoarece pattern-urile de EEG sunt puternic influențate de condiții neuro-patologice, dereglări metabolice, acțiunea drogurilor și a altor stimulente nervoase.

Ritmul Delta este o caracteristică predominantă în EEG înregistrate în timpul somnului adânc al unui om. În acest stadiu, undele Delta au de obicei amplitudini mari (75 – 200 uV) și prezintă coerență puternică pe toată suprafață scalpului.

Ritmul Theta apare rar la adulți, fiind preponderent la rozătoare. În acest caz, gama de frecvențe pentru acest ritm este mai largă (4 – 12Hz), iar undele au amplitudini mari și o caracteristică de tip “dinți de fierăstrău”. În ipoteză, ritmurile theta la rozătoare servesc drept liant în transferul informației între diverse structuri cerebrale. La om, activitatea în banda Theta poate să apară datorită anumitor stări emoționale sau cognitive. O altă cauză de apariție ar fi încetinirea undelor Alpha din motive patologice.

Ritmul Alpha predomină în timpul perioadei conștiente, în special în zonele posterioare ale capului. Pentru evidențierea corectă a acestei benzi, subiectul trebuie să stea cu ochii închiși și să fie relaxat. Astfel, ritmurile alpha caracterizează starea de calm și liniște, putând fi întrerupte sau atenuate de momente de atenție (mai ales vizuală) și efort mental. Deși ritmurile nu au o bandă de frecvență similară cu Alpha, topografia și semnificația fiziologică ale acestora sunt diferite. Acestea sunt corelate cu funcțiile cortexului motor și apar cu predispoziție în zona centrală a capului. Acestea sunt blocate de funcțiile motoare.

Unda Beta este caracteristică stărilor de alertă și atenție, după cum a fost evidențiat în diverse studii umane sau animale.

Unda Gamma este corelată cu procesarea de informație (e.g. recunoașterea unui stimul senzorial) și cu declanșarea unor mișcări voluntare. În concluzie, poate fi asimilat modelul conform căruia cele mai lente ritmuri corticale corespund stării de latență cerebrală, în timp ce ritmurile rapide sunt asociate procesării de informație [Blinowska & Durka].

Figura 2.1. Sursa: [2]. Ritmurile EEG caracteristice, de sus în jos: δ(0.5-4 Hz), θ(4-8 Hz), α(8-13 Hz), β(13-30 Hz). Forma de undă din partea inferioară reprezintă Înregistrarea EEG aferentă unei crize epileptice; de observat că amplitudinea este cu un ordin de mărime mai mare

Sistemul 10-20 de poziționare al electrozilor

Sistemul 10/20 este o metodă internațională recunoscută pentru descrierea localizării punctelor optime pentru plasarea electrozilor pe scalp, în vederea realizării unui EEG, și se bazează pe relația dintre plasamentul unui electrod și zona din cortexul cerebral pe care o măsoară. Numerele ‚10’ și ‚20’ se referă la distanța dintre doi electrozi adiacenți, deoarece aceștia pot fi plasați la o distanță de 10%sau 20 % unul față de celălalt, pe axa longitudinală sau transversală a scalpului [13].

Figura 2.2. a) (stânga) Sursa:[13] Vedere laterală în secțiune a zonelor cerebrale cu identificarea electrozilor;

b) (dreapta) Sursa:[13] Vedere de sus a zonelor cerebrale cu identificarea pozițiilor electrozilor.

Fiecare zonă cerebrală de interes este identificată cu ajutorul unei combinații de caractere. Astfel, lobii sunt identificați cu ajutorul literelor, iar poziția în emisfere cu ajutorul numerelor.

Tabel 2.1. Sursa:[13] Convenție de notații pentru identificarea zonelor de amplasament a electrozilor EEG. (*) Nu există lob central, dar notația cu litera ‚C’ se pastrează pentru o identificare mai ușoară

Poziția marcată cu ‚z’ (zero) se referă la un electrod plasat de-alungul liniei mediene. Numerele pare (2, 4, 6, 8) desemnează electrozi plasați pe emisfera dreaptă, iar numerele impare (1, 3, 5, 7) se referă la electrozi plasați în emisfera stângă.

Pentru trasarea segmentelor de-a lungul cărora se vor stabili punctele optime de plasare a electrozilor se folosesc patru puncte de referință. Pentru axa longitudinală, segmentul are ca limite naison (punctul dintre frunte și nas) și inion (cel mai îndepărtat punct al scalpului posterior, observabil, cu ajutorul unei proeminențe). Pentru axa transversală, cele două puncte de interes simetrice sunt localizate în zona pre-auriculară, în partea anterioară urechii.

Tipuri de conectivitate cerebrală

Conceptul de conectivitate cerebrală se referă la un model de legături anatomice (“conectivitate anatomică”), dependențe statistice (“conectivitate funcțională”) sau interacțiuni cauzale (“conectivitate efectivă”) între diferite părți ale unui sistem nervos. Aceste unități pot corespunde unor neuroni individuali, populații de neuroni sau regiuni cerebrale separate din punct de vedere anatomic.

Modelele de conectivitate au la bază legături structurale cum ar fi sinapsele, căile de fibre nervoase, seturi de relații cauzale sau statistici măsurate ca deplasare de informație, coerență sau cross-corelație.

Activitatea neurală, și, prin extensie, codul neural, sunt constrânse de conectivitate, care devine, astfel, crucială în elucidarea modului de procesare al informației de către neuroni și structurile neuronale.

Conectivitatea anatomică

Conectivitatea anatomică se referă la o rețea de legături fizice sau structurale (sinaptice) ce leagă seturi de neuroni sau elemente neuronale. Acestor conexiuni le sunt asociate atribute structurale încapsulate în parametrii ce descriu legătura biofizică (e.g. puterea legăturii sinaptice).

Figura 3.1. Sursa: [8]. Model de graf de conectivitate anatomică (graf neorientat).

Modelul fizic al legăturilor anatomice este relativ stabil pe perioade scurte de timp (secunde – minute). Pentru perioade mai îndelungate de timp (ore către zile), structura de inteconectivitate devine obiectul unor modificări morfologice seminificative. Din punct de vedere al metodelor de observație, este de luat în considerare că numai un studiu invaziv poate releva cu certitudine conexiuni axonale directe. În contrast, metode de tip neinvaziv (e.g. Diffusion-weighted Imaging Techniques) DTI, au o rezoluție spațială insuficientă, dar sunt folositoare pentru evidențierea schimbărilor temporale în căile de fibre ale întregului creier.

Conectivitatea funcțională

Acest tip de conectivitate este, în contrast cu modelul prezentat anterior, un concept statistic. În general, conectivitatea funcțională evidențiază deviații de la independența statistică dintre regiuni cerebrale distribuite și chiar separate spațial (e.g. regiuni satelitare).

Figura 3.2. Sursa: [8]. Model de graf de conectivitate funcțională (graf neorientat).

Dependența statistică poate fi estimată prin măsurarea gradelor de corelație, covariație, a coerenței spectrale sau a blocării de faza (phase – locking). Conectivitatea funcțională este evaluată, în cele mai multe dintre cazuri, în contextul întregului set de elemente ale unui sistem neuronal, indiferent dacă părți ale acestui sistem sunt conectate fizic sau nu. Spre deosebire de legăturile structurale, conectivitatea funcțională este puternic dependentă de factorul timp, deoarece relațiile statistice dintre neuroni variază rapid și produc schimbări observabile în modelul abordat, cu granularitate de zeci sau sute de milisecunde. Este important de notat că conectivitatea funcțională nu face referire la efecte direcționale specifice și nici la un anumit model structural subzistent.

Conectivitatea efectivă

Poate fi privită ca o sinergie între conectivitatea funcțională și cea structurală, deoarece, aceasta are rolul de a descrie rețele de efecte direcționale provocate de anumiți neuroni și recepționate de alții. În principiu, efectele cauzale pot fi analizate ca perturbații sistematice ale modelului său, din perspectiva cauză-efect, prin analiză în domeniu timp a semnalelor eșantionate. Anumite tehnici pentru extragerea conectivității efective necesită specificarea unui model bazat pe anumiți parametri structurali, în timp ce alte metode de analiză studiază sistemul că un model general, fără necesitatea specificării unor restricții structurale particulare.

Figura 3.3. Sursa: [8]. Model de graf de conectivitate efectivă (graf orientat).

În mod formal, modelele de conectivitate cerebrală pot fi reprezentate sub formă grafică sau sub formă matriceală. Conectivitatea de tip structural este reprezentată de un graf direcționat, cu ponderi pe legături. Ponderile legăturilor reprezintă densitatea sau eficacitatea conexiunilor. Alternativ, aceste ponderi pot fi de tip binar, și pot sugera prezența sau absența unei legături dintre două elemente structurale separate. În contrast, conectivitatea funcțională este reprezentată prin matrici simetrice de elemente, care codifică dependența statistică sau relativă dintre două elemente ale sistemului (neuroni, regiuni cerebrale, voxeli etc.). În practică, se folosește un nivel de încredere pentru stabilirea existenței unei legături funcționale între elementele structurale ale creierului. Prin acesta, valorile elementelor matricei sunt cuantificate binar, fapt ce permite o vizualizare mai ușoară rezultatelor precum și posibilitatea transpunerii reprezentării într-un graf nedirecționat. În contrast, conectivitatea efectivă are asociată o matrice nesimetrică de elemente. În acest caz, fiecare element al structurii efective are o influență ponderată de la toate celelalte elemente, iar acesta influențează ponderat orice alt element. Acest tip de perspectivă duce la generarea unui graf direcționat.

Modele de conectivitate cerebrală

Analiza modelelor de conectivități structurale ale creierului (e.g. matricile de conectivitate largă a cortexului cerebral) permit cuantificarea unei game variate de caracteristici ale rețelelor [Sporns et al., 2004]. Rezultatele au demonstrat că cortexul cerebral este format din clustere de regiuni corticale dense și cuplate reciproc care sunt interconectate global. Aceste modele de conectivitate nu pot fi descrise că fiind complet subordonate unor reguli de coerență, dar nici aleatoare. În fapt, acestea combină aspecte ale acestor extreme. Rețelele corticale de dimensiuni mari preiau anumite proprietăți de la rețelele de dimensiuni inferioare, printre care valori mari pentru coeficienții de cuplare, și sunt supuse unor seturi specifice de constrângeri structurale. Studiile au arătat că influența structurală a regiunilor cerebrale independente permite identificarea și clasificarea nucleelor de conectivitate, situate preponderent în zonele prefrontale și parietale ale cortexului. Deși descrierea completă a rețelelor structurale încă nu a fost realizată [Sporns et al., 2005], utilizarea tehnicilor neinvazive de imagistică prin difuzie promite rezultate bune în atingerea acestui scop.

Rețele cerebrale funcționale și efective

Studii ale modelelor de conectivitate funcțională (bazate pe coerență și corelații) între diferite regiuni corticale arată că rețelele cerebrale funcționale prezintă caracteristici small-world [Achard et al., 2006], relevând organizarea structurală a conexiunilor anatomice. Analiza teoretică a graficelor de conectivitate cerebrală funcțională ajută la identificarea nucleelor funcționale. Aceste nuclee funcționale sunt puternic conectate și reprezintă puncte centrale de integrare și transport al informației cerebrale.

Noțiunea de integrare funcțională definește abilitatea creierul de creea centri delocalizați de procesare a informației specializate, reprezentați de regiuni cerebrale distribuite la nivelul întregului cortex. Acest tip de integrare prezintă diferit noțiunea de conectivitate, întrucât între doi centri intercorelați statistic, forma de cale de transfer al informației are un alt înțeles față de de cel prezentat de conectivitatea anatomică. Astfel, evidențierea unei conexiuni funcționale între două regiuni cerebrale ajută la creearea unei cai de comunicație formată din noduri de transfer și legăturile direcționate dintre acestea. Aceste căi reprezintă secvențe de asocieri statistice între regiunile cerebrale măsurate, dar nu asigură maparea pe structura anatomică de transfer a informației. Drept urmare, acest tip de conectivitate și studiul nivelului de integrare funcțională nu permite evidențierea caracteristicilor unei rețele cerebrale. Caracteristicile unei rețele cerebrale pot fi evidențiate și interpretate corect numai cu ajutorul unor măsuri de conectivitate anatomică. Acest tip de conectivitate, combinată cu studiul integrării funcționale, permite analiza conectivității cerebrale efective [Rubinov & Sporns, 2010].

Conectivitatea cerebrală efectivă a fost studiată folosind diverse tehnici. Modelarea covariată a permis identificarea unor diferențe semnificative în conectivitatea efectivă dintre diferite regiuni cerebrale, atunci când a fost analizată în cadrul unor acțiuni cognitive. Acest aspect reflectă dependența acestor modele de factorul timp sau de acțiunea cognitivă efectuată.

Conectivitatea Granger a fost aplicată cu succes în secvențe temporale de EEG-uri și fMRI-uri și a contribuit la furnizarea informațiilor despre interacțiuni direcționate între elemente neurale în timpul unor acțiuni cognitive sau comportamentale [Brovelli et al., 2004]. Combinația de stimulare magnetică transcraniala cu tehnicile de imagistică funcțională permite folosirea unor perturbații locale în stimularea unor rețele cerebrale în timp ce acestea sunt angrenate în realizarea anumitor acțiuni (E.g. combinația de TMS și electroencefalografia de înalta densitate au relevat o reducere puternică a conectivității efective în timpul somnului non-REM în comparație cu starea de conștiență – când persoana este trează) [Massimini et al, 2005].

Analiza semnalelor și a metodelor

Modelul autoregresiv

În practică, procesul autoregresiv simplu, de ordin 1, poate fi definit ca o serie de timp a cărei valoare curentă depinde numai de eroarea asociată și de valoarea precedentă.

În ecuația (1) μ reprezintă puterea influenței valorii precedente (μ ∈(0 ,1)), iar reprezintă termenii eroare sau inovarea curentă. Acești termeni au rolul de a conferi proprietatea de variabilitate unei serii de timp și, de obicei, reprezintă secvențe de variabile aleatoare cu varianță finită și medie 0 [Horvath & Johnston].

Studiile au arătat că utilizarea modelului autoregresiv simplu nu este, în general, utilă pentru procesele măsurabile, astfel este necesară aplicarea modelelor de ordin superior.

Seria de timp autoregresivă de ordin p (AR(p)) este descrisă de ecuația:

unde și este o variabilă aleatoare de medie 0 și varianță finită.

Analiza modelului multivariat autoregresiv

Majoritatea proceselor măsurabile utilizate pentru achiziția de informație pot fi analizate pentru evidențierea unei interdependențe între două măsurători corespunzătoare unor puncte spațiale diferite sau pentru evidențierea corelației dintre două mărimi fizice diferite, utilizând modele autoregresive (AR). Aceste modele au la bază seriile de timp.

Modelul AR caracterizează dependențele existente într-un set de date, mai precis, influența (din trecut) pe care o variabilă o are asupra altei variabile din setul de date. Cu alte cuvinte, dacă considerăm o serie de timp univariată, ca în cazul modelulului autoregresiv univariat, elementele acesteia dețin informații legate de procesul care a generat seria de timp. Pentru a descrie acest proces se poate modela valoarea curentă a variabilei ca o sumă liniară ponderată a valorilor precedente. Ordinul modelului este dat de numărul de elemente precedente folosite și de ponderile ce caracterizează seria de timp. Astfel, AR se deosebește de tehnicile de regresie care evaluează corelația instantanee dintre date [Korhonen et. al, 1996].

Modelul multivariat autoregresiv (MVAR) reprezintă o extindere a modelului descris mai sus pentru serii de timp multiple astfel încât vectorul valorilor curente a tuturor variabilelor este modelat ca o sumă liniară a valorilor precedente. Dacă considerăm M serii de timp generate de M variabile dintr-un set de date, și p ordinul modelului atunci MVAR estimează următoarea valoare sub formă unei serii de timp de dimensiune M, X(n), ca o combinație liniară a p vectori de elemente precedenți.

Astfel, modelul MVAR este definit prin ecuația [Winterhalder et al., 2005]:

unde x(n) = [x(1), x(2),…, x(M)] este elementul n al unei serii de timp de dimensiune M, sunt matricele MxM ale parametrilor (ponderilor) și și este zgomot aditiv Gaussian cu medie zero, covarianță R și o variație exprimată prin matricea Σ.

Modelul MVAR poate fi clasificat în funcție de tipul matricei coeficienților astfel:

Modelul autoregresiv multivariat invariant în timp (tiMVAR)

În acest caz nu variază în timp modelul autoregresiv multivariat invariant în timp este numit tiMVAR (time-invariant) iar calculul coeficienților este următorul:

unde

unde coeficienții reprezintă efectul interacțiunii liniare dintre si .

Modelul autoregresiv multivariat variant în timp (tvMVAR)

În acest caz variază în timp, iar coeficienții se calculează după următorul algoritm:

Pentru descrierea estimării algoritmului am folosit următorul model reprezentativ. Considerăm K traiectorii ale unui proces stochastic de dimensiune M. Mai mult, fie a m-a componentă (m=1,…,M) al traiectoriei k (k = 1,…,K) în timp .

Întreaga traiectorie va fi potrivită modelului MVAR. Pentru adaptarea estimativă (folosind algoritmul RLS) a parametrilor modelului A, am definit indicele de performanță care este minimizat de:

unde este predicția instantanee a erorii matricei.

Aceasta reprezintă diferențele dintre răspunsurile dorite și cele primite de estimarea RLS. Matricea colectează ultima estimare p a seriilor de timp denotă suma tuturor componentelor pătrate ale matricei . Utilizarea exponențială a factorului de greutate cu este folosită pentru a asigura că datele din trecut sunt uitate în ordine pentru a permite posibilitatea urmăririi următoarelor variații într-un mediu non staționar. Parametrii matricei pentru care se minimizează, este descrisă de ec. (10) [Moller et. al., 2000]:

După mai multe transpoziții incluzând și folosirea multiplă a matricii inversate, va fi obținută următoarea soluție:

unde reprezintă rândul k al matricei . este numită secvență de control și este definită astfel:

Matricea identitate este adăugata pentru a se obține valori de start stabile pentru calcularea recursivă.

Ecuația (18) are două soluții, în funcție de alegerea adaptării de factor c. Pentru 0 < c < 1 covarianța empirică cu creșterea exponențială a greutății este folosită conform formulei:

Aici efectul ușor al unui trecut bogat, respectând axele de timp, este redus când c crește. Deci 0 < c < 1 ar trebui ales pentru potrivirea instantanee a seriilor de timp non staționare. În cazul special când c = 0 semnificația empirică a covarianței, definită de formula:

este folosită. Mai mult, algoritmul cu c = 0 este apropiat particular de modelul staționar al seriilor de timp.

Estimarea RLS este controlată de valoarea de adaptare a factorului c. Modelând datele EEG, vom obține o medie satisfăcătoare între viteza de adaptare și variația estimării pentru 0.01 < c < 0.05 [Moeller et al, 2001].

Metode de analiză a coerenței semnalelelor

În formalismul studiului semnalelor discrete măsurabile este necesară aplicarea unui model pentru descrierea sau aproximarea evoluției în timp a acelor semnale. Aceste modele pot duce la evidențierea unor caracteristici importante ale semnalelor măsurate, caracteristici inexistente la serii de timp pur-aleatoare.

Conceptul modelelor autoregresive a apărut din necesitatea evidențierii unei legi de variație în timp a unui semnal (non)staționar, precum și a studiului influenței dintre diverse semnale sau mărimi fizice.

Metode neparametrice

Coerența directă

Pentru modelul de coerență simplă, se consideră analiza unui sistem format din două serii de timp complexe x1[n] si x2[n] staționare:

x1[n] X1(Ѡ), x2[n] X2(Ѡ), 0<= Ѡ <= 2π

Se definesc seriile aleatoare

Respectiv

Se definește astfel coerența pătratică a două semnale:

unde este densitatea cross-spectrala, iar sunt densități spectrale ale semnalelor analizate.

Un mod de estimare al acestui tip de coerență este folosind metoda generalizată a lui Welch cu estimarea puterii spectrale:

unde * reprezintă forma complex conjugată a unui semnal, și și sunt transformatele Fourier rapide (TFR) ale segmentelor cu indexul i din n semnale și cu k = 1, . . . , N [Wang & Tang, 2004].

Reducând forma analizată și stabilind n = 1 segment analizat, gradul de coerență între două semnale considerate staționare (de dimensiune suficient de scăzută încât să se presupună staționaritatea) se reprezintă simplificat.

Este ușor de calculat că valoarea ecuației precedente este unitară, ceea ce duce la concluzia că (22) nu se poate aplica pentru o singură fereastră de semnal. Totuși, pentru a putea fi evaluată și aceasta situație, densitatea spectrală de putere dintre doua procese, , este calculată ca fiind transformata Fourier a corelației dintre semnalele x și y. Astfel,

unde este funcția de corelație dintre semnalele x și y.

Pe baza formulei (23) se va realiza ulterior o formă rapidă de determinare a gradului de coerență instantanee între diverse semnale cerebrale măsurate. Avantajul major al formulei precedente este simplitatea sa și eficiența în calcul, permițând analiza în timp real al unui număr foarte mare de canale EEG.

Coerența parțială

Conform definiției, sistemele multivariate permit studiul și particularizarea unui număr mai mare de semnale, în vederea stabilirii existentei sau absenței unei influențe între aceste semnale. În acest scop, pe baza extragerii influențelor de tip linear din toate procesele analizate, analiza coerenței bivariate este extinsă către algoritmul de coerență parțială.

Se consideră două procese de interes în analiză: procesul X și procesul Z. Informația despre coerența dintre acestea este obținută din:

Ec. (24) reprezintă valoarea absolută normalizată a cross-spectrului parțial , obținut luând în considerare procesul adițional sistemului Z:

unde reprezintă densitatea spectrală de putere dintre două semnale și densitați spectrale de putere proprii ale semnalelor X, Y, Z. este argumentul funcției.

Se observă că este dependentă de lărgimea benzii de frecvență în care se află semnalele. Din aceasta cauză, în cazul benzilor înguste de frecvență este greu de estimat o relație între fazele semnalelor. În plus, studierea a doua semnale care se influențează reciproc devine dificila.

Metode parametrice

Spre deosebire de metodele neparametrice, metodele parametrice ce vor fi analizate permit detectarea influențelor direcționate dintre procese în sisteme multivariate. Conceptul care stă la baza metodelor parametrice este noțiunea de cauzalitate introdusă de Granger.

Cauzalitatea Granger se bazează pe relația cauză-efect. Astfel, se spune că un proces P1 influențează Granger procesul P2 dacă informația conținută în valorile precedente ale lui P1 ajută la o mai bună predicție a valorilor lui P2.

În cazul acestei analize se va lua în considerare noțiunea de cauzalitate Granger liniară, deoarece orice proces autoregresiv n-dimensional de ordin p studiat este liniar prin definiția:

unde reprezintă vectorul stării la momentul t al sistemului, eroarea sau zgomot de tip Gaussian. […] reprezinta matricea de ponderi ce poartă informația despre influența cauzală dintre procese. Influența valorilor precedente ale unui proces individual asupra valorilor curente este modelată de elementele de pe diagonala matricei ponderilor [Erla et al., 2009].

Matricea de covarianță a erorii conține informații despre interacțiunile spontane și, astfel, despre influențele non-cauzale dintre procese. Schimbările elementelor diagonale din matricea de covarianță vor permite analiza cauzalității Granger, deoarece elementele eroare oferă informații despre interacțiuni spontane ce nu pot fi obținute din valorile precedente ale proceselor.

Indexul de cauzalitate Granger

Algoritmul presupune evaluarea și identificarea influențelor direcționate dintre componentele unui sistem n-dimensional. Pentru identificarea interacțiunii sunt folosite modele MVAR de dimensiune n și (n-1) ale procesului dinamic care este influențat de . Subsistemul de dimensiune (n-1) este extras cu regula . Pe baza celor doua sisteme se obțin matricele de covarianță = var(, respectiv = var(.

Formula pentru cuantificarea cauzalității Granger lineare variante în timp este:

Deoarece varianța valorilor reziduale a sistemului de ordin n este mai mică decât varianța reziduală a sistemului de ordin (n-1), este de așteptat ca indicele să aibă valori pozitive.

Pentru aplicarea GCI variant în timp in cazul semnalelor non-staționare se folosește o tehnică de adaptare a parametrilor bazata pe un model de filtrare adaptivă, algoritmul RLS (recursive least square). Această abordare permite vizualizarea în domeniul timp a interacțiunilor dintre componentele sistemului precum, urmărind evoluția acestora.

Partial directed coherence

Spre deosebire de indexul de cauzalitate Granger, PDC a fost introdusă pentru studiul și evaluarea interacțiunilor cauzale în domeniul frecvență, permițând diferențierea dintre influențele directe sau indirecte (realizate prin noduri auxiliare). Algoritmul PDC de detecție a interacțiunilor dintre componentele si este:

unde reprezintă componenta matricei ponderilor de pe linia i, coloana j, în domeniul frecvență:

Valoarea reprezintă o valoare normată între 0 și 1. Pentru detectarea existenței unei conectivități este necesar ca indicele să fie nenul.

Prelucrarea semnalelor EEG prin algoritmi paraleli

Descrierea algoritmului de analiză pentru seturi mari de date

În exemplele enunțate anterior pentru studiul coerenței dintre semnale, algoritmii folosiți au fost aplicați pe un număr restrâns de canale (până în 5 canale) și un număr mic de eșantioane.

Deoarece semnalele simulate au fost generate în așa fel încât să se pună evidenția coerența în mod puternic, numărul mic de canale a putut fi justificat. În schimb, în situațiile reale, un număr mic de canale de EEG nu oferă nivelul de precizie necesar pentru extragerea de informație specializată utilă.

După cum a fost menționat anterior, în studiile clinice se folosesc, în general, peste 20 de electrozi pentru măsurarea EEG. Se poate demonstra că timpul destinat evaluării coerenței dintre semnale cu ajutorul metodelor precizate anterior crește pătratic cu dimensionalitatea sistemului. Astfel, dacă pentru un sistem format din 5 procese a câte 1024 de eșantioane, durata de evaluare a coerenței durează 60 de secunde, la creșterea numărului de canale la 20, timpul necesar calculului depășește 10 minute.

În cazul bazelor de date reale analizate în această lucrare, numărul canalelor este de 32 și 64, cu 8064 și 20000 de esantioante per canal. Este evident faptul că aplicarea unui algoritm de tip GCI pe sisteme atât de mari este foarte dificilă, din două considerente principale: necesarul de memorie al sistemului de calcul și durata procesului de calcul.

Spre deosebire de metodele neparametrice studiate, PDC și GCI se bazează pe calculul matricilor de coeficienți MVAR, fapt ce limitează considerabil dimensionalitatea sistemului analizat.

În general, în situațiile reale, este necesară extragerea de informații dintr-un sistem în timp real. Mai mult decât atât, se presupune că, acolo unde este necesară o astfel de prelucrare, resursele de calcul pot fi limitate, făcând imposibilă o analiză cu GCI a unui sistem EEG format dintr-un număr mare de canale.

O posibilitate de rezolvare a problemelor menționate este trierea canalelor dintr-un sistem cu dimensionalitate mare, prin metode cu o complexitate redusă. Astfel, în lucrarea de față se propun o serie de pași pentru reducerea considerabilă a numărului de procese printr-o metodă neparametrică, care să reducă timpul de calcul considerabil.

În cele ce urmează se va prezenta algoritmul de analiză pentru situațiile descrise. Principalii pași ai analizei sunt:

Aplicarea algoritmului de coerență directă pentru fiecare pereche de semnale din sistemul studiat.

Împărțirea pe benzi de frecvență a nivelelor de coerență.

Sortarea și ordonarea datelor obținute.

Determinarea canalelor de interes pe baza datelor sortate

Aplicarea algoritmilor de analiză complexă pe canalele considerate la pasul anterior ca fiind „de interes”.

Stabilirea pragurilor de confidență și evaluarea automată a conectivității direcționate dintre semnale.

Aplicarea algoritmului de coerență directă. Acest pas este cel care reduce considerabil timpul de analiză al sistemului studiat. Deoarece algoritmul de coerență directă nu ține cont de direcția interacțiunii, timpul de calcul este direct proporțional cu

(31)

unde N este numărul de canale al sistemului, și se prespupune a fi evoluția timpului de calcul în cazul metodelor cu detecție direcționată a conectivității.

În baza analizei inițiale se va obține o matrice de dimensiune

(32)

Unde q este numărul de perioade de timp in care a fost împărțit semnalul, iar este numărul de puncte de pe axa frecvențelor. Împărțirea pe perioade de timp a semnalului este necesară din cauza caracterului ne-staționar al semnalelor EEG. Astfel, analiza trebuie realizată pe eșantioane de semnal de-ajuns de scurte pentru a se considera ipoteza de staționaritate.

Împărțirea pe ritmuri EEG. Pentru ușurarea procesului de sortare și pentru reducerea influenței erorilor, se va realiza o mediere a valorilor din benzile de frecvență a ritmurilor EEG. Astfel, matricea de analiză se va reduce la:

(33)

unde este numărul de benzi de frecvență.

Sortarea și ordonarea datelor obținute se realizează cu reducerea dimensionalității matricii, prin extragerea maximului de pe fiecare bandă de frecvență, din toate eșantioanele de semnal q. De exemplu, pentru setul de date corespunzător canelelor x și y, este găsită o valoare care reprezintă maximul de pe matricea asociată acestora:

(34)

unde reprezintă matricea de conectivitate asociată canalelor X si Y.

Cu matricea obținută anterior se va trece la o matrice de dimensiune () x 3, unde prima coloană reprezintă nivelul de conectivitate calculat anterior, iar celelălalte două coloane sunt indexul canalelor considerate.

Determinarea canalelor de interes se realizează prin extragerea din matricea obținută anterior a liniilor cu o valoare de conectivitate peste o valoare prag. Dacă valoarea prag stabilită este prea permisivă, se va limita numărul de perechi de canale considerate.

Aplicarea algoritmilor de analiză se realizează ca în exemplele menționate anterior. În studiul lucrării curente, se va utiliza Granger Causality Index pentru analiza canalelor de interes.

Stabilirea pragurilor de confidență si evaluarea automată a conectivității. Acest pas are rolul de dezambiguizare a datelor obținute la pasul anterior. Astfel, se va stabili un nivel de confidență și pe baza acestuia, un algoritm de detecție automată a conectivității. Pentru diferențierea conectivității de trecerile scurte peste pragul de siguranță, se va impune o perioadă prestabilită de rezistență deasupra pragului pentru validarea conectivității.

Accelerarea procesării prin paralelizarea algoritmilor

Ritmul de îmbunătățiri în viteză de procesare secvențială a sistemelor de calcul avansate a început să încetinească în ultimii ani. Deși soluțiile multi-scalare și vectoriale continuă să avanseze, ritmul este prea lent pentru adaptarea în prelucrarea statistică a datelor necesare în aplicații complexe cum ar fi prelucrarea informațiilor cerebrale. Pentru aceste tipuri de date, se pretează procesarea pe sisteme de calcul masiv-paralele, tocmai datorită arhitecturii specifice de tip multi-nucleu.

Procesoare de tip GPGPU și NVIDIA CUDA

GPU sau graphics processing unit este un tip de processor dedicat pentru a manipula elemente grafice (pixeli, texeli – texture element) cu scopul de a ușura afișarea acestora pe ecran. Acestea au fost introduse la începutul anilor ’90 pentru a permite îmbunătățirea performanțelor grafice pe stațiile de lucru.

Acceleratoarele grafice moderne sunt folosite în toate gamele de aplicații, de la stații de lucru avansate până la telefoane mobile și tablete media. Structura masiv paralelă le face mult mai eficiente decât procesoarele de uz general în aplicațiile unde prelucrarea blocurilor de date se poate face în paralel (e.g. decodarea de conținut video, simulări științifice, randare 3D etc.) .

Arhitectura funcțională a plăcilor video este relativ simplă și este reprezentată de o matrice de nuclee de procesare (Shader Core), cărora le sunt distribuite “sarcini” de către un thread dispatch unit (unitate de ordonare și punere în procesare a firelor de execuție), asemănător ca la CPU multi-nucleu. În prezent, plăcile video moderne au ajuns să aibă și până la 2880 de nuclee de procesare pe același cip (Nvidia GTX TITAN).

Acronimul GPGPU provine de la General Purpose Graphics Processing Unit și denumește adaptarea recentă a acceleratoarelor grafice pentru a suporta prelucrări de date comandate de procesorul principal al sistemului. Inițial, nucleele computaționale ale unui procesor video erau foarte rudimentare, destinate unor operații foarte simple. În timp, odată cu integrarea hardware a posibilităților de randare 3D, anti-aliere etc; funcționalitatea acestora a crescut considerabil și converg către setul de instrucțiuni matematice ale procesoarelor de uz general.

Cu toate acestea, un procesor grafic duce lipsă de elementele fundamentale necesare sistemelor de calcul de uz general. Astfel, acesta nu poate fi folosit de sine stătător, ci mai degrabă, ca un coprocesor matematic pentru accelerarea aplicațiilor dedicate. Deși realizarea unei aplicații care să beneficieze de avantajele procesorului video este mai complicată, în anumite domenii (e.g. cercetare, simulări moleculare, predicții și corelații statistice etc.) utilizarea sa aduce o îmbunătățire considerabilă a vitezei de execuție, permițând dezvoltarea programelor complexe de analiză în timp real.

Din punct de vedere al performanțelor în ceea ce privește un singur fir de execuție, GPU-urile nu pot lua locul CPU-urilor din cauza nucleelor de calcul rudimentare. De altfel, vitezele la care lucrează procesoarele video sunt mult inferioare celor de la procesoarele de uz general, de obicei maxim 1 GHz (mai mari pentru GPU integrate), fiind inadecvate rulării algoritmilor secvențiali (e.g. compilatoare). Cu toate acestea, în cazul procesării masiv-paralele, dezvoltarea noii paradigme de calcul a dus la creșterea considerabilă a puterii de calcul totale ale unui procesor video, mai ales datorită numărului foarte mare de nuclee de procesare.

În figura următoare este reprezentată schema simplificată a unui sistem heterogen de procesare, format din procesorul principal, memoria fizică a sistemului, procesorul grafic și memoria sa dedicată. Cu cifre sunt indicați principalii pași din execuția unui program secvențial-paralel într-un sistem ca cel descris anterior.

Figura 5.1. Sursa: [15] Structura funcțională a unui sistem heterogen de procesare CPU-GPU.

La baza unui nucleu de procesare video se află o structura vectorială ce permite rularea în regim SIMD, fapt ce favorizează lucrul cu vectori, matrici, texturi etc.

CUDA (Compute Unified Device Architecture) este o platformă de integrare software-hardware dezvoltată de NVIDIA și implementată în toată gama de acceleratoare grafice ale companiei. Această tehnologie oferă programatorului acces direct la setul de instrucțiuni virtual și la memoria de performanță ale plăcilor video.

În ultima decadă, procesoarele grafice ale companiei NVIDIA au luat locul procesoarelor convenționale în clusterele de super-calculatoare, mai ales datorită posibilității lor deosebite de a accelera procesarea în cazul simulatoarelor și a algoritmilor de analiză statistică.

Paralelizarea prelucrării

Algoritmii de prelucrare a semnalelor EEG au la bază calcule statistice, operații pe matrici și transformări punctuale pe datele prelucrate. Se va încerca implementarea unor algoritmi de prelucrare masiv-paralelă pe plăcile video, pentru accelerarea.

În vederea accelerării vitezei de prelucrare a semnalelor, se va propune o metodă de paralelizare a calculului cross-densității spectrale pentru un număr mare de canale. Datorită caracterului nestaționar al semnalelor EEG, semnalele utilizate vor fi împărțite în eșantioane de dimensiuni reduse pentru a se putea considera cvasi-staționaritatea lor.

Astfel, structura datelor de intrare va avea forma unei matrici tridimensionale N*m*q, unde N – numărul de canale analizate, m – numărul de bucăți în care sunt împărțite semnalele și q – numărul de eșantioane per bucată de semnal.

Secvența de lucru se va desfășura în trei etape:

Citirea și structurarea datelor de intrare: după citirea datelor din fișier, acestea vor fi împărțite în m segmente egale. Este important ca datele să fie segmentate egal, deoarece dimensiuni diferite ale vectorilor de intrare cauzează timp suplimentari de preparare a prelucrării.

Aplicarea TFR (Transformata Fourier Rapidă): Datele inițiale de intrare vor fi realiniate într-un vector, dimensiunea setului de date rămânând același. Algoritmul nativ de FFT implementat în librăria CuFFT, proprietara Nvidia, va executa prelucrarea în paralel pe tot setul de date, transferul către placă video realizându-se o singură dată pentru datele de intrare și încă o dată pentru datele de ieșire. Observație: Este important de luat în considerare dimensiunea memoriei fizice instalate pe placă grafică. Dacă setul de date de intrare este prea mare este necesară o segmentare suplimentară a datelor, pentru evitarea erorilor de tip „memorie insuficientă”.

Calculul paralel al cross-densității spectrale: După obținerea vectorului de ieșire cu eșantioanele FFT ale semnalelor, se va calcula cross-densitatea spectrală de putere pentru fiecare pereche de segmente de date în parte. În acest caz, se vor realiza transferuri multiple către și dinspre placă grafică, operații efectuate pentru fiecare pereche de segmente. Datorită lipsei de simetrie în stabilirea perechilor de date al algoritmului bidirecțional, paralelizarea se va realiza la nivel de calcul punctual al cross-densității spectrale și nu pe tot setul de date.

Pentru evaluarea performanțelor și comparația cu tipul de procesare secvențială s-au elaborat coduri pentru prelucrarea paralelă și secvențială a metodei descrise anterior. Aceste programe vor fi descrise, pentru simplitate, în pseudocod, în secțiunea următoare:

Pseudocod secvential

citeste canale[nr_canale][dimensiune_canal]

pentru i=1..nr_canale

pentru j=1.. nr_esantioane_canal /dimensiune_segment_canal

fft_segment[i][j]=calculeaza_fft_secvential(i,j, dim_segment_canal)

sfarsit bucla

sfarsit bucla

pentru i = 1..nr_segmente_fft

pentru j = (i+1)..nr_segmente_fft

pentru k = 1..dimensiune_segment

cpsd[i][j][k] = fft_segment[i][k] * conjugat(fft_segment[j][k])

cpsd[i][j][k] = val_absoluta(cpsd[i][j][k])

cpsd[i][j][k] = transforma_in_decibeli(cpsd[i][j][k])

sfarsit bucla

sfarsit bucla

sfarsit bucla

Pseudocod paralel

citeste canale[nr_canale][dimensiune_canal]

segmente[nr_canale][nr_segmente_canal][dimensiune_segment] = segmentare(canale)

transfera_in_memoria_video(segmente)

fft_segmente = fft_paralel_pe_segmente(segmente)

pentru i = 1..nr_segmente_fft

pentru j = (i+1)..nr_segmente_fft

transfera_in_memoria_video(fft_segmente[i])

transfera_in_memoria_video(fft_segmente[k])

cpsd[i][k]=mul_vect_GPU(fft_segmente[i], fft_segmente[k])

cpsd[i][k]=transforma_in_dB_GPU(cpsd[i][k])

transfera_in_memoria_fizica(cpsd[i][k])

sfarsit bucla

sfarsit bucla

În cazul programului pseudo-cod destinat implementării în CUDA, funcțiile terminate în particula „GPU” menționează execuția pe procesorul grafic și formează partea paralelizabilă a programului. În rest, toate celelalte instrucțiuni se execută pe procesorul sistemului și formează partea secvențială a aplicației.

Descrierea semnalelor

Descrierea semnalelor simulate

Semnalele utilizate în lucrare au fost obținute cu ajutorul următoarelor modele MVAR:

Modelul Chavez:

Modelul Gourevitch

Modelul Shelter01

Figura 6.1. Model pentru evaluarea metodei hibride

Descrierea datelor reale

Baza de date „EEG Motor Movement/Imagery Dataset” – PhysioBank

Baza de date din care au fost extrase datele pentru analiză, „EEG Motor Movement/Imagery Dataset” reprezintă un set de date format din aproximativ 1500 de înregistrări EEG multi-canal cu durate de unul sau două minute. Acest set de date a fost obținut cu ajutorul a 109 voluntari.

Fiecare subiect a realizat 14 seturi de activități formate din acțiuni motorii atunci când a fost supus unor stimuli vizuali. Măsurarea semnalelor s-a efectuat cu 64 de electrozi EEG, folosing sistemul BCI2000 și o rată de eșantionare de 160 Hz.

Cele 14 înregistrări experimentale s-au bazat pe activitățile:

– Înregistrare-referință pentru evaluarea semnalelor (cu ochii deschiși) – un minut.

– Înregistrare-referință pentru evaluarea semnalelor (cu ochii închiși) – un minut.

– Experimentul 1.Un stimul vizual sub formă unui punct apare din partea stânga sau din partea dreapta pe un ecran. La apariția țintei de interes, subiectul încleștează și descleștează pumnul corespunzător direcției din care a apărut stimulul. După dispariția acestuia, subiectul se relaxează.

– Experimentul 2. Un stimul vizual sub formă unui punct apare din partea stânga sau din partea dreapta pe un ecran. La apariția țintei de interes, subiectul își imaginează că încleștează și descleștează pumnul corespunzător direcției din care a apărut stimulul. După dispariția acestuia, subiectul se relaxează.

– Experimentul 3. Un stimul apare din partea de sus sau de jos a ecranului. Subiectul încleștează și descleștează ambii pumni (dacă țintă apare de sus) sau ambele picioare (dacă țintă apare de jos). După dispariția stimulului, subiectul se relaxează.

– Experimentul 4. Un stimul apare din partea de sus sau de jos a ecranului. Subiectul își imaginează că încleștează și descleștează ambii pumni (dacă țintă apare de sus) sau ambele picioare (dacă țintă apare de jos). După dispariția stimulului, subiectul se relaxează.

– Experimentul 1

– Experimentul 2

– Experimentul 3

– Experimentul 4

– Experimentul 1

– Experimentul 2

– Experimentul 3

– Experimentul 4

Montajul. Semnalele au fost măsurate cu ajutorul a 64 de electrozi poziționați conform sistemului 10-10 (excepție făcând pozițiile Nz, F9, F10, FT9, FT10, A1, A2, TP9, TP10, P9 și P10).

În lucrarea de față, pentru evidențierea anumitor caracteristici de performanță a algoritmilor de interes, s-au extras măsurătorile subiectului cu numărul 1 din baza de date. Pe acest subiect, au fost analizate experimental 6 experimente din lista prezentată, după cum urmează:

3 încercări pentru experimentul 1 (încercările 3, 7 și 11).

3 încercări pentru experimentul 2 (încercările 4, 8 și 12).

Baza de date „DEAP”

Acest set de date a fost realizat în vederea analizei emoțiilor cu ajutorul corelațiilor dintre semnalele EEG, semnale fiziologice și scurte filme video vizualizate de către subiecți. Din totalitatea măsurătorilor efectuate pe această baza de date, în lucrarea de față se vor studia semnalele EEG, măsurate cu ajutorul a 32 de electrozi plasați pe suprafață scalpului.

Măsurătorile au fost realizate cu o frecvența de eșantionare de 512 Hz, ca apoi să se realizeze o reducere la 128 Hz pentru datele prelucrate. Semnalele măsurate au fost filtrate cu un filtru trece-bandă (4 – 45 Hz). Un avantaj considerabil al acestei baze de date este eliminarea artefactelor datorate clipirii subiecților prin măsurarea pe canale auxiliare a EOG (electro-oculografie). Rezultatele EOG au ajutat la prelucrarea și filtrarea integrală a semnalelor EEG de impulsurile datorate mișcării oculare.

În cadrul prezentului studiu se va alege aleator un subiect și se vor studia algoritmii de analiză propuși pentru evidențierea conectivității cerebrale. Prin aplicarea algoritmului hibrid de analiză se va încerca reducerea numărului de canale analizat, pentru că apoi să poată fi testată conectivitatea direcționată pe un set mai restrâns de date (maxim 5-6 canale).

Rezultate și discuții

Performanțele implementării paralele

Analiza performanțelor pentru variantele secvențială și paralelă a calculului cross-densității spectrale se va realiza în două moduri:

Aplicația 1. Pentru un sistem cu un număr fix de 64 de canale se va varia dimensiunea temporală a sistemului. La o rată de eșantionare teoretică de 250 Hz, testele se vor baza pe dimensiuni ale sistemului de la 10 la 400 de secunde.

Aplicația 2. În acest caz, se menține fixată lungimea semnalelor (50000 de eșantioane – 200 de secunde) și se variază numărul de canale în 6 puncte (16, 32, 48, 64, 96 și 128 de canale).

Pentru o mai bună comparație, măsurarea duratelor celor două aplicații începe numai în momentul execuției propriu-zise a algoritmilor descriși anterior în pseudocod și nu va cuprinde și simularea/citirea datelor din fișier. Este de remarcat că, în cazul algoritmilor rulați pe procesoarele CUDA, timpii de execuție nu sunt direct proporționali cu volumul computațional, ci depind puternic de organizarea datelor și de transferul dintre cele două memorii (memoria dedicată plăcii video și memoria sistemului). Măsurarea duratei de execuție pe GPU va cuprinde și acești timpi de transfer. Rezultatele duratelor de timp pentru execuțiile cu datele de intrare descrise anterior se regăsesc în tabelele 2 și 3.

Tabel 7.1. Datele obținute în urma rulării experimentale din aplicația 1

Tabel 7.2. Datele obținute în urma rulării experimentale din aplicația 2

Figura 7.1. Comparație grafică a imbunătațirii vitezei de execuție

În cazul ambelor aplicații se poate observă o eficientizare considerabilă a timpului de lucru în cazul metodei de prelucrare paralelă față de prelucrarea secvențială cu ajutorul Matlab.

Prima aplicație a avut că scop evidențierea evoluției timpului de lucru pentru programul paralel, în condițiile în care numărul de transferuri dintre memorii rămâne același (nu se sechimba numărul de bucle din partea secvențială a codului). Astfel, în mod evident, odată cu creșterea dimensiunilor datelor de intrare pentru prelucrarea TFR și a algoritmului de cpsd se observă îmbunătățirea performanțelor acestui tip de procesare (de la 8x – în cazul a 2500 de eșantioane per semnal – la 87x – în cazul a 100000 de eșantioane per semnal). Practic, durata programului începe să intre pe un trend ascendent după trecerea de 20000 de eșantioane per semnal, ceea ce confirmă overhead-ul necesar lansării unui program în CUDA, precum și a latențelor crescute de transfer între memorii.

În al doilea exemplu, se modifică modifică dimensiunea părții secvențiale ale programului și, astfel, a numărului de transferuri de date dintre memorii. Deși durata de calcul ajunge la aproape 8 secunde pentru 128 de canale, îmbunătățirea față de Matlab rămâne de to aprox 80x.

Analiza performanțelor în cazul implementării clasice. Comparație PDC-GCI

Acest capitol abordează analiză performanțelor algoritmilor prezentați anterior (cap. 3). Rezultatele s-au obținut pe sisteme de semnale reale și simulate, generate pe baza modelului autoregresiv.

Rezultate obținute pe baza semnalelor simulate

Se vor prezenta în ordinea rezultatelor obținute în urmă analizei sistemelor autoregresive propuse în capitolul 4. Următorii algoritmi de analiză vor fi prezentați:

PDC-I (PDC cu coeficienți invariabili în timp)

GCI

Procesele analizate au fost generate cu un număr de 2148 de eșantioane. Primele 100 de eșantioane de semnal sunt înlăturate pentru a permite stabilizarea conectivităților matematice dintre procese, prin uitarea factorului eroare de la începutul generării semnalelor. Pentru fiecare dintre sistemele analizate se consideră o frecvența de eșantionare de 256 Hz. Astfel, în cazul PDC-I, unde reprezentarea se realizează în domeniul frecvența, frecvența maximă este de 128 Hz.

Se consideră că sistemele simulate au o durata de 8 secunde, durata reprezentată pe axa temporală în cazul analizei cu GCI. Graficele reprezentate conțin N*N subgrafice reprezentative pentru influențele direcționate dintre procesele prezentate. Fiecare subgrafic este denumit cu ch x->ch y pentru evidențierea sensului de influență. Pe diagonalele principale ale graficelor, în cazul PDC-I, este reprezentată puterea spectrală a semnalului corespunzător.

În cazul algoritmului GCI, valorile de pe diagonală principala sunt nule, deoarece calculul valorii interacției unui singur canal nu are sens.

Deoarece algoritmul de coerență directă nu ține cont de sensul interacțiunii, reprezentarea interacțiunii se realizează între canale, făcând abstracție de direcție. Se vor reprezenta un număr de grafice. La fel ca și în cazul PDC-I, analiza se realizează în frecvența, frecvența maximă din spectru fiind egală cu fs/2 (jumătatea frecvenței de eșantionare a semnalului). Valorile pentru coerență directă au fost calculate folosind un număr de fs puncte în care s-a efectuat transformarea de TFR (Transformată Fourier Rapidă). O metodă simplă de comparație a coerenței dintre semnalele analizate este medierea valorilor din spectru. Deși, în practică este necesară împărțirea pe benzi de frecvența, pentru sistemele abordate, pentru simplificarea observațiilor, s-a ales medierea întregului spectru. Reprezentarea mediei se realizează cu o linie roșie, orizontală de-alungul întregului grafic.

Metodele propuse constituie avantaje și dezavantaje în funcție de situația analizată. În funcție de modelul liniar al proceselor, se va incerca diferențierea pe cazuri între algoritmii de analiză a conectivității și a cauzalității. La baza construcției modelelor stă metoda arborilor direcționați.

Metoda arborilor direcționați este principiu de analiză a sistemelor multivariate pentru ușurarea examinării conectivităților dintre noduri (procese) pe baza săgeților (interacțiunile dintre procese detectate de tehnicile de analiză aplicate).

Pentru evaluarea performanțelor metodelor se vor folosi pentru început 4 tipuri de sisteme simulate:

Un sistem autoregresiv din 4 canale

Un sistem format din 3 procese de zgomot alb

Un sistem autoregresiv cu parametri de influență variabili in timp

Exemplul 1

În primul exemplu se vor analiza rezultatele obtinute de metodele studiate pe baza unui sistem variat autoregresiv liniar. Scopul exemplului este de a determina performanțele algoritmilor si nu de a caracteriza procesele simulate. Sistemul de la baza acestei analize este descris de graful direcționat din figura de mai jos, unde nodurile reprezinta canalele simulate, iar sagețile, influența dintre acestea. E.g. o săgeata între nodul x1 si x2 sugerează influența de la x1 la x2, iar lipsa săgeții reprezintă absența conectivității direcționate dintre aceste canale. Este de remarcat că fiecare legătura are aceeași pondere in graf (1 sau 0), deoarece scopul exemplului este de a releva prezenta sau absenta conectivității dintre canale. Evoluția proceselor este descrisă de sistemul:

Figura 7.2. GCI pentru modelul multivariat autoregresiv invariant în timp din exemplul 1. Conectivitatea dintre canalele de interes este evidențiată corect prin valorile nenule a GCI.

Figura 7.3. Partial directed coherence.

Se observă că atat PDC, cat si GCI releva rezultate corecte pentru situația analizată. Se observă legaturile dintre canalele 1->2, 2->3 si 2->4. Datorită faptului că semnalele analizate sunt pur staționare, PDC oferă informații utile. Din cauza naturilor sale cauzale, valorile obținute cu cei doi algoritmi nu sunt întotdeauna relevante decat pentru evidențierea conectivității, nu si a puterii de interacțiune dintre canale. În schimb, niciunul dintre algoritmi nu a putut identifica influența indirect ă dintre canalele 1->3 si 1->4.

Exemplul 2

În exemplul anterior a fost testată capabilitatea algoritmilor de a detecta existența conectivității între procesele unui sistem multivariat. În acest exemplu se va testa corectitudinea afisării evaluării informațiilor in cazul a trei procese formate din zgomot alb gaussian de medie 0. Scopul acestui test este de a evidentia situațiile în care algoritmii analizați pot oferi „alarme false” ().

Figura 7.4. Indexul de cauzalitate Granger pentru un sistem 3-dimensional de zgomot alb din exemplul 2. Se observă influențe detectate eronat în cazul ch1->ch2 si ch2->ch3.

Figura 7.5. Partial Directed Coherence pentru un sistem 3-dimensional de zgomot alb din exemplul 2. Nu sunt detectate conectivitații false intre canale.

Exemplul 3

Datorită caracterului non-staționar al semnalelor EEG, este necesară observarea capacității de analiză a algoritmilor studiați pe seturi de date al caror spectru variază în timp. Pentru evidențierea acestor caracteristici, urmatorul test se bazeaza pe un sistem multivariat al carui comportament se schimbă dupa t = T/2 (T = durata totala a proceselor). În acest fel, se simuleaza comportamentul de tip ERP (event-related potential), care presupune stimularea subiectului cu scopul de a observa schimbări în caracteristicile semnalelor studiate.

Inițial, în prima jumatate a semnalului:

La t = T/2:

Figura 7.6. Granger Causality Index. Conectivitatile inter-canale sunt bine evidențiate înainte si după simularea apariției stimulului. Se poate observa momentul apariției conectivității direcționate între canalele ch3->ch2, precum si pentru ch4->ch1. În aceasta situație este posibila observarea cresterii coeficientului in cazul ch2->ch1 unde coeficienții de conectivitate au fost crescuți. Desi greu de evidențiat, se observă o crestere minoră a influenței de la ch1->ch3.

Figura 7.7. Partial Directed Coherence pentru sistemul din exemplul 3. Se observă apariția unui influențe detectate eronat in cazul ch4->ch2.

Modelul Chavez

Acest sistem este format din doua canale si are ordinul p = 1. Interacțiunea este realizată într-un singur sens, de la canalul 1 la canalul 2 (x2<-x1).

a)

b) c)

Figura 7.8. a) Coerența directă pentru modelul Chavez. b) Granger Causality Index (stanga) c) Partial Directed Coherence (dreapta)

Cei trei algoritmi de analiză dau rezultate bune pe modelul Chavez. PDC-I și GCI detectează corect interacțiunea direcționată de la canalul 1 la canalul 2. Algoritmul coerenței directe dă o valoare orientativă pentru modelul propus. Media sub 0.5 (marcată cu linia roșie) se datorează interacțiunii într-un singur sens. Analiza pe porțiuni a MSC arată un grad de conectivitate mai mare la frecvențe medii (50-80 Hz) decât PDC-I.

Modelul Gourevitch

Acest model se diferențiaza de primul prin creșterea ordinului cu o unitate și prin bidirecționalitatea influențelor dintre canale.

Figura 7.9. a) Coerența directă pentru modelul Gourevitch. b) Granger Causality Index (stanga) c) PDC-I (dreapta)

Analiza pe modelul Gourevitch arată influențe directe puternice pentru fiecare algoritm. În cazul PDC-I și GCI se observă reprezentarea conectivității bidirecționale pe diagonala secundară a graficelor. Mai mult decât atât, se observă creșterea valorilor maxime ale coeficienților de conectivitate atât în cazul GCI, cât și în cazul PDC-I.

Datorită influenței bidirecționate, analiza coerenței directe dintre cele două canale dă un grad ridicat de certitudine, media de conectivitate ridicându-se la peste 0.5 din maxim. Mai mult, în acest caz, spre deosebire de modelul Chavez se poate face o delimitare pe benzi de frecvența a spectrului obținut. Din această analiză poate fi evidențiată corespondență cu modelul teoretic. Amplitudinea medie din spectrul coerenței directe se situează în jurul valorii de 0.8, în concordanță cu media coeficienților din matricea ponderilor asociată sistemului simulat.

Modelul Shelter01

Urmǎtorul model analizat este Shelter01, constituit din 5 procese. Ordinul sistemului autoregresiv este 3.

Figura 7.10. a) Coerenṭa directǎ pentru modelul Shelter01. b) GCI pentru modelul Shelter01. c) PDC-I pentru modelul Shelter01.

Pentru modelul Shelter01 se observă o ambiguizare a rezultatelor algoritmului de coerenṭă directă. Astfel, dacă la modelul Gourevitch, puteau fi evidențiate cu certitudine canalele de interes, aici este necesară o analiză mai amănunțită pentru a putea fi extrase informații folositoare în legătură cu interacṭia dintre canale. Cu toate acestea, deși pragul de 0.5 nu este depășit nici măcar în cazul unde influența este bidirecțională (ch4 <-> ch5), prin ordonarea mediilor realizate pe fiecare spectru se pot înlătura canalele care nu prezintă conectivitate (a căror medie spectrală este mult inferioară celorlalte). Astfel, se poate evidenția corect lipsa de interacțiune dintre canalele ch3<->ch4, ch2<->ch5 și ch1<-> ch5.

Analiza cu GCI oferă informații corecte despre interacțiunile direcționate dintre canale. Se observă interacțiuni puternice dinspre canalul 1 înspre canalele 2 și 4, confirmând modelul teoretic pe baza căruia au fost generate semnalele.

Analiza modelului complex Shelter01 oferă informații noi referitoare la tipul de analiză. În acest caz, algoritmul cu cele mai folositoare rezultate este tot Granger Causality Index, deoarece oferă valori ridicate ale coeficienṭilor de conectivitate, crescând certitudinea interacțiunilor. Analiza cu PDC invariant oferă rezultate corecte, iar interacțiunea se poate observa în cadrul jumătății inferioare ale spectrului (până în 60 Hz).

Deși cele mai corecte rezultate au fost oferite de GCI și PDC, creșterea ordinului sistemului la 3 și a dimensionalitatii la 5 canale a mărit semnificativ timpul computațional. Se observă că timpul de lucru depășește cadrul pentru un răspuns real-time pentru acest model (pe anumite sisteme). Algoritmul de coerenṭă directă a fost efectuat pe sistemul de calcul folosit, în sub o secundă, cu o eficienṭă vizibil crescută față de PDC și GCI.

Model pentru evaluarea metodei hibride

In urma analizei pentru seturi mari de date cu metoda propusa in capitolul 3, canalele obṭinute ca fiind cele mai semnificative au fost 1, 2, 3, 4, 5.

După analiza vizuală a graficului de mai sus, se pot observa coerenṭele puternice descrise de modelul teoretic din capitolul 4 în cazul grupului de semnale 1-2-3-4-5. În interacțiunile cu canalele 1 și 2 se pot remarca valori mari ale coerenței în jumătatea inferioară a spectrului (până în 20-25 Hz), în timp ce, pentru canalul 5, amplitudini ridicate ale coeficientului de coerenṭă se află în a două jumătate a spectrului de frecvenṭe (după 32 Hz).

Se poate observa că GCI confirmă interacțiunile din grupul semnalelor 1, 2, 3, 4, 5 – considerate a fi de interes. Marcajul cu linie neagră punctată din fiecare grafic arată decizia de influență dintre canale. Astfel, remarcăm identificarea corectă a conectivitǎṭilor dintre canale.

Beneficiul algoritmului propus, pe date simulate, este reducerea de două ori a dimensionalitatii sistemului de intrare (de la 10 canale la 5 canale) pentru analiza GCI.

Figura 7.11. a) Coerenṭǎ directǎ pentru modelul de evaluare a metodei hibride pentru semnalele determinate de algoritm ca fiind de interes. b) GCI pentru modelul de evaluare a metodei hibride

Rezultate obtinute pe baza semnalelor reale

Pentru bazele de date reale descrise în capitolul 6 s-a realizat selecționarea canalelor de interes conform algoritmului bazat pe coerenṭa directă, propus în capitolul 5.1. Deoarece baza de date DEAP a fost realizată pe încercări unice ale experimentelor, analiză acesteia s-a realizat cu stabilirea numărului de trials la 1.

Graficele de coerenṭă directă prezintă date tridimensionale, amplitudinea din spectrele de coerentă fiind reprezentată prin codul culorilor, iar pe cele două axe: pe axa orizontală – domeniul frecvența, iar pe axa verticală – domeniul timp.

Analiza bazei de date PhysioBank

Figura 7.12 a) Graficele algoritmului de coerenṭǎ directa pentru semnalele considerate a fi de interes, extrase din baza de date PhysioBank b) Graficele GCI pentru semnalele considerate a fi de interes, extrase din baza de date PhysioBank

Setul de date analizat în cazul acestei baze de date are dimensiunea 3 x 64 x 1920 (număr de încercări x număr de canale x eșantioane/canal). Durata secțiunilor analizate este de 12 secunde, la o frecvențǎ de eșantionare de 160 Hz. Reducerea la 12 secunde a porțiunilor de semnal analizate (dintr-un total de aproximativ 120 de secunde) a fost realizată pentru o mai bună evidentiere a performanțelor algoritmului. Analiză acestui set de date are rezultate slabe în cazul selecției canalelor de interes. În figura 7.12 a) se poate observa saturarea consistentă a spectrului de coerenṭă directă în proporție ridicată în cazul interacțiunii dintre canalele 23-24, 23-25 și 23-26. În subgraficele de coerenṭă dintre canalele 25-26 și 35-36, după frecvența de 45 de Hz, unde algoritmul propus nu a realizat normalizarea valorilor, se poate observa o saturare totală a amplitudinilor datorită erorilor de calcul din cazul frecventelor ridicate.

În figura 7.12 b), indicele GCI are valori reduse, iar detecția de conectivitate are un grad scăzut de certitudine pentru toate canalele, excepție făcând detecția realizată în cazul interacțiunii dintre canalele ch1->ch2, unde valoarea mai ridicată a indicelui GCI (aproximativ 0.2) oferă mai multă siguranță asupra corectitudinii interpretării. Rezultatele analizei cu Granger Causality Index evidențiază performanța scăzută a algoritmului de selecție a canalelor de interes.

Analiza bazei de date „DEAP”

Figura 7.13 a) Graficele algoritmului de coerenṭǎ directa pentru semnalele considerate a fi de interes, extrase din baza de date DEAP. b) Graficele analizei cu Granger Causality Index pentru semnalele considerate a fi de interes, extrase din baza de date DEAP

Deoarece baza de date DEAP a fost realizată pe încercări unice ale experimentelor, analiza acesteia s-a realizat cu stabilirea numărului de trials la 1. Setul de date de intrare are dimensionalitatea (număr canale) x (număr eșantioane/canal).

Aplicarea alogritmului hibrid pe acest set de date a înregistrat rezultate vizibil mai bune decât în cazul bazei de date PysioBank. Se poate observă în figura 7.13 a) o împrăștiere mai largă a valorilor din spectrul coerenței directe. Acest fapt ajută la o mai bună interpretare a valorilor de conectivitate bidirecțională dintre canale. În intervalul 0 – 4 Hz se poate observa efectul filtrării din prelucrarea bazei de date și a algoritmului propus. Reducerea la 0 a valorilor (din spectrul de coerenṭă directă) s-a realizat din cauza valorilor reduse din spectrul benzii de frecvența 0 – 4 Hz a semnalelor analizate, pentru a elimina erorile de calcul.

Selecția realizată pe baza coerenței directe a restrâns dimensionalitatea sistemului la 5 canale analizate, canalele 1, 2, 26, 29 și 31. Acestea corespund senzorilor plasați în pozițiile Fp1, AF3, FC2, AF4 și Fz.

În figura 7.13 b), analiza canalelor de interes arată multiple conectivitǎṭi direcționate cu valori ale coeficientului GCI mai ridicate decât în cazul bazei de date PhysioBank, oferind un grad de certitudine crescut asupra metodei de selecție a canalelor de interes. Se poate observă detecția unor influențe repetate mai puternice, de la canalul 31 (Fz) către canalele 26 (FC2) și 29 (AF4), precum și de la canalul 1 (Fp1) la canalele 26 (FC2) și 29 (AF4).

Concluzii si perspective

Din rezultatele prezentate în capitolul 4 și 6, pe baza modelelor propuse și a proceselor simulate se observă faptul că algoritmii de estimare a conectivității reușesc să determine în mod corect interacțiunile dintre două sau mai multe procese analizate.

Indiferent de ordinul sistemului multivariat, metodele parametrice utilizate în analiză au evidențiat un grad mare de certitudine în decizia influențelor și a direcțiilor acestora din sistemele simulate. În cazul metodei de coerenṭă simplă, rezultatele obținute tind să se ambiguizeze odată cu creșterea ordinului și dimensiunii sistemului analizat. Deși gradul de certitudine scade în situațiile descrise anterior și în cazul modelelor cu conectivitate unidirecțională, gradul de coerentă rămâne ridicat și se poate realiza diferențierea între procesele ce prezintă influență reciprocă și procesele independente.

Studiul celor două metode parametrice le caracterizează ca fiind potrivite în funcție de domeniul de aplicabilitate. Astfel, pentru analiza proceselor staționare, cele mai bune rezultate sunt obținute cu metodă PDC-I (Baccala). Adaptată corespunzător, această metodă poate oferi informații folositoare în legătură cu ritmurile EEG studiate. În schimb, pentru analiza semnalelor fiziologice și cerebrale, cu un caracter nestaṭionar, pentru detecția influențelor pe perioade temporale crescute, PDC-I duce la ambiguizarea rezultatelor, cea mai potrivită soluție fiind utilizarea metodei GCI.

Metoda hibridă descrisă în capitolul 5.1 pentru selectarea canalelor EEG de interes a arătat rezultate foarte bune pentru sistemele simulate, putând evidenția conectivitǎṭ corecte între canalele sistemului, fără să introducă în analiză niciunul dintre procesele independente. Cu toate acestea, metoda propusă analizează cele mai puternice nivele de conectivitate, luate în ordine, având probleme în identificarea grupurilor independente. Practic, s-a identificat corect grupul de neuroni (simulaṭi), cu cea mai puternicǎ interacṭiune, lăsându-se deoparte cel de-al doilea grup al căror neuroni prezentau interacțiuni mai slabe.

Aplicarea metodei hibride pe semnale reale a reușit clasificarea seturilor de date și reducerea dimensionalitatii la 5-6 canale, propice analizei cu GCI. Rezultatele obținute au fost mai slabe, datorită caracterului mai ambiguu al semnalelor reale, iar, în cazul bazei de date Physiobank s-au observat erori considerabile în calculul coerenței directe, erori ce au dus la o clasificare și ordonare greșită a canalelor de interes. Rezultatul mai bun din cazul bazei de date DEAP se datorează prelucrării semnalelor cu un filtru trece-bandă și împotriva artefactelor EOG. Deși algoritmul de coerenṭă directă a dat rezultate bune pe semnalele simulate în capitolele 4 și 6, este necesară o adaptare a acestuia pentru a putea fi folosit cu succes pe date cu un caracter cvasi-staționar, cum sunt semnalele reale.

Aplicația de paralelizare a algoritmilor de prelucrare a informațiilor cerebrale a dat rezultate superioare calculului secvențial utilizat în majoritatea programelor ce stau la baza acestei lucrări. Deși exemplul analizat a fost de o complexitate redusă, s-au observat îmbunătățiri de aproape 80 de ori în ceea ce privește timpul de execuție față de programul analog rulat în mediul Matlab. În aplicațiile reale, datorită imposibilității paralelizării algoritmilor în anumite situații, îmbunătățirile ar trebui să se reducă. Un exemplu este calcularea indicilor matricii ponderilor pe baza algoritmului RLS, unde caracterul regresiv al metodei îngreunează paralelizarea calculelor.

Lucrarea de față a realizat o analiză a metodelor de prelucrare a informațiilor cerebrale și a posibilităților de îmbunătățire a acestora, cu scopul de a studia posibilitatea aplicării metodelor parametrice în situațiile de zi cu zi, în viitorul apropiat. Reducerea timpului de lucru este necesar pentru a conferi echipamentelor viitorului posibilitatea de a interacționa în timp real cu utilizatorul. Deși metoda de selecție a semnalelor EEG de interes a dat rezultate mai slabe, aceasta poate fi îmbunătățită printr-o adaptare dinamică a segmentelor de semnal analizate, precum și prin o mai bună filtrare a erorilor datorate amplitudinilor reduse din spectrul de frecvența al semnalelor.

Cele menționate anterior în legătură cu domeniul de aplicabilitate a PDC și GCI, împreună cu paralelizarea algoritmilor și execuția pe coprocesoare matematice dedicate cum sunt acceleratoarele grafice, pot impulsiona răspândirea către uz general a sistemelor EEG, sporind ritmul de evoluṭie al acestei tehnologii emergente.

Bibliografie

[1] Achard, S., Salvador, R., Whitcher, B., Suckling, J., Bullmore, E., “A resilient, low-frequency, small-world human brain functional network with highly connected association cortical hubs” in Journal of Neuroscience, Nr. 26/2006, pp. 63-72.

[2] Blinowska, K. & Durka, P., “Electroencephalography (EEG)” in Wiley Encyclopedia of Electrical and Electronics Engineering, John Wiley & Sons, Varsovia, vol. 6, pp. 394-402.

[3] Erla, S., Faes, L, Tranquillini, E., Orrico, D., Nollo, G., “Multivariate Autoregressive Model with Instantaneous Effect to Improve Brain Connectivity Estimations” in International Journal of Bioelectromagnetism, Vol. II, Nr. 2/2009, pp. 74-79.

[4] Horvath, Z. & Johnston, R., AR(1) Time Series Proccess Econometrics 7590, Univeristy of Utah, Department of Mathematics, accesat la data de 20 iunie 2014, de pe http://www.math.utah.edu/~zhorvath/ar1.pdf.

[5] Massimini, M., Ferrarelli, F., Huber, R., Esser, S.K., Singh, H., Tononi, G., “Breakdown of cortical effective connectivity during sleep” in Science, Nr. 309/2005, pp. 2228-2232.

[6] Moller, E., Schack, B., Arnold, M., Witte, H., “Instantaneous multivariate EEG coherence analysis by means of adapative high-dimensional autoregressive models” in Journal of Neuroscience Methods, Nr. 105/2001, pp/. 143-158.

[7] Rubinov, M. & Sporns, O. , “Complex network measures of brain connectivity: Uses and interpretations”, in Neuro Image, Nr. 52/ 2010, pp. 1059-1069.

[8] Sporns, O., “Brain Conectivity” in Scholarpedia , 2(10):4695, accesat la data de 22 iunie 2014, de pe http://www.scholarpedia.org/article/Brain_connectivity.

[9] Sporns, O., Chialvo, D.R., Kaier, M., Hilgetag, C.C., “Organization, development and function of complex brain networks” in Cognitive Sciences, Vol. 8, Nr. 9/ 2004, pp.418-425.https://bccn2009.org/~triesch/courses/cogs1/readings/sporns04.pdf

[10] Sporns, O., Tononi, G., Kötter, R., “The human connectome: A structural description of the human brain” in PLOS Computational Biology. Nr. 1/2005, pp.245-251.

[11] Koelstra, S., Muhl, C., Soleymani, M., Lee, J-S., Yazdani, A., Ebrahimi, T., Pun, T., Nijholt, A., Patras, I., “DEAP: A Database for Emotion Analysis using Psychological Signals” in IEEE Transactions on Affective Computing, Vol. III, Nr. 1/2012, pp. 18-31.

[12] Teplan, M., “Fundamentals of EEG Measurment” in Measurement Science Review, Nr. 2/2002, Sectiunea 2, pp.1-11.

[13] Trans Cranial Technologies, “10/20 System Positioning”, tDCS manuals, accesat la data de 20 iunie 2014, de pe http://www.trans-cranial.com/local/manuals/10_20_pos_man_v1_0_pdf.

[14] Wang, S.Y. & Tang, M.X., “Exact Confidence Interval for Magnitude-Squared Coherence Estimates” in IEEE Signal Processing Letters, Vol. II, Nr. 3/ 2004, pp.326-329.

[15] CUDA, „Example of CUDA processing flow”, accesat la data de 03 iulie 2014, de pe http://en.wikipedia.org/wiki/CUDA.

[16] Winterhalder, M., Schelter, B., Hessem W., Scwab, K., Leistritz, L., Klan, D., Bauer, R., Timmer, J., Witte, H., “Comparison of linear signal processing techniques to infer directed interactions in multivariate neural systems”, ȋn Signal Processing , Nr. 85/ 2005, pp. 2137-2160 .

ANEXE

Bibliografie

[1] Achard, S., Salvador, R., Whitcher, B., Suckling, J., Bullmore, E., “A resilient, low-frequency, small-world human brain functional network with highly connected association cortical hubs” in Journal of Neuroscience, Nr. 26/2006, pp. 63-72.

[2] Blinowska, K. & Durka, P., “Electroencephalography (EEG)” in Wiley Encyclopedia of Electrical and Electronics Engineering, John Wiley & Sons, Varsovia, vol. 6, pp. 394-402.

[3] Erla, S., Faes, L, Tranquillini, E., Orrico, D., Nollo, G., “Multivariate Autoregressive Model with Instantaneous Effect to Improve Brain Connectivity Estimations” in International Journal of Bioelectromagnetism, Vol. II, Nr. 2/2009, pp. 74-79.

[4] Horvath, Z. & Johnston, R., AR(1) Time Series Proccess Econometrics 7590, Univeristy of Utah, Department of Mathematics, accesat la data de 20 iunie 2014, de pe http://www.math.utah.edu/~zhorvath/ar1.pdf.

[5] Massimini, M., Ferrarelli, F., Huber, R., Esser, S.K., Singh, H., Tononi, G., “Breakdown of cortical effective connectivity during sleep” in Science, Nr. 309/2005, pp. 2228-2232.

[6] Moller, E., Schack, B., Arnold, M., Witte, H., “Instantaneous multivariate EEG coherence analysis by means of adapative high-dimensional autoregressive models” in Journal of Neuroscience Methods, Nr. 105/2001, pp/. 143-158.

[7] Rubinov, M. & Sporns, O. , “Complex network measures of brain connectivity: Uses and interpretations”, in Neuro Image, Nr. 52/ 2010, pp. 1059-1069.

[8] Sporns, O., “Brain Conectivity” in Scholarpedia , 2(10):4695, accesat la data de 22 iunie 2014, de pe http://www.scholarpedia.org/article/Brain_connectivity.

[9] Sporns, O., Chialvo, D.R., Kaier, M., Hilgetag, C.C., “Organization, development and function of complex brain networks” in Cognitive Sciences, Vol. 8, Nr. 9/ 2004, pp.418-425.https://bccn2009.org/~triesch/courses/cogs1/readings/sporns04.pdf

[10] Sporns, O., Tononi, G., Kötter, R., “The human connectome: A structural description of the human brain” in PLOS Computational Biology. Nr. 1/2005, pp.245-251.

[11] Koelstra, S., Muhl, C., Soleymani, M., Lee, J-S., Yazdani, A., Ebrahimi, T., Pun, T., Nijholt, A., Patras, I., “DEAP: A Database for Emotion Analysis using Psychological Signals” in IEEE Transactions on Affective Computing, Vol. III, Nr. 1/2012, pp. 18-31.

[12] Teplan, M., “Fundamentals of EEG Measurment” in Measurement Science Review, Nr. 2/2002, Sectiunea 2, pp.1-11.

[13] Trans Cranial Technologies, “10/20 System Positioning”, tDCS manuals, accesat la data de 20 iunie 2014, de pe http://www.trans-cranial.com/local/manuals/10_20_pos_man_v1_0_pdf.

[14] Wang, S.Y. & Tang, M.X., “Exact Confidence Interval for Magnitude-Squared Coherence Estimates” in IEEE Signal Processing Letters, Vol. II, Nr. 3/ 2004, pp.326-329.

[15] CUDA, „Example of CUDA processing flow”, accesat la data de 03 iulie 2014, de pe http://en.wikipedia.org/wiki/CUDA.

[16] Winterhalder, M., Schelter, B., Hessem W., Scwab, K., Leistritz, L., Klan, D., Bauer, R., Timmer, J., Witte, H., “Comparison of linear signal processing techniques to infer directed interactions in multivariate neural systems”, ȋn Signal Processing , Nr. 85/ 2005, pp. 2137-2160 .

ANEXE

Script SimulatedModels.m

close all

Chavez

GenerateSignals

PlotThemAll

Gourevitch

GenerateSignals

PlotThemAll

Shelter01

GenerateSignals

PlotThemAll

Script GenerateSignals.m

regrOrd = size(weighted,3);

for k=1:nTrials

vector=zeros(channelNo,maxLen);

for i = 5:maxLen

for ch = 1:channelNo

%chInf – influence channel

for infCh = 1:channelNo

for order = 1: regrOrd

% Add an influence from channel [infCh] from [order]

% samples ago

vector(ch,i) = vector(ch,i) + vector(infCh, i – order) * weighted(infCh, ch, order);

end

end

vector(ch,i) = vector(ch,i) + normrnd(0,1);

end

end

%let's discard the first points of this IIR models–> we don't need to see

%the random effect of the series generation (discard first 100 points)

out(k,:,:) = vector(:,discarded+1:end);

end

signals = out;

Script Chavez.m

clear all

algName = cell(3,1);

algName{1} = 'Magnitude squared coherence';

algName{2} = 'Granger Causality Index';

algName{3} = 'Partial Directed Coherence';

modelName = 'Chavez';

channelNo = 2;

nTrials = 1;

fs = 256;

maxLen = 2148;

discarded = 100;

esLen = maxLen – discarded;

out = zeros(nTrials, channelNo, esLen);

weighted = zeros(2,2,1);

% Ch1 influences channels [2]

weighted(1,1,1) = 0.6;

% Ch2 influences channels [3 4]

weighted(2,1,1) = 0.5;

weighted(2,2,1) = 0.6;

Script Gourevitch.m

clear all

algName = cell(3,1);

algName{1} = 'Magnitude squared coherence';

algName{2} = 'Granger Causality Index';

algName{3} = 'Partial Directed Coherence';

modelName = 'Gourevitch';

channelNo = 2;

nTrials = 1;

fs = 256;

maxLen = 2148;

discarded = 100;

esLen = maxLen – discarded;

out = zeros(nTrials, channelNo, esLen);

weighted = zeros(2,2,1);

% Ch1 influences channels [1 2]

weighted(1,1,1) = 1.3434;

weighted(1,1,2) = -0.9025;

weighted(1,2,1) = -0.8;

% Ch2 influences channels [1 2]

weighted(2,1,1) = -0.9;

weighted(2,2,1) = -1.05;

weighted(2,2,2) = -0.85;

Script Shelter01.m

clear all

algName = cell(3,1);

algName{1} = 'Magnitude squared coherence';

algName{2} = 'Granger Causality Index';

algName{3} = 'Partial Directed Coherence';

modelName = 'Shelter01';

channelNo = 5;

nTrials = 1;

fs = 256;

maxLen = 2148;

discarded = 100;

esLen = maxLen – discarded;

out = zeros(nTrials, channelNo, esLen);

weighted = zeros(2,2,1);

% Ch1 influences channels [1 2]

weighted(1,1,1) = 1.34;

weighted(1,1,2) = -0.9025;

weighted(1,2,2) = 0.5;

weighted(1,3,3) = -0.4;

weighted(1,4,2) = -0.5;

% Ch4 influences channels [4 5]

weighted(4,4,1) = 0.35;

weighted(4,5,1) = -0.25;

weighted(5,4,1) = 0.35;

weighted(5,5,1) = 0.35;

Script PlotThemAll.m

figure('name',[modelName ' – ' algName{1}]);

counter = 1;

for stCh = 1:channelNo

for toCh = stCh+1:channelNo

mcxy = mean(mscohere(squeeze(signals(1,stCh,:)),squeeze(signals(1,toCh,:)),[],[],fs,fs));

subplot(5,ceil(sum(1:channelNo-1)/5), counter),

mscohere(squeeze(signals(1,stCh,:)),squeeze(signals(1,toCh,:)),[],[],fs,fs),

axis([0 fs/2 0 1]),

axis 'auto y',

title(['ch', int2str(stCh) ,'<->ch', int2str(toCh)]);

hold on

plot(ones(fs/2,1)*mcxy,'-r', 'LineWidth', 2);

counter = counter + 1;

end

end

gc = GCI_calc(4,signals,0.01);

figure('name',[modelName ' – ' algName{2}]);

minimum = min(min(min(gc)));

for i=0:channelNo^2-1

stCh = mod(i,channelNo)+1;

toCh = fix(i/channelNo)+1;

subplot(channelNo,channelNo,i+1),

plot(gc(:,toCh, stCh),'-b', 'LineWidth', 2),

axis([0 esLen minimum 1]),

title(['ch', int2str(stCh) ,'->ch', int2str(toCh)]);

end

signals = squeeze(signals(1,:,:));

pdc_analysis_TVA

pdc = c.pdc;

figure('name',[modelName ' – ' algName{3}]);

for i=0:channelNo^2-1

stCh = mod(i,channelNo)+1;

toCh = fix(i/channelNo)+1;

subplot(channelNo,channelNo,i+1),

plot(squeeze(pdc(toCh, stCh,:)),'-b', 'LineWidth', 2),

axis([0 128 0 1]),

title(['ch', int2str(stCh) ,'->ch', int2str(toCh)]);

end

Script simulare Hibrida

close all

% clear all

clc

signals = generateFromMatrix10ch(2148,10,5,100);

process_data(signals, 0.01, 128);

Script generateFromMatrix10ch.m

function out = generateFromMatrix10ch(channelSize,nChannels,nTrials,discarded);

out = zeros(nTrials, nChannels, channelSize – discarded);

weighted = zeros(10,10,5);

% Ch1 influences channels [1 3]

weighted(1,1,1) = 0.3;

weighted(1,3,1) = 0.6;

% Ch2 influences channels [1 2]

weighted(2,1,4) = 0.65;

weighted(2,2,1) = 0.6;

% Ch4 influences channels [2 4 4]

weighted(4,2,5) = 0.6;

weighted(4,4,1) = 1.2;

weighted(4,4,2) = -0.7;

%Ch5 influences channels [2 3]

weighted(5,2,1) = 0.8;

weighted(5,3,5) = -0.9;

%Ch7 influences channels [2 6 9]

weighted(7,6,4) = -0.5;

weighted(7,9,4) = 0.1;

%Channels 3, 6, 8, 9 and 10 don't have influences on the system

%Channels 8, 9 and 10 are completely isolated (not influenced)

regrOrd = size(weighted,3);

for k=1:nTrials

vector=zeros(nChannels,channelSize);

recursiveWeghts = weighted;

for i = 6:channelSize

for ch = 1:nChannels

%chInf – influence channel

for infCh = 1:nChannels

for order = 1: regrOrd

% Add an influence from channel [infCh] from [order]

% samples ago

vector(ch,i) = vector(ch,i) + vector(infCh, i – order) * recursiveWeghts(infCh, ch, order);

end

end

vector(ch,i) = vector(ch,i) + normrnd(0,1);

end

end

%let's discard the first points of this IIR models–> we don't need to see

%the random effect of the series generation (discard first 100 points)

out(k,:,:) = vector(:,discarded+1:end);

end

Script process_data.m

function process_data(channels, adapCon, sampling_rate)

trials = size(channels,1);

channelNo = size(channels,2);

esLen = size(channels,3);

discarded = 100;

order = 10;

plotStartValue = 1;

esLength = 1/sampling_rate;

duration_in_sec = esLen/sampling_rate;

% recompute even if nothing's changed (0 – false, 1 – true)

forceCompute = 1;

% if computation is necessary (e.g. variables sizes have changed)

compute = 0;

% let's not waste computation time if we want to do post-processing

% search for all the variables

variables = who;

if size(variables,1) < 10 || forceCompute == 1

compute = 1;

% channels = generateFromMatrix(esLen,channelNo,trials,discarded);

else

if ~isequal(size(channels),[trials channelNo esLen])

compute = 1;

% channels = generateFromMatrix(esLen,channelNo,trials,discarded);

end

end

% writeToFile(channels);

% get values for a spectrogram

ordered_ch = signalCorrelationMsCohere(channels, sampling_rate);

channels = channels(:, ordered_ch, :);

channelNo = size(channels,2);

if compute == 1

gc = GCI_calc(order,channels,adapCon);

end

% plot the raw GC results

figure('name','GCI plots');

for i=0:channelNo^2-1

stCh = mod(i,channelNo)+1;

toCh = fix(i/channelNo)+1;

subplot(channelNo,channelNo,i+1),

hold on

axis([0 esLen-plotStartValue -0.5 0.5]);

plot(gc(:,toCh, stCh)),

title(['ch', int2str(stCh) ,'->ch', int2str(toCh)]);

end

% plot the signals

figure('name','signal plots');

for i=1:channelNo

subplot(channelNo,1,i),plot(squeeze(channels(1,i,:)));

end

confRanges = zeros(trials, size(gc,1), size(gc,2), size(gc,3));

ordered = confRanges;

conf_high = zeros(esLen, channelNo, channelNo);

conf_low = conf_high;

gcThreshArrays = zeros(trials,size(gc,1),size(gc,2),size(gc,3));

upperThresh = zeros(size(gc,1),size(gc,2),size(gc,3));

% Compute thresholds and confidence range

if compute == 1

t = cputime;

progressBar = waitbar(0,'Step 1/3 Computing randomized thresholds');

for k = 1:channelNo

for i = 1:trials

permuted = channels;

% do the channel's permutation

vector = squeeze(permuted(i,k,:));

permIndexes = randperm(numel(vector));

vector = vector(permIndexes);

permuted(i,k,:) = vector;

alteredGc = GCI_calc(order, permuted(i,:,:), adapCon);

gcThreshArrays(i,:,k,:) = alteredGc(:,k,:);

ETR = (cputime – t)/((k-1)*trials+i-1)*((trials*channelNo) – ((k-1)*trials+i-1));

waitbar(((k-1)*trials+i-1)/(trials*channelNo),progressBar,['Step 1/3 Computing randomized thresholds' ' ETR: ' int2str(ETR) ' seconds']);

% disp(['For channel ',int2str(k),'-> Calculating trial ',int2str(i),' from ',int2str(trials),'..']);

end

end

close(progressBar)

for j = 1:channelNo

gcThreshArrays(:,:,k,j) = sort(gcThreshArrays(:,:,k,j));

end

if trials > 1

aux = channels;

ETR = 0;

progressBar = waitbar(0,'Step 2/3 Computing confidence ranges');

t = cputime;

for i = 1:trials

trialRandomized = randi(trials, 1, trials);

aux = channels(trialRandomized, :, :);

confRanges(i, :, : , 🙂 = GCI_calc(order, aux, adapCon);

ETR = (cputime – t)/(trials-1)*(trials-i);

waitbar(i/trials,progressBar, ['Step 2/3 Computing confidence ranges ' ' ETR: ' int2str(ETR) ' seconds']);

end

close(progressBar)

disp('Sorting channel quantas..');

progressBar = waitbar(0, 'Step 3/3 Sorting confidence ranges');

t = cputime;

for i = 1:esLen

for j = 1:channelNo

for k = 1:channelNo

ordered(:, i, j, k) = sort(confRanges(:,i, j, k));

end

end

ETR = (cputime – t)/(esLen -1)*(esLen -i);

waitbar(i/(esLen ),progressBar, ['Step 3/3 Sorting confidence ranges ' ' ETR: ' int2str(ETR) ' seconds']);

end

close(progressBar)

upperThresh = squeeze(gcThreshArrays(ceil(0.975 * trials),:,:,:));

% set the confidence levels for all graphs

for i=0:channelNo^2 – 1

toCh = mod(i,channelNo)+1;

stCh = fix(i/channelNo)+1;

conf_high(:, toCh, stCh) = ordered(ceil(0.95 * trials), :, toCh, stCh);

conf_low(:, toCh, stCh) = ordered(ceil(0.05 * trials), :, toCh, stCh);

end

end

end

[influence_start, influence_end] = disp_results(gc, sampling_rate, conf_low, upperThresh, plotStartValue);

influence_lines = zeros(esLen, channelNo, channelNo);

for j = 1:channelNo

for k = 1:channelNo

counter = 1;

valid = 0;

for i = 1:esLen

if (influence_start(j, k, counter)) <= i && (influence_end(j, k, counter) >= i)

influence_lines(i, j, k) = 1;

valid = 1;

else if valid==1

counter = counter + 1;

valid = 0;

end

end

end

end

end

xaxis = esLength:esLength:duration_in_sec;

xaxis = xaxis.';

figure('name','Thresholds and GCI');

% plot the gcThreshold

for i = 0:channelNo^2-1

toCh = mod(i, channelNo) + 1;

stCh = fix(i/channelNo) + 1;

subplot(channelNo, channelNo,i+1),

plot(xaxis(plotStartValue:end),gc(plotStartValue:end, toCh, stCh),'-b');

axis([1 duration_in_sec -1 1.1]);

hold on

plot(xaxis(plotStartValue:end),upperThresh(plotStartValue:end, toCh, stCh),'-r');

plot(xaxis(plotStartValue:end),conf_low(plotStartValue:end, toCh, stCh), '-k', 'LineWidth', 1);

% plot(conf_high(plotStartValue:end, toCh, stCh), '-g', 'LineWidth', 1);

plot(xaxis(plotStartValue:end),influence_lines(plotStartValue:end, toCh, stCh), '–k', 'LineWidth', 2);

title(['ch', int2str(ordered_ch(stCh)) ,'->ch', int2str(ordered_ch(toCh))]);

if stCh == channelNo && toCh == 1

xlabel('Timp (sec)');

ylabel('GCI coeff');

end

end

legend('Granger Causality Index', 'Prag de validare', 'Prag de confidenta', 'Validare conectivitate');

Script signalCorrelationMsCohere.m

function ordered_ch = signalCorrelationMsCohere(channels, sampling_rate)

channelNo = size(channels,2);

trials = size(channels,1);

sample_no = size(channels,3);

section_size = sampling_rate*2;

signal_no_of_sections = sample_no / section_size;

signal_no_of_sections = floor(signal_no_of_sections);

sample_no = signal_no_of_sections*section_size;

signalRanges=[4; 8; 12; 30; 45];

meanCoh = zeros(channelNo, channelNo, signal_no_of_sections, size(signalRanges,1));

msCoh = zeros(channelNo, channelNo, signal_no_of_sections, sampling_rate);

progressBar = waitbar(0, 'Step 0/3 Sorting channels');

tic;

no_of_cycles = trials * sum(1:(channelNo-1));

aux_counter = 1;

for trialNo = 1:trials

startToCh = 2;

for stCh = 1:channelNo

for toCh = startToCh:channelNo

for samples = 0:section_size:sample_no-1

lastCoh = squeeze(msCoh(stCh, toCh, samples/section_size + 1, :));

lastMeanCoh = squeeze(meanCoh(stCh, toCh, samples/section_size + 1, :));

ch1 = squeeze(channels(trialNo,stCh,samples+1:samples+section_size));

ch2 = squeeze(channels(trialNo,toCh,samples+1:samples+section_size));

msc = mscohere(ch1,ch2,sampling_rate,[],sampling_rate*2-1,sampling_rate);

diviz = sampling_rate/(2*size(msc,1));

H=sigwin.hamming(section_size);

win=generate(H);

ch1 = ch1.*win;

fftCh1 = abs(fft(ch1,section_size));

ch2 = ch2.*win;

fftCh2 = abs(fft(ch2,section_size));

meanFFT = mean(fftCh1)*mean(fftCh2);

for i = 1:floor(45/diviz)

numitor = fftCh1(i)*fftCh2(i);

if numitor < (meanFFT * 0.1)

msc(i) = 0;

end

end

msc = msc./trials;

% Split the mscohere vector in the five evaluated EEG

% rhytms

aux(1) = mean(msc(1:floor(4/diviz)));

aux(2) = mean(msc(ceil(4/diviz):floor(8/diviz)));

aux(3) = mean(msc(ceil(8/diviz):floor(12/diviz)));

aux(4) = mean(msc(ceil(12/diviz):floor(30/diviz)));

aux(5) = mean(msc(ceil(30/diviz):floor(45/diviz)));

meanCoh(stCh, toCh, samples/section_size + 1, 🙂 = lastMeanCoh + aux.';

msCoh(stCh, toCh, samples/section_size + 1, 🙂 = lastCoh + msc;

end

ETR = toc / aux_counter * (no_of_cycles – aux_counter);

waitbar(aux_counter / no_of_cycles ,progressBar, ['Step 0/3 Sorting channels ' ' ETR: ' int2str(ETR) ' seconds'])

aux_counter = aux_counter + 1;

end

startToCh = startToCh + 1;

end

end

close(progressBar);

counter = 1;

startToCh = 2;

coherences = zeros(channelNo, channelNo);

for stCh = 1:channelNo

for toCh = startToCh:channelNo

coherences(stCh, toCh) = max(max(meanCoh(stCh, toCh, :, :)));

end

startToCh = startToCh + 1;

end

counter = 1;

startToCh = 2;

aux = coherences;

for stCh = 1:channelNo

for toCh = startToCh:channelNo

[Y, I] = max(aux);

[Y2, I2] = max(Y);

aux(I(I2), I2) = -1;

sortedCoherences(counter,:) = [Y2 I(I2) I2];

counter = counter+1;

end

startToCh = startToCh + 1;

end

disp('Signal coherences:\n');

coherences

disp('Sorted signal coherences:\n');

sortedCoherences

counter = 1;

max_ch = 6;

i = 1;

ordered_ch = 0;

while counter<max_ch

if(~ismember(sortedCoherences(i,2),ordered_ch))

ordered_ch(counter) = sortedCoherences(i,2);

counter = counter + 1;

end

if(~ismember(sortedCoherences(i,3),ordered_ch))

ordered_ch(counter) = sortedCoherences(i,3);

counter = counter + 1;

end

i = i + 1;

end

plotMsCohereForHybrid;

Script plotMsCohereForHybrid.m

figure

channelNo = size(ordered_ch,2);

ordered_ch = sort(ordered_ch);

counter = 1;

sampleLen = 1/sampling_rate;

x = 0.5:0.5:sampling_rate/2;

signal_duration_in_sec = sample_no/sampling_rate;

y = 0:(section_size/sampling_rate):(sample_no/sampling_rate);

for stCh = 1:channelNo

for toCh = stCh+1:channelNo

subplot(5,floor(sum(1:channelNo-1)/5), counter),

imagesc(x, y, squeeze(msCoh(ordered_ch(stCh),ordered_ch(toCh),:,:)));

axis([0 45 0 sample_no/sampling_rate]);

colormap jet;

title(['ch', int2str(ordered_ch(stCh)) ,'<->ch', int2str(ordered_ch(toCh))]);

hold on

if sum(1:channelNo-1)-1 == counter

ylabel('Timp (sec)');

xlabel('Frecventa (Hz)');

end

counter = counter + 1;

end

end

Script disp_results.m

function [influence_start, influence_end] = disp_results(gc, sampling_rate, conf_low, threshold, start_evaluation_value);

% Number of samples

esLen = size(gc,1);

channel_number = size(gc,2);

% Time for signal to stay above threshold in order to be validated as true

% coherence

validation_period = sampling_rate;

% The influence decisions will be made based solely on the comparison with

% the computed threshold.

influence_start = zeros(channel_number,channel_number,ceil(esLen/validation_period));

influence_end = influence_start;

for i = 1:channel_number

for j = 1:channel_number

counter = 1;

for k = 1:esLen

if influence_start(i,j,counter) == 0

if gc(k, i, j) > threshold(k, i, j) && gc(k, i, j) > conf_low(k, i, j) && gc(k, i, j) > 0

influence_start(i, j, counter) = k;

end

else if (influence_end(i, j, counter) == 0) && (gc(k, i, j) < threshold(k, i, j) || gc(k, i, j) < conf_low(k, i, j) || gc(k, i, j) < 0)

if(k – influence_start(i, j, counter)) < validation_period

%RESET start value if coherence is too short

influence_start(i, j, counter) = 0;

else

influence_end(i, j, counter) = k;

counter = counter+1;

continue;

end

end

end

end

if (influence_start(i, j, counter) > 0 && influence_end(i, j, counter) == 0)

influence_end(i, j, counter) = esLen;

end

disp(['Influence DETECTED -> from ch' int2str(j) ' to ch' int2str(i) ' between ' int2str(influence_start(i, j) + 1) ' and ' int2str(influence_end(i, j) + 1)]);

end

end

end

Script Main_DEAP.m

close all

% clear all

clc

signal_structure = load('C:\Users\Tudor\Desktop\Licenta\DEAP database\s01.mat');

% The DEAP database

% sizes(40 x 40 x 8064) –> video/trial x channel x data

signals = signal_structure.data;

% Reduce the size first

signals = signals(11, 1:32, :);

process_data(signals, 0.01, 128);

close all

Script Main_Physiobank.m

close all

% clear all

clc

sampling_rate = 160;

[s1, x] = plotATM('C:\Users\Tudor\Desktop\Licenta\Surse DB\S001R03_edfm');

[s2, x] = plotATM('C:\Users\Tudor\Desktop\Licenta\Surse DB\S001R07_edfm');

[s3, x] = plotATM('C:\Users\Tudor\Desktop\Licenta\Surse DB\S001R11_edfm');

signals = zeros(3,64,size(s1,2));

signals(1,:,:) = s1(1:64,:);

signals(2,:,:) = s2(1:64,:);

signals(3,:,:) = s3(1:64,:);

doCpsd(signals);

writeToFile(signals);

process_data(signals, 0.01, sampling_rate);

Script Main_SimulatedForCuda.m

close all

clc

sampling_rate = 160;

signals = rand(16,50000);

doCpsd(signals);

writeToFile(signals);

Script doCpsd.m

function doCpsd(channels);

batch = size(channels,1);

signalSize = size(channels,2);

evalSize = 200;

fftArray = zeros(batch,signalSize);

for i = 1:batch

for j = 1:signalSize/evalSize

fftArray(i,(j-1)*evalSize + 1:j*evalSize) = fft(squeeze(channels(i,(j-1)*evalSize + 1:j*evalSize)));

end

end

s = zeros(sum(1:(batch-1)),signalSize);

% a = load('C:\Users\Tudor\Desktop\Licenta\Cod sursa\cudaFFT.m');

counter = 1;

tic

for i = 1:size(channels,1)

for k = i+1:size(channels,1)

aux1 = squeeze(fftArray(i,:));

aux2 = squeeze(fftArray(k,:));

for j = 1:signalSize

s(counter,j) = abs(aux1(j)*conj(aux2(j)));

s(counter,j) = 20 * log10(s(counter,j));

end

counter = counter + 1;

end

end

c = toc

%

% figure

% plot(a(189,:));

% title('From CUDA');

%

% figure

% plot(squeeze(s(189,:)));

% title('From MATLAB');

Script CUDA_Matlab_comp.m

close all

app1_cuda = [2.048 1.58 2.06 1.77 2.64 3.49];

app1_matlab = [7.47 14.92 30.12 61.5 153 302];

x1 = [2500 5000 10000 20000 50000 100000];

for i = 1:6

improvement_1(i) = 1/(app1_cuda(i)/app1_matlab(i));

end

app2_cuda = [1.01 1.49 1.84 2.77 4.85 7.76];

app2_matlab = [8.76 36.7 84 153 348 621];

x2 = [16 32 48 64 96 128];

for i = 1:6

improvement_2(i) = 1/(app2_cuda(i)/app2_matlab(i));

end

figure('name','Castig CUDA vs. Matlab');

hold on

plot(1:6,improvement_1, '-ok', 'linewidth',2);

plot(1:6,improvement_2, '-ob', 'linewidth',2);

ylabel('T(cuda) / T(matlab)');

xlabel('Indice simulare');

legend('Aplicatia 1','Aplicatia 2');

Script rulare exemple

clear all

close all

channelNo = 2;

nTrials = 3;

fs = 128;

maxLen = 1124;

discarded = 100;

esLen = maxLen – discarded;

figure

signals = LinearIntStat(maxLen,channelNo,nTrials,discarded);

counter = 1;

for stCh = 1:channelNo

for toCh = stCh+1:channelNo

mcxy = mean(mscohere(squeeze(signals(1,stCh,:)),squeeze(signals(1,toCh,:)),[],[],fs,fs));

subplot(2,3, counter),

mscohere(squeeze(signals(1,stCh,:)),squeeze(signals(1,toCh,:)),[],[],fs,fs),

axis([0 fs/2 0 1]),

title(['ch', int2str(stCh) ,'<->ch', int2str(toCh)]);

hold on

plot(ones(fs/2,1)*mcxy,'-r', 'LineWidth', 2);

counter = counter + 1;

end

end

gc = GCI_calc(4,signals,0.01);

figure('name','GCI plots – EX1');

minimum = min(min(min(gc)));

for i=0:channelNo^2-1

stCh = mod(i,channelNo)+1;

toCh = fix(i/channelNo)+1;

subplot(channelNo,channelNo,i+1),

plot(gc(:,toCh, stCh),'-b', 'LineWidth', 2),

axis([0 esLen minimum 1]),

title(['ch', int2str(stCh) ,'->ch', int2str(toCh)]);

end

signals = squeeze(signals(1,:,:));

pdc_analysis_TVA

pdc = c.pdc;

figure('name','PDC plots – EX1');

for i=0:channelNo^2-1

stCh = mod(i,channelNo)+1;

toCh = fix(i/channelNo)+1;

subplot(channelNo,channelNo,i+1),

plot(squeeze(pdc(toCh, stCh,:)),'-b', 'LineWidth', 2),

axis([0 128 0 1]),

title(['ch', int2str(stCh) ,'->ch', int2str(toCh)]);

end

%EX2

channelNo = 3;

signals = RandomLin(maxLen,channelNo,nTrials,discarded);

figure('name','GCI plots – EX2');

minimum = min(min(min(gc)));

for i=0:channelNo^2-1

stCh = mod(i,channelNo)+1;

toCh = fix(i/channelNo)+1;

subplot(channelNo,channelNo,i+1),

plot(gc(:,toCh, stCh),'-b', 'LineWidth', 2),

axis([0 esLen minimum 1]),

title(['ch', int2str(stCh) ,'->ch', int2str(toCh)]);

end

signals = squeeze(signals(1,:,:));

pdc_analysis_TVA

pdc = c.pdc;

figure('name','PDC plots – EX2');

for i=0:channelNo^2-1

stCh = mod(i,channelNo)+1;

toCh = fix(i/channelNo)+1;

subplot(channelNo,channelNo,i+1),

plot(squeeze(pdc(toCh, stCh,:)),'-b', 'LineWidth', 2),

axis([0 128 0 1]),

title(['ch', int2str(stCh) ,'->ch', int2str(toCh)]);

end

channelNo = 4;

signals = generateFromMatrix(maxLen,channelNo,nTrials,discarded);

gc = GCI_calc(4,signals,0.01);

figure('name','GCI plots – EX3');

minimum = min(min(min(gc)));

for i=0:channelNo^2-1

stCh = mod(i,channelNo)+1;

toCh = fix(i/channelNo)+1;

subplot(channelNo,channelNo,i+1),

plot(gc(:,toCh, stCh),'-b', 'LineWidth', 2),

axis([0 esLen minimum 1]),

% title(['ch', int2str(stCh) ,'->ch', int2str(toCh)]);

end

signals = squeeze(signals(1,:,:));

pdc_analysis_TVA

pdc = c.pdc;

figure('name','PDC plots – EX3');

for i=0:channelNo^2-1

stCh = mod(i,channelNo)+1;

toCh = fix(i/channelNo)+1;

subplot(channelNo,channelNo,i+1),

plot(squeeze(pdc(toCh, stCh,:)),'-b', 'LineWidth', 2),

axis([0 128 0 1]),

title(['ch', int2str(stCh) ,'->ch', int2str(toCh)]);

end

Scripturi generare exemple

function out = RandomLin(channelSize,nChannels,nTrials,discarded);

out = zeros(nTrials, nChannels, channelSize – discarded);

for k=1:nTrials

vector=zeros(nChannels,channelSize);

for i = 1:channelSize

for ch = 1:nChannels

vector(ch,i) = vector(ch,i) + normrnd(0,1);

end

end

%let's discard the first points of this IIR models–> we don't need to see

%the random effect of the series generation (discard first 100 points)

out(k,:,:) = vector(:,discarded+1:end);

end

function out = LinearIntStat(channelSize,nChannels,nTrials,discarded);

out = zeros(nTrials, nChannels, channelSize – discarded);

weighted = zeros(4,4,4);

% Ch1 influences channels [2]

weighted(1,2,2) = 0.65;

% Ch2 influences channels [3 4]

weighted(2,3,4) = 0.9;

weighted(2,4,1) = 0.8;

regrOrd = size(weighted,3);

for k=1:nTrials

vector=zeros(nChannels,channelSize);

for i = 5:channelSize

for ch = 1:nChannels

%chInf – influence channel

for infCh = 1:nChannels

for order = 1: regrOrd

% Add an influence from channel [infCh] from [order]

% samples ago

vector(ch,i) = vector(ch,i) + vector(infCh, i – order) * weighted(infCh, ch, order);

end

end

vector(ch,i) = vector(ch,i) + normrnd(0,1);

end

end

%let's discard the first points of this IIR models–> we don't need to see

%the random effect of the series generation (discard first 100 points)

out(k,:,:) = vector(:,discarded+1:end);

end

function out = generateFromMatrix(channelSize,nChannels,nTrials,discarded);

out = zeros(nTrials, nChannels, channelSize – discarded);

weighted = zeros(4,4,5);

% Ch1 influences channels [1 3]

weighted(1,1,1) = 0.3;

weighted(1,3,1) = 0.6;

% Ch2 influences channels [1 2]

weighted(2,1,4) = 0.65;

% Ch3 influences channels [3]

weighted(3,3,3) = 0.5;

% Ch4 influences channels [2 4 4]

weighted(4,2,5) = 0.6;

regrOrd = size(weighted,3);

for k=1:nTrials

vector=zeros(nChannels,channelSize);

recursiveWeghts = weighted;

for i = 6:channelSize

for ch = 1:nChannels

%chInf – influence channel

for infCh = 1:nChannels

for order = 1: regrOrd

% Add an influence from channel [infCh] from [order]

% samples ago

vector(ch,i) = vector(ch,i) + vector(infCh, i – order) * recursiveWeghts(infCh, ch, order);

end

end

vector(ch,i) = vector(ch,i) + normrnd(0,1);

end

%Simulate some kind of triggers at the middle of the signal

if i==round(channelSize/2)

recursiveWeghts(4,1,2) = 0.6;

recursiveWeghts(4,2,5) = 0;

recursiveWeghts(2,2,1) = 0;

recursiveWeghts(3,2,3) = 0.8;

weighted(3,3,3) = 0;

end

end

%let's discard the first points of this IIR models–> we don't need to see

%the random effect of the series generation (discard first 100 points)

out(k,:,:) = vector(:,discarded+1:end);

end

Script pdc_analysis_TVA.m

Fisier sursa CUDA mass_sort.cu

// includes, system

#include <stdlib.h>

#include <stdio.h>

#include <string.h>

#include <math.h>

// includes, project

#include <cuda_runtime.h>

#include <cufft.h>

#include <helper_functions.h>

#include <helper_cuda.h>

#include <time.h>

#define BATCH 16

#define NO_OF_SIGNALS_IN_FILE 16

#define SIGNAL_SIZE 50000

#define EVAL_SIZE 200

cufftHandle plan;

// Complex data type

typedef float2 Complex;

int lines = 0, i = 0, j = 0;

float *c;

__global__ void cpsd(float *c, const Complex *a, const Complex *b, int maximumIdx)

{

int id = blockIdx.x*blockDim.x+threadIdx.x;

float real = (a[id].x * b[id].x + a[id].y * b[id].y);

float imag = (a[id].y * b[id].x – a[id].x * b[id].y);

float aMod = sqrtf(a[id].x * a[id].x + a[id].y * a[id].y);

float bMod = sqrtf(b[id].x * b[id].x + b[id].y * b[id].y);

c[id] = real * real + imag * imag;

c[id] = sqrtf(c[id]);

c[id] = fabsf(c[id]);

c[id] = 20.0f * log10(c[id]);

}

cudaError_t mscohere(float *c, const Complex *a, const Complex *b, size_t size)

{

Complex *dev_a = 0;

Complex *dev_b = 0;

float *dev_c = 0;

cudaError_t cudaStatus;

int blockSize, gridSize;

// Number of threads in each thread block

blockSize = 1024;

// Number of thread blocks in grid

gridSize = (int)ceil((float)size/blockSize);

// Choose which GPU to run on, change this on a multi-GPU system.

cudaStatus = cudaSetDevice(0);

// Allocate GPU buffers for three vectors (two input, one output) .

cudaStatus = cudaMalloc((void**)&dev_c, size * sizeof(float));

cudaStatus = cudaMalloc((void**)&dev_a, size * sizeof(Complex));

cudaStatus = cudaMalloc((void**)&dev_b, size * sizeof(Complex));

// Copy input vectors from host memory to GPU buffers.

cudaStatus = cudaMemcpy(dev_a, a, size * sizeof(Complex), cudaMemcpyHostToDevice);

cudaStatus = cudaMemcpy(dev_b, b, size * sizeof(Complex), cudaMemcpyHostToDevice);

// Launch a kernel on the GPU with one thread for each element.

cpsd<<<gridSize, blockSize>>>(dev_c, dev_a, dev_b, SIGNAL_SIZE);

// cudaThreadSynchronize waits for the kernel to finish, and returns

// any errors encountered during the launch.

cudaStatus = cudaThreadSynchronize();

// Copy output vector from GPU buffer to host memory.

cudaStatus = cudaMemcpy(c, dev_c, size * sizeof(float), cudaMemcpyDeviceToHost);

Error:

cudaFree(dev_c);

cudaFree(dev_a);

cudaFree(dev_b);

return cudaStatus;

}

void readSignalsFromFile(cufftComplex (*signal)[SIGNAL_SIZE]){

int ch = 0;

FILE *f = fopen("C:\\Users\\Tudor\\Desktop\\Licenta\\Cod sursa\\signalOutput50000_16.txt","r");

while(!feof(f))

{

ch = fgetc(f);

if(ch == '\n')

lines++;

}

fclose(f);

f = fopen("C:\\Users\\Tudor\\Desktop\\Licenta\\Cod sursa\\signalOutput50000_16.txt","r");

for(i = 0; i < lines; i++)

{

for (j = 0 ; j < NO_OF_SIGNALS_IN_FILE; j++)

{

fscanf(f,"%g", &signal[j][i].x);

signal[j][i].y = 0;

}

printf("\n");

}

fclose(f);

}

void write_FFT_to_file(float signal[BATCH][BATCH][SIGNAL_SIZE]) {

FILE *cudaFFTFile = fopen("C:\\Users\\Tudor\\Desktop\\Licenta\\Cod sursa\\cudaFFT.m","w");

for(int i = 0; i < BATCH; i++) {

for(int k = i + 1; k < BATCH; k++) {

// fprintf(cudaFFTFile, "Signal %i –> %i \n", i, k);

for(int j = 0; j < SIGNAL_SIZE; j++) {

fprintf(cudaFFTFile," %6.2f",signal[i][k][j]);

}

fprintf(cudaFFTFile, "\n");

}

}

fclose(cudaFFTFile);

}

int main()

{

auto c = new float[BATCH][BATCH][SIGNAL_SIZE];

auto signal = new cufftComplex[NO_OF_SIGNALS_IN_FILE][SIGNAL_SIZE];

auto real_fft = new float[BATCH][SIGNAL_SIZE];

readSignalsFromFile(signal);

clock_t start = clock();

Complex* d_signal;

cudaMalloc((void**)&d_signal, sizeof(cufftComplex)*SIGNAL_SIZE*BATCH);

cudaMemcpy(d_signal, signal, sizeof(cufftComplex)*SIGNAL_SIZE*BATCH, cudaMemcpyHostToDevice);

/*Create a 1D FFT plan.*/

checkCudaErrors(cufftPlan1d(&plan, EVAL_SIZE, CUFFT_C2C, BATCH*SIGNAL_SIZE/EVAL_SIZE));

/*Use the CUFFT plan to transform the signal in place.*/

checkCudaErrors(cufftExecC2C(plan, (cufftComplex *)d_signal, (cufftComplex *)d_signal, CUFFT_FORWARD));

checkCudaErrors(cudaMemcpy(signal, d_signal, sizeof(cufftComplex)*SIGNAL_SIZE*BATCH, cudaMemcpyDeviceToHost));

/*Destroy the CUFFT plan.*/

cufftDestroy(plan);

cudaFree(d_signal);

for(int i = 0; i < BATCH; i++) {

for(int k = i + 1; k < BATCH; k++) {

cudaError_t cudaStatus = mscohere(c[i][k], signal[i], signal[k], SIGNAL_SIZE);

if (cudaStatus != cudaSuccess) {

fprintf(stderr, "cudaMemcpy failed!");

}

}

}

clock_t end = clock();

float seconds = (float)(end – start) / CLOCKS_PER_SEC;

printf("%f", seconds);

seconds = (float)(end – start) / CLOCKS_PER_SEC;

write_FFT_to_file(c);

}