Approccio basato su Compressive Sensing e modellistica...

DIPARTIMENTO DI INGEGNERIA DELL’INFORMAZIONE

Corso di Laurea Magistrale in Bioingegneria

Approccio basato su Compressive

Sensing e modellistica Taylor-Fourier

per stimare la risposta emodinamica

da segnali fNIRS

Laureando

Marco Cescon

Relatore

Prof. Giovanni Sparacino Correlatori

Ing. Sabrina Brigadoi

Ing. Guglielmo Frigo

– ANNO ACCADEMICO 2015/2016 –

- 19 Aprile 2016 -

INDICE

INDICE ................................................................................ I

ABSTRACT ............................................................................ 1

Capitolo 1

LA SPETTROSCOPIA NEL VICINO INFRAROSSO (NIRS) ..... 3

1.1 ASPETTI GENERALI E STORICI ................................................ 3

1.2 I SETTORI APPLICATIVI DELLA NIRS ...................................... 5

1.3 PUNTI DI FORZA E CRITICITA’ DELLA NIRS ...................... 6

Bassa invasività ................................................................................................................ 6

Complementarità/Intercambiabilità con altre tecniche ............................................ 7

Trasportabilità ................................................................................................................. 8

Praticità e facilità di acquisizione .................................................................................. 8

Basso costo di acquisto e manutenzione ..................................................................... 9

Variabilità dell’SNR ........................................................................................................ 9

Bassa risoluzione spaziale .............................................................................................. 9

Elevato rumore di fondo ............................................................................................. 10

1.4 LA RISPOSTA EMODINAMICA (HRF) .................................... 10

1.5 IL RUMORE ...................................................................................... 13

Il rumore strumentale................................................................................................... 13

L’errore sperimentale ................................................................................................... 14

i

INDICE

2

Le componenti fisiologiche ........................................................................................ 15

1. Componente cardiaca ..................................................................................... 15

2. Componente respiratoria ................................................................................ 15

3. Onda di Mayer ................................................................................................. 15

4. Oscillazioni a basse frequenze ....................................................................... 16

1.6 SCOPO DELLA TESI E OUTLINE ............................................ 16

Capitolo 2

ALGORITMI PER LA RIMOZIONE DEL RUMORE NEL

SEGNALE NIRS ................................................................. 17

2.1 PANORAMICA GENERALE ....................................................... 17

2.2 APPROCCIO DI FILTRAGGIO PASSABANDA .................... 18

2.3 MEDIA CONVENZIONALE ....................................................... 18

2.4 METODI CHE SFRUTTANO IL CANALE DI

RIFERIMENTO ........................................................................................ 20

2.4.1 METODO DI SOTTRAZIONE ............................................................................. 20

2.4.2 METODO CHE INCLUDE IL FILTRAGGIO CON METODO BAYESIANO

(REMCOBA) .................................................................................................................. 21

Capitolo 3

STIMA DELLA RISPOSTA EMODINAMICA MEDIANTE

CS-TFT ............................................................................ 25

3.1 CREAZIONE DI UN DATASET SIMULATO ......................... 26

3.1.1 DATASET IDEALE .............................................................................................. 26

3.1.2 DATASET REALISTICO ....................................................................................... 33

Capitolo 4

ALGORITMO CS-TFT ........................................................ 41

ii

4.1 PRE-PROCESSING ......................................................................... 41

4.2 MEDIA DELLE HRF ...................................................................... 53

4.3 ALCUNI DETTAGLI IMPLEMENTATIVI................................ 54

Capitolo 5

RISULTATI ......................................................................... 57

5.1 SCENARIO COMPLETAMENTE SIMULATO ....................... 57

5.1.1 DETTAGLI IMPLEMENTATIVI ........................................................................... 58

5.1.2 STIMA DELLE HRF ............................................................................................. 60

5.1.3 STIMA DEL RUMORE ........................................................................................... 60

5.1.4 CONFRONTO COL FILTRAGGIO PASSABANDA ................................................ 65

5.2 SCENARIO SEMI-SIMULATO .................................................... 71

5.2.1 DETTAGLI IMPLEMENTATIVI ........................................................................... 72

5.2.2 STIME OTTENUTE E CONFRONTO CON FILTRAGGIO BP E METODO

REMCOBA .................................................................................................................... 72

5.3 DISCUSSIONE ................................................................................. 89

5.4 RMSE................................................................................................... 89

Capitolo 6

CONCLUSIONI E SVILUPPI FUTURI ................................... 93

6.1 SVILUPPI FUTURI .......................................................................... 93

APPENDICE ....................................................................... 97

BIBLIOGRAFIA ................................................................. 109

iii

1

ABSTRACT

La spettroscopia nel vicino infrarosso, definita anche con l’acronimo “NIRS”,

ovvero Near-Infrared Spectroscopy, è una tecnica non invasiva per il monitoraggio

dei cambiamenti emodinamici cerebrali. Si basa principalmente sull’utilizzo

della luce nel rosso e vicino infrarosso e misura le variazioni nell’assorbimento

ottico della luce nel tempo. Dalla differenza tra l’intensità di luce immessa dalla

sorgente e quella misurata dal detettore, è possibile ricavare le variazioni di

concentrazione sia dell’emoglobina ossigenata che deossigenata che

avvengono nella zona indagata, ed esaminare gli effetti metabolici dell’ossigeno

durante l’attività cerebrale.

Il segnale acquisito, tuttavia, risulta composto oltre che dal segnale utile,

ovvero la risposta emodinamica cerebrale, anche da più componenti rumorose

legate a variazioni fisiologiche, quali in particolare le oscillazioni

emodinamiche dovute al battito cardiaco, alla respirazione, e alla pressione

sanguigna, queste ultime a bassa e bassissima frequenza. In letteratura vari

algoritmi sono stati proposti per ridurre questo rumore fisiologico nel segnale

acquisito e migliorare di conseguenza la stima della risposta emodinamica

cerebrale. Tuttavia, non sempre questi metodi sono risultati efficaci (un

semplice filtraggio passa banda, ad esempio, non rimuove completamente il

rumore fisiologico in quanto una parte di esso ha le stesse frequenze della

risposta emodinamica cerebrale).

In questo lavoro è stato sviluppato un nuovo algoritmo per la riduzione del

rumore fisiologico e la stima della risposta emodinamica cerebrale. Tale

algoritmo, definito con l’acronimo CS-TFT, fa uso del campionamento

compresso e implementa la trasformata di Taylor Fourier per ottimizzare la

ABSTRACT

2

stima del rumore fisiologico, che può essere poi sottratta al segnale acquisito.

Tale stima viene effettuata su un segnale acquisito da un canale di riferimento

che contiene le soli componenti rumorose, presenti anche nel canale standard,

ma che non contiene alcuna risposta cerebrale. Dopo la riduzione del rumore

fisiologico, dal canale standard si può stimare la risposta emodinamica

cerebrale. L’algoritmo è stato testato sia su un dataset completamente simulato

che su uno semi-simulato, al fine di testarne il funzionamento in casistiche

differenti. I risultati sono stati poi confrontati con i risultati ottenuti con due

metodi pubblicati in letteratura (il filtraggio passa-banda, e un approccio

basato sul filtro bayesiano).

3

CAPITOLO 1

LA SPETTROSCOPIA NEL VICINO

INFRAROSSO (NIRS)

1.1 ASPETTI GENERALI E STORICI

La spettroscopia nel vicino infrarosso, abbreviata in inglese con l’acronimo

“NIRS”, ovvero Near-Infrared Spectroscopy, è una tecnica non invasiva per il

monitoraggio dei cambiamenti emodinamici cerebrali. Utilizza la luce nel rosso

e vicino infrarosso e misura le variazioni nell’assorbimento ottico della luce nel

tempo. Dalla differenza tra l’intensità di luce immessa dalla sorgente e quella

misurata dal detettore, è possibile ricavare le variazioni di concentrazione sia

dell’emoglobina ossigenata che deossigenata che avvengono nella zona

indagata, ed esaminare gli effetti metabolici dell’ossigeno durante l’attività

cerebrale. La NIRS utilizza lunghezze d’onda comprese nella banda spettrale

che costituisce il rosso e il vicino infrarosso (650 – 950 nm), che risultano

completamente innocue per l’organismo [1][2].

I fotoni emessi dalla sorgente non seguono un percorso rettilineo

corrispondente ad una linea retta tra la sorgente ed il detettore, bensì

percorrono un tragitto più lungo. Essi infatti vengono in parte assorbiti, ma in

maniera preponderante subiscono dispersione (scattering), che fa si che il loro

andamento somigli ad un random walk. L’assorbimento dei fotoni da parte dei

LA SPETTROSCOPIA NEL VICINO INFRAROSSO (NIRS)

4

tessuti come pelle ed ossa è particolarmente basso a queste lunghezze d’onda,

mentre diventa particolarmente rilevante invece grazie all’emoglobina, che

presenta uno spettro di assorbimento differente a seconda che essa sia

ossigenata (HbO2) o meno (HbR).

L’applicazione della NIRS al cervello è attribuita in particolare a Jöbsis [3], che

verso la fine degli anni ‘70 compì importanti studi negli Stati Uniti su questo

nuovo tipo di tecnica di imaging. Jöbsis infatti, è considerato il padre della

NIRS, in quanto dimostrò per primo la possibilità di monitorare non

invasivamente le variazioni di concentrazioni di emoglobina (ossigenata e

deossigenata) nel cervello con l’utilizzo della luce a particolari lunghezze

d’onda [4].

Andando a ritroso, già nei primi anni del XIX secolo fu utilizzata la tecnologia

dell’onda continua per eseguire le prime diagnosi non invasive sul tessuto

umano, tra cui busto e testa, tramite la trasmissione della luce attraverso il

corpo [4]. Con il termine “onda continua” s’intende che la strumentazione si

basa unicamente sulla misurazione dell’intensità della luce, dove viene misurata

l’intensità della luce risultante dopo che questa ha attraversato i tessuti. Tale

tecnica si differenzia in particolare rispetto a quelle basate sul dominio in

frequenza, che misurano il tempo di volo dei singoli fotoni.

Hoppe-Seyler descrive nel 1862 lo spettro dell’ossiemoglobina mentre l’inglese

Stoke, due anni più tardi, aggiunse lo spettro della deossiemoglobina

scoprendo di conseguenza l’importanza dell’emoglobina stessa per il trasporto

dell’ossigeno. Nel 1876 il tedesco Von Vierordt analizzò i tessuti misurando i

cambiamenti spettrali della luce penetrando il tessuto quando la circolazione

sanguigna era occlusa, mentre nel 1894 Hüfner, determinò

spettroscopicamente le quantità assolute e relative di entrambe le tipologie di

emoglobina in vitro. Da quel momento e per alcuni decenni non venne

condotta alcuna rilevante ricerca sul campo fino agli anni ‘30, quando lo studio

sulla determinazione spettroscopica dell’ossigenazione del tessuto venne

ripreso e portato avanti da numerosi ricercatori, fra i quali il tedesco Nicolai,

che continuò gli studi di Von Vierordt, mentre Matthes e Gross

determinarono per la prima volta lo spettro dell’ossiemoglobina e

deossiemoglobina nel tessuto umano usando due lunghezze d’onda, la prima

nella zona rossa mentre la seconda nella banda infrarossa [4].


5

Un importante passo in avanti fu soprattutto la scoperta della legge Beer-

Lambert [5] del matematico francese Bouguer, nel 1729: essa fu estesa poi nel

1852 dal tedesco Beer per quantificare le concentrazioni, ma non era

applicabile ai tessuti biologici in quanto non si basava sui fenomeni privi di

scattering: risulta pertanto relativamente recente (1988) la legge Beer-Lambert

modificata e sviluppata dall’inglese Delpy, che nel suo studio tenne conto

anche dell’effetto diffusione. Le prime ricerche durante gli anni ‘80 e ‘90 si

concentrarono sull’uso della luce nel vicino infrarosso sia in neonati che adulti,

principalmente per monitorare i cambiamenti emodinamici e di flusso in

seguito ad eventi clinici di interesse. Ad inizi anni ‘90 la risonanza magnetica

funzionale (fMRI) fece la sua comparsa, fornendo intere immagini cerebrali del

segnale dipendente dal livello dell’ossigeno nel sangue.

Tornando al lavoro svolto da Jöbsis, la sua scoperta portò alla successiva

progettazione e costruzione di numerosi strumenti NIRS basati sulla

misurazione dell’intensità della luce. Negli stessi anni in cui la fMRI faceva la

sua comparsa, anche la NIRS veniva applicata per la prima volta in campo

funzionale (fNIRS). L’anno considerato più importante per lo sviluppo e la

nascita di tale nuova tecnica è sicuramente il 1993: quattro gruppi di ricerca

[6][7][8][9], dimostrarono la possibilità di analizzare l’attività emodinamica

cerebrale in seguito a particolari stimoli con l’utilizzo della fNIRS. [4][10].

1.2 I SETTORI APPLICATIVI DELLA NIRS

Negli ultimi tre decenni questa tecnologia ha continuato a crescere: ad oggi

infatti gli studi con la fNIRS sono stati condotti in svariati ambiti, dal campo

psicologico a quello fisiologico e nello studio delle patologie cerebrali.

Importante è considerare che la fNIRS può essere utilizzata per studiare solo

le aree corticali, perché la luce non riesce a raggiungere zone cerebrali più

profonde, più adatte ad essere studiate con la fMRI. Gli ambiti applicativi della

NIRS sono molteplici, e va menzionato il fatto che oltre ad essere utilizzata

per monitorare il cervello, essa è stata impiegata in particolare anche sul

tessuto muscolare, anche se con risultati meno accurati, in quanto il cromoforo


6

associato (citocromo ossidasi) possiede un grado di assorbimento dei raggi

infrarossi di circa un ordine di grandezza inferiore rispetto all’emoglobina,

rendendo quindi più difficile la lettura finale. Ulteriori studi con la NIRS sono

stati inoltre realizzati nell’ambito della chirurgia neurovascolare, in cui si è

dimostrata la capacità di individuare eventi ischemici dal segnale NIRS o nel

monitoraggio di pazienti che presentavano gravi traumi cranici, oltre all’analisi

sul decorso di alcune malattie quali ad esempio l’Alzheimer. Sempre in ambito

chirurgico, la NIRS è stata valutata anche come metodo nella diagnosi di

arteriosclerosi delle arterie carotidee e vertebrali o per la misurazione non

invasiva del pH [11].

In ambito clinico, è utilizzata per il monitoraggio della saturazione cerebrale

dei neonati, soprattutto in terapia intensiva neonatale. Ha fornito inoltre un

importante apporto nello sviluppo della riabilitazione motoria. E’ soprattutto

in ambito cognitivo che la fNIRS ha avuto grande applicazione, in particolare

in popolazioni in età evolutiva e pediatrica, dove altre tecniche di monitoraggio

dello sviluppo cognitivo cerebrale sono di più difficile utilizzo.

1.3 PUNTI DI FORZA E CRITICITA’ DELLA NIRS

L’interesse verso la NIRS è determinato in particolare da una serie di aspetti

che la rendono vantaggiosa rispetto ad altri tecniche di monitoraggio cerebrale.

Vediamone i principali:

Bassa invasività

La tecnica NIRS misura i cambiamenti di assorbimento ottico nel

tempo. Questo viene fatto in modo non invasivo disponendo dei

sensori ottici direttamente sul cuoio capelluto del soggetto. I sensori

inviano e ricevono fasci di luce a lunghezze d’onda nel rosso e vicino

infrarosso, registrando la differenza di intensità tra la luce inviata e

quella misurata dopo aver attraversato il cranio ed aver raggiunto il

cervello. Grazie al basso assorbimento ottico dei tessuti biologici a


7

queste lunghezze d’onda, la luce riesce a penetrare alcuni centimetri,

raggiungendo la corteccia cerebrale. E’ una tecnica completamente

non invasiva in quanto non necessita di particolari somministrazioni

di farmaci o traccianti, come avviene ad esempio nella tomografia ad

emissioni di positroni (PET) con i traccianti radioattivi. [12]

Figura 1.1: Esempio di disposizione dei probe durante un’acquisizione fNIRS effettuata all’università di Zurigo. Le fibre ottiche sono disposte su una cuffia a contatto col cuoio capelluto. Le fibre in grigio sono le sorgenti, mentre quelle nere i detettori [4]

Complementarità/Intercambiabilità con altre tecniche

La tecnica NIRS può essere considerata una valida sostituta dell’fMRI

in soggetti e situazioni dove quest’ultima può difficilmente essere

utilizzata (ad esempio neonati, pazienti con pacemaker, etc.). Il

confronto può essere esteso anche all’elettroencefalografia (EEG), in

quanto le prestazioni offerte da queste tre tipologie di tecniche


8

possono essere paragonabili: la NIRS ad esempio ha una risoluzione

temporale particolarmente buona (l’ordine è di qualche millisecondo)

se confrontata con l’fMRI (pari invece a circa 2 – 3 secondi) ma

inferiore rispetto all’EEG. A livello spaziale invece, la situazione si

inverte: la NIRS permette una risoluzione spaziale comparabile con

quella offerta dall’EEG (entrambe dell’ordine dei centimetri) ma

nettamente inferiore alla fMRI (i cui voxel sono dell’ordine dei

millimetri, solitamente 3x3x3 mm) [10].

Trasportabilità

Generalmente, gli strumenti e le attrezzature necessarie a svolgere

un’acquisizione NIRS sono più facilmente trasportabili in quanto più

piccoli e di dimensioni ridotte rispetto ad una normale risonanza

magnetica che invece non ha alcuna possibilità di movimentazione, e

risultano particolarmente vantaggiosi soprattutto in casistiche

particolari e delicate, quali ad esempio la non deambulazione del

paziente o la sua impossibilità di spostamento (ad esempio se il

paziente è in terapia intensiva). In questi casi, la strumentazione NIRS

può essere facilmente posizionata ai piedi del letto senza particolari

difficoltà, ed impiegata inoltre anche su pazienti portatori di protesi

metalliche, che non possono essere sottoposti ad fMRI.

Praticità e facilità di acquisizione

Le procedure sperimentali NIRS sono condotte in un ambiente più

controllabile a livello sperimentale e quindi più ecologico, rispetto ai

rumorosi scanner chiusi usati in fMRI. La NIRS inoltre può tollerare il

movimento del soggetto in proporzioni maggiori rispetto alla fMRI,

purché la sonda NIRS rimanga fissa e stabile sul cuoio capelluto del

soggetto. Ciò ha permesso alla NIRS di essere impiegata anche sui

soggetti che più difficilmente si adattavano a tale situazione, quali ad

esempio bambini o pazienti anziani affetti da particolari patologie.

Questo ha anche permesso alla NIRS di essere usata per studi

richiedenti soggetti in movimento, per esempio per studiare gli effetti

di un esercizio, o la posizione su segnali cerebrali.


9

Basso costo di acquisto e manutenzione

La popolarità crescente di questa tecnica è in parte dovuta al basso

costo dell’attrezzatura se comparata alla risonanza magnetica o alla

PET. Gli strumenti ottici infatti, oltre ad essere più piccoli e più

trasportabili sono anche meno costosi da acquistare e richiedono

inoltre minori costi di manutenzione (contrariamente ad esempio ad

un normale scanner MRI).

La NIRS tuttavia presenta anche alcuni aspetti critici o addirittura di

svantaggio nel confronto con altre tecniche di imaging diagnostico. Le

riportiamo di seguito.

Variabilità dell’SNR

Il segnale acquisito, ed in particolare l’SNR (Signal to noise ratio) ha una

variabilità che si differenzia molto da soggetto a soggetto, in quanto

influenzato da diversi fattori, quali ad esempio il colore dei capelli ma

soprattutto lo spessore del cranio, che varia molto sia l’assorbimento

che lo scattering del segnale. La luce emessa infatti, attraversa in

sequenza il cuoio capelluto, il cranio e liquido cerebro-spinale, e il suo

assorbimento è influenzato in particolar modo da questi 3 strati, che

variano in spessore e proprietà ottiche molto in base all’età, al sesso, e

ad altre caratteristiche genetiche. Oltre a ciò, anche il colore dei capelli

ha una particolare influenza sull’assorbimento di luce (le tonalità scure

infatti conducono ad un assorbimento maggiore rispetto a quelle

bionde o più chiare) [12].

Poiché la NIRS è altamente sensibile alle fonti esterne di luce,

l’acquisizione viene solitamente effettuata in una stanza semi buia, in

modo da non influenzare o contaminare con luce esterna il fascio

luminoso misurato dal detettore [10].

Bassa risoluzione spaziale

Come accennato in precedenza, nel confronto con le altre

metodologie (in particolare con la fMRI), la risoluzione spaziale

offerta dalla NIRS è alquanto limitata.


10

Nelle acquisizioni fNIRS, solitamente la sorgente e ildetettoresono

