Guida Matlab per il corso di Elaborazione di segnali biologici

X

−

y(n) = a y + u

k n−k n

k=1 2

Supponiamo u gaussiano a media nulla e varianza σ .

n u

Operazioni preliminari

Prima di partire con l’identificazione di un modello, bisogna sempre fare delle operazioni

preliminari. 15

load dati_lab % contiene s

1 s = s - mean ( s ) ; % sottraggo la media

2

3 N = length ( s ) ;

4 Fs = 1000;

5 Ts = 1/ Fs ; % periodo di campionamento

6 t = 0: Ts :( N -1) * Ts % vettore dei tempi

7 Sottrarre la media è importante perché garantisce che la realizzazione del processo su cui

andremo a lavorare sia a media nulla, che è un prerequisito fondamentale.

Identificazione dei parametri

Identificare un modello AR(p) di ordine p dato, vuol dire identificare (p + 1) parametri:

2

a , ..., a e σ .

1 p u

th = ar (s ,p , ’ yw ’) ; % stima i parametri tramite Yule - Walker

1 Estrazione dei parametri:

th . a % coeff del polinomio al den della FDT

1 th . noisevariance % stima di Jmin ( varianza del rumore u)

2

3 fprintf ( ’\n - - - Modello AR (% d ) - - -\ n ’ , p ) ;

4 fprintf ( ’ Coefficienti a_k :\ n ’) ;

5 disp ( th . a )

6 fprintf ( ’ Varianza rumore : %.4 f \ n ’ , th . NoiseVariance ) ;

7 Dalla teoria sappiamo che per ogni modello AR il coefficiente a = 1 sempre, per cui il

0

primo valore del vettore sarà sempre 1.

th.a

Errore di predizione − −

L’errore di predizione e = s ŝ = s y si calcola come l’uscita del filtro inverso. Se la

n n n n n

stima è buona deve approssimare una realizzazione di rumore bianco (ovvero con campioni

scorrelati tra loro).

err = filter ( th .a , 1 , s ) ;

1 err_norm = err / sqrt ( th . noisevariance ) ; % err normalizzato

2

3 plot (t , err_norm )

4 title ([ ’ Errore di predizione normalizzato con modello AR ( ’

5 num2str ( p ) ’) ’ ])

Se il grafico di associato ad un certo ordine p appare come una funzione che

err norm

oscilla molto velocemente attorno al valor medio, senza pattern e senza trend

crescente o decrescente, allora approssima dei campioni scorrelati. Allora il relativo

modello AR di ordine p è sufficientemente adatto a descrivere i dati.

16

Comunque a volte è difficile valutare quale ordine p è migliore rispetto ad un altro guardando

solamente il grafico dell’errore di predizione normalizzato. Ci sono altri metodi indicativi.

5. Scelta dell’ordine ottimo p

Se l’ordine p non è noto, si può trovare tramite due criteri di parsimonia.

I due criteri fanno uso di due funzioni, definite come segue, ma sono anche già calcolate nel

momento di identificazione di un modello tramite la funzione ar.

AIC ( p ) = log ( Jmin ) + 2*( p +1) / N

1 FPE ( p ) = (( N +p -1) /( N -p -1) ) * Jmin

2 th . Report . Fit . FPE % valore FPE (p)

3 Se sto confrontando vari ordini p per capire quale modello AR(p) approssimerebbe meglio i

dati, l’ordine ottimo è quello per cui si verifica min(AIC) o min(FPE).

ordini = 1:25;

1 AIC = zeros ( length ( ordini ) ,1) ;

2 for i =1: length ( ordini )

3 p = ordini ( i ) ;

4 th = ar (s ,p , ’ yw ’) ;

5 J ( i ) = th . noisevariance ;

6 AIC ( i ) = log ( J ( i ) ) +2*( p +1) / N ;

7 end

8

9 [v , p_ottimo_aic ] = min ( AIC ) ; % trova AIC (p) e l ’ ordine p

10 th = ar (s , p_ottimo_aic , ’ yw ’) ; % stima dei parametri

11 Poi si può visualizzare l’andamento del valore di AIC al variare dell’ordine p ed evidenziare

l’ordine ottimo scelto.

figure

1 box on

2 hold on

3 plot ( ordini , AIC , ’bo - ’)

