Tema 5: Misura di un segnale costante in rumore
Significato della ddp a priori versus la ddp a posteriori
Modello scalare
dell’osservato:
Yi  S  Wi , i  1,2,
Modello vettoriale
dell’osservato:
y  N  mY , CY 
,N
campioni di rumore
segnale costante aleatorio
Wi  N  0, W2  , IID
S  N  S , S2 
y  [Y1 Y2
y  S1  w
1  [11 1]
W
WN ]T
mY  S 1
Matrice di covarianza degli osservati: CY  E yy T   mY mYT   S2 11T   W2 I
Varianza degli osservati:  Y2  E Yk2   Y2   S2   W2
2
 2

S
T
w  [W1 W2
Vettore valor medio degli osservati:
YN ]T
Segnale costante in rumore: tre realizzzioni
N.B. Vogliamo misurare il valore s
assunto da S nella particolare
realizzazione che viene osservata!!
20
s=12.8
18
16
20
14
y (t )  s  w(t )
16
12
s=9.4
volt
18
14
8
12
volt
10
6
10
4
8
2
6
0
0
20
40
60
80
100
(msec)
120
140
160
180
60
80
100
(msec)
120
140
160
180
200
4
2
20
0
0
20
40
60
80
100
(msec)
120
140
160
180
200
s=7.99
18
16
yi  y (iT )  s  w(iT )
12
,N
volt
i  1, 2,
14
10
8
 2   2  4    1  SNR  0 dB
S
W
S  10
N  10
6
4
2
0
0
20
40
200
File.m: segncost_realiz.m
% Grafico di una realizzazione del segnale osservato
% “tempo-continuo” e tempo-discreto
clc
clear
N=10;
sc=20;
% numero di campioni prelevati
% fattore di sovracampionamento per simulare il "tempo continuo"
etas=10
sig2s=4
sig2w=4
% valor medio del segnale costante
% varianza del segnale costante
% varianza del rumore
s=randn(1,1)*sqrt(sig2s)+etas;
% genera un valore di S con ddp Gaussiana
wtc=randn(N*sc,1)*sqrt(sig2w);
ytc=s+wtc;
% genera il rumore "tempo continuo" con ddp Gaussiana
% calcola il segnale "tempo continuo" totale
y=ytc(1:sc:N*sc);
% estrae gli N campioni osservati
plot(ytc,'r')
xlabel('(msec)')
% traccia il grafico del segnale totale “tempo-continuo“
% scala temporale – esemplificativa %(tempo di campionamento= sc msec)
ylabel('volt')
axis([0 N*sc etas-5*sqrt(sig2s) etas+5*sqrt(sig2s)])
hold on;
plot([0 N*sc],[s s],'k')
% traccia il grafico del livello costante;
hold on;
stem([1:sc:N*sc],y)
% traccia i campioni osservati
Misura di un segnale costante in rumore
Errore della misura
Algoritmo di misura:
Sˆ  f Y1 , Y2 ,
, YN 
  S  Sˆ
Indici di qualità dell’algoritmo: valore
medio e potenza media dell’errore
Significato della ddp a priori:
Ho a disposizione solo fS(s), ovvero
nessun osservato; cosa possiamo
fare?


m  E    E S  Sˆ

2
Sˆ  c    S  c  Sˆ  arg min E  S  c 
c
Sˆ   S
 m  E    E S   S   0

 
P  E  2  E S  2   2
   S
S
 




P  E  2   E S  Sˆ

 Sˆ   S
Significato della ddp
a posteriori:
Avendo a disposizione N
osservati, cosa possiamo fare per
migliorare l’accuratezza della
misura?
2
Misura del segnale costante in rumore
Algoritmo A: media campionaria
Algoritmo B:
minimo errore quadratico medio
N  1 N 
1
ˆ
Sb 
Y

S

i

1   N  N i 1  1   N
N
1
1
Sˆa   Yi  1T y
N i 1
N
 
Valore medio
dell’errore
di misura:
m a  mb  0
Potenza media
statistica
dell’errore
di misura:
P a 
 W2
N
1 T
1
1 y
S
N
1  N
2
N
1
 2 , 