disposti sullo scalpo ad una distanza di circa 2-3 cm l’un dall’altro. La

luce può raggiungere in questo modo una profondità di circa 3 cm

sotto la corteccia cerebrale.

Ciò non consente di “ispezionare” i tessuti cerebrali più profondi, in

quanto l’SNR risultante avrebbe valori troppo bassi. La risoluzione

spaziale può essere tuttavia migliorata attraverso l’uso di array ad alta

densità, già utilizzati ad esempio nella tomografia ottica diffusa

(DOT). [13]

Elevato rumore di fondo

Il segnale acquisito non solo contiene la risposta cerebrale di interesse,

ma anche quello che viene definito rumore fisiologico. Come

vedremo nel paragrafo 1.5, il segnale acquisito infatti, è costituito da

un insieme di componenti, quali principalmente le risposte evocate

cerebrali, a cui si aggiungono una serie di variazioni emodinamiche

dovute al battito cardiaco, alla respirazione, e ad altre componenti

fisiologiche. Di conseguenza, si rende necessaria una procedura di

elaborazione del segnale NIRS stesso, che conduca all’estrazione delle

risposte cerebrali, per evitare di catalogare variazioni emodinamiche

dovute alla fisiologia o al rumore come risposte cerebrali.

1.4 LA RISPOSTA EMODINAMICA (HRF)

La risposta emodinamica (HRF, Hemodynamic Response Function) è la

componente di principale interesse contenuta nel segnale fNIRS, in quanto

rappresenta l’attivazione cerebrale che avviene in seguito alla presentazione al

soggetto di particolari stimoli. Quando al soggetto viene chiesto di eseguire un

particolare compito, i neuroni responsabili iniziano a lavorare, incrementando

il loro bisogno metabolico. Il cervello risponde a questa richiesta inducendo un

aumento del flusso sanguigno cerebrale verso le particolari aree del cervello

interessate. Per un principio di conservazione (il principio di coupling neuro-

vascolare) l’aumento di flusso sanguigno avviene in quantità molto maggiore


11

rispetto al reale consumo di ossigeno da parte dei neuroni. Pertanto, quello che

si misura è un aumento dell’emoglobina ossigenata, e una diminuzione di

quella deossigenata, al contrario di quello che ci si potrebbe aspettare. La

fNIRS riesce a misurare in maniera distinta, a differenza ad esempio dell’fMRI,

sia la variazione di HbO2 che di HbR, che hanno un andamento simile a

quanto raffigurato nella Figura 1.2:

Figura 1.2: Andamento della concentrazione di emoglobina ossigenata (rosso) e deossigenata (verde) in seguito alla presentazione di uno stimolo al tempo 0.

Solitamente, la HRF ha un andamento lento, con una durata di circa 10 – 15

secondi, molto maggiore pertanto rispetto alle risposte elettriche neuronali

misurate con l’EEG (dell’ordine dei millisecondi).

I modelli parametrici impiegati per stimare o simulare la HRF sono vari, ma il

più usato è il modello canonico dato dalla somma di due funzioni gamma-

variate secondo la formula:

( ) ( )

( ) (1.1)


12

dove

( )

(

)

( )

( ) (1.2)

dove regola l’ampiezza, mentre e rispettivamente la larghezza e l’istante

iniziale della risposta. Il parametro infine, regola il rapporto tra la risposta e il

suo undershoot (l’ampiezza della seconda funzione gamma che modella

l’undershoot della risposta). In figura 1.3 è raffigurato l’andamento atteso della

risposta emodinamica, secondo il modello canonico appena menzionato.

Figura 1.3: Risposta emodinamica realizzata mediante modello canonico (1.1) con i seguenti

parametri: , , , , , ,

Da anni ormai i ricercatori studiano ed implementano varie tecniche e

algoritmi con l’obiettivo di rimuovere la componente rumorosa fisiologica,

migliorando pertanto l’SNR.


13

1.5 IL RUMORE

Abbiamo già menzionato in precedenza la difficoltà di estrarre dal segnale

NIRS le sole risposte emodinamiche cerebrali, in quanto il segnale misurato

contiene anche alcune componenti fisiologiche di disturbo che alterano le

HRF stesse, perché possono averne la stessa frequenza e fase. Le misurazioni

NIRS sono caratterizzate principalmente da tre tipologie di rumore: il rumore

strumentale, l’errore sperimentale e le fluttuazioni fisiologiche. I primi due tipi

di rumore possono essere facilmente rimossi, tramite dei filtri e un’accurata

predisposizione del disegno sperimentale.

Vediamo in dettaglio le componenti fondamentali.

Il rumore strumentale

E’ costituito solitamente dal rumore di misura dovuto alle componenti

hardware, nonché alla luce di fondo presente nella stanza, e

solitamente presenta uno spettro di frequenze abbastanza uniforme.

Dal momento che le frequenze di campionamento degli strumenti

NIRS sono solitamente molto più elevate rispetto a quelle dei segnali

emodinamici di interesse, si stima pertanto che le frequenze alte siano

date soprattutto dal rumore strumentale. Molte di queste possono

essere separate dalle componenti fisiologiche mediante un semplice

filtro passabasso, o modificando l’intensità del fascio o il guadagno

dell’amplificazione dell’intensità ricevuta [12].


14

Figura 1.4: Simulazione del rumore di misura (finestra di 1000 s)

L’errore sperimentale

Come descritto in precedenza, la NIRS è poco sensibile ai movimenti

della testa: se l’accoppiamento optodo-testa è buono, il segnale non

presenta particolari artefatti anche in caso di movimenti bruschi del

soggetto. Tuttavia, non sempre risulta possibile garantire che non vi

sia alcun distacco dell’optodo dalla cute, causando visibili variazioni

sul segnale acquisito simili a volte a “salti” o impulsi, che possono

compromettere la baseline del segnale o modificare l’ampiezza della

HRF. La rimozione di simili artefatti risulta a volte difficile. Algoritmi

che analizzano le tempistiche delle fluttuazioni del segnale, andando a

cercare brusche variazioni in ampiezza e/o deviazione standard,

vengono solitamente utilizzati per riconoscere gli artefatti da

movimento, che possono poi essere corretti utilizzando varie tecniche

proposte in letteratura. Molti errori sperimentali sono evitabili grazie

alla competenza dell’operatore nel preparare la strumentazione


15

necessaria e il soggetto al fine da ridurre le fluttuazioni indesiderate

[12].

Le componenti fisiologiche

Una delle più insidiose fonti di rumore è data dai segnali fisiologici

apprezzabili anche a livello superficiale del cuoio capelluto e della

corteccia cerebrale sottostante. Tali componenti fisiologiche sono

principalmente causate dal battito cardiaco, dalla respirazione e dalla

cosiddetta onda di Mayer. Vi sono poi ulteriori oscillazioni definite a

bassa frequenza.

Queste quattro categorie sono illustrate sotto in dettaglio.

1. Componente cardiaca. Il battito cardiaco è assimilabile ad una

sinusoide avente una frequenza compresa tra 60 e 80 pulsazioni al

minuto e un’ampiezza che varia nel range tra 0 e 500 nM. Questa

componente non è cruciale nella stima delle HRF. Questo tipo di

oscillazioni è facilmente eliminabile mediante un filtro passa-

basso, in quanto la sua frequenza (~1 Hz) è ben lontana da quella

di interesse della HRF.

2. Componente respiratoria. Rappresentabile sempre attraverso

una sinusoide, ha anch’essa un’ampiezza variabile tra 0 e 500 nM,

ma frequenza molto più bassa, ovvero intorno ai 0.2 Hz, molto

simile alla frequenza della risposta emodinamica. Il filtraggio

pertanto non risulta in questo caso un metodo efficace per

rimuovere questo tipo di artefatto, in quanto rimuoverebbe anche

la risposta emodinamica, compresa nella stessa banda di

frequenze.

3. Onda di Mayer. Definita anche onda “vasomotore”, è dovuta

probabilmente alle contrazioni dei muscoli presenti sulle pareti dei

vasi sanguigni. Ha una frequenza piuttosto bassa (0.03 – 0.18 Hz)

che si sovrappone a quella della HRF, ma ampiezza maggiore di

quest’ultima, che può variare in base all’età e al sesso del soggetto.


16

Anche in questo caso, risulta impossibile rimuovere tale tipo di

componente mediante un filtraggio, in quanto verrebbe eliminata

anche la risposta emodinamica, dal momento che i due segnali

hanno frequenza simile.

4. Oscillazioni a basse frequenze. Si tratta di sinusoidi molto

lente, aventi una frequenza media di circa 0.004 Hz e ampiezza

variabile che può seguire alcuni trend, che si modificano

tipicamente in base alle variazioni di pressione nei vasi sanguigni.

1.6 SCOPO DELLA TESI E OUTLINE

L’obiettivo che ci prefissiamo, è quello di implementare l’algoritmo CS-TFT

idealizzato da Frigo et al. [13] su un segnale fNIRS, al fine di valutarne

l’efficienza nella ricostruzione del rumore presente in modo da poterlo poi

rimuovere, in un secondo momento, riuscendo di conseguenza ad ottenere

una stima delle HRF estratte. Al fine di testare la vera efficienza dell’algoritmo,

si cercherà di implementarlo su un particolare database apposito, in modo da

valutarne concretamente il risultato confrontandolo poi con lo stesso ottenuto

da altri algoritmi già presenti in letteratura e menzionati in precedenza.

17

CAPITOLO 2

ALGORITMI PER LA RIMOZIONE

DEL RUMORE NEL SEGNALE

NIRS

2.1 PANORAMICA GENERALE

Il numero totale di metodi sino ad oggi sperimentati è elevato e una loro

quantificazione precisa risulta difficile, ma riportiamo qui i principali e più

utilizzati, in particolare quelli univariati, ovvero a variante unica.

I metodi a una variante comprendono alcuni dei primissimi e più semplici

algoritmi, quali in particolare il filtraggio passabanda e la media convenzionale

(definita con l’acronimo CA, ovvero Conventional Averaging). Vi sono poi

ulteriori metodologie implementate più di recente che, sempre considerando

un singolo segnale alla volta, hanno un funzionamento più complesso, e che

sfrutta un particolare tipo di canale, il canale di riferimento, che contiene solo

il rumore fisiologico di fondo e non la risposta cerebrale, e che può quindi

essere utilizzato per una più accurata riduzione del rumore fisiologico nel

canale standard. Tra questi vi sono la regressione ai minimi quadrati [15]

ReMCoBA, che utilizza un metodo Bayesiano [16], e l’approccio parametrico

basato sul filtro di Kalman [10].

ALGORITMI PER LA RIMOZIONE DEL RUMORE NEL SEGNALE NIRS

18

Vediamoli in dettaglio:

2.2 APPROCCIO DI FILTRAGGIO PASSABANDA

Uno dei più semplici metodi utilizzati per filtrare i segnali NIRS è il comune

filtraggio passabanda. Esso permette la rimozione di tutta quella parte di

segnale le cui frequenze non sono di particolare interesse per l’analisi, ovvero

le alte frequenze e le bassissime frequenze. Questo semplice metodo sfrutta

solo i segnali del canale standard permettendo di rimuovere tutte le oscillazioni

che presentano sia frequenze piuttosto elevate, come ad esempio quelle

cardiache, sia quelle molto lente, inferiori allo 0.01 Hz. Questo metodo è stato

ampiamente utilizzato in molti studi, non solo come metodo stesso, ma anche

come passo preliminare di pre-processing prima di applicare al segnale NIRS

algoritmi di filtraggio più avanzati. Tuttavia, il filtraggio non permette di

rimuovere completamente il rumore fisiologico di fondo. Ciò avviene in

particolare alle basse frequenze, in quanto vi è completa sovrapposizione tra il

rumore fisiologico e le HRF, poiché hanno frequenze molto simili tra loro

[REF].

2.3 MEDIA CONVENZIONALE

Il CA è un altro metodo piuttosto semplice, utilizzato per stimare le risposte

emodinamiche medie dal segnale acquisito. Il funzionamento si basa su una

prima suddivisione in trial del segnale, tutti della stessa durata. Una volta

ottenuta la serie di trial, se ne calcola la media aritmetica secondo la formula:

( )

∑ ( ( ) ( )) ( )

∑ ( )

(2.1)


19

dove è il segnale acquisito, è la HRF del trial i-esimo(con i che va da 1 a

N, dove N è il numero totale di trial), e è la componente rumorosa

(contenente sia il rumore fisiologico che di misura).

Dato che il segnale totale acquisito può essere considerato formato dalla

somma della risposta emodinamica più una componente di rumore, si può

ipotizzare che il rumore abbia media nulla, ovvero:

( ) ( )

∑ ( )

( ) (2.2)

Si ipotizza pertanto che il valore atteso del segnale acquisito tenda ad essere

quanto più simile possibile alle vere risposte emodinamiche, in luce della media

nulla del rumore e di un numero di HRF inserite che risulti più elevato

possibile. Tuttavia, uno dei principali aspetti critici della media convenzionale è

che essa richiede un alto numero di trial per generare una stima affidabile. Il

problema, pertanto, è la lunghezza dell’acquisizione NIRS necessaria per

ottenere un numero di trial sufficienti. Un altro aspetto negativo è legato al

fatto che, sebbene il rumore di misura possa essere considerato con valore

atteso tendente a zero, questo non sempre è vero per il rumore fisiologico. In

tal modo, verrebbe meno una delle ipotesi di partenza di questo metodo,

inficiandone la validità. Potremmo dunque dire che la media convenzionale

risulta un metodo piuttosto semplicistico per la stima delle HRF, in quanto

presenta delle limitazioni piuttosto importanti specie qualora il rumore

fisiologico di fondo e le risposte emodinamiche si presentino in fase [16].

Nell’algoritmo qui implementato, la media convenzionale è stata utilizzata

nella fase finale del codice, dopo aver ridotto il rumore fisiologico, come

tecnica di stima della HRF media partendo da risposte emodinamiche “pulite”

in ogni trial.


20

2.4 METODI CHE SFRUTTANO IL CANALE DI RIFERIMENTO

2.4.1 METODO DI SOTTRAZIONE

Per ottimizzare la riduzione del rumore fisiologico dal segnale acquisito in un

canale NIRS standard, Saager e Berger, [15], proposero un’interessante

soluzione. Riducendo la distanza tra sorgente e detettore (canale short-

separation o SS), il fascio luminoso non raggiunge la corteccia e il segnale

NIRS misurato contiene quindi solo la rappresentazione del rumore fisiologico

extra-cerebrale. Questo segnale può essere utilizzato per ridurre il rumore

fisiologico presente nei canali NIRS standard, in quanto viene acquisito

simultaneamente e, non contenendo alcuna informazione sull’attività cerebrale,

non si corre il rischio di eliminare parte del segnale utile, come nei metodi

descritti in precedenza. Il metodo proposto da Saager e Berger (2005), che

prevedeva la sottrazione del segnale SS fittato al segnale del canale standard

con dei minimi quadrati lineari, è stato il precursore di tutti i successivi

algoritmi basati sull’idea di sfruttare il canale SS per ridurre il rumore

fisiologico.

I canali standard (o LS), che contengono sia rumore fisiologico che la risposta

cerebrale, e quelli SS, che contengono solo il rumore fisiologico, possono

infatti essere descritti nel seguente modo:

(2.3)

(2.4)

dove sono le risposte emodinamiche, rappresenta il rumore

fisiologico e il rumore di misura, che sarà diverso nel canale LS e SS,

ma derivante dalla stessa distribuzione. I segnali LS e SS vengono acquisiti in

contemporanea, per cui si può supporre che il rumore fisiologico presente in

entrambi i segnali sia molto simile.


21

La risposta HRF si stima quindi sottraendo il canale SS al canale LS, ovvero:

( ) ( ) ( ) (2.5)

( ) è un generico canale convenzionale LS, mentre ( ) è il canale di

riferimento (SS), moltiplicato per un fattore di scala , che minimizza lo scarto

quadratico medio tra ( ) e ( ).

Figura 2.1: Tipica disposizione di sorgente e detettore. Solo il canale costituito da una distanza sorgente-detettore di circa 3 cm, ovvero LS (detector 2), permette un attraversamento del tessuto cerebrale da parte della luce infrarossa [20]

2.4.2 METODO CHE INCLUDE IL FILTRAGGIO CON METODO

BAYESIANO (REMCOBA)

“ReMCoBA”(abbreviazione di “Reference-channel Modelling Corrected Bayesian

Approach”) è un metodo di riduzione del rumore fisiologico che sfrutta il canale

SS [16]. È considerato un metodo non parametrico nella stima della risposta

emodinamica. Brevemente, l’algoritmo consiste in due passi principali: il primo

passo consiste nella stima, trial per trial, del rumore fisiologico, identificato dal

segnale SS fittando un modello dato da una somma di sinusoidi a diversa


22

frequenza, e nella sua sottrazione dal segnale LS; il secondo passo consistente

nel filtrare con un approccio bayesiano i trial corretti del canale LS, per ridurre

ulteriormente il rumore residuo.

Più dettagliatamente, il segnale LS, ( ), si può considerare una sommatoria

costituita da tre componenti, la HRF, indicata con ( ), il rumore fisiologico,

( ), e il rumore di misura, ( ):

( ) ( ) ( ) ( ) (2.6)

Nel dettaglio, il rumore ( ) può essere modellizzato mediante una

sommatoria di sinusoidi così strutturata:

( ) ∑ ( ) ( ) ( )

(2.7)

dove ( ) è l’errore del modello, e sono le ampiezze delle sinusoidi,

mentre è la frequenza e l’intercetta (offset).

Il canale SS si può considerare composto dallo stesso rumore fisiologico del

canale LS, ( ), scalato però di una costante che è dovuta al differente

percorso effettuato dai fotoni nei due canali. Al rumore fisiologico nel canale

SS si aggiunge anche il rumore di misura casuale, definito :

( ) ( )

( ) (2.8)

Il coefficiente è determinato con un fit ai minimi quadrati tra i segnali LS e

SS, utilizzando i primi 30 secondi di acquisizione, dove al soggetto non è

ancora stato chiesto di eseguire alcun compito. In tal modo si può supporre

che non sia presente risposta emodinamica cerebrale in questo periodo

temporale, per ridurre il rischio che il coefficiente s la fitti e sia quindi

sovrastimato.


23

Il segnale SS viene utilizzato per stimare i parametri del modello del rumore

fisiologico ( ): i valori dei parametri , e vengono calcolati trial per trial

con un fit ai minimi quadrati mentre e sono ricavati dalla densità di

potenza spettrale, anch’essa calcolata in ogni trial. Il valore più comunemente

usato per è 3, a rappresentare il rumore dovuto alla respirazione, all’onda di

Mayer e alle onde con frequenza molto bassa.

Una volta stimato il modello del rumore fisiologico, esso viene sottratto, trial

per trial, al segnale LS:

( ) ( ) ∑[ ( ) ( )]

(2.9)

Il segnale corretto ( ) viene poi filtrato, trial per trial, mediante un

approccio Bayesiano non parametrico. ( ) può essere modellizzato

mediante la somma di 2 segnali: ( ), che è costituito dalle HRF e ( ) che

rappresenta invece il rumore di misura:

( ) ( ) + ( ) (2.10)

La stessa equazione, scritta in forma matriciale, si presenta come:

(2.11)

dove e contengono gli n campioni di un trial.

Un modello stazionario autoregressivo di ordine 4 è stato usato come

descrizione a priori di v, mentre l’informazione a priori su u è stata fornita

modellando la sua “smoothness” come la realizzazione di un processo

stocastico ottenuto dalla cascata di due integratori a rumore bianco gaussiano

( ) e varianza . In queste ipotesi, si può ottenere una stima lineare della

risposta emodinamica, mediante la formula della stima a minima varianza

d’errore, ovvero:


24

( ) (2.12)

con

( ) (2.13)

( ) (2.14)

dove è la matrice di covarianza del rumore di processo, A è una matrice

quadrata di Toeplitz contenente i coefficienti del modello autoregressivo e

è la varianza del rumore di processo stimata dai dati calcolati in un intervallo

che solitamente va da 1.5 secondi prima dello stimolo a 2.5 dopo [16]. è la

matrice di covarianza di u, è la matrice quadrata n-dimensionale

triangolare inferiore di Toeplitz che contiene i coefficienti delle derivate

seconde di , ovvero che ha sulla prima colonna i valori [1, -2, 1, 0, 0, …, 0]T e

con determinato, trial per trial, con il criterio di discrepanza di

Twomey (1965).

Dopo aver filtrato i singoli trial, la risposta emodinamica media viene calcolata

mediando i trial con la media convenzionale. Un filtro di Savitzky-Golay del

terzo ordine e una finestra di 25 campioni viene infine applicato alla HRF

media per ridurre il rumore residuo alle alte frequenze.

Il metodo ReMCoBA risulta particolarmente buono nella stima delle HRF sia