4 plot ( p_ottimo , AIC ( p_ottimo ) , ’r * ’)

5 xlabel ( ’p ’)

6 title ( ’ AIC ( p ) ’)

7 Prendere come modello finale non ci dà la certezza che questo sia il

AIC(p ottimo)

modello più adatto, soprattutto per N elevato. Cosa guardo?

- Ordini superiori? No, avranno (teoricamente) valori di AIC più grandi.

- Ordini inferiori? Se la curva di AIC al variare di p è molto piatta, è probabile che anche

un ordine inferiore potrebbe essere sufficiente a descrivere bene il processo.

17

6. Validazione di un modello AR(p)

2

Il modello identificato (p, a , ..., a e σ ) è adatto a descrivere il processo fisico s che stiamo

1 2 u

esaminando? Abbiamo visto che l’errore di predizione si calcola come l’uscita del filtro

inverso e se la stima è buona deve approssimare una realizzazione di rumore bianco. Nella

pratica, l’errore approssima una realizzazione di un rumore bianco se la stima dei coefficienti

di correlazione è: (

−1−k

N

P 1 ,k = 0

e e

n n+k

n=0

ρ

ˆ (k) = =

e −1

N

P piccolo, altrove

2

e

n

n=0

Test di Bianchezza (Anderson)

Passaggi: ̸

1. L’ipotesi nulla è H = ”e è bianco”. Sotto H , lo stimatore ρ̂ (k) (per k = 0) segue

0 n 0 e 1

2 ,

approssimativamente una distribuzione normale con media m = 0 e varianza σ = N

dove N è il numero di campioni.

th = ar (s , p_ottimo , ’ yw ’) ;

1 err = filter ( th .a ,1 , s ) ;

2

2. Stimare i coefficienti di correlazione ρ̂ (k) per k = 1, . . . , k , tramite una function:

e max

k_max = 25;

1 function ro = stima_coeff (y , k_max ) % y sara ’ err

2 ro = zeros ( k_max ,1) ;

3 for k =1: k_max

4 ro ( k ) = sum ( y (1: end - k ) .* y (1+ k : end ) ) ;

5 end

6 ro_zero = sum ( y .^2) ;

7 ro = ro / ro_zero ;

8 end

9

3. Dato N numero di campioni e fissato il livello di significatività α, trovo gli estremi

−

dell’intervallo di confidenza [−β, β] in cui cade (1 α)% dei valori di ρ

ˆ (k)

e

N = length ( err ) ;

1 sigma = sqrt (1/ N ) ;

2 alpha = 0.05;

3

4 bb = norminv ([ alpha /2 1 - alpha /2] , 0 , sigma ) ; % vettore con

5 -b ,+ b

beta = bb (2) ;

6

4. Calcolo la percentuale dei valori di ρ

ˆ (k) che cadono fuori dall’intervallo

e

fuori = length ( find ( abs ( ro ) > beta ) ) ;

1 perc = fuori / k_max *100;

2 18

disp ([ ’I valori fuori dalla soglia sono ’ , num2str ( perc ) , ’% ’

3 ]) ;

5. Regola di decisione: se al più una frazione α delle stime cade fuori dalla fascia [−β, β],

allora accetto l’ipotesi H

0

if perc > alpha *100

1 disp ( ’ Ipotesi nulla RIFIUTATA : errore NON bianco ’) ;

2 else

3 disp ( ’ Ipotesi nulla NON RIFIUTATA : errore bianco ’) ;

4 end

5

Dopo aver ricavato l’intervallo di confidenza [−β.β] si può visualizzare il grafico che contiene

il valore della stima dei coefficienti di autocorrelazione e i limiti delle soglie dell’intervallo di

confidenza.

figure

1 box on

2 hold on

3 plot (1: k_max , ro , ’o - - ’) % coeff . ro stimati

4 plot ([0 k_max ] , [0 0] , ’r - - ’) % linea dello zero (m)

5 plot ([0 k_max ] , [ bb (1) bb (1) ] , ’k - - ’) % limite inf (-b)

6 plot ([0 k_max ] , [ bb (2) bb (2) ] , ’k - - ’) % limite sup (+ b)

7 title ( ’ Stima dei coefficienti \ rho per il test di Anderson ’)

8 xlabel ( ’k ’) , ylabel ( ’ Valori normalizzati ’)

9 ≈

Se quasi tutti i puntini stanno tra le linee nere, il modello è valido (errore rumore

bianco). Se molti puntini stanno fuori, allora il modello non spiega tutto: residuo ha

ancora struttura 19

Lab 05: Stima Spettrale

Ci sono due diversi metodi per la stima dello spettro di un segnale. Ill primo è il metodo

basato sulla trasformata di Fourier, che a sua volta sottintende due metodi. Il secondo è

basato sull’utilizzo dei modelli ARMA.

Dal grafico della stima spettrale posso identificare quali sono le frequenze principali nel

segnale x.

Operazioni preliminari

load dati_lab % contiene x

1 x = x - mean ( x ) ; % sottraggo la media

2

3 N = length ( x ) ;

4 Fs = 1000;

5 Ts = 1/ Fs ; % periodo di campionamento

6 t = 0: Ts :( N -1) * Ts % vettore dei tempi

7 1. Metodi FT-based

Metodo diretto o Periodogramma

1 2

|F

P̂ (ω) = T [x(n)]|

P ER N

Dato un segnale x con frequenza di campionamento F s è possibile calcolare la stima del suo

spettro tramite il metodo del periodogramma.

FTx = fft (x , N ) ; % Trasformata di Fourier

1 P_per = ( abs ( FTx ) .^2) / N ; % Periodogramma ( DSP )

2 f_FT = (0: Fs / N : Fs - Fs / N ) ; % Vettore frequenze

3 P_per = P_per (1: N /2) ; % Taglio la seconda meta ’

4 f_FT = f_FT (1: N /2) ;

5

6 plot ( f_FT , P_per )

7 xlabel ( ’ Frequenze ( Hz ) ’)

8 title ( ’ Densita ’ spettrale di potenza ( Periodogramma ) ’)

9 20

La stima è abbastanza precisa (la distorsione del metodo del periodogramma nella stima

della DSP è minima)

Metodo indiretto

P̂ (ω) = F T [

R̂ (m)], dove R̂ (m) è la stima dell’autocorrelazione.

BT x x

M = round ( N /2) ; % identico al k_max

1 autocorr = zeros (1 ,2* M +1) ; % vettore [-M ,0 , M]

2 autocorr ( M +1: end ) = stima_autocorr (x , M ) ; % autocorrelazione

3 per k =0: M

autocorr (1: M ) = fliplr ( autocorr ( M +2: end ) ) ; % completa i valori

4 negativi ( funzione pari )

5 P_bt = fft ( autocorr ,2* M +1) ; % DSP

6 f_FT = 0: Fs /(2* M +1) : Fs - Fs /(2* M +1) ; % vettore frequenze

7 P_bt = P_bt (1:( M +1) ) ; % taglio la seconda

8 meta ’

f_FT = f_FT (1:( M +1) ) ;

9

10 plot ( f_FT , abs ( P_bt ) )

11 xlabel ( ’ Frequenze ( Hz ) ’)

12 title ( ’ Densita ’ spettrale di potenza ( metodo indiretto ) ’)

13 Il modulo è usato perché il vettore è complesso a causa dello shift di M della trasfor-

abs()

mata.

2. Metodi Parametrici basati sui modelli ARMA

2 2

|H(ω)|

P̂ (ω) = σ

AR u

Dato un segnale x con frequenza di campionamento F s è possibile calcolare la stima dello

spettro del modello AR(p) identificato sul segnale.

p = 30; % ordine

1 th = ar (x ,p , ’ yw ’) ; % stima AR (p) con Y -W

2 [H , F ] = freqz (1 , th .a ,N , Fs ) ; % risposta in frequenza

3 P_ar = ( abs ( H ) .^2) * th . noisevariance ; % DSP

4

5 plot (F , P_ar )

6 xlabel ( ’ Frequenze ( Hz ) ’)

7 title ( ’ Densita ’ spettrale di potenza del modello AR ( p ) ’)

8 Si può anche creare una funzione che restituisce direttamente lo spettro una volta ottenuti

i parametri del modello:

function [ P_ar , F ] = ar_dsp (a , sigmau2 , NP , Fs )

1 [H , F ] = freqz (1 ,a , NP , Fs ) ;

2 21

P_ar = ( abs ( H ) .^2) * sigmau2 ;

3 end

4 La DSP stimata a partire dal modello AR è più regola