1  N 1 1
N
S
W
 N  W2
P 

  P
1  N N
b
a
Nota bene: entrambi sono combinazione lineare di v.a. Gaussiane,
quindi entrambe sono v.a. Gaussiane
Algoritmo B: calcolo dell’errore quadratico medio

S Y  N S Y , S2 Y
Potenza media
statistica
dell’errore
di misura:
S Y

 Sˆb  E S Y  S Y
1
1
   1T y 
S
N
1  N

P b  E S  Sˆb


2
SY
 
 W2
N
   E  S     E  S    
2
2

 EY ES Y S   S Y

SY

2
2
S ,Y
SY
W

2

E




  P a

Y
SY
N

 
2
2
1
 2 , 
1

1
N
S
Nota bene: P   P  P
b
a
a
W
Misura della matrice di covarianza
 
1
CY  E (y  mY )(y  mY )T    S2 11T   W2 I   S2   W2  
11T 
1 
1  
Misura della matrice di covarianza del vettore degli
osservati, utilizzando K realizzazioni del vettore y
1 K
1 K
T
T
ˆ
ˆ Ym
ˆ Y , dove m
ˆ Y   yi
CY   y i y i  m
K i 1
K i 1
Misura_mat_cov =

I    Y2C

i j
 1

Ci , j   
1  

i j
Coeff. di correlazione
tra Yi e Yj
Valore effettivo:
 S  10,