su dati simulati che su dati reali, ed ha il grande vantaggio di non richiedere

modelli a priori sulle risposte emodinamiche stesse, essendo un metodo non

parametrico nella stima della HRF. Tuttavia, l’accuratezza del metodo nella

ricostruzione del rumore fisiologico sul singolo trial non è ottimale, poiché lo

spettro di frequenze è costruito su un numero finito e limitato di campioni.

25

CAPITOLO 3

STIMA DELLA RISPOSTA

EMODINAMICA MEDIANTE

CS-TFT

In questo capitolo viene presentato l’algoritmo di stima della risposta

emodinamica che consta di tre passaggi fondamentali: 1) la definizione delle

componenti spettrali del rumore fisiologico mediante una tecnica di super-

risoluzione basata sul paradigma del campionamento compresso (compressive

sensing, CS); 2) la sottrazione del modello del rumore fisiologico stimato

ottenuto mediante una trasformata di Taylor-Fourier (TFT) dal segnale LS; 3)

la stima della risposta emodinamica, ottenuta mediando i profili delle risposte

riconducibili al medesimo stimolo o compito funzionale presentato al

soggetto. Di seguito verranno anche descritti i dataset utilizzati per testare

l’algoritmo proposto.

STIMA DELLA RISPOSTA EMODINAMICA MEDIANTE CS-TFT

26

3.1 CREAZIONE DI UN DATASET SIMULATO

L’algoritmo che qui proponiamo è stato testato su due dataset, che simulano,

in maniera più o meno veritiera, lo scenario che si può presentare durante

un’acquisizione di segnale fNIRS. I due dataset differiscono per il grado di

complessità e verosimiglianza: il primo, detto anche ideale, è totalmente

simulato e ricavato da un modello matematico dei segnali in esame (la risposta

emodinamica e le diverse componenti rumorose); il secondo, detto realistico,

utilizza del rumore fisiologico acquisito, e quindi reale, e vi aggiunge delle

risposte emodinamiche simulate. Nel primo caso, tutte le grandezze in gioco

sono controllabili ed è possibile verificare se le diverse componenti rumorose

vengono stimate con sufficiente accuratezza. Anche nel secondo caso, la

risposta emodinamica è nota a priori ed è possibile valutare l’accuratezza

dell’algoritmo nella sua stima; tuttavia, il rumore è reale, e quindi lo scenario

molto più simile a quello di una reale acquisizione NIRS.

3.1.1 DATASET IDEALE

Nel cosiddetto dataset ideale, sia la risposta emodinamica sia le componenti

spurie dovute al rumore fisiologico sono generate a partire da un modello

matematico che è stato derivato dalla pratica sperimentale.

In particolare, il modello per la generazione della risposta emodinamica è lo

stesso già menzionato nel capitolo 1, ovvero:

( ) ( )

( ) (3.1)

dove le due funzioni sono funzioni gamma variante:

( )

(

)

( )

( ) (3.2)


27

il cui istante iniziale dalla funzione gradino:

( ) { ( )

(3.3)

In questo modello, introdotto in [1], i parametri , e determinano

l’ampiezza, la durata e il tempo di inizio della risposta, mentre il parametro

definisce il rapporto tra il picco della risposta e l’undershoot seguente. In

tabella x sono riportati i valori medi e relative deviazioni standard delle

distribuzioni gaussiane da cui sono derivati i valori dei parametri usati per

creare il dataset in esame, in modo tale da lasciare una certa variabilità sia tra

trial che tra soggetti. Il coefficiente viene impostato pari a coerentemente

con i valori nominali presentati in [17]. L’intervallo inter-stimolo è stato posto

pari a 12 s.

Nella versione simulata, le HRF sono state inserite ad intervalli di tempo non

regolari, in modo da simulare una reale acquisizione NIRS, dove le risposte

emodinamiche, per esempio dovute ad un compito di finger tapping, non hanno

cadenza costante per evitare eventuali sincronizzazioni di fase tra la risposta

emodinamica e quella fisiologica.

Figura 3.1: Esempio di risposte emodinamiche aggiunte al segnale completamente simulato


28

Il rumore fisiologico invece, definito con la sigla , è stato implementato

attraverso una sommatoria di 5 sinusoidi, ciascuna con ampiezza, frequenza e

fase diverse, secondo la formula:

∑ ( )

(3.4)

dove , e sono rispettivamente l’ampiezza, la frequenza e la fase della

sinusoide i-esima, mentre è il periodo di campionamento.

I valori attribuiti alle diverse componenti rumorose sono riportati in Tabella

3.1 [16].

Frequency (Hz) Amplitude (nM)

Very low freq.

Low freq.

Vasomotor

Respiratory

Cardiac

Tabella 3.1: Frequenza e ampiezza delle componenti sinusoidali del modello di rumore fisiologico.

Le singole componenti del rumore fisiologico utilizzate nella simulazione sono

raggruppate in un’unica sommatoria che contiene il rumore fisiologico dovuto

al battito cardiaco, alla respirazione, alla pressione sanguigna (quest’ultima

definita onda di Mayer) e alle oscillazioni alle basse e bassissime frequenze.

Alcuni esempi di queste oscillazioni sono rappresentati nelle figure seguenti:


29

Figura 3.2: Rumore fisiologico simulato: andamento dell’onda di Mayer.

Figura 3.3: Rumore fisiologico simulato: andamento dell’onda lenta.

Figura 3.4: Rumore fisiologico simulato: andamento dell’onda molto lenta.


30

Figura 3.5: Rumore fisiologico simulato: andamento del rumore fisiologico dovuto alla respirazione.

Figura 3.6: Rumore fisiologico simulato: andamento del rumore fisiologico dovuto al battito cardiaco.


31

Figure 3.7: Rumore fisiologico simulato: grafico del rumore fisiologico simulato dato dalla sommatoria delle 5 componenti analizzate separatamente in precedenza (onda di Mayer, onda lenta e molto lenta, rumore cardiaco e respiratorio).

Figura 3.8: Esempio di segnale Long.


32

Figura 3.9: Esempio di segnale Short.

Infine, rumore bianco additivo e Gaussiano (con ampiezza 400 nM e

deviazione standard di 180 nM) è stato simulato per rappresentare il rumore

strumentale. La somma di rumore strumentale, rumore fisiologico e risposte

emodinamiche da origine ai canali LS, mentre un paio di canali SS per soggetto

sono stati simulati semplicemente sommando il rumore strumentale a quello

fisiologico, senza aggiungervi le HRF.

L’algoritmo è stato implementato inizialmente su 3 soli soggetti, al fine da

poterne estrarre i primi risultati ed ottenere le prime indicazioni. In

corrispondenza di ciascuna risposta evocata è stato implementato

parallelamente un segnale di trigger, costituito da impulsi aventi 5 diversi

possibili valori pari a 1, 2, 3, 4 o 5: i trigger pari a 1 e 2 indicano lo start di una

risposta emodinamica: il trigger 1 simula una HRF ottenuta mediante finger

tapping dell’indice della mano destra, ed è stata aggiunta solo in canali simulati

sopra l’emisfero sinistro. Viceversa, il trigger 2, indica una HRF osservabile nei

canali dell’emisfero destro, a seguito del finger tapping effettuato con l’indice

della mano sinistra. I trigger pari a 3, 4 o 5 non sono stati utilizzati.


33

Per ciascuno dei due trigger considerati, vi sono un totale di 10 canali: 5 a

risposta “piena” (360 ± 20 nM per il trigger 1 e 420 ± 20 nM per il trigger 2), 2

a risposta “ridotta” (180 ± 20nM per il trigger 1 e 210 ± 20 nM per il trigger

2), 1 a risposta nulla (dove non sono state inserite HRF), e due canali di

riferimento SS. Le numerazioni sono raffigurate nello schema seguente:

Numerazione canali LS

(piena, o “full”) LS

(ridotta o “half”) LS

(0 nM) SS

TRIGGER = 1 1-3-5-11-19 7-9 17 13-15

TRIGGER = 2 21-27-29-31-39 23-25 37 33-35

Tabella 3.2: Schema della numerazione dei canali utilizzato su ciascun soggetto

3.1.2 DATASET REALISTICO

Nel dataset realistico, sono stati simulati due tipi di risposte emodinamiche,

differenti nella forma, nell’ampiezza e nelle latenze. La prima tipologia prevede

una ampiezza di nM e una latenza di s, laddove la

seconda tipologia prevede una ampiezza di nM e una latenza di

s. La distanza tra una risposta e la successiva (indipendentemente

dalla loro tipologia) è definita come una variabile aleatoria uniforme

nell’intervallo s.

Un esempio di risposta emodinamica simulata è fornito in Figura 3.10:


34

Figura 3.10: Esempio di risposta emodinamica simulata con i seguenti parametri: , , , , , , [10].

Il rumore fisiologico è invece tratto direttamente da acquisizioni di segnali

fNIRS in condizioni di resting state, dove si suppone non siano presenti risposte

cerebrali evocate. Nei canali SS, nessuna HRF è stata aggiunta al segnale

misurato.

Il fatto di conoscere a priori la risposta emodinamica consente di valutare

l’effettiva capacità dell’algoritmo di rimuovere le componenti spurie e

ricostruire fedelmente l’andamento della risposta emodinamica. A tal fine, nel

capitolo successivo vengono definiti degli indici di performance in grado di

quantificare la bontà della procedura di de-noising e l’accuratezza con cui è

possibile stimare ampiezza e latenza della HRF stimata.

I dati reali di resting state sono stati acquisiti mediante un continuous wave NIRS

system (CW6, TechEn, Medford, MA, USA) ad una frequenza di

campionamento pari a 50 Hz, presso la Optics Division del Martinos Center

for Biomedical Imaging, Massachussets General Hospital, Boston, MA, USA.

La configurazione delle sorgenti e detettori è fornita nella Figura 3.11 [10]: i

nodi identificati da un numero rosso rappresentano le sorgenti, i nodi


35

identificati da un numero blu i detettorie le linee verdi indicano le coppie

sorgente-detettore che formano un canale di misurazione (52 canali totali). Si

noti la differenza tra le distanze sorgente-detettore: nei canali standard LS è

pari a circa 3 cm (visibile ad esempio tra la sorgente n° 1 e il detettore n° 1),

nei canali di riferimento SS è invece pari a 0.8 cm (si osservi la linea che

congiunge la sorgente n° 1 e il detettore n° 17). L’array di sorgenti e detettori

era posizionato sulla testa del soggetto in modo da coprire le aree motorie e

frontali bilateralmente.

Figura 3.11: Schema di configurazione delle sorgenti e dei detettori.

Lo schema riassuntivo di tutti i canali è rappresentato in figura 3.12, dove si

distinguono 36 canali Long e 16 Short.


36

Figura 3.12: Rappresentazione dei canali secondo la distanza sorgente-detettore: in verde i Long Channels, in rosso gli Short Channels

I dati raccolti sono relativi a un campione di sette soggetti sani (numerati

secondo la sequenza 9 - 12 - 13 - 14 - 15 - 16 - 17) a cui è stato chiesto di

sedere su di una sedia, in posizione comoda e rilassata e rimanere fermi con gli

occhi chiusi per circa 8 min. Per tutti i soggetti, ad eccezione del soggetto n°

15, l’acquisizione è stata ripetuta più di una volta, in sessioni distinte di resting

state. In tal modo, si sono ottenuti 14 dataset distinti di resting state.

Al fine di simulare uno scenario il più possibile prossimo a quello reale, ai

campioni dei canali LS relativi a periodi di resting state viene sommata una

risposta emodinamica simulata. Infatti, durante i periodi di resting state, è lecito

assumere che le variazioni di concentrazione di emoglobina ossigenata siano

principalmente dovute al rumore fisiologico e non ad attività cerebrale evocata,

data l’assenza di compiti specifici richiesti al soggetto. La risposta


37

emodinamica (HRF) viene generata da un’apposita funzione che implementa il

modello presentato in [16]. Nel caso in esame, ai segnali acquisiti sono

sommate 28 risposte HRF della durata di circa 15 secondi ciascuna. Un

esempio pratico è fornito dalla Figura 3.13 che mostra come le risposte non

siano identiche tra loro ma esibiscano una variabilità che riguarda non solo

l’ampiezza, ma anche la latenza del picco e la durata della risposta.

Figura 3.13: Sequenza di HRF sommate a ciascun canale LS acquisito in resting state.

I canali sono distinti in 3 gruppi: quelli dell’emisfero sinistro (sorgenti da 1 a

6), quelli dell’emisfero destro (sorgenti da 7 a 12) e quelli della parte frontale

(sorgenti da 13 a 15). Osservando la figura 3.11, si nota che la parte sinistra e

destra hanno pertanto lo stesso numero di canali (15), mentre i 6 rimanenti

sono attribuiti alla parte frontale.


38

Per quanto riguarda quest’ultimo, i parametri inseriti nel modello per la

simulazione delle HRF sono i seguenti:

Parametro Α Β p

Valore 1600

10 1.3 0.003

2 0.003

0.5 0.0001

1.8 0.0001

0.5

0.0001 2

Nell’emisfero destro sono stati mantenuti gli stessi parametri del precedente ad

eccezione di α, che presenta invece un valore più basso:

Parametro Α Β p

Valore 1600

10 1.3 0.003

2 0.003

0.5 0.0001

1.8 0.0001

0.5

0.0001 2

Essendo α il parametro che determina l’ampiezza delle risposte HRF, l’altezza

dei picchi è diversa a seconda che si tratti di un canale appartenente

all’emisfero destro o sinistro: nel primo infatti è stata mantenuta più bassa,

ovvero circa 75 nM, mentre nel sinistro più elevata, pari a circa 250 nM. A

ciascun parametro, è stata inoltre sommata una componente casuale a

distribuzione normale, al fine di simulare anche la variabilità fisiologica delle

HRF, tipica del caso reale dove le risposte non risultano mai identiche.

Nei canali della parte frontale invece, non sono state inserite risposte

emodinamiche: in tal modo possiamo testare la capacità dell’algoritmo di non

stimare HRF quando queste non sono presenti, ovvero di non generare

attivazioni spurie. Inoltre, si può testare la capacità dell’algoritmo nel ridurre il

rumore fisiologico di fondo.

Nelle figure 3.14 e 3.15 è rappresentata la media delle HRF inserite, sia per

l’emisfero sinistro che destro:


39

Figura 3.14: Media delle HRF inserite nei canali dell’emisfero sinistro.

Figura 3.15: Media delle HRF inserite nei canali dell’emisfero destro.

41

CAPITOLO 4

ALGORITMO CS-TFT

4.1 PRE-PROCESSING

Una volta caricati i dati relativi al soggetto e alle condizioni di test desiderate

[righe 29-89], si procede alla conversione dell’unità di misura dei campioni

dalla intensità luminosa rilevata dal sensore NIRS alle nano-moli (nM),

utilizzando la legge modificata di Beer-Lambert. Così facendo, la misura

riguarda non più i fotoni luminosi recepiti dal sensore, ma la variazione di

concentrazione dell’emoglobina ossigenata (HbO2) nella zona analizzata.

Ogni canale di acquisizione è formato da una coppia sorgente-detettore. A

seconda della distanza tra questi dispositivi è possibile distinguere i canali short

separation (~0.8 cm), o brevemente canali SS, dai canali long separation (~ 3 cm),

o brevemente canali LS. Nei primi, il segnale acquisito contiene solo rumore

fisiologico e non attività emodinamica cerebrale. Nei secondi, al contrario, è

presente anche la componente informativa, cioè la HRF. Per discriminare gli

uni dagli altri, si è misurata la distanza sorgente-detettore, si è calcolato il

valore medio tra massima e minima distanza e si è utilizzato questo valore

come soglia: i canali con una distanza sorgente-detettore superiore alla soglia

sono classificati come LS, viceversa i canali con una distanza sorgente-

detettore inferiore alla soglia sono classificati come SS [righe 141-157].

ALGORITMO CS-TFT

42

Ogni canale LS è stato associato al canale SS maggiormente correlato.

L’intento è quello di individuare il canale SS i cui campioni siano

un’approssimazione quanto più possibile verosimile del rumore fisiologico

presente nel canale LS in esame. Dal punto di vista implementativo, tale

correlazione viene calcolata attraverso la somma delle differenze in valore

assoluto dei due segnali. Ad ogni canale LS viene associato il canale SS con lo

scostamento minimo [righe 168-174].

I segnali relativi ai canali SS e LS vengono opportunamente filtrati [righe 231-

244] e il segnale SS viene poi scalato per renderlo di ampiezza simile a quello

LS [righe 288-302]. In particolare, viene applicato un filtro passa-banda

nell’intervallo 0.01-0.55 Hz, per eliminare le alte frequenze e i trend alle basse

frequenze, pur preservando la HRF, che ha un’estensione spettrale centrata

intorno a 0.08 Hz.

L’effetto del filtro passa-banda è quello infatti di ridurre significativamente la

distorsione introdotta da una parte del rumore fisiologico, specialmente quello

indotto dal battito cardiaco, che ha una frequenza media di circa 1.1 Hz

nell’adulto. A titolo di esempio, le Figure 4.1 – 4.2 mostra l’effettiva riduzione

del rumore alle alte e bassissime frequenze prodotta dal filtro: con la rimozione

della continua e delle frequenze più basse, il segnale risulta mediamente

centrato attorno allo 0; con il taglio delle frequenze più elevate, anche le

oscillazioni spurie si riducono.

Una volta filtrati i segnali, si rende necessario scalare il canale SS su quello LS

per compensare il diverso percorso intrapreso dai fotoni nei due canali a

lunghezza diversa. Nel caso in esame, la cosiddetta procedura di scaling viene

realizzata mediante un algoritmo di fit polinomiale lineare che consente di

definire non solo il fattore di scala ma anche un eventuale offset additivo.

È bene sottolineare che per ogni coppia di canali LS e SS, la procedura di

scaling viene effettuata su una ridotta porzione di segnale acquisito ad inizio

esperimento e prima della presentazione al soggetto di alcuno stimolo; in

entrambi i canali è lecito quindi attendersi l’assenza di risposta cerebrale.

Altrimenti, ci sarebbe il rischio che l’algoritmo di fit compensi eccessivamente

l’ampiezza del canale SS, fittando anche la HRF, e vanifichi tutti i passaggi

successivi.

ALGORITMO CS-TFT

43

Figura 4.1: Confronto tra segnale originale e segnale filtrato con filtro passa-banda

Figura 4.2: Confronto tra una porzione del segnale originale e la stessa del segnale filtrato

con filtro passa-banda

ALGORITMO CS-TFT

44

In Figura 4.3, si può notare che la fase di fit permette di sovrapporre quasi

perfettamente i canali LS e SS.

Figura 4.3: Confronto tra il canale LS e quello SS sia scalato che non scalato

Una volta filtrato e scalato, il segnale di ogni canale viene suddiviso in trials,

basandosi sul segnale di trigger disponibile. In corrispondenza di ogni trigger,

viene selezionata la porzione di segnale che va da 1 secondo prima del trigger a

13 secondi dopo il trigger, per una durata complessiva di 14 secondi.

Giunti a questo punto, l’algoritmo prevede quattro ulteriori fasi di

elaborazione [14]: l’identificazione delle frequenze delle componenti di rumore

fisiologico, la definizione di un modello di Taylor-Fourier dello stesso rumore,

la sottrazione del rumore e la stima della risposta HRF media. Tale

procedimento è schematizzato in Figura 4.4 ed è illustrato in dettaglio nelle

sezioni che seguono.

ALGORITMO CS-TFT

45

Figura 4.4: Schema dell’algoritmo CS-TFT per la stima della risposta emodinamica [5].

La fase di stima spettrale si basa su una tecnica di super-risoluzione, che mira a

superare i limiti della trasformata discreta di Fourier (DFT) causati dal

troncamento a N coefficienti, e a ridurre di almeno un ordine di grandezza

l’incertezza sulla frequenza. La tecnica in questione sfrutta i principi della

teoria matematica del campionamento compresso: pur lavorando su un

ristretto insieme di campioni, si è in grado di fornire stime con

un’accuratezza pari a quella che si avrebbe lavorando con

campioni. In questo contesto, il parametro viene identificato con il nome di

fattore di super-risoluzione [18].

Evidentemente, passare da a campioni, con è possibile solo sotto

opportune ipotesi matematiche. Matematicamente, questo passaggio fornisce

una soluzione solo se il segnale in esame è sparso nel dominio della frequenza,

ovvero esibisce un numero di coefficienti non nulli o comunque significativi

che è molto inferiore alla sua cardinalità. In tal caso, la teoria del

campionamento compresso garantisce che è possibile arrivare ad una

soluzione risolvendo un problema di minimizzazione di norma . Il problema

non prevede una soluzione in forma chiusa, ma può essere risolto mediante

tecniche di ottimizzazione convessa o tecniche iterative cosiddette greedy, come

ALGORITMO CS-TFT

46

nel caso in esame dove si applica una versione dell’algoritmo Orthogonal

Matching Pursuit [19].

Grazie all’applicazione di questa tecnica di super-risoluzione, pur lavorando su

campioni, si riesce a fornire una stima della DFT su una griglia di

campioni, risultanti cioè da un fattore . La scelta di questi

parametri è dettata dal fatto che, campionando a 50 Hz, ciò consente di

passare da una risoluzione di mHz a mHz. In particolare, la scelta di

un numero dispari di campioni garantisce che i coefficienti della DFT del

segnale si dispongano su una griglia frequenziale con simmetria centrale.

Inoltre, la scelta di un fattore di interpolazione che sia co-primo consente di

minimizzare la probabilità di casi di mal condizionamento numerico nel

calcolo dei nuclei di Dirichlet.

Nel caso in esame, sono stati acquisiti circa 8 minuti si segnale su ciascun

soggetto, ed ogni canale consta di un numero medio di 24239 campioni.

L’algoritmo viene applicato al segnale SS associato al canale LS in questione.

Applicando la DFT su porzioni consecutive di campioni, si

ottengono diversi segmenti, che matematicamente possono essere

rappresentati da un vettore colonna , con una leggera sovrapposizione (pari a

5 campioni) tra un segmento e il successivo.

La DFT del vettore delle misure risulta dalla seguente formula matriciale:

(4.1)

dove è la matrice degli esponenziali complessi della DFT. L’idea alla

base della tecnica di super-risoluzione è quella di sostituire la matrice con

una matrice che renda conto sia del passaggio dal dominio del tempo a quello

della frequenza sia della proiezione su una griglia composta da un numero

superiore di punti. Allo scopo l’equazione (4.1) viene riformulata come:

(4.2)

ALGORITMO CS-TFT

47

dove è la matrice di dimensione le cui colonne contengono dei

nuclei Dirichlet centrati nelle frequenze della griglia più risoluta, i.e. in

con , mentre il vettore di dimensione

assume valori non nulli o comunque significativi solo in corrispondenza delle

attuali frequenze delle componenti rumorose.

A titolo di esempio, in Figura 4.5 è fornita una rappresentazione

tridimensionale della matrice D per e . In un certo senso, la

matrice può essere interpretata come una sorta di dizionario, da cui estrarre

singoli vettori colonna che sommati assieme approssimino lo spettro del

segnale in esame. I singoli elementi della matrice possono essere calcolati dalla

seguente espressione:

( ) [ (

)]

[ (

)]

(

)

(4.3)

dove il primo termine rappresenta la versione campionata del nucleo di

Dirichlet, mentre l’esponenziale complesso consente di traslarlo nel punto

desiderato della griglia più risoluta.

ALGORITMO CS-TFT

48

Figura 4.5: Rappresentazione tridimensionale della matrice D nel caso di e .

Per poter applicare la teoria del campionamento compresso è necessario

ipotizzare che il vettore sia sparso, ossia abbia solo pochi coefficienti non

nulli, ipotesi plausibile nel caso in esame dal momento che ci attendiamo che il

rumore fisiologico consti di al massimo 5 componenti spettrali.

In tale contesto, il sistema si presenta sotto-determinato perché prevede un

numero di incognite molto maggiore del numero di dati e quindi ammette

un numero infinito di soluzioni. Ciò nonostante, sfruttando l’ipotesi di sparsità

del vettore è possibile definire la soluzione del problema come il vettore

con minima norma :

(4.4)

ALGORITMO CS-TFT

49

A tal proposito, è bene sottolineare che la norma è definita come la distanza

in modulo tra i due vettori, mentre la norma quantifica il numero di

coefficienti non nulli della differenza dei due vettori. Per questo motivo, la

norma viene usata per definire il livello di sparsità del segnale, ossia quale sia

il massimo livello di compressione dell’informazione in esso contenuta. Dal

punto di vista formale, la miglior soluzione al problema formulato

nell’equazione 4.4 prevedrebbe la minimizzazione di una norma , tuttavia dal

punto di vista matematico un simile modo di procedere risulta improponibile

in pratica: infatti, il problema presenta complessità non polinomiale e porge

infinite soluzioni. Una plausibile alternativa è fornita dalla norma che pur

non garantendo la massima sparsità restituisce una versione ragionevolmente

compressa del segnale iniziale e dispone di numerosi algoritmi ottimizzati per

la soluzione del problema. Per esempio, in questo caso l’equazione 4.4 viene

risolta mediante l’approccio iterativo noto come Orthogonal Marching Pursuit [19]

che si arresta quando l’energia dei residui della soluzione corrente è inferiore a

una determinata soglia o quando la soluzione corrente ha raggiunto il livello di

sparsità desiderato.

Le frequenze delle componenti significative vengono stimate a partire dagli

indici non nulli del vettore soluzione :

‖ ‖

(4.5)

Stimate le frequenze delle componenti di rumore fisiologico del segnale, è

necessario stabilire con accuratezza la loro ampiezza e fase iniziale. Qualsiasi

algoritmo che adotti una trasformata discreta di Fourier, basandosi su un

modello statico di segnale, fornisce solo una stima mediata sull’intervallo di

osservazione considerato. Nel caso in esame, se il segnale non si mantiene

stazionario nell’arco della finestra di campioni, le stime di ampiezza e fase

iniziale cercheranno di fornire una media dei valori assunti ma non

garantiranno una accuratezza assoluta, in quanto le ipotesi della trasformata

non sono soddisfatte in toto.

ALGORITMO CS-TFT

50

Dal punto di vista dell’analisi funzionale, la trasformata DFT può essere

rappresentata dal suo nucleo o kernel, che viene formulato come segue:

(4.6)

dove è il passo di campionamento. È immediato constatare come se si cerca

di rappresentare un segnale tempo variante solo con questi nuclei,

l’informazione relativa alle dinamiche del segnale venga persa o comunque

distorta e influisca sulla accuratezza delle stime delle componenti statiche.

Per superare questo limite, si ricorre alla cosiddetta trasformata di Taylor-

Fourier (nel seguito indicata con l’acronimo TFT), che svolge un’espansione in

serie di Taylor-Mclaurin della DFT canonica. Analogamente all’equazione 4.6,

il nucleo della TFT può essere espresso come:

( )

(4.7)

dove denota il massimo ordine di derivazione adottato nell’espansione.

Grazie a questa trasformata, si ottiene una miglior accuratezza nella

ricostruzione del rumore fisiologico dovuta al fatto che ciascuna componente

del canale SS viene rappresentata mediante K vettori di coefficienti TFT

linearmente indipendenti tra loro e ciascuno rappresentante peculiari

caratteristiche del segnale. In particolare, i termini di ordine rendono

conto della parte stazionaria del segnale, mentre i termini di ordine superiore

sono in grado di rappresentare al meglio eventuali variazioni, siano esse

regolari e periodiche o transitorie e limitate nel tempo. Un buon compromesso

tra complessità computazionale e accuratezza delle stime è stato individuato in

, ossia in un’espansione in serie di Taylor che si arresti al secondo

ordine di derivazione. Gli eventuali andamenti lineari o rapide oscillazioni

dovute a funzionamenti non lineari dei sensori o a piccoli spostamenti dei

sensori vengono perlopiù proiettati sui termini di ordine superiore e

ALGORITMO CS-TFT

51

consentono una stima più accurata della reale ampiezza delle specifiche

sorgenti di rumore.

Per facilitare il calcolo della TFT, i nuclei definiti dall’equazione 4.7 vengono

usati come colonne di una matrice dove per ogni punto della griglia super-

risoluta vengono definiti vettori. Il vettore dei coefficienti TFT è il

risultato della stima ai minimi quadrati:

( ) (4.8)

dove è il vettore delle misure della partizione di segnale SS in esame e

indica la matrice trasposta hermitiana di B [riga 366]. Il passaggio dai

coefficienti alla stima del rumore fisiologico è immediato e altro non è che la

combinazione lineare dei nuclei TFT opportunamente pesati per il

vettore dei coefficienti [righe 350-354]:

[∑ ∑

] (4.9)

Questo modello del rumore fisiologico può essere sottratto dal rispettivo

segmento nel canale LS. Tale sottrazione non garantisce però la soppressione

di tutto il rumore presente nel canale. Pertanto, si individuano tutti i trial

relativi al medesimo tipo di stimolo e li si media tra di loro, per stimare la

risposta emodinamica media. Tutti i trial sono stati corretti per la baseline,

sottraendo la media del segnale in un intervallo centrato intorno all’istante di

presentazione dello stimolo (da -1 a 1 s), in modo che in tutti i trial la HRF

partisse da 0.

ALGORITMO CS-TFT

52

Figura 4.6: Long Signal raffigurato prima (blu) e dopo la sottrazione del rumore ricostruito mediante algoritmo CS-TFT (verde)

L’operazione di media aritmetica fornisce una stima media della risposta

emodinamica del soggetto al dato stimolo, il cui livello di regolarità dipende dal

numero di trial su cui si effettua la media. Non volendo aumentare a dismisura

questo numero, così da evitare acquisizioni troppo lunghe che possono

indurre fenomeni di stanchezza o abituazione nel soggetto, eventuali

oscillazioni residue sono eliminate mediante un filtraggio di Savitzky-Golay di

ordine 3 e lunghezza pari a 25 campioni.

ALGORITMO CS-TFT

53

Figura 4.7: Media delle HRF dell’emisfero destro prima (rosso) e dopo il filtraggio mediante il filtro di Savitzky-Golay (blu).

4.2 MEDIA DELLE HRF

Su ciascun canale poi, è stata calcolata la media delle HRF ricostruite mediante

algoritmo CS-TFT, mantenendo separata la distinzione ancora una volta tra i

canali dell’emisfero destro e sinistro [righe 472-474 e 477-479]. Infine, le HRF

medie appartenenti a canali dello stesso emisfero sono state mediate tra di

loro, per ottenere una singola HRF media rappresentativa di una delle 3 aree di

interesse: area motoria destra [righe 540-542], area motoria sinistra [righe 536-

538]e area frontale [righe 490-496]. Tali medie sono state confrontate sia a

livello grafico che numerico con le medie delle HRF vere ottenute nelle stesse

aree, e allo stesso tempo anche con quelle ottenute sia mediante ReMCoBA

ALGORITMO CS-TFT

54

che con il metodo di filtraggio passabanda senza uso del canale SS (BP W/O

SS).

4.3 ALCUNI DETTAGLI IMPLEMENTATIVI

Come descritto in particolare nel capitolo 3, per quanto riguarda lo scenario

semi-simulato, sono state aggiunte 28 risposte emodinamiche simulate, ad

intervalli di circa 15 secondi. Il numero 28 non è scelto tuttavia a caso, in

quanto si è deciso di aggiungere le HRF fino al termine del segnale di resting

state (eccezion fatta per i primi 30 secondi di segnale, necessari per lo scaling,

mantenuti originali e su cui nessuna risposta è stata aggiunta), in modo tale che

la media risultante fosse calcolata sul numero massimo di risposte

emodinamiche inseribili.

Le caratteristiche peculiari del codice sono contraddistinte dai seguenti

parametri in particolare:

– N = numero di campioni della griglia originale = 501

– P = fattore di interpolazione della griglia super-risoluta = 11

– Numero di campioni di overlap = 5

– Numero di ordine massimo della serie di Taylor = 2

Tali valori non sono stati modificati particolarmente nello sviluppo del codice,

ma sono rimasti inalterati durante tutta la fase di implementazione fin qui

realizzata, in quanto, come primo studio, si è cercato principalmente di creare

un primo algoritmo in grado di filtrare il segnale NIRS e stimare la risposta

emodinamica, con risultati soddisfacenti rispetto agli altri metodi e tempi

computazionali non eccessivi. I parametri elencati, pertanto, sono da

considerarsi come valori di “default”, ma nessuno può attualmente garantire

che, cambiando anche uno solo di essi, si raggiungano risultati sicuramente

peggiori.

ALGORITMO CS-TFT

55

Valutando le prestazioni del codice dal punto di vista del costo

computazionale, emerge in particolare che i comandi Matlab più onerosi

dell’intero programma sono esattamente due, eseguiti in sequenza uno dopo

l’altro. Essi si riferiscono in particolare alla creazione della matrice super-

risoluta D di dimensioni NxH (dove N = 501 e H = 11*501 = 5511), e alla

norma calcolata sulla diagonale della matrice DT*D, necessaria alla costruzione

della griglia super-risoluta, ovvero:

D = creaMatriceD(N,H) (4.10)

NormACols = sqrt(diag(D'*D)) (4.11)

Attualmente, questi comandi vengono calcolati ogni volta che l’algoritmo

esamina un canale diverso dal precedente, ovvero per un totale di 36 volte nel

nostro caso specifico, pari al numero di canali, per ogni soggetto. I tempi di

tale calcolo, considerato come riferimento un computer di medie prestazioni,

risultano circa di 5 e 7 secondi ciascuno. La creazione della matrice D su

ciascun canale è stata inizialmente mantenuta in virtù di possibili modifiche

alle dimensioni della stessa, con l’ipotesi iniziale che N e H potessero appunto

variare a seconda della tipologia di canale o dell’SNR del segnale analizzato.

Tuttavia, tale ipotesi non è stata ancora attualmente considerata o

implementata, e i valori di N e H risultano inalterati sia in tutti i canali, che in

ciascuno dei 14 data set su cui l’algoritmo è stato testato. Sempre considerando

la matrice D, un ulteriore miglioramento in campo computazionale

deriverebbe dal fatto che la gamma di frequenze attualmente implementata

nella matrice è eccessivamente ampia rispetto alla banda poi effettivamente

considerata: il range su cui varia la matrice super-risoluta va infatti da 0 Hz a

7.8125 Hz, mentre le effettive frequenze di interesse, tra cui quella della

risposta emodinamica, variano in un intervallo molto più ristretto, ovvero da

0.01 Hz a 0.55 Hz (come si nota anche dal filtraggio passabanda che viene

implementato sul segnale prima dell’avvio dell’algoritmo).

57

CAPITOLO 5

RISULTATI

5.1 SCENARIO COMPLETAMENTE SIMULATO

L’algoritmo CS-TFT è relativamente recente e per questo in questo capitolo

viene testato su segnali creati “artificialmente” mediante funzioni random e

componenti sinusoidali, al fine di riprodurre un segnale NIRS realistico.

I segnali NIRS analizzati e confrontati nell’implementazione dell’algoritmo,

sono essenzialmente due: il Long Signal (LS) e il Short Signal (SS). Il primo è

una somma di varie componenti, quali il rumore fisiologico (dato da una

somma di sinusoidi), il rumore di misura, e le HRF, mentre il secondo

comprende solamente il rumore fisiologico e di misura (non include pertanto

le risposte emodinamiche). La struttura dei due canali, come già anticipato nel

capitolo 2, è la seguente:

(5.1)

(5.2)

dove sono le risposte emodinamiche, il rumore fisiologico e

il rumore di misura.

RISULTATI

58

Figura 5.1: Esempio di segnale LS (sopra) e SS (sotto) utilizzati nello scenario completamente simulato.

L’algoritmo è stato inizialmente testato scansionando un solo canale alla volta,

dove l’utente poteva inserire il numero del canale Long e Short su cui

effettuare il calcolo, oltre ovviamente al valore di trigger.

5.1.1 DETTAGLI IMPLEMENTATIVI

Essendo le risposte emodinamiche inserite a cadenza non costante, il segnale è

stato partizionato secondo un numero di campioni sempre diverso, dove

ciascun segmento di segnale poteva contenere, in certi casi, anche più di una

singola onda HRF, specie se molto ravvicinate tra loro. Le partizioni sono

state realizzate facendo riferimento al segnale di trigger, dove lo start

RISULTATI

59

corrispondeva ad un valore di trigger pari a 1, mentre lo stop al successivo

trigger sempre uguale ad 1, ma con la condizione che tra i due valori di trigger

di start e stop e pari a 1 non ce ne fossero altri di diverso valore, quali ad

esempio 3, 4 o 5. Al fine di rendere poi, al momento del calcolo della griglia

super-risoluta, le porzioni tutte di uguale lunghezza, si è fatto uso dello zero-

padding, prendendo come riferimento la porzione più lunga, e adattando le

altre alla stessa lunghezza mediante aggiunta di zeri, ininfluenti alla stima.

Figura 5.2: confronto tra la HRF vera (blu) e quella ricostruita mediante algoritmo CS-TFT

(rossa).

RISULTATI

60

5.1.2 STIMA DELLE HRF

Il risultato illustrato in Figura 5.2, che confronta l’onda HRF stimata con

quella vera: vi è la totale corrispondenza temporale dei rispettivi picchi su

ciascuna risposta emodinamica. Meno ottimale risulta invece la corrispondenza

in ampiezza, in quanto la HRF stimata risulta, nei picchi, sempre inferiore

rispetto a quella vera: in alcuni casi, si nota dal grafico una differenza di anche

70 nM. Nella parte inferiore dell’onda invece, ovvero per valori negativi di

HbO2, l’onda CS-TFT sottostima sempre quella vera, ma questa volta con una

differenza minore e a volte quasi trascurabile. Potremmo pertanto ipotizzare

una cattiva ricostruzione dell’intercetta, probabilmente dovuta allo scaling tra i

due segnali (quello LS e quello SS) effettuato subito dopo il filtraggio, in

quanto l’onda vera sembra essere leggermente “shiftata” verso l’alto rispetto a

quella ricostruita, cosa che avviene su tutti i picchi analizzati. Un’altra ipotesi è

che la costante sottostima sia dovuta ad una rimozione di parte del segnale

utile durante le operazioni di sottrazione del rumore fisiologico.

Essendo in una condizione ideale di totale simulazione, importante risulta

anche l’analisi sull’efficacia della rimozione del rumore fisiologico. Da figura

5.2, si desume come tale sottrazione non sia sufficiente a rimuovere

definitivamente le onde fisiologiche presenti nel segnale di resting state, cioè

dove nessuna risposta emodinamica era presente. Il grafico denota pertanto,

che una volta finita la HRF, il segnale non diventa mai nullo, ma sembra

oscillare attorno allo zero con una sinusoide piuttosto regolare.

5.1.3 STIMA DEL RUMORE

L’obiettivo dell’algoritmo, come precedentemente menzionato, è la

ricostruzione del rumore di fondo, al fine da poterlo poi rimuovere, in un

secondo momento, dal canale Long. Confrontando il rumore stimato con

quello vero lungo tutto il segnale (Figura 5.3), emerge che l’algoritmo fornisce

in uscita un segnale che sembra fittare molto bene il rumore vero, evitando di

fittare le alte frequenze, frutto del rumore di misura. Infatti, l’idea di utilizzare

RISULTATI

61

un modello che fitti solo particolari frequenze presenti nel segnale SS (cioè

quelle fisiologiche), per poi sottrarlo al segnale LS, invece di sottrarre

direttamente il segnale SS, deriva proprio dalla necessità di evitare di introdurre

ulteriore rumore nel canale LS, soprattutto rumore alle alte frequenze.

L’algoritmo qui sviluppato quindi raggiunge bene l’obiettivo, stimando solo la

parte del segnale SS di interesse.

Figura 5.3: Confronto tra il rumore vero (rosso) e quello ricostruito mediante algoritmo (blu).

RISULTATI

62

Figura 5.4: Dettaglio del confronto tra il rumore vero (blu) e quello ricostruito mediante algoritmo (rosso), raffigurato in un range temporale di 100 secondi

In virtù di tale risultato, si è deciso di osservare i residui, dati dalla sottrazione

tra il rumore vero e quello stimato (Figura 5.5):

RISULTATI

63

Figura 5.5: Grafico dei residui del rumore stimato rispetto a quello vero. La media di essi, raffigurata in rosso, è pari a 0.16773

Ad eccezione di alcuni picchi ben visibili, i residui denotano un’ampiezza

abbastanza costante, la cui media totale è approssimabile allo zero (~0.17).

Osservando anche in questo caso una porzione di tali valori (Figura 5.6), si

osserva ancora una volta un andamento regolare, simile ad una sinusoide.

RISULTATI

64

Figura 5.6:. Residui del rumore stimato rispetto a quello vero, raffigurati in un range temporale di 100 secondi

Osservando il particolare trend dell’onda residua, si è deciso di calcolarne la

trasformata di Fourier discreta (DFT) al fine di verificare se vi fosse una

frequenza dominante, ed in caso affermativo confrontarla con la stessa

frequenza delle principali componenti del rumore fisiologico inserito

manualmente.

Come si evince dal grafico della trasformata (Figura 5.7), i residui hanno un

andamento che può propriamente essere assimilabile ad una sinusoide, in

quanto costituiti da una frequenza principale pari a circa 0.23 Hz. L’unica

componente di rumore fisiologico che ha frequenza simile è quella respiratoria

(0.2 0.03 Hz).

RISULTATI

65

Figura 5.7: Grafico della DFT dei residui.

Potremmo dire pertanto che l’algoritmo fatica in questo caso a stimare almeno

una componente del rumore fisiologico, che sembrerebbe principalmente

quella causata dalla respirazione.

5.1.4 CONFRONTO COL FILTRAGGIO PASSABANDA

Con lo scopo di analizzare i primi risultati forniti dal CS-TFT, si è deciso di

confrontarli con gli stessi dati dalla media convenzionale realizzata a seguito

del filtraggio BP già menzionato nel capitolo 2. Come enunciato in precedenza

, è stata mantenuta la separazione nella ricostruzione e nel confronto delle

HRF distinte sia per gli emisferi destro e sinistro, che sia per ampiezza, ovvero

piene (chiamate “Full” nei grafici sottostanti) e ridotte (“Half”).

RISULTATI

66

A)

B)

C)

