ANALISI DI SEGNALI BIOMEDICI: STIMA DELLE PERIODICITA`

ANALISI DI SEGNALI BIOMEDICI:
STIMA DELLE PERIODICITA’
Se x(t) è un segnale (deterministico) limitato su [0, N0-1], nullo fuori da tale intervallo,
l’autocorrelazione (AC) Rxx(m) è definita come:
m= ritardo (“lag”)
Rxx(m)
Rxx(m) rappresenta il grado di similitudine medio fra x(k) e x(k+m) e dà una misura della
memoria del sistema.
Proprietà dell’autocorrelazione:
Rxx(-m) = Rxx(m)
| Rxx(m)| | Rxx(0)|
Quanto detto è relativo al caso di segnali TD. Analogamente per segnali TC.
AUTOCORRELAZIONE - TD
L’autocorrelazione (AC) Rxx [m] indica le relazione lineare
fra due punti del segnale x distanti m.
Problema:
I segnali reali sono quasi-periodici (di durata finita), perciò i
picchi della Rxx [m] non sono molto pronunciati, con difficoltà
nella loro individuazione.
Es: Segnale vocale: Segnali periodici (vocalici) hanno valori
elevati dell’autocorrelazione, mentre i suoni non vocalici
(consonanti) o vocalici ma fortemente disfonici sono
caratterizzati da bassi valori di Rxx [m].
ESEMPIO: SEGNALE VOCALICO
PERIODO T ?
PERIODO T
/a/ post-surgical
0.6
0.4
0.4
Normalised amplitude [arb.units]
Normalised amplitude [arb. units]
/a/ pre-surgical
0.6
0.2
0
-0.2
-0.4
-0.6
0.2
0
-0.2
-0.4
-0.6
2
4
6
Time [s]
8
10
x 10
1
2
3
4
-3
PATOLOGICO PRE-INTERVENTO
5
6
7
8
Time [s]
POST-INTERVENTO
T = 1/F0 - F0 = frequenza fondamentale
9
10
11
-3
x 10
ESEMPIO
L’autocorrelazione (AC) può essere usata come primo test di verifica sulla periodicità del segnale
vocale: suoni non vocalici (consonanti) o disfonici (dovuti a patologie o malfunzionamento) hanno
valori bassi di AC. Negli esempi Fs=44.1kHz. A seconda del tipo di segnale si analizzano frame di
dimensione 5mst 20ms  220 N  880  100 lag  400  2.5mst 10ms
Valore basso di autocorrelazione
non vocalico
vocalico
Valore elevato di autocorrelazione
lag
ESEMPI
Neonato sano
Tenore
lag
Paralisi cordale
AUTOCORRELAZIONE - AC
Alcune proprietà dell’AC:
(*=trasposto coniugato)
Matrice di AC costruita
con M+1 valori dell’AC:
E’ una matrice che viene utilizzata spesso nell’analisi dei sistemi dinamici. E’
una matrice hermitiana (rxx(k)=rxx(-k)) e Toeplitz (tutti gli elementi lungo ogni
diagonale sono uguali).
AUTOCOVARIANZA AC - TD
Se il segnale è a media non nulla, l’ autocorrelazione può essere mascherata da
tale valore.
L’autocovarianza Cxx(m) è la misura della memoria del sistema, relativamente
ai suoi scostamenti rispetto al valore medio :
AC e autocorrelazione sono legate dalla relazione:
2
R x (m
x )  C x (m
x )  x
Quanto detto è relativo al caso di segnali discreti (TD). Analogamente per segnali
continui (TC).
ESEMPIO
x(n) = n per 1n3
x(0)=0; x(n)=0 per n>3
 N0=4
Valore medio:
Autocorrelazione:
ESEMPIO (CONT.)
AC:
R x (m)  C x (m)  x 2
Cx(0)=3.5-(1.5)2=1.25 Cx(1)=2- (1.5)2=0.25 CX(2)=0.75-(1.5)2=-0.5 Cx(m)=-2.25 , m3
ES3 MATLAB confrcorr.m
3 segnali random (scorrelati) a media nulla di 64, 512, 4096 punti risp. La AC dovrebbe
essere uguale a zero. I segnali con pochi dati (blu e rosso) non consentono di ottenere
una media nulla.
2
C x (m
x )  R x (m
x ) x
v. funzioni Matlab: showcorr.m, showdtcorr.m, showcorr_sine, covstat.m
ESEMPIO
La funzione di AC ha grande importanza nell’analisi dei segnali biomedici.
ES.: Valutare l’indipendenza di un atto respiratorio dagli altri. Da punto di vista clinico, è
utile per stabilire le connessioni neuronali nel midollo allungato che generano il ritmo
respiratorio.
L’approccio comunemente utilizzato è quello di costruire un modello matematico di rete
neuronale che simuli il funzionamento dei neuroni reali.
Il modello deve produrre un’oscillazione simile al ritmo respiratorio, ed avere altre
caratteristiche fisiologiche.
Attualmente, i modelli hanno la proprietà che ogni atto respiratorio è indipendente dagli
altri. Ma il segnale reale verifica questa condizione?
Per rispondere a questa domanda, si definisce:
v(n) = volume d’aria inspirato all’n-mo atto respiratorio;
N = numero di atti respiratori (N=120 in questo esempio);
Si calcola la AC e si osserva il risultato ottenuto:
ATTI RESPIRATORI
Max per lag=0
AC >0 fino a lag=25
AC <0 per lag>25
Questo contraddice l’ipotesi di indipendenza  i modelli attuali non tengono conto di elementi fisiologici importanti
ES2 MATLAB ranproc2.mat
Flusso respiratorio di 8
ratti.
Ogni
registrazione
comprende circa 30-40
atti respiratori.
Fs= 75Hz.
I segnali sono simili, ma
non uguali.
Ci
sono
differenze
anche all’interno della
singola
registrazione
(realizzazione).
flusso respiratorio di 2 ratti a riposo (n.1 e n.2) - Fs=75 Hz
Example2.m
1
Caso n.2
25
0
-1
0
500
1000
1500
1
0
-1
0
500
1000
1500
0
-1
0
500
1000
1500
0
-1
0
500
1000
1500
2000
0
500
1000
1500
2000
0
500
1000
1500
2000
0
500
1000
1500
2000
0
500
1000
1500
2000
1
0
1
0
-1
2000
1
0
-1
2000
1
1
-1
2000
Caso n.4
20
Caso n.6
15
tempo [s]
Caso n.8
10
Caso n.3
5
Caso n.5
0
Caso n.1
0
Caso n.7
ampiezza
Grafici ottenuti con
MATLAB
1
0
-1
PSD (periodogramma)
ES2/6 MATLAB - ranproc2.mat
Flusso respiratorio
Freq.respiratoria
1° registrazione
Picco nella PSD (Power Spectral Density): frequenza respiratoria media del
soggetto. Un confronto fra questi valori su tutti i casi può dare un’indicazione sulla
ipotesi di processo stocastico per questo esempio.
ES.6 MATLAB – PLOTRP2.M
Grafico della derivata del segnale vs. il segnale: grafici sovrapposti
0.3
0.2
0.1
0
-0.1
-0.2
-0.3
-0.4
-0.8
-0.6
-0.4
-0.2
0
0.2
0.4
0.6
0.8
plotrp2.m - ranproc2.mat
PSD del 4° segnale: scala lineare
20
15
10
5
0
0
5
10
15
20
25
20
25
PSD del 4° segnale: scala logaritmica
5
10
0
10
-5
10
-10
10
0
5
10
15
Freq. (Hz)
MATLAB ES6 - ranproc3.mat
Segnale EMG, Fs= 500Hz: contrazione del muscolo della lingua estrinseco
genio-glosso. Registrazioni di 8 contrazioni diverse (realizzazioni) dello stesso
soggetto, ciascuna contenente 256 dati.
PSD della 5° registrazione
Non ci sono picchi evidenti nella PSD, ma tutti gli spettri sono simili, indicando la
provenienza dallo stesso processo stocastico.
ES.6 MATLAB - PLOTRP3.M
4
4
x 10 Grafico della derivata del segnale vs. il segnale: grafici sovrapposti
3
2
1
0
-1
-2
-3
-3
-2
-1
0
1
2
3
4
x 10
ranproc3.mat
7
6
PSD del 4° segnale: scala lineare
x 10
4
2
0
0
50
100
150
200
250
200
250
PSD del 4° segnale: scala logaritmica
8
10
6
10
4
10
0
50
100
150
Freq. (Hz)
ES4 MATLAB hrv1.mat
Frequenza
cardiaca
istantanea in
battiti al
minuto (cioè:
60/durata
battito) in
soggetto a
riposo.
“memoria”
del sistema
ES.4 – HRV1COV.M
Freq. cardiaca istantanea
Heart rate
100
90
80
70
0
50
100
150
time
AC normalizzata
200
250
300
0
50
100
150
lag
200
250
300
1
0.5
0
-0.5
ES.4 MATLAB hrv.mat
A sinistra: 10 registrazioni di frequenza cardiaca dallo stesso individuo. A
destra: AC delle 10 registrazioni.
HRVCOV.M - HRV.MAT
AC per le prime due registrazioni
40
20
0
-20
0
50
100
150
200
lag
AC per le 10 registrazioni
0
50
100
250
300
250
300
40
20
0
-20
150
lag
200
SEGNALE: ENERGIA E POTENZA
Per un segnale deterministico x(t), la potenza
istantanea all’istante t è (*=trasposto coniugato):
La potenza istantanea all’istante t data
dall’interazione di due segnali x(t) e y(t) è:
La potenza media su un intervallo [t0, t0+T]
è:
L’ energia di un segnale x(t) è
l’integrale della potenza nel tempo :
L’energia data dall’interazione di due
segnali x(t) e y(t) è exy():
Due segnali si dicono scorrelati se
exy() =0
DENSITA’ SPETTRALE
La densità spettrale di energia (Energy Spectral
Density, ESD), o spettro di energia, del segnale x(t) è:
La cross-densità spettrale di energia (Energy
Spectral Density, ESD), o spettro di energia, fra due
segnali x(t) e y(t) è:
Il teorema di Parseval afferma che l’energia totale del segnale è indipendente dalla scelta
della sua rappresentazione nel tempo o in frequenza (conservazione dell’energia):
Poiché il prodotto di due trasformate di Fourier è la trasformata di Fourier della convoluzione
delle due funzioni nel tempo, il cross-spettro di energia è (teorema di Wiener-Khintchine):
PROCESSO STOCASTICO
I valori futuri di un segnale aleatorio (stocastico) non possono essere
predetti esattamente. I segnali aleatori sono tali per loro natura
(meccanismo interno) o rappresentano meccanismi che non conosciamo
esattamente.
Esiste variabilità fra i risultati di esperimenti diversi ed internamente al
singolo risultato.
Es.: un gruppo di individui a cui viene misurata la pressione costituisce
l’insieme degli “esperimenti” e il singolo individuo è una “realizzazione”.
Il valore di una singola realizzazione ad uno specifico istante di tempo è una
variabile aleatoria, cioè una variabile il cui valore dipende da un evento
casuale. Ogni realizzazione è diversa dalle altre.
Scopo: ottenere rappresentazioni (modelli) di realizzazioni di processi
stocastici da cui ricavare le proprietà ed i parametri del processo
stocastico.
Quasi tutti i segnali biomedici possono essere visti come realizzazioni di
processi stocastici, e spesso utilizzare valori medi invece della singola
realizzazione non è il principale obbiettivo dell’analisi.
PROCESSO STOCASTICO
Insieme di successioni (TC o TD) ciascuna derivante da un diverso
esperimento:
x(n;i)
i = i-ma successione (esperimento, o “realizzazione”)
n = indice temporale
Fissato i, si utilizza la notazione semplificata x(n).
Un processo stocastico è stazionario se il suo valore medio è costante
e l’AC dipende solo dalla differenza temporale m=n2-n1.
Due processi stocastici stazionari x(n) e y(n) sono caratterizzati da
cross-correlazione rxy(m) e cross-covarianza cxy(m):
ε = valor medio; * = complesso coniugato
ESEMPIO
L’uscita del processo stocastico è influenzata da componenti di rumore
aleatorie in ogni realizzazione, pertanto i vari segnali in uscita pi(t) non
descrivono esattamente il processo reale.
SEGNALI STOCASTICI
Per segnali stocastici, al posto della
media si considera il valore atteso o
valore medio o momento del 1°
ordine, cioè il valore a cui converge
la media di un’insieme di osservazioni
(realizzazioni). (p=densità di
probabilità)
Il valore atteso di x2 è il momento del
2° ordine.
La varianza è la deviazione
quadratica della variabile aleatoria dal
suo valore medio.
Bias = differenza fra il valore vero dei
parametri e il valore atteso della loro
stima. Uno stimatore è consistente se
bias e varianza tendono a zero
all’aumentare del numero delle
osservazioni.
POWER SPECTRAL DENSITY
Power Spectral Density (PSD) = estensione statistica della ESD a
segnali stocastici.
Il teorema (equazione) di Parseval afferma che l’energia di un segnale x(t) è
distribuita fra le sue componenti in frequenza in modo tale che l’energia ad
ogni frequenza f è proporzionale al quadrato dell’ampiezza di X(f), |X(f)|2 .
Il teorema di Parseval costituisce uno strumento efficiente per analizzare la
distribuzione dell’energia in un segnale.
La Densità Spettrale di Potenza (PSD, Power Spectral Density) di un
segnale (di lunghezza N) è rappresentata con il grafico di |X(f)|2 in funzione di
f, diviso per il n. N di campioni di segnale:
Pxx (f ) 
X( f )
N
2
PSD
La Power Spectral Density (PSD) è definita come la trasformata di
Fourier (discreta) della sequenza di autocorrelazione rxx(m):
rxx(m)=Ex(n+m)x*(n)
(T = periodo di campionamento):
Questa equazione mette in relazione le proprietà di “memoria”
(autocorrelazione) con la funzione di densità spettrale di potenza di un
sistema, e afferma che la potenza ad ogni frequenza riflette l’ampiezza
della componente sinusoidale a quella frequenza nella funzione di
autocorrelazione.
PSD PER RUMORE BIANCO
Rumore bianco w(n): processo a media nulla scorrelato per tutti i valori di m,
tranne che per m=0, che corrisponde alla varianza w. L’AC è:
PSD costante  f
(m)=successione delta discreta.
w(n) ha componenti a tutte le frequenze: da qui il nome di “rumore bianco”.
EEG: il problema degli artefatti /1
artefatti muscolari
Realizzazioni dello stesso processo stocastico
Dipartimento di Elettronica Informatica e Sistemistica - UNIVERSITA’ di BOLOGNA
EEG di superficie /1
La corteccia cerebrale contiene
diversi tipi di cellule nervose che
possono venire raggruppate in due
classi principali:
Neuroni piramidali
Interneuroni
L’EEG registra variazioni del campo
elettrico generato da gruppi di neuroni
piramidali, mentre la MEG registra
variazioni del campo magnetico indotto
dal variare del campo elettrico generato
dagli stessi neuroni.
Dipartimento di Elettronica Informatica e Sistemistica - UNIVERSITA’ di BOLOGNA
EEG di superficie /1
basi fisiologiche della generazione del segnale:
Il segnale EEG è il risultato dell’attività elettrica assonale e denditica dei neuroni
piramidali corticali – serve un’attività sincrona di un numero elevato di neuroni (105)
per generare un segnale apprezzabile sulla superficie dello scalpo
Ampiezza segnale: dell’ordine delle decine di μV
Bande di frequenza di interesse: δ, θ, α, β, γ (range: 0.5-oltre 30Hz)
sistemi per l’acquisizione del segnale
segnale registrato attraverso elettrodi di superficie.
numero elettrodi: 19, 32, 64, 128,...
modi per registrarlo: tanti tipi di elettrodo, tra cui le cuffie EEG
Frequenza di campionamento: 200 : 1000 samples/sec
Digitalizzazione a 12 bit
impedenza di contatto elettrodo-cute da mantenere sotto i 5 Kohm
Dipartimento di Elettronica Informatica e Sistemistica - UNIVERSITA’ di BOLOGNA
EEG sistema 10/20
Dipartimento di Elettronica Informatica e Sistemistica - UNIVERSITA’ di BOLOGNA
EEG di superficie /2
gli elettrodi vengono posizionati secondo il sistema di riferimento
internazionale 10-20.
Ciascun elettrodo e’ definito rispetto:
- all’ area cerebrale sottostante (F= frontale, P= parietale, C=‘centrale’
per il vertice, T= temporale, O= occipitale);
- alla linea mediana (numero pari per gli elett. destri, dispari per gli elett.
sinistri, e z per gli elett. mediani).
Ad es, F3 indica un elettrodo frontale sinistro, Cz un elettrodo centrale
mediano.
Dipartimento di Elettronica Informatica e Sistemistica - UNIVERSITA’ di BOLOGNA
EEG: ritmi di base /1
Passando dallo stato
di veglia a quello di
sonno e coma, le
onde EEG diventano
progressivamente
piu’ ampie e dalle
componenti frequenziali
più basse
(regola generale, con le
dovute eccezioni..)
Dipartimento di Elettronica Informatica e Sistemistica - UNIVERSITA’ di BOLOGNA
EEG: ritmi di base /2
banda
Hz
delta
0,5-3
teta
3-7
Sonno profondo, memoria episodica
(corteccia del cingolo anteriore ippocampo)
alfa
8-13
Rilassamento mentale / oscillazione “idling”
(occipitale / occhi chiusi, somatosensoriale
(ritmo μ))
beta
14-30
Attenzione, concentrazione, aree corticali
attivate, attivazione aree motorie
>30
Attenzione, concentrazione, aree corticali
attivate, processi di integrazione della
percezione
gamma
Funzione localizzazione
Condizioni patologiche (coma)
Dipartimento di Elettronica Informatica e Sistemistica - UNIVERSITA’ di BOLOGNA
EEG: ritmi di base /3
Onde delta
Onde theta
Onde alpha
Onde beta
Dipartimento di Elettronica Informatica e Sistemistica - UNIVERSITA’ di BOLOGNA
ESEMPIO: EEG
Attività elettrica fra le membrane dei neuroni corticali del cervello. Si distinguono 3
principali tipi di attività, che possono variare circa 15 volte durante il sonno notturno
regolare:
Ritmo  = stato di veglia: 8-12 Hz, bassa ampiezza
Ritmo  = sonno leggero: 4-8 Hz
Ritmo  = sonno profondo: <4 Hz, ampiezza elevata (fase del sogno: simile al ritmo
, ma con rapidi movimenti oculari (REM=Rapid Eye Movements))
Si registra l’EEG durante la notte.
Si assegna un “indice” di sonno
all’EEG ogni 30 sec., sulla base
del ritmo dominante in quei 30
sec.
Caso patologico: la transizione
da un ritmo ad un altro avviene in
meno di 30 sec.
Problema:
distinguere i
3 ritmi
ES5 MATLAB – eegsim1.m
20s. di segnale simulato– transizione
rapida da  ad  (dopo 10s. circa)
Segnale EEG simulato, campionato a 50
Hz per circa 4 min. Simula i 3 ritmi.
Ogni ritmo dura almeno 15s.
E’ un segnale semplificato: le transizioni
sono più rapide del caso reale, non si
considera lo stato REM, le variazioni in
ampiezza sono minori del caso reale.
Ipotesi di ergodicità: si dispone di
un’unica realizzazione del processo
stocastico.
Spettrogramma: grafico tempofrequenza dell’intensità del segnale.
Ritmo  (0-50s), seguito da ritmo 
(50-80s) e poi  (80-120s).
120s-150s: si torna al ritmo , poi si
procede con ritmi a freq. più bassa
v. eegsim1.m, tfar.m, montage.mat