8.0342 3.9898 4.1356 4.1333
8
2 2  4
3.9898 8.0060 3.9277 4.0460
4 8 4 4
[  1]
N  4, K  5000
4.1356 3.9277 8.1985 4.1425
4
4
8
4
4.1333 4.0460 4.1425 8.0706
4
4
4
8
S
W
4
4
4
File.m: segncost_matcov.m
% Misura della matrice di covarianza degli osservati di segnale costante in rumore
clc
clear
N=4;
% numero di campioni
K=5000; % numero di esperimenti (prove) Monte Carlo
etas=10
sig2s=4
sig2w=4
% valor medio del segnale costante
% varianza del segnale costante
% varianza del rumore
s=randn(1,K)*sqrt(sig2s)+etas;
w=randn(N,K)*sqrt(sig2w);
%
%
%
%
genera ne valori di segnale costante con ddp Gaussiana
genera ne vettori di rumore con ddp Gaussiana,
organizzati in una matrice
calcola ne vettori osservati, organizzati in una matrice
valvettmedio=mean(y.');
y=y-repmat(valvettmedio.',1,K);
valmatcov=y*y.'/ne;
%
%
%
%
calcola il vettor valor medio degli osservati
sottrae dagli osservati il vettore valor medio
calcola la matrice di covarianza degli osservati,
sfruttando la loro organizzazione in matrice
valmatcov
% stampa la valutazione della matrice di covarianza
% degli osservati
y=repmat(s,N,1)+w;
Simulazione stimatori e istogramma degli errori
Wi  N  0, W2  , IID
y  S1  w
S  N  S , S2 
2
N
 2 , 

1  N
S
1
Sˆa  1T y
N
Sˆb   Sˆa 
1
S
1  N
W
P a 
 W2
N
, P b   P a
- Generare K=50000 realizzazioni del vettore osservato Y di N=2 campioni,
per S  10,  2   2  4
[istruzioni utili: randn, repmat]
S
W
- Implementare i due stimatori, calcolare e visualizzare l’istogramma
normalizzato dei loro errori di stima e confrontarli con gli andamenti
teorici (ddp Gaussiane a valor medio nullo e varianza P )
[istruzioni utili: mean, hist, bar]
Istogrammi (normalizzati) degli errori
stimatore a
S  10 ,  2   2  4 (  1)
0.4
S
0.35
N 2
0.3
W
K=50000 realizzazioni
ddp
0.25
0.2
0.15
stimatore b
0.1
0.4
0.05
0.35
0
-6
-4
-2
0
errore
2
4
0.3
6
ddp
0.25
1
Sˆa  1T y
N
0.2
0.15
0.1

1
Sˆb  1T y 
S
N
1  N
0.05
0
-6
-4
-2
0
errore
2
4
6
File.m: stimesegncost_istog.m
% Calcolo istogrammi degli errori degli algoritmi di misura del segnale costante
clc
clear
N=2;
% numero di campioni
K=50000; % numero di esperimenti (prove) Monte Carlo
etas=10
sig2s=4
sig2w=4
% valor medio del segnale costante
% varianza del segnale costante
% varianza del rumore
g=sig2s/sig2w;
% gamma
s=randn(1,K)*sqrt(sig2s)+etas; % genera K valori di segnale costante con ddp Gaussiana
w=randn(N,K)*sqrt(sig2w);
% genera K vettori di rumore con ddp Gaussiana,
% organizzati in una matrice
y=repmat(s,N,1)+w;
% calcola K vettori osservati, organizzati in una matrice
sa=mean(y);
% calcola K misure di tipo a
sb=g*N/(g*N+1)*sa+etas/(g*N+1); % calcola K misure di tipo b
esa=sa-s;
esb=sb-s;
% calcola i K errori delle misure di tipo a
% calcola i K errori delle misure di tipo b
[ista,aa]=hist(esa,50);
[istb,ab]=hist(esb,50);
% calcola l'istogramma degli errori delle misure di tipo a
% calcola l'istogramma degli errori delle misure di tipo b
figure
bar(aa,ista/K/(aa(2)-aa(1)))
xlabel('errore')
ylabel('ddp')
title('a')
figure
bar(ab,istb/K/(ab(2)-ab(1)))
xlabel('errore')
ylabel('ddp')
title('b')
% esegue i grafici degli istogrammi normalizzati
Curve di prestazioni degli algoritmi di misura
P  E  2  versus N
versus 
e
- Ripetere la simulazione per un numero di campioni N variabile da 2 a 50,
con 1000 realizzazioni, per S  10,  2   2  4, quindi valutare e visualizzare
S
W
la potenza media degli errori P di entrambi gli algoritmi al variare di N
[istruzioni utili: randn, repmat, for, mean, std, plot]
- Ripetere la simulazione per    2 /  2 variabile, per S  10,  2  4
S
W
W
 2 [1, 20] ; valutare e visualizzare la potenza media degli errori P di entrambi
S
gli algoritmi al variare di 
[istruzioni utili: randn, repmat, for, mean, std, plot]
Prestazioni degli algoritmi in funzione di N
S  10 ,  2   2  4
S
W
20
stimatore a
stimatore b
18
K=1000 realizzazioni
16
14
medie
12
2
stimatore a
stimatore b
10
1.8
8
1.6
6
1.4
2
0
0
5
10
15
20
25
30
num. campioni
35
40
45
• Il valore medio dell’errore è nullo
per entrambi gli stimatori
50
potenze errore
4
1.2
P versus N
1
0.8
0.6
0.4
0.2
0
0
5
10
15
20
25
30
num. campioni
35
40
45
50
• La precisione della misura migliora con il numero N di campioni osservati
• Il vantaggio dovuto alla conoscenza a priori si riduce all’aumentare di N
Prestazioni degli algoritmi in funzione di 
N  2, S  10 ,  W2  4
20
stimatore a
stimatore b
18
K=1000 realizzazioni
16
   2 / 2
14
S
W
  S2   W2
medie
12
10
8
sigw2=4
2.2
6
4
2
2
1.8
0.5
1
1.5
2
2.5
3
gamma
3.5
4
4.5
5
potenze errore
0
• Il valore medio dell’errore è nullo
per entrambi gli stimatori
stimatore a
stimatore b
1.6
1.4
1.2
P versus 
1
0.8
0.6
• Al crescere di  per
a priori si riduce
0
0.5
1
1.5
2
2.5
gamma
3
3.5
4
4.5
5
 2 ed N fissati, l’utilità delle informazioni
W
File.m: stimesegncost_N.m
% Calcolo della media e della varianza dell’errore dei due algoritmi al variare del numero di campioni N
clc
clear
Nmax=50; % numero massimo di campioni
ne=1000; % numero di esperimenti (prove) Monte Carlo
etas=10 % valor medio del segnale costante
sig2s=4 % varianza del segnale costante
sig2w=4 % varianza del rumore
g=sig2s/sig2w;
% gamma
s=randn(1,ne)*sqrt(sig2s)+etas; % genera ne valori di segnale costante con d.d.p. Gaussiana
w=randn(Nmax,ne)*sqrt(sig2w);
% genera ne vettori di rumore con d.d.p. Gaussiana,
% organizzati in una matrice
y=repmat(s,Nmax,1)+w;
% calcola ne vettori osservati, organizzati in una matrice
for indice=1:Nmax/2,
N=indice*2;
% per numero di campioni N variabile a passi di 2:
sa=mean(y(1:N,:));
% calcola ne stime di tipo a con numero di campioni
sb=g*N/(g*N+1)*sa+etas/(g*N+1); % calcola ne stime di tipo b con numero di campioni N
esa=sa-s;
esb=sb-s;
msa(indice)=mean(sa);
vesa(indice)=mean(esa.^2);
msb(indice)=mean(sb);
vesb(indice)=mean(esb.^2);
end
%
%
%
%
%
%
calcola
calcola
calcola
calcola
calcola
calcola
a=2*[1:Nmax/2];
figure
plot(a,msa,'r*',a,msb,'b*')
xlabel('num. campioni')
ylabel('medie')
legend('stimatore a','stimatore b')
axis([0 Nmax 0 2*etas])
figure
plot(a,vesa,'r*',a,vesb,'b*')
hold on;
plot(a,sig2w./a,'r',a,sig2w*g./(g*a+1),'b')
xlabel('num. campioni')
ylabel('varianze')
legend('stimatore a','stimatore b')
gli ne errori delle stime di tipo a
gli ne errori delle stime di tipo b
media delle stime di tipo a
potenza errore delle stime di tipo a
media delle stime di tipo b
potenza errore delle stime di tipo b
% esegue i grafici
% anche delle varianze teoriche
File.m: stimesegncost_gamma.m
% Calcolo della media e della varianza dell’errore dei due algoritmi al variare del rapporto gamma
clc
clear
N=2;
% numero di campioni
ne=10000; % numero di esperimenti (prove) Monte Carlo
etas=10 % valor medio del segnale costante
sig2si=1 % varianza iniziale del segnale costante
sig2w=4 % varianza del rumore
for indice=1:20,
sig2s=sig2si*indice;
g=sig2s/sig2w;
% per rapporto gamma variabile:
% gamma
s=randn(1,ne)*sqrt(sig2s)+etas; % genera ne valori di segnale costante con d.d.p. Gaussiana
w=randn(N,ne)*sqrt(sig2w);
% genera ne vettori di rumore con d.d.p. Gaussiana,
% organizzati in una matrice
y=repmat(s,N,1)+w;
% calcola ne vettori osservati, organizzati in una matrice
sa=mean(y);
% calcola ne stime di tipo a
sb=g*N/(g*N+1)*sa+etas/(g*N+1); % calcola ne stime di tipo b
esa=sa-s;
esb=sb-s;
msa(indice)=mean(sa);
vesa(indice)=mean(esa.^2);;
msb(indice)=mean(sb);
vesb(indice)=mean(esb.^2);;
end
%
%
%
%
%
%
calcola
calcola
calcola
calcola
calcola
calcola
gli ne errori delle stime di tipo a
gli ne errori delle stime di tipo b
media delle stime di tipo a
potenza errore delle stime di tipo a
media delle stime di tipo b
potenza errore delle stime di tipo b
a=sig2si*[1:20]/sig2w;
% esegue i grafici
figure
plot(a,msa,'r*',a,msb,'b*')
xlabel('gamma')
ylabel('medie')
legend('stimatore a','stimatore b')
axis([min(a) max(a) 0 2*etas])
figure
plot(a,vesa,'r*',a,vesb,'b*')
hold on;
plot([a(1) a(20)],[sig2w/N sig2w/N],'r',a,sig2w*a./(a*N+1),'b') % anche delle varianze teoriche
xlabel('gamma')
ylabel('varianze')
legend('stimatore a','stimatore b')
title('sigw^2=4')