D)

Figura 5.8: HRF dell’emisfero sinistro: A) HRF a risposta piena a confronto: media delle HRF stimate mediante algoritmo CS-TFT

(rossa), media delle HRF vere simulate (nera), HRF medie di ogni soggetto ottenute mediante algoritmo CS-TFT (verdi).

B) Media delle HRF a risposta piena a confronto: HRF stimate mediante algoritmo CS-TFT (rossa), filtraggio passabanda BP (blu), confrontata con la media delle HRF vere (nero).

C) HRF a risposta ridotta a confronto: media delle HRF stimate mediante algoritmo CS-TFT (rossa), media delle HRF vere simulate (nera), HRF medie di ogni soggetto ottenute mediante algoritmo CS-TFT (verdi).

D) Media delle HRF a risposta ridotta a confronto: HRF stimate mediante algoritmo CS-TFT (rossa), filtraggio passabanda BP (blu), confrontata con la media delle HRF vere (nero).

Come si vede in Figura 5.8 in entrambi i casi, sia per le HRF piene (riquadri A

e B) che ridotte (quadri C e D), l’algoritmo CS-TFT conduce a risultati migliori

RISULTATI

67

rispetto al BP. Accurata sembra essere la stima in latenza, mentre vi è una

leggera sottostima, da parte del CS-TFT rispetto alla HRF vera, per quanto

riguarda l’ampiezza del picco.

A)

B)

C)

D)

Figura 5.9: HRF dell’emisfero destro: A) HRF a risposta piena a confronto: media delle HRF stimate mediante algoritmo CS-TFT

(rossa), media delle HRF vere simulate (nera), HRF medie di ogni soggetto ottenute mediante algoritmo CS-TFT (verdi).

B) Media delle HRF a risposta piena a confronto: HRF stimate mediante algoritmo CS-TFT (rossa), filtraggio passabanda BP (blu), confrontata con la media delle HRF vere (nero).

C) HRF a risposta ridotta a confronto: media delle HRF stimate mediante algoritmo CS-TFT (rossa), media delle HRF vere simulate (nera), HRF medie di ogni soggetto ottenute mediante algoritmo CS-TFT (verdi).

D) Media delle HRF a risposta ridotta a confronto: HRF stimate mediante algoritmo CS-TFT (rossa), filtraggio passabanda BP (blu), confrontata con la media delle HRF vere (nero).

RISULTATI

68

Per quanto riguarda la stima delle HRF nell’emisfero destro (Figura 5.9), si

osserva un fit ottimo per quanto riguarda il CS-TFT rispetto alle HRF vere,

mentre particolarmente instabile risulta la parte discendente dell’onda,

costituita da evidenti oscillazioni presenti sia sulle HRF piene che ridotte.

Per confrontare le performance dei metodi in maniera quantitativa si sono

calcolati degli indici di errore di stima. Questo è stato possibile poiché ci si

trova in uno scenario simulato dove le HRF vere sono note. In tal modo è

stato possibile quantificare l’accuratezza di ciascun metodo nella stima

dell’andamento temporale delle risposte emodinamiche e nella stima di alcune

metriche derivanti da essa, come l’ampiezza e la latenza del picco massimo. In

particolare la stima dell’errore percentuale in ampiezza, è calcolata mediante la

formula:

| ( ) ( )|

| ( )| (5.3)

dove è il segnale vero (ovvero le HRF simulate note), mentre è la

HRF stimata da uno degli algoritmi. Come si desume dalla formula, l’errore è

stimato in percentuale e considera il rapporto tra la differenza in ampiezza dei

rispettivi picchi (val_picco) e il valore di picco del segnale vero.

E_amp % (LEFT side, FULL HRF)

CS-TFT BP

[W/O SS]

E_amp % (LEFT side, HALF HRF)

CS-TFT BP

[W/O SS]

Subj1 10,04 3,74 Subj1 9,47 6,02

Subj2 25,05 50,36 Subj2 8,10 50,34

Subj3 16,46 5,10 Subj3 17,02 14,53

Tabella 5.1: Tabella errore in ampiezza emisfero sinistro. Tabella contenente i valori di errore in ampiezza delle risposte emodinamiche stimate rispetto a quelle vere, sia per le piene che per le ridotte, considerando solo quelle stimate nei canali dell’emisfero sinistro, e calcolate distintamente su ciascun volontario.

RISULTATI

69

E_amp % (RIGHT side, FULL HRF)

CS-TFT BP

[W/O SS]

E_amp % (RIGHT side, HALF HRF)

CS-TFT BP

[W/O SS]

Subj1 12,26 23,89 Subj1 13,10 36,36

Subj2 12,57 8,49 Subj2 66,34 46,97

Subj3 7,45 19,28 Subj3 10,23 43,43

Tabella 5.2:. Tabella errore in ampiezza emisfero destro. Tabella contenente i valori di errore in ampiezza delle risposte emodinamiche stimate rispetto a quelle vere, sia per le piene che per le ridotte, considerando solo quelle stimate nei canali dell’emisfero destro, e calcolate distintamente su ciascun volontario.

L’errore percentuale nella stima della latenza del picco, misurando pertanto

l’eventuale “sfasamento temporale” tra la HRF vera e quella stimata è calcolato

come:

| ( ) ( )|

| ( )| (5.4)

dove t_picco indica l’istante temporale (in secondi) del picco massimo.

E_lat % (LEFT side, FULL HRF)

CS-TFT BP

[W/O SS]

E_lat % (LEFT side, HALF HRF)

CS-TFT BP

[W/O SS]

Subj1 0,50 5,58 Subj1 0,00 10,21

Subj2 6,07 9,10 Subj2 6,22 13,69

Subj3 12,75 9,19 Subj3 7,66 3,83

Tabella 5.3: Tabella errore in latenza emisfero sinistro. Tabella contenente i valori di errore

in latenza delle risposte emodinamiche stimate rispetto a quelle vere, sia per le piene che per

le ridotte, considerando solo quelle stimate nei canali dell’emisfero sinistro, e calcolate

distintamente su ciascun volontario.

RISULTATI

70

E_lat %

(RIGHT side, FULL HRF)

CS-TFT BP

[W/O SS]

E_lat %

(RIGHT side, HALF HRF)

CS-TFT BP

[W/O SS]

Subj1 0,45 0,46 Subj1 0,00 1,16

Subj2 6,95 8,80 Subj2 25,47 27,79

Subj3 2,29 2,32 Subj3 2,32 2,32

Tabella 5.4: Tabella errore in latenza emisfero destro. Tabella contenente i valori di errore in

latenza delle risposte emodinamiche stimate rispetto a quelle vere, sia per le piene che per le

ridotte, considerando solo quelle stimate nei canali dell’emisfero destro, e calcolate


Infine, al fine di determinare con accuratezza l’errore commesso nella stima

dell’andamento temporale della HRF, si è calcolato un altro indice di errore,

l’errore totale, mediante la formula:

(5.5)

In termini pratici, l’errore totale considera la differenza, campione per

campione, tra le due risposte emodinamiche, e ne calcola la norma quadratica.

E_tot % (LEFT side, FULL HRF)

CS-TFT BP

[W/O SS]

E_tot % (LEFT side, HALF HRF)

CS-TFT BP

[W/O SS]

Subj1 1,21 10,70 Subj1 1,47 34,30

Subj2 11,33 26,06 Subj2 14,62 43,81

Subj3 11,07 16,25 Subj3 11,88 31,76

Tabella 5.5: Tabella errore totale emisfero sinistro. Tabella contenente i valori di errore

totale delle risposte emodinamiche stimate rispetto a quelle vere, sia per le piene che per le

ridotte, considerando solo quelle stimate nei canali dell’emisfero sinistro, e calcolate


RISULTATI

71

E_tot %

(RIGHT side, FULL HRF)

CS-TFT BP

[W/O SS]

E_tot %

(RIGHT side, HALF HRF)

CS-TFT BP

[W/O SS]

Subj1 2,73 9,75 Subj1 3,60 25,18

Subj2 14,65 14,99 Subj2 122,37 77,12

Subj3 6,59 16,25 Subj3 56,28 101,54

Tabella 5.6: Tabella errore totale emisfero destro. Tabella contenente i valori di errore totale

delle risposte emodinamiche stimate rispetto a quelle vere, sia per le piene che per le ridotte,

considerando solo quelle stimate nei canali dell’emisfero destro, e calcolate distintamente su

ciascun volontario.

In generale, dalle Tabelle 5.1, 5.2, 5.3, 5.4, 5.5 e 5.6, si osserva che le

percentuali di errore date dal BP sono sempre maggiori di quelle offerte dal

CS-TFT. La differenza minima si osserva nell’errore in latenza, in quanto i

risultati sono simili, mentre più marcato risulta lo scarto nella stima del picco.

Dopo questo primo test, si è passati a testare l’algoritmo in una situazione più

vicina a quella reale, implementandolo in uno scenario semi-simulato.

5.2 SCENARIO SEMI-SIMULATO

In tale scenario, il canale standard è stato realizzato sommando delle HRF

simulate a dati reali di resting state acquisiti su alcuni volontari. Rispetto al caso

completamente simulato realizzato nella prima fase dunque, il rumore

fisiologico che in precedenza era stato creato in simulazione mediante una

sommatoria di 5 sinusoidi a frequenze diverse, ora viene invece realmente

acquisito.

Come enunciato nel paragrafo 3.1.2, i segnali utilizzati per testare l’algoritmo

sono in particolare 14, provenienti da 7 soggetti sani e acquisiti in sessioni

distinte di resting state.

RISULTATI

72

5.2.1 DETTAGLI IMPLEMENTATIVI

Lo schema di funzionamento dell’algoritmo qui implementato, non di molto si

discosta rispetto a quello adottato nella prima fase, eccezion fatta per le

partizioni di segnale che qui non contengono alcun zero-padding, in quanto le

HRF sono state inserite ad intervalli regolari. Si è in tal modo provveduto a

segmentare tutto il segnale acquisito in parti uguali facendo corrispondere

l’inizio di ogni segmento con l’inizio di ciascun trial. L’obiettivo rimane il

medesimo, ovvero quello di stimare il rumore fisiologico a partire dal canale di

riferimento (SS), sottrarlo al segnale del canale standard (LS) per ridurre il

rumore fisiologico presente ed infine stimare la HRF media. Con lo scopo di

analizzare l’algoritmo in una visione più ampia e dettagliata, si è deciso di

estendere lo studio confrontando le prestazioni fin qui raggiunte con quelle

ottenute mediante algoritmi già implementati in passato, quali ad esempio la

media convenzionale realizzata a seguito del filtraggio passabanda (BP) e la

stima ReMCoBA, descritte entrambe nel capitolo 2.

5.2.2 STIME OTTENUTE E CONFRONTO CON FILTRAGGIO BP E

METODO REMCOBA

Nei grafici sottostanti sono raffigurate, separatamente per emisfero destro e

sinistro, le risposte emodinamiche medie stimate sia con l’algoritmo CS-TFT,

che con i metodi usati per il confronto, ovvero ReMCoBA e filtraggio

passabanda (BP W/O SS). Sono riportate inoltre in verde le medie delle HRF

di ciascuno dei 14 soggetti considerati stimate utilizzando il campionamento

compresso di Taylor Fourier, per dare un’idea della variabilità inter-individuale

nelle stime ottenute con questo algoritmo.

RISULTATI

73

A)

B)

Figura 5.10: Emisfero sinistro, HRF a confronto: A) Media delle HRF stimate mediante algoritmo CS-TFT (rossa), media delle HRF vere

simulate (nera), HRF medie di ogni soggetto ottenute mediante algoritmo CS-TFT (verdi).

B) HRF stimate mediante algoritmo CS-TFT (rossa), filtraggio passabanda BP (blu), ReMCoBA (magenta) confrontata con la media delle HRF vere (nero).

Osservando la Figura 5.10, appare chiara fin da subito la forte differenza nella

stima in ampiezza del picco della risposta emodinamica. Nessuno dei metodi

qui a confronto riesce a fittare il picco della HRF vera, che si discosta

notevolmente con uno scarto di circa 50 nM per il CS-TFT e addirittura di 100

nM per l’algoritmo ReMCoBA. L’algoritmo CS-TFT e il filtraggio passabanda

sembrano condurre a risultati simili per quanto riguarda la stima del picco, ma

diversi nella fase discendente dell’onda, in quanto il filtraggio comporta una

deriva molto maggiore (circa 50 nM) sia rispetto all’ HRF vera che a

ReMCoBA e CS-TFT. La latenza del picco sembra tuttavia essere abbastanza

buona.

RISULTATI

74

A)

B)

Figura 5.11: Emisfero destro, HRF a confronto: A) Media delle HRF stimate mediante algoritmo CS-TFT (rossa), media delle HRF vere

simulate (nera), HRF medie di ogni soggetto ottenute mediante algoritmo CS-TFT (verdi).

B) HRF stimate mediante algoritmo CS-TFT (rossa), filtraggio passabanda BP (blu), ReMCoBA (magenta) confrontata con la media delle HRF vere (nero).

Sul lato opposto, ovvero nell’emisfero destro (Figura 5.11), la situazione

appare alquanto diversa, in quanto il fit tra HRF vera e stimata mediante CS-

TFT appare ottimo per quanto riguarda la stima del picco. Tuttavia, la HRF

stimata non riesce a tornare a baseline, rimanendo ad un valore di circa 15 – 20

nM, per entrambi i metodi implementati. Anche nella fase ascendente vi è uno

scostamento, dove entrambe le HRF ricostruite non fittano i dati reali

risultando “in anticipo” rispetto all’onda simulata di circa 1 secondo.

L’algoritmo ReMCoBA tuttavia ha una peggiore stima in ampiezza del picco,

dal momento che, come nel caso dell’emisfero sinistro, presenta un andamento

più smooth: il picco medio stimato è infatti inferiore rispetto al vero di circa 20

nM. Poiché l’ampiezza delle risposte inserite nell’emisfero destro risulta circa

un terzo di quella dell’emisfero sinistro, il rapporto segnale-rumore nei canali

dell’emisfero destro è molto più basso. Questo è il motivo principale

dell’aumento di variabilità inter-individuale tra le stime medie dei diversi

soggetti, che risultano anche divergere maggiormente dalla HRF vera.

Oltre alla stima dell’errore calcolata separatamente per gli emisferi destro e

sinistro come nello scenario completamente simulato, nei canali dove non

RISULTATI

75

erano state aggiunte HRF simulate (canali posizionati frontalmente), si è

calcolato l’RMSE (cioè la radice quadrata dell’errore quadratico medio) come

metrica quantitativa, per confrontare la capacità degli algoritmi di stimare

l’assenza di HRF e, di conseguenza, la loro abilità nell’evitare di stimare bias.

Per ognuna delle metriche calcolate, una serie di t-test a due code a campioni

appaiati sono stati calcolati per confrontare le performance degli algoritmi

anche da un punto di vista statistico.

In Tabella 5.7 sono riportati, per l’emisfero sinistro, i valori medi dell’errore in

ampiezza di ogni soggetto e in Figura X gli scatter plot (dove il CS-TFT è

confrontato sia con il filtraggio passa-banda [BP W/O SS] che con l’algoritmo

ReMCoBA).

E_amp % (LEFT side)

CS-TFT BP

[W/O SS] ReMCoBA

Subj9 rs1 33,25 31,65 49,29

Subj9 rs2 21,18 18,98 38,54

Subj12 rs1 27,42 24,31 47,14

Subj12 rs2 26,85 41,33 29,19

Subj13 rs1 28,22 14,88 37,32

Subj13 rs2 17,27 20,03 34,37

Subj14 rs1 35,78 25,84 50,46

Subj14 rs2 29,63 54,88 53,91

Subj14 rs3 28,70 16,07 36,58

Subj15 rs1 37,80 43,71 33,09

Subj16 rs1 13,13 24,04 29,65

Subj16 rs2 13,23 24,17 14,47

Subj17 rs1 18,57 28,83 21,50

Subj17 rs2 17,84 16,78 30,34

MEAN 24,92 27,54 36,13

STD 17,09 17,26 14,62

Tabella 5.7: Tabella errore in ampiezza emisfero sinistro. Tabella contenente i valori di errore in ampiezza delle risposte emodinamiche stimate rispetto a quelle vere, considerando solo quelle stimate nei canali dell’emisfero sinistro, e calcolate distintamente su ciascun volontario.

RISULTATI

76

T-TEST E_amp (LEFT)

CS-TFT vs. BP

CS-TFT vs.

ReMCoBA

BP vs.

ReMCoBA

P-VALUE 0,1960 0,0001 0,0193

Tabella 5.8: Risultati dei t-test, calcolati confrontando ciascuno dei 3 algoritmi con gli altri

due. Per ciascun confronto, sono riportati i valori del p. Un p < 0.05 indica differenze

statisticamente significative tra gli algoritmi confrontati.

E’ evidente come la stima media in ampiezza dell’algoritmo ReMCoBA risulta

peggiore rispetto ad entrambi gli altri metodi (36,13 % contro il 27,53 % del

filtraggio BP W/O SS e i 24,92 % del CS-TFT). Il motivo potrebbe essere

riconducibile ad un’eccessiva riduzione del rumore fisiologico da parte di

questo algoritmo, che inevitabilmente comporta una rimozione di parte del

segnale utile. Un riscontro in questo verso si ha dalla miglior performance del

metodo BP W/O SS rispetto a ReMCoBA, che, al contrario di ReMCoBA,

non prevede nessuna sottrazione del rumore fisiologico proveniente dal canale

SS. La miglior performance del metodo CS-TFT potrebbe indicare una miglior

capacità ed accuratezza di questo algoritmo, rispetto a ReMCoBA, nello

stimare il rumore fisiologico dal canale SS.

A)

B)

Figura 5.12: A) Scatter plot dell’errore in ampiezza per confrontare la performance dell’algoritmo CS-

TFT (asse delle x) e di quello di filtraggio passabanda (asse delle y) per l’emisfero sinistro. Ogni asterisco rappresenta un soggetto di Tabella 5.7.

B) Scatter plot dell’errore in ampiezza per confrontare la performance dell’algoritmo CS-TFT (asse delle x) e di quello di Bayes (asse delle y) per l’emisfero sinistro. Ogni cerchio rappresenta un soggetto di Tabella 5.7.

RISULTATI

77

Lo scatter plot conferma quanto appena descritto: nel confronto tra BP W/O

SS e CS-TFT i risultati sono disposti equamente ai lati della diagonale (7 valori

contro 7), mentre lo stesso non si può dire nel secondo grafico, dove l’errore

in stima dell’ampiezza da parte di ReMCoBA è maggiore in 13 soggetti su 14

rispetto all’algoritmo CS-TFT.

Le Tabelle 5.9 e 5.10 riportano gli analoghi risultati ottenuti per l’emisfero

destro.

E_amp %

(RIGHT side) CS-TFT

BP

[W/O SS] ReMCoBA

Subj9 rs1 80,31 86,55 73,76

Subj9 rs2 80,15 76,00 59,76

Subj12 rs1 60,29 102,36 54,63

Subj12 rs2 71,75 261,95 58,74

Subj13 rs1 40,27 145,92 29,24

Subj13 rs2 47,15 72,98 35,30

Subj14 rs1 27,57 255,51 24,90

Subj14 rs2 48,12 71,81 44,71

Subj14 rs3 70,82 152,22 38,73

Subj15 rs1 95,78 87,97 66,81

Subj16 rs1 65,91 43,22 48,55

Subj16 rs2 76,00 42,91 74,74

Subj17 rs1 88,16 196,89 89,82

Subj17 rs2 65,87 127,93 44,11

MEAN 65,58 123,16 53,13

STD 56,42 67,50 42,03

Tabella 5.9: Tabella errore in ampiezza emisfero destro. Tabella contenente i valori di errore in ampiezza delle risposte emodinamiche stimate rispetto a quelle vere, considerando solo quelle stimate nei canali dell’emisfero destro, e calcolate distintamente su ciascun volontario.

RISULTATI

78

T-TEST E_amp (RIGHT)

CS-TFT vs. BP

CS-TFT vs.

ReMCoBA

BP vs.

ReMCoBA

P-VALUE 0,0083 0,0003 0,0024




Sul lato destro, a più basso SNR, notiamo una situazione inversa rispetto al

sinistro per quanto riguarda il confronto tra CS-TFT e ReMCoBA, con il

primo che commette una maggiore percentuale di errore di stima del picco

rispetto al secondo. Confrontando la media dell’errore col grafico di Figura

5.11 B), tale risultato appare inveritiero, in quanto il grafico denota un fitting

ottimo del picco HRF da parte del CS-TFT, decisamente migliore rispetto a

ReMCoBA.

Tuttavia, va detto che il confronto di tale grafico con i risultati in Tabella 5.9

sarebbe improprio, in quanto basato su misurazioni di diverso tipo: dal

momento che l’errore in ampiezza è infatti calcolato in modulo, non si

percepisce la tipologia della stima (se ad esempio vi è una sotto o sovrastima),

cosa di cui invece l’andamento dell’HRF nel grafico tiene in considerazione.

Osservando il grafico sembrerebbe pertanto che il CS-TFT ha una maggior

accuratezza del ReMCoBA rispetto alla stima sul picco, ma i valori della tabella

denotano che quest’ultimo ha una maggior precisione.

Il metodo BP W/O SS ha invece una performance nettamente inferiore agli

altri due metodi, mettendo in evidenza come la sottrazione del rumore

fisiologico sia di particolare importanza soprattutto nelle situazioni di dati con

basso SNR.

RISULTATI

79

A)

B)

Figura 5.13: A) Scatter plot dell’errore in ampiezza per confrontare la performance dell’algoritmo CS-

TFT (asse delle x) e di quello di filtraggio passabanda (asse delle y) per l’emisfero destro. Ogni asterisco rappresenta un soggetto di Tabella 5.9.

B) Scatter plot dell’errore in ampiezza per confrontare la performance dell’algoritmo CS-TFT (asse delle x) e di quello di Bayes (asse delle y) per l’emisfero destro. Ogni cerchio rappresenta un soggetto di Tabella 5.9.

Le Tabelle 5.11 e 5.13 riportano gli errori percentuali per latenza,

rispettivamente nell’emisfero sinistro e destro. I valori dei p-value sono

riportati, per gli stessi due casi, in Tabella 5.12 e 5.14.

RISULTATI

80

E_lat %

(LEFT side) CS-TFT

BP

[W/O SS] ReMCoBA

Subj9 rs1 6,20 5,88 6,19

Subj9 rs2 7,18 5,62 9,67

Subj12 rs1 45,73 29,28 48,20

Subj12 rs2 13,82 29,04 27,08

Subj13 rs1 11,88 13,60 14,11

Subj13 rs2 7,96 13,38 15,26

Subj14 rs1 9,14 12,08 6,90

Subj14 rs2 6,20 7,67 10,35

Subj14 rs3 4,24 19,44 10,17

Subj15 rs1 10,07 28,07 18,92

Subj16 rs1 4,44 7,57 10,77

Subj16 rs2 8,81 5,68 9,10

Subj17 rs1 7,80 10,29 8,38

Subj17 rs2 11,02 12,65 11,91

MEAN 11,04 14,30 14,79

STD 11,97 14,51 16,70

Tabella 5.11: Tabella errore in latenza emisfero sinistro. Tabella contenente i valori di errore in latenza delle risposte emodinamiche stimate rispetto a quelle vere, considerando solo quelle stimate nei canali dell’emisfero sinistro, e calcolate distintamente su ciascun volontario.

T-TEST E_lat (LEFT)

CS-TFT vs. BP

CS-TFT vs.

ReMCoBA

BP vs

ReMCoBA

P-VALUE 0,0915 0,0025 0,3976

Tabella 5.12: Risultati dei t-test, calcolati confrontando ciascuno dei 3 algoritmi con gli altri due. Per ciascun confronto, sono riportati i valori del p. Un p < 0.05 indica differenze statisticamente significative tra gli algoritmi confrontati.

RISULTATI

81

Emerge dai risultati che in latenza l’algoritmo da noi proposto si comporta

bene, fornendo risultati migliori, anche se non con ampio margine di

differenza, sia rispetto all’approccio ReMCoBA che al semplice filtraggio BP

W/O SS (11,04% del CS-TFT a fronte di 14,3% del BP W/O SS e 14,79% del

ReMCoBA). Lo stesso risultato è visibile dallo scatter plot, dove in entrambi i

grafici i valori percentuali dell’errore propendono in tutti e due i casi a favore

del CS-TFT (4-10 per il CS-BP e 1-12 per il CS-ReMCoBA).

A)

B)

Figura 5.14:

A). Scatter plot dell’errore in latenza per confrontare la performance dell’algoritmo CS-TFT

(asse delle x) e di quello di filtraggio passabanda (asse delle y) per l’emisfero sinistro. Ogni

asterisco rappresenta un soggetto di Tabella 5.11.

B) Scatter plot dell’errore in latenza per confrontare la performance dell’algoritmo CS-TFT

(asse delle x) e di quello di Bayes (asse delle y) per l’emisfero sinistro. Ogni cerchio

rappresenta un soggetto di Tabella 5.11.

RISULTATI

82

E_lat %

(RIGHT side) CS-TFT

BP

[W/O SS] ReMCoBA

Subj9 rs1 46,22 57,88 59,20

Subj9 rs2 33,15 38,49 33,92

Subj12 rs1 48,91 42,75 70,18

Subj12 rs2 84,49 124,79 122,02

Subj13 rs1 36,15 47,50 21,24

Subj13 rs2 85,16 103,86 88,80

Subj14 rs1 69,29 179,46 91,45

Subj14 rs2 38,73 138,91 103,65

Subj14 rs3 55,14 153,65 92,02

Subj15 rs1 43,19 81,31 49,35

Subj16 rs1 21,88 76,01 26,03

Subj16 rs2 37,15 23,04 18,73

Subj17 rs1 25,87 12,52 28,66

Subj17 rs2 32,09 56,86 63,30

MEAN 46,96 81,22 62,04

STD 50,61 47,33 52,61

Tabella 5.13: Tabella errore in latenza emisfero destro. Tabella contenente i valori di errore in latenza delle risposte emodinamiche stimate rispetto a quelle vere, considerando solo quelle stimate nei canali dell’emisfero destro, e calcolate distintamente su ciascun volontario.

T-TEST E_lat (RIGHT)

CS-TFT vs. BP

CS-TFT vs.

ReMCoBA

BP vs.

ReMCoBA

P-VALUE 0,0047 0,0128 0,0210


RISULTATI

83

Sul lato destro a più basso SNR le differenze tra i vari metodi implementati

sono maggiori e i valori ottenuti ben distanti tra loro: il vantaggio, ancora una

volta, è osservabile particolarmente per l’algoritmo qui proposto (46,96% di

errore di stima della latenza), che offre una percentuale di errore minore di

circa il 15% rispetto all’approccio ReMCoBA, e del 34% rispetto al BP W/O

SS. Nello scatter plot emerge che in soli 3 casi su 14 i valori di errore in latenza

ottenuti da ReMCoBA e BP W/O SS sono stati inferiori a quelli dati dal CS-

TFT.

A)

B)

Figura 5.15: A) Scatter plot dell’errore in latenza per confrontare per l’emisfero destro la performance

dell’algoritmo CS-TFT (asse delle x) e di quello di filtraggio passabanda (asse delle y) per l’emisfero destro. Ogni asterisco rappresenta un soggetto di Tabella 5.13.

B) Scatter plot dell’errore in latenza per confrontare per l’emisfero destro la performance dell’algoritmo CS-TFT (asse delle x) e di quello di Bayes (asse delle y) per l’emisfero destro. Ogni cerchio rappresenta un soggetto di Tabella 5.13.

Infine, come per lo scenario completamente simulato, si è calcolata anche la

percentuale di errore. Le Tabelle 5.15 e 5.17 riportano gli errori percentuali

calcolati sulla stima dell’intera curva, rispettivamente nell’emisfero sinistro e

destro. I valori dei p-value sono riportati, per gli stessi due casi, in Tabella 5.16

e 5.18.

RISULTATI

84

E_tot %

(LEFT side) CS-TFT

BP

[W/O SS] ReMCoBA

Subj9 rs1 27,87 31,98 35,32

Subj9 rs2 11,32 11,30 18,77

Subj12 rs1 93,08 104,67 131,41

Subj12 rs2 81,11 177,95 49,20

Subj13 rs1 19,60 38,20 20,36

Subj13 rs2 14,45 35,54 20,26

Subj14 rs1 34,27 129,14 46,27

Subj14 rs2 24,46 85,15 40,58

Subj14 rs3 17,02 63,19 22,79

Subj15 rs1 93,95 96,40 63,64

Subj16 rs1 11,11 23,33 15,88

Subj16 rs2 16,87 23,99 17,96

Subj17 rs1 19,97 39,22 17,58

Subj17 rs2 22,85 50,07 23,23

MEAN 34,85 65,01 37,38

STD 36,84 45,33 28,87

Tabella 5.15: Tabella errore totale emisfero sinistro. Tabella contenente i valori di errore

totale delle risposte emodinamiche stimate rispetto a quelle vere, considerando solo quelle

stimate nei canali dell’emisfero sinistro, e calcolate distintamente su ciascun volontario.

T-TEST E_tot (LEFT)

CS-TFT vs. BP

CS-TFT vs.

ReMCoBA

BP vs.

ReMCoBA

P-VALUE 0,0021 0,2980 0,0104




RISULTATI

85

Dalla Tabella 5.15 emerge come il metodo BP W/O SS commetta l’errore

maggiore (quasi doppio rispetto agli altri metodi, 65% contro i 34% e 37% di

CS-TFT e ReMCoBA), sottolineando ulteriormente l’importanza di utilizzare il

segnale proveniente dal canale SS per ridurre il rumore fisiologico e stimare la

HRF. L’algoritmo CS-TFT commette in media un errore leggermente inferiore

rispetto a ReMCoBA, anche se non ci sono differenze statisticamente

significative tra i risultati dei due algoritmi. Lo scatter plot dei risultati di questi

due algoritmi infatti, denota che i valori sono abbastanza bilanciati attorno alla

diagonale, anche se il loro numero propende a favore del CS-TFT (2 contro 8,

con i restanti disposti intorno alla diagonale).

A)

B)

Figura 5.16: A) Scatter plot dell’errore totale per confrontare la performance dell’algoritmo CS-TFT (asse

delle x) e di quello di filtraggio passabanda (asse delle y) per l’emisfero sinistro. Ogni asterisco rappresenta un soggetto di Tabella 5.15.

B) Scatter plot dell’errore totale per confrontare la performance dell’algoritmo CS-TFT (asse delle x) e di quello di Bayes (asse delle y) per l’emisfero sinistro. Ogni cerchio rappresenta un soggetto di Tabella 5.15.

RISULTATI

86

E_tot % (RIGHT side)

CS-TFT BP

[W/O SS] ReMCoBA

Subj9 rs1 278,23 260,30 252,53

Subj9 rs2 864,14 822,30 513,62

Subj12 rs1 317,94 433,11 265,74

Subj12 rs2 443,18 2654,46 421,13

Subj13 rs1 86,35 647,12 63,75

Subj13 rs2 192,60 363,23 181,30

Subj14 rs1 181,89 1444,77 205,60

Subj14 rs2 314,62 621,07 284,99

Subj14 rs3 362,88 1010,02 183,08

Subj15 rs1 448,53 655,96 335,73

Subj16 rs1 224,31 274,90 204,34

Subj16 rs2 268,68 224,14 245,85

Subj17 rs1 450,04 1022,67 404,89

Subj17 rs2 325,40 839,40 174,59

MEAN 339,91 805,25 266,94

STD 461,75 593,74 297,75

Tabella 5.17: Tabella errore totale emisfero destro. Tabella contenente i valori di errore

totale delle risposte emodinamiche stimate rispetto a quelle vere, considerando solo quelle

stimate nei canali dell’emisfero destro, e calcolate distintamente su ciascun volontario.

T-TEST E_tot (RIGHT)

CS-TFT vs. BP

CS-TFT vs.

ReMCoBA

BP vs.

ReMCoBA

P-VALUE 0,0073 0,0077 0,0026


due. Per ciascun confronto, sono riportati i valori del p. Un p < 0.05 indica

differenze statisticamente significative tra gli algoritmi confrontati.

RISULTATI

87

Come per le altre tipologie di errore in precedenza considerate, anche in

questo caso il lato destro, a più basso SNR, presenta valori di errore medio

totale molto maggiori e più eterogenei. ReMCoBA questa volta appare in netto

“vantaggio” rispetto a CS-TFT, con uno scarto in percentuale del 73%.

Nettamente superiore invece l’errore percentuale medio commesso da BP

W/O SS, che supera di circa il triplo ReMCoBA (805,25% contro 266,94%) e

di più del doppio CS-TFT (805,25% contro 339,91%).

Nello scatter plot che raffigura CS-TFT e ReMCoBA a confronto, per la prima

volta rispetto ai grafici precedenti emerge che gli errori propendono tutti a

sfavore di CS-TFT (13 contro 1), anche se buona parte di questi, per la

precisione 9, sono situati a ridosso della diagonale.

A)

B)

Figura 5.17:

A) Scatter plot dell’errore totale per confrontare la performance dell’algoritmo CS-TFT (asse

delle x) e di quello di filtraggio passabanda (asse delle y) per l’emisfero destro. Ogni

asterisco rappresenta un soggetto di Tabella 5.17.

B) Scatter plot dell’errore totale per confrontare la performance dell’algoritmo CS-TFT (asse

delle x) e di quello di Bayes (asse delle y) per l’emisfero destro. Ogni cerchio rappresenta

un soggetto di Tabella 5.17.

La media totale e la relativa deviazione standard per le tre metriche di errore

considerate, per ciascun algoritmo implementato, e per i due emisferi, è

raffigurata nei grafici a barre sottostanti. Si nota sin da una prima osservazione

RISULTATI

88

che le percentuali di errore si differenziano molto a seconda dell’emisfero

considerato: con HRF di ampiezza più alta (emisfero sinistro, circa 250 nM) e

quindi dati con SNR più alto, le prestazioni migliori sono ottenute

dall’algoritmo qui descritto, ovvero CS-TFT, che permette di ottenere una

stima migliore delle HRF sia in ampiezza che in latenza. Come dai grafici e

dalle tabelle precedenti, emerge come la stima dell’ampiezza e della latenza del

picco di ReMCoBA sia in media peggiore del più semplice BP W/O SS. In

quanto ad errore totale invece, CS-TFT e ReMCoBA sembrano avere circa gli

stessi risultati, con una miglior stima ancora una volta leggermente a favore di

CS-TFT e in entrambi i casi nettamente migliore di BP W/O SS.

A)

B)

Figura 5.18: A) Emisfero sinistro, grafico della media dell’errore con deviazione standard: percentuale media di errore in ampiezza (sinistra), latenza (centro) e totale (destra) per ciascun algoritmo: CS-TFT (blu), BP W/O SS (azzurro), e ReMCoBA (giallo). B) Emisfero destro, grafico della media dell’errore con deviazione standard: percentuale media di errore in ampiezza (sinistra), latenza (centro) e totale (destra) per ciascun algoritmo: CS-TFT (blu), BP W/O SS (azzurro), e ReMCoBA (giallo).

Osservando l’emisfero destro tuttavia, a più basso SNR, ci si accorge

immediatamente che ReMCoBA e CS-TFT offrono prestazioni molto più

simili tra loro rispetto a quanto notato sul lato sinistro, mentre il metodo BP

W/O SS si discosta di molto, e sembra non “competere” con le prestazioni

offerte dai due metodi precedentemente menzionati. Considerando infatti i

soli ReMCoBA e CS-TFT, le differenze in percentuali di errore possono dirsi

trascurabili, specie se paragonate ai valori ottenuti da BP W/O SS: una

prestazione migliore di ReMCoBA rispetto a CS-TFT per quanto riguarda

RISULTATI

89

l’ampiezza (53% di ReMCoBA contro il 65% di CS-TFT), mentre si nota

un’inversione di posizioni nella stima dell’errore in latenza (46% contro 62%).

Nell’errore totale, invece, ReMCoBA commette un errore inferiore (266%

contro 340%).

5.3 DISCUSSIONE

A conclusione delle analisi effettuate distintamente su entrambi gli emisferi,

potremmo pertanto dire che l’algoritmo di CS-TFT e ReMCoBA offrono

entrambi interessanti soluzioni per la riduzione del rumore fisiologico e la

stima della HRF: entrambi hanno una resa decisamente superiore al normale

filtraggio passabanda BP W/O SS. In particolare, il metodo del

campionamento compresso sembra avere i migliori risultati in assoluto quando

il segnale contiene delle risposte emodinamiche di ampiezza più elevata (~250

nM), mentre lo stesso non si può dire se queste hanno invece il picco intorno

nell’intorno di 75nM, ovvero con dati che presentano un più basso SNR.

Un’ipotesi potrebbe essere dunque quella di utilizzare il metodo CS-TFT nella

stima del rumore fisiologico, e di applicare in cascata un filtraggio bayesiano

per migliorare la stima della HRF, come in ReMCoBA.

5.4 RMSE

Per i canali situati in zona frontale, a cui non è stata aggiunta alcuna risposta

emodinamica simulata, si è deciso, come termine di confronto quantitativo tra

algoritmi, di calcolare la radice quadrata dell’errore quadratico medio (o Root

Mean Square Error, RMSE), ovvero:

√∑ (

)

(5.6)

RISULTATI

90

dove N è il numero di campioni contenuti in un trial.

A)

B)

Figura 5.19: A) Parte frontale, resting state (RS) confronto: RS stimato mediante algoritmo CS-TFT (rosso), vero (nero), e resting state di ogni soggetto ottenuto mediante algoritmo CS-TFT (verdi). B) Parte frontale, media dei resting state (RS) a confronto: RS stimato mediante algoritmo CS-TFT (rosso), filtraggio passabanda BP (blu), ReMCoBA (magenta) confrontate con la media del resting state vero (nero).

In tal caso, il segnale risulta completamente nullo, mentre la stima

dovrebbe essere il più vicina possibile allo zero. Tuttavia, come si può notare

in Figura 5.19, essa presenta delle oscillazioni causate dal rumore fisiologico di

fondo che persistono nonostante l’uso del segnale proveniente dal canale SS. Il

rumore fisiologico residuo oscilla intorno alla linea della 0: questo risultato è

molto promettente, in quanto evidenzia come, in una situazione reale, la

possibilità di stimare una HRF dove questa non è presente, sia molto ridotta. I

risultati ottenuti sono visibili nella tabella sottostante:

RISULTATI

91

RMSE

FRONTAL CS-TFT

BP

[W/O SS] ReMCoBA

Subj9 rs1 40,23 31,92 26,57

Subj9 rs2 18,98 26,64 15,83

Subj12 rs1 52,54 45,08 45,49

Subj12 rs2 46,36 81,20 42,34

Subj13 rs1 37,21 49,82 29,45

Subj13 rs2 31,17 47,11 35,00

Subj14 rs1 41,68 100,49 43,49

Subj14 rs2 29,16 65,52 30,36

Subj14 rs3 42,76 71,10 44,84

Subj15 rs1 32,19 51,81 33,23

Subj16 rs1 26,31 36,46 24,84

Subj16 rs2 40,74 29,92 29,81

Subj17 rs1 45,76 71,07 49,33

Subj17 rs2 40,75 44,67 34,93

MEAN 37,56 53,77 34,68

STD 9,00 21,44 9,45

Tabella 5.19: Tabella RMSE. Tabella contenente i valori di RMSE calcolate distintamente su ciascun volontario sui segnali di resting state dei canali frontali.

T-TEST RMSE

CS-TFT vs. BP

CS-TFT vs.

ReMCoBA

BP vs.

ReMCoBA

P-VALUE 0,0042 0,0368 0,0003


RISULTATI

92

Come emerge dai valori della tabella, i risultati di ReMCoBA e CS-TFT

sembrano simili. Tuttavia, i T-test mostrano una differenza significativa tra i

due algoritmi, con una miglior performance di ReMCoBA rispetto al CS-TFT.

È probabile che il filtro bayesiano applicato in ReMCoBA riesca a ridurre

ulteriormente il rumore fisiologico rimasto e quindi a prevalere nei risultati sul

CS-TFT. La soluzione proposta prima, di applicare il filtro bayesiano in

cascata al CS-TFT, rimane valida anche qui. Il solo filtraggio BP W/O SS non

porta a risultati altrettanto soddisfacenti, sottolineando ulteriormente

l’importanza di ridurre il rumore fisiologico con l’uso del canale SS.

93

CAPITOLO 6

CONCLUSIONI E SVILUPPI

FUTURI

6.1 SVILUPPI FUTURI

Come accennato nei capitoli precedenti, l’algoritmo Compressive Sensing

mediante trasformata di Taylor-Fourier presentato in questa tesi è ancora agli

albori per quanto riguarda il suo utilizzo nella riduzione del rumore fisiologico

e stima della risposta emodinamica in segnali acquisiti con la NIRS. Essendo

questo uno dei primi tentativi di sperimentazione, e data la poca storicità del

metodo, è indubbio che i primi risultati forniti siano da considerarsi provvisori

e atti ad essere assolutamente oggetto di uno studio più approfondito che

accerti, a seguito di molteplici prove e sperimentazioni, l’effettiva validità della

nuova soluzione proposta.

Osservando il codice implementato, molti sono ancora i punti su cui si può

lavorare cercando eventuali soluzioni migliori. L’algoritmo è stato applicato sia

nello scenario totalmente simulato che in quello semi-simulato senza alcuna

modifica: questo, da un lato mostra la stabilità e generalità dell’algoritmo,

dall’altro il fatto che potenzialmente è migliorabile, in quanto i parametri erano

stati settati basandosi sullo scenario totalmente simulato. Inoltre, va aggiunto

che non sono stati eseguiti test particolarmente intensivi sul codice, variando

CONCLUSIONI E SVILUPPI FUTURI

94

ad esempio numerose volte i parametri fondamentali. Non sono stati trattati

segnali particolari (ad esempio segnali ricchi di artefatti da movimento, dove

l’algoritmo potrebbe andare in crisi e offrire risultati completamente

discordanti rispetto agli attuali), e non si è cercato, per ora, di migliorare il

costo computazionale dell’algoritmo. In questa tesi l’algoritmo è stato testato

su un numero ridotto di soggetti (7 soggetti per un totale di 14 dataset): per

una migliore attendibilità ed accuratezza della performance dell’algoritmo è

auspicabile testare l’algoritmo su un numero maggiore di soggetti. Inoltre, la

validazione dell’algoritmo su dati totalmente reali è essenziale per una sua

futura applicabilità in contesti reali.

Si potrebbe pertanto ridurre la dimensione di D al solo range di frequenze di

interesse, riducendo in tal modo il tempo di calcolo di circa il 90%.

Il costo computazionale totale risulta attualmente inaccettabile qualora si

volesse far lavorare l’algoritmo in real-time, in quanto occorrono circa 20

minuti per ottenere i risultati di un singolo soggetto con 36 canali (sempre

considerando un computer di medie prestazioni). E’ ovvio pertanto che

qualora l’obiettivo fosse raggiungere una rielaborazione del segnale in tempo

reale, ulteriori sono ancora i progressi da realizzare: uno di questi, in prima

ipotesi, potrebbe essere quello di implementare il parallelismo del codice, oltre

indubbiamente a semplificare le parti di programma costituite dai comandi più

onerosi, come accennato poc’anzi.

Un interessante test per l’algoritmo presentato potrebbe essere dedicato alla

stima della media data da un numero di HRF sempre minore, sia per

osservarne il cambiamento della stima che per stabilire magari un numero

minimo di risposte sotto il quale la media risulta completamente distorta e

insignificante. Confrontare le performance dell’algoritmo proposto e di quelli

di confronto potrebbe fornire importanti indicazioni sull’abilità dell’algoritmo

proposto di convergere ad una stima stabile con un numero minore di trial

disponibili, essenziale, per esempio, in disegni sperimentali che coinvolgono

bambini o particolari popolazioni cliniche, che devono per necessità essere di

breve durata.

Un’interessante ed ulteriore prova potrebbe essere un’analisi di sensibilità, che

faccia variare tali parametri, uno o più alla volta, su un vasto campione di dati

CONCLUSIONI E SVILUPPI FUTURI

95

simulati, con l’obiettivo di raggiungere la miglior stima dell’errore sia in

ampiezza che in latenza. In tal modo si potrebbe riuscire anche a determinare,

ad esempio, se tali parametri dipendono molto dai dati considerati (e quindi

devono essere adattati e “personalizzati” di volta in volta su ciascun differente

segnale), o se invece possono essere mantenuti costanti in ogni segnale

analizzato. In particolare, tra i parametri considerati, l’ordine a cui viene

troncata la serie di Taylor, potrebbe essere quello più interessante da

modificare e allo stesso tempo forse il più determinante: non sappiamo

attualmente se facendo lavorare l’algoritmo con ordini maggiori si abbia una

ricostruzione ancora migliore del rumore fisiologico, dove anche le frequenze

molto più vicine tra loro vengano rilevate in maniera chiara e distinta, o se

invece si abbia una sorta di “saturazione” dove anche implementando ordini

nettamente maggiori, i miglioramenti visibili risultino alquanto trascurabili

rispetto agli attuali.

Finora i termini derivativi di ordine superiore sono stati usati solo per cercare

di riprodurre al meglio le variazioni nel tempo del rumore fisiologico, ma si

potrebbe iniziare a studiarle singolarmente. La derivata inoltre, potrebbe essere

utilizzata anche nelle valutazioni sull’ampiezza del segnale, al fine di

determinare se siamo in presenza di un particolare artefatto o se si tratta invece

di solo rumore fisiologico.

97

APPENDICE

% ********************************************** 1 % Algoritmo CS-TFT implementato su segnali fNIRS 2 % ********************************************** 3 % 4 % ------------- Marco Cescon ------------------- 5 % ------------- matr. 1039022 ------------------ 6 % --------- Anno accademico 2015/16 ------------ 7 % 8 % ********************************************** 9 % "Approccio basato su Compressive Sensing e 10 % modellistica Taylor-Fourier per stimare la 11 % risposta emodinamica da segnali fNIRS" 12 % ********************************************** 13 % 14 % Scenario semi-simulato 15 % 14 segnali, 52 canali per ciascun segnale 16 17 function stimatore_fNIRS(soggetto,rs) 18 close all 19 clc 20 21 if nargin ~= 2 22 display('Use: stimatore_fNIRS_human(soggetto,rs,trigger,Fs)') 23 return 24 end 25 26 27 % Conversione dei dati da densità ottiche a nanomoli 28 switch soggetto 29 case 9 30 if rs == 1 31 load('20120917_6009_RS1.nirs','-mat'); 32 elseif rs == 2 33 load('20120917_6009_RS2.nirs','-mat'); 34 else disp('ERROR!! invalid RS') 35 return 36 end 37 case 12 38

APPENDICE

98

if rs == 1 39 load('20120921_6012_RS1.nirs','-mat'); 40 elseif rs == 2 41 load('20120921_6012_RS2.nirs','-mat'); 42 else disp('ERROR!! invalid RS') 43 return 44 end 45 case 13 46 if rs == 1 47 load('20120925_6013_RS1.nirs','-mat'); 48 elseif rs == 2 49 load('20120925_6013_RS2.nirs','-mat'); 50 else disp('ERROR!! invalid RS') 51 return 52 end 53 case 14 54 if rs == 1 55 load('20120925_6014_RS1.nirs','-mat'); 56 elseif rs == 2 57 load('20120925_6014_RS2.nirs','-mat'); 58 elseif rs == 3 59 load('20120925_6014_RS3.nirs','-mat'); 60 else disp('ERROR!! invalid RS') 61 return 62 end 63 case 15 64 if rs == 1 65 load('20120925_6015_RS1.nirs','-mat'); 66 else disp('ERROR!! invalid RS') 67 return 68 end 69 case 16 70 if rs == 1 71 load('20120927_6016_RS1.nirs','-mat'); 72 elseif rs == 2 73 load('20120927_6016_RS2.nirs','-mat'); 74 else disp('ERROR!! invalid RS') 75 return 76 end 77 case 17 78 if rs == 1 79 load('20120927_6017_RS1.nirs','-mat'); 80 elseif rs == 2 81 load('20120927_6017_RS2.nirs','-mat'); 82 else disp('ERROR!! invalid RS') 83 return 84 end 85 otherwise 86 disp('ERROR!! Invalid subject') 87 return 88 end 89 90 coeff = 1e9; % Conversione in nM 91

APPENDICE

99

dod = hmrIntensity2OD(d); 92 conc = hmrOD2Conc(dod,SD,[6 6]); 93 concHbO = squeeze(conc(:,1,:)); 94 concHbO = coeff*concHbO; 95 96 nHRF = 28; 97 t_camp = t(2)-t(1); 98 % Valori inseriti di default 99 durata_hrf = 15; 100 distance = 40; 101 102 display('----------------------------------------------------------------') 103 pause(.5) 104 clc 105 106 disp(['Numero HRF: ' num2str(nHRF)]) 107 disp(['Durata HRF: ' num2str(durata_hrf)]) 108 109 %------------------------------------------------------------------------- 110 NumChannels = size(concHbO,2); 111 112 channel_LEFT = [1 2 3 4 5 6]; % Sorgenti emisfero SX 113 channel_RIGHT = [7 8 9 10 11 12]; % Sorgenti emisfero DX 114 channel_MIDDLE = [13 14 15]; % Sorgenti frontali 115 116 errorAmpRic_LEFT = []; 117 errorAmpRic_RIGHT = []; 118 errorLatRic_LEFT = []; 119 errorLatRic_RIGHT = []; 120 errorTotRic_LEFT = []; 121 errorTotRic_RIGHT = []; 122 123 errorAmpFiltr_LEFT = []; 124 errorAmpFiltr_RIGHT = []; 125 errorLatFiltr_LEFT = []; 126 errorLatFiltr_RIGHT = []; 127 errorTotFiltr_LEFT = []; 128 errorTotFiltr_RIGHT = []; 129 130 HRFvera_LEFT = []; 131 HRFric_LEFT = []; 132 HRFfiltr_LEFT = []; 133 HRFvera_RIGHT = []; 134 HRFric_RIGHT = []; 135 HRFfiltr_RIGHT = []; 136 137 rmse_ric = []; 138 rmse_filtr = []; 139 140 for i = 1:NumChannels 141 sorg = SD.MeasList(i,1); 142 detet = SD.MeasList(i,2); 143 144

APPENDICE

100

distX = abs(SD.SrcPos(sorg,1) - SD.DetPos(detet,1)); 145 distY = abs(SD.SrcPos(sorg,2) - SD.DetPos(detet,2)); 146 147 dist(i) = norm([distX distY]); 148 end 149 150 % I valori in uscita delle distanze sono solo 2, pertando è sufficiente 151 % fare la media tra il valore minimo e massimo: LS sono i valori sopra la 152 % media, SS invece quelli sotto 153 meanDist = mean([min(dist) max(dist)]); 154 155 SS = find(dist<meanDist); 156 LS = find(dist>meanDist); 157 158 numCh_LS = length(LS); 159 numCh_SS = length(SS); 160 161 clear meanDist valMinDiff indMinDiff diff_Ch 162 163 diff_Ch = zeros(numCh_SS,1); % Alloco la matrice delle differenze 164 abbinam_LS_SS = zeros(numCh_LS,2); % Alloco la matrice degli abbinamenti LS - 165 SS 166 167 for i = 1:numCh_LS 168 for j = 1:numCh_SS 169 diff_Ch(j) = sum(abs(concHbO(:,LS(i)) - concHbO(:,SS(j)) )); 170 end 171 [valMinDiff,indMinDiff] = min(diff_Ch); 172 abbinam_LS_SS(i,:) = [LS(i) SS(indMinDiff)]; % Inserisco gli abbinamenti 173 LS-SS 174 end 175 176 clear valMinDiff diff_Ch indMinDiff 177 display(' Start calculating...') 178 pause(.5) 179 180 181 % ************************************************************************* 182 % ************************************************************************* 183 % ************************************************************************* 184 for iter = 1:numCh_LS 185 cha_long = abbinam_LS_SS(iter,1); 186 cha_short = abbinam_LS_SS(iter,2); 187 src_long = SD.MeasList(cha_long,1); % Sorgente del canale lungo 188 189 clc 190 fprintf('\n\n') 191 display(['Iteration ' num2str(iter) ' of ' num2str(numCh_LS)]) 192 fprintf('\n\n') 193 display(['SELECTED SUBJECT: ' num2str(soggetto) ', RS: ' num2str(rs)]) 194 display(['SELECTED LONG CHANNEL: ' num2str(cha_long) ' (Src=' 195 num2str(SD.MeasList(cha_long,1)) ', Det=' num2str(SD.MeasList(cha_long,2)) 196 ')']) 197

APPENDICE

101

display(['SELECTED SHORT CHANNEL: ' num2str(cha_short) ' (Src=' 198 num2str(SD.MeasList(cha_short,1)) ', Det=' num2str(SD.MeasList(cha_short,2)) 199 ')']) 200 %--------------------------------------------------------------------------201 - 202 203 if ismember(src_long,channel_LEFT) 204 hr_amp = 1600; % Emisfero SX: ampiezza effettiva di 250 nmol 205 titleText = ['Channel: LEFT, HRF amp= ' num2str(hr_amp) ' nmol']; 206 disp(['Ampiezza HRF: ' num2str(hr_amp)]) 207 elseif ismember(src_long,channel_RIGHT) 208 hr_amp = 500; % Emisfero DX: ampiezza effettiva di 75 nmol 209 titleText = ['Channel: RIGHT, HRF amp= ' num2str(hr_amp) ' nmol']; 210 disp(['Ampiezza HRF: ' num2str(hr_amp)]) 211 else hr_amp = 0; % Fronte: non aggiungo alcuna HRF 212 titleText = ['Channel: MIDDLE, HRF amp= NO HRF inserted']; 213 disp(['Ampiezza HRF: no HRF']) 214 end 215 216 % Lo calcolo comunque per tutti i canali per mantenere poi il trigger 217 218 [raw.data,vett_hrf,raw.trigger,t_hrf,hrf_avg]=aggiungiHRF_marco(t,concHbO(:,cha219 _long),durata_hrf,nHRF,distance,hr_amp); 220 if ismember(src_long,channel_MIDDLE) 221 raw.data = concHbO(:,cha_long); % Annullo tutte le risposte evocate 222 della fronte 223 end 224 225 sig_long = raw.data; % segnale long separation 226 sig_short = concHbO(:,cha_short); 227 sig_trig = raw.trigger; % segnale di trigger corrispondente 228 % ************************************************************************ 229 230 %% Filtraggio passa-banda delle tracce ----------------------------------- 231 % % filtraggio passa-banda del segnale raw nell'intervallo 0.005 - 0.55 Hz 232 load LP_filter.mat % variabili: filtro - FIR_Hann_lp, ordine - No 233 load HP_filter.mat % variabili: filtro - FIR_Hann_hp, ordine - No 234 235 % segnale long separation 236 sig_lonf = filtfilt(FIR_Hann_lp.Numerator,1,sig_long); % filtraggio LP 237 sig_lonf = filtfilt(FIR_Hann_hp.Numerator,1,sig_lonf); % filtraggio HP 238 sig_lonf = sig_lonf(No+1:end); % tronco l'andata a regime del filtro 239 240 % segnale short separation 241 sig_shorf = filtfilt(FIR_Hann_lp.Numerator,1,sig_short); % filtraggio LP 242 sig_shorf = filtfilt(FIR_Hann_hp.Numerator,1,sig_shorf); % filtraggio HP 243 sig_shorf = sig_shorf(No+1:end); % tronco l'andata a regime del filtro, 244 elimino i primi 100 campioni 245 246 sig_trig = sig_trig(No+1:end); % segnale di trigger corrispondente 247 vett_hrf_bis = vett_hrf(No+1:end); % segnale con le HRF 248 % ************************************************************************ 249 250

APPENDICE

102

%% Estrazione del trigger di interesse ------------------------------------ 251 index_start = 1; % Per circa 30 secondi non ci sono impulsi HRF, lo 252 scaling lo calcoliamo qui 253 index_stop= min(find(sig_trig==1))-1; % Trovo l'indice di stop pari al 254 primo trigger meno 1 255 % ************************************************************************* 256 257 % Trovo gli indici di start di ciascuna porzione di 258 % segnale relativa al trigger desiderato ---------------------------------- 259 iStart = find(sig_trig==1); % Start = trigger 260 numPorz = length(iStart); % Numero di HRF (porzioni segnale) 261 262 HRF_start = -1; % La porzione di HRF la considero da -1 a 13 sec 263 HRF_stop = 13; % La porzione di HRF la considero da -1 a 13 sec 264 tHRF = HRF_start:t_camp:HRF_stop; % Considero un asse dei tempi da -1 a 13 265 ind_start = iStart-abs(HRF_start)/t_camp; 266 ind_stop = iStart+HRF_stop/t_camp; 267 durataPorz = ind_stop(1)-ind_start(1)+1; 268 269 % Matrice con tutte le risposte evocate vere e solo filtrate di un singolo 270 canale 271 totalHRF_vere = zeros(durataPorz,numPorz); 272 totalHRF_filtr = zeros(durataPorz,numPorz); 273 274 for i = 1:numPorz 275 totalHRF_vere(:,i) = vett_hrf_bis(ind_start(i):ind_stop(i)); 276 totalHRF_filtr(:,i) = sig_lonf(ind_start(i):ind_stop(i)); 277 end 278 279 if ismember(src_long,channel_MIDDLE) 280 % Annullo se non è stata aggiunta HRF al canale 281 mediaHRFvera_singleCh = zeros(durataPorz,1); 282 else 283 mediaHRFvera_singleCh = mean(totalHRF_vere,2); 284 end 285 clear iStart iStop durate_trig nCampAftPeak nCampBefPeak 286 287 % ************************************************************************* 288 % Scaling ----------------------------------------------------------------- 289 % Effettuato sui primi 30 secondi di segnale, dove forzatamente non ci sono 290 % impulsi 291 [xData, yData] = prepareCurveData( sig_shorf(index_start:index_stop), 292 sig_lonf(index_start:index_stop) ); 293 ft = fittype( 'poly1' ); % Tipologia di fitting polinomiale 294 295 % Trovo le costanti p1 e p2. 296 fitresult = fit(xData,yData, ft); 297 p1 = fitresult.p1; % P1: y = p1*x + p2 298 p2 = fitresult.p2; % P2: y = p1*x + p2 299 300 sig_shorf1 = p1*sig_shorf + p2; % Scaling del segnale short -------------- 301 clear index_start index_stop 302 % ************************************************************************ 303

APPENDICE

103

304 % Super risoluzione ------------------------------------------------------ 305 Ns = length(sig_shorf1); 306 Ts = t(2)-t(1); % Periodo di campionamento 307 Fs = 1/Ts; % Freq di campionamento 308 N = 501; 309 P = 11; % fattore interpolazione 310 H = N*P; % granularità griglia risultante 311 Qh = Fs/H; % nuovo quanto frequenziale 312 h_axis = 0:Qh:Qh*(H-1); % asso spettrale super-risoluto 313 tot_iter = ceil(Ns/N); 314 max_iter = 20; % massimo numero di componenti da stimare 315 316 %% Ciclo stima CSTFM 317 % stimo le frequenze sul segnale short separation 318 frq_cstfm = zeros(tot_iter,max_iter*2); % inizializzo vettore stime 319 sig_noise = zeros(N,tot_iter); 320 overlap = 5; % Serve per ricostruire bene il segnale anche nei bordi 321 simm = 1; % simmetria del supporto 322 323 %--------------------------------- 324 % Sono le due funzioni più onerose computazionalmente di tutto il 325 % programma: 326 D = creaMatriceD(N,H); % Circa 5 secondi 327 NormACols = sqrt(diag(D'*D)); % NormACols è un vettore colonna di lunghezza 328 H=5011, circa 7 secondi 329 %--------------------------------- 330 331 for k = 1:tot_iter 332 if k == 1 333 y_run = sig_shorf1(1:N); 334 else 335 iStart = N*(k-1)-overlap; 336 iStop = iStart+N-1; 337 if iStop > Ns 338 y_run = [sig_shorf1(iStart:end); zeros(iStop-Ns,1)]; 339 else 340 y_run = sig_shorf1(iStart:iStop); 341 end 342 end 343 344 v = fft(y_run); 345 suppo = RunOMP_IMEKO(v,D,max_iter,NormACols,simm); 346 347 if suppo(1) == 1 348 suppo(1) = 2; 349 end 350 351 frq_cstfm(k,1:length(suppo)) = h_axis(suppo); 352 353 % Definizione template 354 % ricostruzione del segnale mediante TFT costruita sulla base del 355 % supporto di cui al passo prima 356

APPENDICE

104

f1 = frq_cstfm(k,1); % frequenza più bassa rintracciata 357 (fondamentale), in Hz 358 ih = frq_cstfm(k,[2:max_iter])./f1; % interarmoniche della 359 fondamentale, in Hz 360 361 B=matrixTaylorFourier_ih(f1,Fs,N,3*ones(max_iter,1),ih,0); 362 % generazione della matrice TFT in base a fondamentale e interarmoniche 363 364 % coefficienti della TFT 365 stima = pinv(B'*B)*B'*y_run; 366 367 % Porzione di segnale ricostruito 368 z(:,k) = B*stima; % vettore ricostruito 369 370 %% Sottrazione template 371 sig_noise(:,k) = y_run - real(z(:,k)); 372 end 373 374 ric = real(z); % Segnale ricostruito 375 sig_end = ric(N-overlap+1:end,1); 376 sig_tot_ric = ric(1:N-overlap,1); 377 378 for i = 2:tot_iter 379 sig_start = ric(1:overlap,i); 380 mean_sig = mean([sig_start sig_end],2); 381 382 sig_tot_ric = [sig_tot_ric; mean_sig; ric(overlap+1:end,i);]; 383 sig_end = ric(N-overlap+1:end,i); 384 end 385 386 % Segnale ricostruito generale 387 sig_tot_ric = sig_tot_ric(1:Ns); 388 389 % Calcolo la risposta evocata 390 rispEvocate = sig_lonf-sig_tot_ric; 391 392 rispEvocateSing = zeros(durataPorz,numPorz); 393 meanHRF_singlCh_filtr = zeros(durataPorz,numPorz); 394 for i = 1:numPorz 395 % Faccio la baseline su ciascuna risposta evocata (sottraggo la 396 % media della stessa considerato tra -1 e 1 sec) 397 HRF_ric = rispEvocate(ind_start(i):ind_stop(i)); 398 HRF_filtr = totalHRF_filtr(:,i); 399 400 baseLineRic = mean(HRF_ric(1:(2/t_camp))); 401 baseLineFiltr = mean(HRF_filtr(1:(2/t_camp))); 402 403 rispEvocateSing(:,i) = HRF_ric - baseLineRic; 404 meanHRF_singlCh_filtr(:,i) = HRF_filtr - baseLineFiltr; 405 end 406 407 % CALCOLO LA MEDIA DELLE RISPOSTE EVOCATE: 408

APPENDICE

105

% Media degli impulsi calcolata a partire dalle porzioni intere di segnale 409 corrispondenti al trigger 410 mediaPorz_ric = mean(rispEvocateSing,2); 411 mediaPorz_filtr = mean(meanHRF_singlCh_filtr,2); 412 413 % Filtro mediano 414 mediaPorz_medianRic = zeros(1,durataPorz); 415 mediaPorz_medianFiltr = zeros(1,durataPorz); 416 417 % Metà banda, la banda intera del filtro è 80 campioni 418 halfWindow = 40; 419 420 for j = 1:durataPorz 421 i= max([j-halfWindow,1]):min([j+halfWindow,durataPorz]); 422 mediaPorz_medianRic(j) = median(mediaPorz_ric(i)); 423 mediaPorz_medianFiltr(j) = median(mediaPorz_filtr(i)); 424 end 425 426 % Filtro Savitzky - Golay 427 mediaPorz_sgolayRic = sgolayfilt(mediaPorz_medianRic,3,51); 428 mediaPorz_sgolayFiltr = sgolayfilt(mediaPorz_medianFiltr,3,51); 429 430 if iter==1 % Inizializzazione delle matrici contenente tutti i canali 431 mediaHRF_VeraChannel = zeros(durataPorz,numCh_LS); 432 mediaHRF_RicChannel = zeros(durataPorz,numCh_LS); 433 mediaHRF_FiltrChannel = zeros(durataPorz,numCh_LS); 434 rispEvocate_channel = zeros(Ns,numCh_LS); 435 end 436 437 rispEvocate_channel(:,iter) = rispEvocate; 438 mediaHRF_RicChannel(:,iter) = mediaPorz_sgolayRic; 439 mediaHRF_FiltrChannel(:,iter) = mediaPorz_sgolayFiltr; 440 mediaHRF_VeraChannel(:,iter) = mediaHRFvera_singleCh; 441 442 if ~ismember(src_long,channel_MIDDLE) 443 [valPeakRic,indPeakRic] = max(mediaPorz_sgolayRic); 444 [valPeakFiltr,indPeakFiltr] = max(mediaPorz_sgolayFiltr); 445 [valPeakVera,indPeakVera] = max(mediaHRFvera_singleCh); 446 447 t_indPeakVera = tHRF(indPeakVera); 448 t_indPeakRic = tHRF(indPeakRic); 449 t_indPeakFiltr = tHRF(indPeakFiltr); 450 451 errorAmpRic = 100*(abs(valPeakVera-valPeakRic)/abs(valPeakVera)); 452 errorLatRic = 100*(abs(t_indPeakVera-t_indPeakRic)/abs(t_indPeakVera)); 453 errorTotRic = 100*((norm(mediaPorz_sgolayRic-454 mediaHRFvera_singleCh') 2)/(norm(mediaHRFvera_singleCh) 2)); 455 456 errorAmpFiltr = 100*(abs(valPeakVera-valPeakFiltr)/abs(valPeakVera)); 457 errorLatFiltr = 100*(abs(t_indPeakVera-458 t_indPeakFiltr)/abs(t_indPeakVera)); 459 errorTotFiltr = 100*((norm(mediaPorz_sgolayFiltr-460 mediaHRFvera_singleCh') 2)/(norm(mediaHRFvera_singleCh) 2)); 461

APPENDICE

106

462 if ismember(src_long,channel_LEFT) 463 errorAmpRic_LEFT = [errorAmpRic_LEFT; errorAmpRic]; 464 errorLatRic_LEFT = [errorLatRic_LEFT; errorLatRic]; 465 errorTotRic_LEFT = [errorTotRic_LEFT; errorTotRic]; 466 467 errorAmpFiltr_LEFT = [errorAmpFiltr_LEFT; errorAmpFiltr]; 468 errorLatFiltr_LEFT = [errorLatFiltr_LEFT; errorLatFiltr]; 469 errorTotFiltr_LEFT = [errorTotFiltr_LEFT; errorTotFiltr]; 470 471 HRFvera_LEFT = [HRFvera_LEFT mediaHRFvera_singleCh]; 472 HRFric_LEFT = [HRFric_LEFT mediaPorz_sgolayRic']; 473 HRFfiltr_LEFT = [HRFfiltr_LEFT mediaPorz_sgolayFiltr']; 474 475 else 476 errorAmpRic_RIGHT = [errorAmpRic_RIGHT; errorAmpRic]; 477 errorLatRic_RIGHT = [errorLatRic_RIGHT; errorLatRic]; 478 errorTotRic_RIGHT = [errorTotRic_RIGHT; errorTotRic]; 479 480 errorAmpFiltr_RIGHT = [errorAmpFiltr_RIGHT; errorAmpFiltr]; 481 errorLatFiltr_RIGHT = [errorLatFiltr_RIGHT; errorLatFiltr]; 482 errorTotFiltr_RIGHT = [errorTotFiltr_RIGHT; errorTotFiltr]; 483 484 HRFvera_RIGHT = [HRFvera_RIGHT mediaHRFvera_singleCh]; 485 HRFric_RIGHT = [HRFric_RIGHT mediaPorz_sgolayRic']; 486 HRFfiltr_RIGHT = [HRFfiltr_RIGHT mediaPorz_sgolayFiltr']; 487 end 488 else 489 rmse_ric_singleCh = sqrt(sum((mediaHRFvera_singleCh-490 mediaPorz_sgolayRic'). 2)/durataPorz); 491 rmse_filtr_singleCh = sqrt(sum((mediaHRFvera_singleCh-492 mediaPorz_sgolayFiltr'). 2)/durataPorz); 493 494 rmse_ric = [rmse_ric; rmse_ric_singleCh]; 495 rmse_filtr = [rmse_filtr; rmse_filtr_singleCh]; 496 end 497 end 498 499 mediaHRF_RicTot = mean(mediaHRF_RicChannel,2); 500 mediaHRF_VeraTot = mean(mediaHRF_VeraChannel,2); 501 mediaHRF_FiltrTot = mean(mediaHRF_FiltrChannel,2); 502 503 % Errore e STD del segnale ricostruito mediante algoritmo 504 MeanErrorAmp_LEFT_Ric = mean(errorAmpRic_LEFT); 505 MeanErrorLat_LEFT_Ric = mean(errorLatRic_LEFT); 506 MeanErrorTot_LEFT_Ric = mean(errorTotRic_LEFT); 507 MeanErrorAmp_RIGHT_Ric = mean(errorAmpRic_RIGHT); 508 MeanErrorLat_RIGHT_Ric = mean(errorLatRic_RIGHT); 509 MeanErrorTot_RIGHT_Ric = mean(errorTotRic_RIGHT); 510 511 StdErrorAmp_LEFT_Ric = std(errorAmpRic_LEFT); 512 StdErrorLat_LEFT_Ric = std(errorLatRic_LEFT); 513 StdErrorTot_LEFT_Ric = std(errorTotRic_LEFT); 514

APPENDICE

107

StdErrorAmp_RIGHT_Ric = std(errorAmpRic_RIGHT); 515 StdErrorLat_RIGHT_Ric = std(errorLatRic_RIGHT); 516 StdErrorTot_RIGHT_Ric = std(errorTotRic_RIGHT); 517 518 % Errore e STD del segnale solamente filtrato 519 MeanErrorAmp_LEFT_Filtr = mean(errorAmpFiltr_LEFT); 520 MeanErrorLat_LEFT_Filtr = mean(errorLatFiltr_LEFT); 521 MeanErrorTot_LEFT_Filtr = mean(errorTotFiltr_LEFT); 522 MeanErrorAmp_RIGHT_Filtr = mean(errorAmpFiltr_RIGHT); 523 MeanErrorLat_RIGHT_Filtr = mean(errorLatFiltr_RIGHT); 524 MeanErrorTot_RIGHT_Filtr = mean(errorTotFiltr_RIGHT); 525 526 StdErrorAmp_LEFT_Filtr = std(errorAmpFiltr_LEFT); 527 StdErrorLat_LEFT_Filtr = std(errorLatFiltr_LEFT); 528 StdErrorTot_LEFT_Filtr = std(errorTotFiltr_LEFT); 529 StdErrorAmp_RIGHT_Filtr = std(errorAmpFiltr_RIGHT); 530 StdErrorLat_RIGHT_Filtr = std(errorLatFiltr_RIGHT); 531 StdErrorTot_RIGHT_Filtr = std(errorTotFiltr_RIGHT); 532 % ----------------------------------------------------------------------- 533 534 % Media delle HRF emisfero SX 535 meanHRFvera_LEFT = mean(HRFvera_LEFT,2); 536 meanHRFric_LEFT = mean(HRFric_LEFT,2); 537 meanHRFfiltr_LEFT = mean(HRFfiltr_LEFT,2); 538 % Media delle HRF emisfero DX 539 meanHRFvera_RIGHT = mean(HRFvera_RIGHT,2); 540 meanHRFric_RIGHT = mean(HRFric_RIGHT,2); 541 meanHRFfiltr_RIGHT = mean(HRFfiltr_RIGHT,2); 542 543 clc 544 fprintf('\n\n\n\n\n\n\n') 545 display('****************** ERRORE *************************') 546 fprintf('\n') 547 display('---Errore sul segnale ricostruito mediante CFT-----') 548 display('---------------- Emisfero SX ----------------------') 549 display(['MEDIA Errore AMPIEZZA = ' num2str(MeanErrorAmp_LEFT_Ric) ' %']) 550 display(['MEDIA Errore LATENZA = ' num2str(MeanErrorLat_LEFT_Ric) ' %']) 551 display(['MEDIA ERRORE TOTALE = ' num2str(MeanErrorTot_LEFT_Ric) ' %']) 552 display('---------------------------------------------------') 553 display('---------------- Emisfero DX ----------------------') 554 display(['MEDIA Errore AMPIEZZA = ' num2str(MeanErrorAmp_RIGHT_Ric) ' %']) 555 display(['MEDIA Errore LATENZA = ' num2str(MeanErrorLat_RIGHT_Ric) ' %']) 556 display(['MEDIA ERRORE TOTALE = ' num2str(MeanErrorTot_RIGHT_Ric) ' %']) 557 display('---------------------------------------------------') 558 fprintf('\n') 559 display('------SD degli errori di ogni singolo canale-------') 560 display('---------------- Emisfero SX ----------------------') 561 display(['SD Errore AMPIEZZA = ' num2str(StdErrorAmp_LEFT_Ric) ]) 562 display(['SD Errore LATENZA = ' num2str(StdErrorLat_LEFT_Ric) ]) 563 display(['SD ERRORE TOTALE = ' num2str(StdErrorTot_LEFT_Ric) ]) 564 display('---------------------------------------------------') 565 display('---------------- Emisfero DX ----------------------') 566 display(['SD Errore AMPIEZZA = ' num2str(StdErrorAmp_RIGHT_Ric) ]) 567

APPENDICE

108

display(['SD Errore LATENZA = ' num2str(StdErrorLat_RIGHT_Ric) ]) 568 display(['SD ERRORE TOTALE = ' num2str(StdErrorTot_RIGHT_Ric) ]) 569 display('---------------------------------------------------') 570 fprintf('\n\n\n\n\n') 571 display('-------Errore sul segnale solo filtrato------------') 572 display('---------------- Emisfero SX ----------------------') 573 display(['MEDIA Errore AMPIEZZA = ' num2str(MeanErrorAmp_LEFT_Filtr) ' %']) 574 display(['MEDIA Errore LATENZA = ' num2str(MeanErrorLat_LEFT_Filtr) ' %']) 575 display(['MEDIA ERRORE TOTALE = ' num2str(MeanErrorTot_LEFT_Filtr) ' %']) 576 display('---------------------------------------------------') 577 display('---------------- Emisfero DX ----------------------') 578 display(['MEDIA Errore AMPIEZZA = ' num2str(MeanErrorAmp_RIGHT_Filtr) ' %']) 579 display(['MEDIA Errore LATENZA = ' num2str(MeanErrorLat_RIGHT_Filtr) ' %']) 580 display(['MEDIA ERRORE TOTALE = ' num2str(MeanErrorTot_RIGHT_Filtr) ' %']) 581 display('---------------------------------------------------') 582 fprintf('\n') 583 display('------SD degli errori di ogni singolo canale-------') 584 display('---------------- Emisfero SX ----------------------') 585 display(['SD Errore AMPIEZZA = ' num2str(StdErrorAmp_LEFT_Filtr) ]) 586 display(['SD Errore LATENZA = ' num2str(StdErrorLat_LEFT_Filtr) ]) 587 display(['SD ERRORE TOTALE = ' num2str(StdErrorTot_LEFT_Filtr) ]) 588 display('---------------------------------------------------') 589 display('---------------- Emisfero DX ----------------------') 590 display(['SD Errore AMPIEZZA = ' num2str(StdErrorAmp_RIGHT_Filtr) ]) 591 display(['SD Errore LATENZA = ' num2str(StdErrorLat_RIGHT_Filtr) ]) 592 display(['SD ERRORE TOTALE = ' num2str(StdErrorTot_RIGHT_Filtr) ]) 593 display('---------------------------------------------------') 594 595 fprintf('\n\n\n\n') 596 display('************** Canali FRONTALI ********************') 597 display(['Media RMSE segnale filtrato: ' num2str(mean(rmse_filtr))]) 598 display(['Media RMSE segnale ricostruito: ' num2str(mean(rmse_ric))]) 599 display('***************************************************') 600 601 % Salvataggio dei dati 602 eval(['save RESULT_subj' num2str(soggetto) '_rs' num2str(rs) ' 603 rispEvocate_channel mediaHRF_RicChannel mediaHRF_VeraChannel 604 mediaHRF_FiltrChannel HRFvera_LEFT HRFvera_RIGHT HRFric_LEFT HRFric_RIGHT 605 HRFfiltr_LEFT HRFfiltr_RIGHT rmse_ric rmse_filtr']) 606

607

109

BIBLIOGRAFIA

[1] GRATTON G., FABIANI M., FRIEDMAN D., FRANCESCHINI M.A.,

FANTINI S., CORBALLIS P., GRATTON E., Rapid Changes of Optical

Parameters in the Human Brain during a Tapping Task, Journal of Cognitive

Neuroscience, 7, 446-456, 1995.

[2] HOSHI Y., TAMURA M., Near-Infrared Optical Detection of Sequential Brain

Activation in the Prefrontal Cortex during Mental Tasks, NeuroImage, 5, 292-297,

1997.

[3] JÖBSIS, F.F., Noninvasive infrared monitoring of cerebral and myocardial sufficiency

and circulatory parameters, Science, 198, 1264-1267, 1977.

[4] SCHOLKMANN F., KLEISER S., JAAKKO METZ A.,

ZIMMERMANN R., MATA PAVIA J., WOLF U., WOLF M., A review on

continuous wave functional near-infrared spectroscopy and imaging instrumentation and

methodology, NeuroImage, Volume 85, 2014, 6-27, 2014.

[5] BOAS D., FRANCESCHINI M., DUNN A., STRANGMAN G.,

Nonivasive Imaging of Cerebral Activation with Di_use Optical Tomography, In Vivo

Optical Imaging of Brain Function, CRC Press, 8, 196-216, 2002.

BIBLIOGRAFIA

110

[6] CHANCE B., ZHUANG Z., UNAH C., ET AL., Cognition-activated low-

frequency modulation of light absorption in human brain., Proc. Natl. Acad. Sci. USA

90 (8), 3770–3774, 1993.

[7] HOSHI Y., TAMURA M., Detection of dynamic changes in cerebral oxygenation

coupled to neuronal function during mental work in man., Neurosci. Lett. 150 (1),

1993.

[8] KATO T., KAMEI A., TAKASHIMA S., ET AL., Human visual cortical

function during photic stimulation monitoring by means of near-infrared spectroscopy. J.

Cereb. Blood Flow Metab. 13 (3), 516–520. 5–8., 1993.

[9] VILLRINGER A., PLANCK J., HOCK C., ET AL., Near infrared

spectroscopy (NIRS): a new tool to study hemodynamic changes during activation of brain

function in human adults. Neurosci. Lett. 154 (1–2), 101–104., 1993.

[10] P. DAL BIANCO, Approaches based on nonlinear Kalman filtering to deal with

physiological noise in the estimation of the hemodynamic response from Functional Near-

Infrared Spectroscopy data, Master Thesis, 2014.

[11] AMBROSIN A., Master Thesis: Kalman filter approach for the estimation of the

hemodynamic response from near-infrared spectroscopy data (fNIRS), Università di

Padova, 2013.

[12] HUPPERT T. J., DIAMOND S. G., FRANCESCHINI M. A., BOAS D.

A., HomER: a review of time-series analysis methods for near infrared spectroscopy of the

brain, Appl Opt., 48(10): D280–D298 1 Aprile 2009.

BIBLIOGRAFIA

111

[13] LLOYD-FOX S., BLASI A., ELWELL C. E., Illuminating the developing

brain: the past, present and future of functional near infrared spectroscopy, Neuroscience

and Biobehavioral Reviews, 34, 269-284, 2010.

[14] FRIGO G, BRIGADOI S., GIORGI G., SPARACINO G.,

NARDUZZI C., Measuring Cerebral Activation From fNIRS Signals: An Approach

Based on Compressive Sensing and Taylor-Fourier Model, IEEE Transactions on

Instrumentation and Measurement, vol. PP, no. 99, 1-9, doi:

10.1109/TIM.2016.2518363.

[15] SAAGER R. B., BERGER A. J., Direct characterization and removal of

interfering absorption trends in two-layer turbid media, Optical Society of America,

Vol. 22, no. 9, 1874 – 1882, , September 2005.

[16] SCARPA F., BRIGADOI S., CUTINI S., SCATTURIN P., ZORZI M.,

DELL'ACQUA R., SPARACINO G., A reference-channelbasedmethodology to

improveestimation of event-related hemodynamic response from fNIRSmeasurements,

NeuroImage, Volume 72, 106-119, 15 May 2013, ISSN 1053-8119.

[17] GLOVER G. H., Deconvolution of Impulse Response in Event-Related BOLD

fMRI, NeuroImage, Volume 9, Issue 4, 416-429, April 1999, ISSN 1053-8119

[18] BERTOCCO M., FRIGO G., NARDUZZI C., On compressed sensing

and super-resolution in DFT-based spectral analysis, Proceedings 19th

IMEKO TC-4 Symposium and 17th IWADC Workshop Advances in

Instrumentation and Sensors Interoperability, 615-620, Barcelona, 18 July

2013.

BIBLIOGRAFIA

112

[19] TROPP J. A., GILBERT A. C., Signal Recovery From Random Measurements

Via Orthogonal Matching Pursuit, in IEEE Transactions on Information Theory,

vol. 53, no. 12, 4655-4666, Dec. 2007, doi: 10.1109/TIT.2007.909108.

[20] http://www.iopscience.iop.org/0967-3334/33/6/925/article

Approccio basato su Compressive Sensing e modellistica...

Documents

Transcript of Approccio basato su Compressive Sensing e modellistica...