Elaborazione numerica dei segnali -...

39
Appunti di Appunti di Elaborazione numerica dei segnali Elaborazione numerica dei segnali Capitolo 6 - Capitolo 6 - Stima spettrale (parte I) Introduzione ............................................................................................... 1 Valutazione della densità spettrale di energia ............................................. 3 Stima della densità spettrale di potenza per segnali aleatori tempo-continui7 Uso del banco di filtri ................................................................................. 9 Risoluzione in frequenza del sistema di misura ................................... 10 Dispersione della stima spettrale ........................................................ 11 Richiami di statistica .......................................................................... 11 Stima della potenza di segnali correlati gaussiani .............................. 15 Campioni correlati ......................................................................... 17 Metodo del periodogramma ...................................................................... 19 Implementazione del periodogramma tramite la DFT ............................... 22 Osservazione: allungamento con zeri della x(n) .................................. 24 Esempio ........................................................................................ 25 Metodo del periodogramma, DFT e banco di filtri .............................. 26 Allungamento della finestra di osservazione ............................................. 29 Polarizzazione della stima spettrale .......................................................... 30 Modifica del periodogramma .................................................................... 33 Media di periodogrammi di sottosequenze: metodo di Bartlett ........... 34 Metodo WOSA (finestre a coseno rialzato sovrapposte) ...................... 35 Smussamento dello spettro: metodo di Blackman-Tukey ..................... 36 Riepilogo e conclusioni ............................................................................ 38 I NTRODUZIONE Affrontiamo il problema dell’analisi spettrale dei segnali e cioè la stima del loro spettro di potenza. Il problema della stima spettrale si può introdurre pensando allo studio dello spettro delle radiofrequenze in una grande città: l’uno accanto all’altro si incontrano i vari canali radio e televisivi, i segnali di telefonia cellulare, i segnali radar e così via. Gli intervalli di frequenza, se considerati con sufficiente larghezza di banda, sono occupati da segnali totalmente indipendenti tra loro, con diversa ampiezza e ciascuno con diverse caratteristiche statistiche. Per analizzare uno spettro di radiofrequenze fatto in questo modo, appare logico utilizzare uno schema di misura del tipo seguente:

Transcript of Elaborazione numerica dei segnali -...

Page 1: Elaborazione numerica dei segnali - users.libero.itusers.libero.it/sandry/download/DSPdownload/DSP_06a.pdf · una frequenza di campionamento f S che rispetti il teorema di campionamento

Appunti diAppunti diElaborazione numerica dei segnaliElaborazione numerica dei segnaliCapitolo 6 - Capitolo 6 - Stima spettrale (parte I)

Introduzione ............................................................................................... 1Valutazione della densità spettrale di energia ............................................. 3Stima della densità spettrale di potenza per segnali aleatori tempo-continui7Uso del banco di filtri ................................................................................. 9

Risoluzione in frequenza del sistema di misura ................................... 10Dispersione della stima spettrale ........................................................ 11Richiami di statistica .......................................................................... 11Stima della potenza di segnali correlati gaussiani .............................. 15

Campioni correlati......................................................................... 17Metodo del periodogramma ...................................................................... 19Implementazione del periodogramma tramite la DFT ............................... 22

Osservazione: allungamento con zeri della x(n).................................. 24Esempio ........................................................................................ 25

Metodo del periodogramma, DFT e banco di filtri.............................. 26Allungamento della finestra di osservazione............................................. 29Polarizzazione della stima spettrale .......................................................... 30Modifica del periodogramma.................................................................... 33

Media di periodogrammi di sottosequenze: metodo di Bartlett ........... 34Metodo WOSA (finestre a coseno rialzato sovrapposte) ...................... 35Smussamento dello spettro: metodo di Blackman-Tukey ..................... 36

Riepilogo e conclusioni ............................................................................ 38

INTRODUZIONE

Affrontiamo il problema dell’analisi spettrale dei segnali e cioè la stima del loro spettro dipotenza.

Il problema della stima spettrale si può introdurre pensando allo studio dello spettro delleradiofrequenze in una grande città: l’uno accanto all’altro si incontrano i vari canali radio etelevisivi, i segnali di telefonia cellulare, i segnali radar e così via. Gli intervalli di frequenza, seconsiderati con sufficiente larghezza di banda, sono occupati da segnali totalmente indipendenti traloro, con diversa ampiezza e ciascuno con diverse caratteristiche statistiche. Per analizzare unospettro di radiofrequenze fatto in questo modo, appare logico utilizzare uno schema di misura deltipo seguente:

Page 2: Elaborazione numerica dei segnali - users.libero.itusers.libero.it/sandry/download/DSPdownload/DSP_06a.pdf · una frequenza di campionamento f S che rispetti il teorema di campionamento

Appunti di “Elaborazione numerica dei segnali” - Capitolo 6

Autore: Sandro Petrizzelli2

Si tratta cioè di utilizzare un ricevitore selettivo, composto da un banco di filtri: ciascun filtroisola il contenuto spettrale in un assegnato intervallo di frequenza. L’uscita di ciascun filtro va in unmisuratore di potenza, che dà una misura della potenza allocata nell’intervallo considerato1.L’insieme delle diverse misure costituisce la misura della densitàspettrale di potenza su tutto lo spettro di frequenze preso inesame.

Il nostro scopo è quello di fornire una misura più precisa possibile di tali potenze,compatibilmente con due vincoli imprescindibili:

• in primo luogo, non abbiamo a disposizione un tempo infinito di osservazione, ma solo untempo di durata finita, che non potrà nemmeno essere troppo lunga2;

• in secondo luogo, in questo contesto ci interessa realizzare un sistema numerico (il cosiddettoanalizzatore di spettro numerico), per cui l’ingresso analogico dovrà essere convertito innumerico (tramite un campionamento a frequenza opportuna) prima di poter essere elaborato.

In altre parole, quindi, possiamo compiere le nostre misure solo su segnali campionati in unintervallo di tempo di durata finita, il che rappresenta il principale limite ai risultati che potremoottenere: potremo ottenere solo delle stime di ciò che stiamo cercando. Lo scopo è ovviamentequello di ottenere delle stime quanto più accurate possibile.

E’ importante notare sin da ora che non facciamo alcuna ipotesi ulteriore3 sui segnali in ingresso.Si parla allora di analisi spettrale non parametrica. Questo tipo di analisi si contrappone allacosiddetta analisi spettrale parametrica, che parte in modo diverso: ad esempio, seconsideriamo lo spettro di un segnale acustico dovuto a vibrazioni oppure ad un segnale vocale,possiamo star certi che i contenuti delle varie bande spettrali siano in qualche modo correlati tra loro,dato appunto che il segnale, nel suo complesso, ha origine uniche. Sulla base di questa

1 La scelta di un banco di filtri si rivela utile, come si vedrà, solo nel caso di un sistema numerico. Al contrario, nel caso di un

sistema di misura analogico, il metodo più praticabile è quello di un ricevitore a supereterodina, nel quale un oscillatore locale, afrequenza variabile, separa i vari intervalli di frequenza portandoli all’interno di una banda prefissata (quella di un filtro selettoreposto a valle del moltiplicatore). In un sistema di questo tipo, uno dei parametri critici è la velocità di spazzolamento dellefrequenze, cioè la velocità con cui l’oscillatore locale varia la propria frequenza di oscillazione in modo da isolare intervalli viavia diversi di frequenza.

2 Come vedremo, la lunghezza finita dei dati a nostra disposizione è la principale limitazione alla qualità della stima della densitàspettrale.

3 E’ ovvio che supponiamo un campionamento effettuato durante un tempo finito e che rispetti il teorema del campionamento, inmodo da evitare qualsiasi problema di aliasing

Page 3: Elaborazione numerica dei segnali - users.libero.itusers.libero.it/sandry/download/DSPdownload/DSP_06a.pdf · una frequenza di campionamento f S che rispetti il teorema di campionamento

Stima spettrale (parte I) - Metodi non parametrici

Autore: Sandro Petrizzelli3

considerazione, l’analisi spettrale parametrica parte dal fissare un modello per la generazione delsegnale in esame e si occupa successivamente di determinare i parametri del modello stesso (da cuiappunto il nome di analisi parametrica)4.

E’ ovvio che i segnali con cui abbiamo a che fare non sono noti a priori, per cui non sonodeterministici, ma aleatori. Date le fluttuazioni casuali in tali segnali, dobbiamo necessariamenteadottare una analisi di tipo statistico, che si basi cioè sulle caratteristiche medie dei segnali. Inparticolare, la funzione di autocorrelazione di un processo casuale è il parametro statisticomedio che useremo per caratterizzare i segnali aleatori nel dominio del tempo, mentre la trasformatadi Fourier di tale funzione, vale a dire la densità spettrale di potenza, fornisce il passaggio daldominio del tempo al dominio della frequenza.

In effetti, l’analisi spettrale parametrica diventa una questione concettualmente molto semplice sesi fa uso del concetto di media di insieme: infatti, se noi abbiamo a disposizione un buon numero direalizzazioni del segnale (ad esempio, nel caso del segnale telefonico, un buon numero diconversazioni), ci basta calcolarne le trasformate discrete di Fourier e mediarne le potenze,componente per componente. Al contrario, raramente si hanno a disposizione numerose repliche delsegnale; spesso, invece, si ha a disposizione una sola replica, per un intervallo di tempo assegnato.Per risalire, allora, allo spettro di potenza, non potremo prescindere da ipotesi come stazionarietàed ergodicità5, che supporremo quindi sempre valide.

VALUTAZIONE DELLA DENSITÀ SPETTRALE DI ENERGIA

Per inquadrare il problema che intendiamo affrontare, cominciamo a ragionare con segnali di tipodeterministico. Supponiamo, in particolare, di disporre di un vettore x(n) che sia frutto delcampionamento, in un intervallo di tempo di durata finita, di un segnale tempo-continuo xa(t). Seavessimo a disposizione il segnale xa(t) per intero, non avremmo problemi a calcolare il suo spettrodi energia. Al contrario, il fatto di disporre solo dei campioni di xa(t) in un intervallo di durata finitaci dice che possiamo solo avvicinarci allo spettro di energia effettiva, ossia possiamo farne solo unastima più o meno precisa.

Per prima cosa, ricordiamo che un segnale tempo-continuo xa(t) ha una energia definita dallanota relazione

∫+∞

∞−

= dt)t(xE2

a

Se questa quantità risulta finita, allora il segnale si dice ad energia finita ed ammette unatrasformata di Fourier, che indichiamo con Xa(f). In base al noto teorema di Parseval, l’energia ditale trasformata coincide con quella di xa(t):

∫∫+∞

∞−

+∞

∞−

== df)f(Xdt)t(xE2

a

2

a

4 Tipico esempio, che sarà ripreso più avanti, è il seguente: una volta scelto uno spettro di potenza per il segnale sotto analisi, si fa

in modo da ricavarlo come uscita di un filtro lineare che sia alimentato in ingresso da un rumore bianco (cioè a spettro uniforme);si tratta allora di determinare i parametri del filtro in modo che lo spettro di potenza in uscita sia quanto più simile possibile allospettro sotto analisi.

5 Ricordiamo che un processo, per essere ergodico, deve essere stazionarietà, mentre invece non è detto il viceversa. Tutti i processiergodici sono stazionari, mentre un processo stazionario non necessariamente è ergodico.

Page 4: Elaborazione numerica dei segnali - users.libero.itusers.libero.it/sandry/download/DSPdownload/DSP_06a.pdf · una frequenza di campionamento f S che rispetti il teorema di campionamento

Appunti di “Elaborazione numerica dei segnali” - Capitolo 6

Autore: Sandro Petrizzelli4

In quest’ultima relazione, si è scritto, in pratica, che l’energia di Xa(f) o di xa(t) si calcola come

area sottesa dalla funzione 2

a )f(X , la quale quindi viene a rappresentare una densità spettrale di

energia (o semplicemente spettro di energia). La indicheremo nel seguito come 2

aXX )f(X)f(S = .

Questa funzione corrisponde inoltre alla trasformata di Fourier di un’altra funzione RXX(τ), cheprende il nome di funzione di autocorrelazione del segnale xa(t) e che è definita nel modoseguente:

→τ+=τ ∫+∞

∞−

dt)t(x)t(x)(R a*aXX ∫

+∞

∞−

τπ− ττ= de)(R)f(S f2jXXXX

Fatti questi richiami, torniamo al nostro problema. Supponiamo di campionare il segnale xa(t) aduna frequenza di campionamento fS che rispetti il teorema di campionamento, in modo da evitarefenomeni di aliasing nel dominio della frequenza. Sia x(n), per -∞<n<∞, la sequenza di campioniottenuta. La trasformata di Fourier di x(n) è notoriamente

∑∑+∞

−∞=

+∞

−∞=

π− −==n

SaSn

fn2j )kff(Xfe)n(x)f(X

Si tratta della ripetizione periodica, a passo fS, dello spettro di xa(t), a meno di un fattore di scalafS. In assenza di aliasing, come stiamo supponendo, X(f) coincide con Xa(f) nel periodofondamentale, ossia per |f|<fS/2.

Calcolando il modulo quadro di X(f), all’interno del periodo fondamentale6, otteniamo la densitàspettrale del segnale campionato:

2

ff )f(Xf)f(X S2

a2

S

2 <=

Quindi, se evitiamo fenomeni di aliasing in frequenza, possiamoricostruire la densità spettrale di energia del segnale xa(t) apartire dalla densità spettrale di energia dei suoi campioni:

2

ff )f(X

f

1)f(X S2

2S

2

a <=

Questa possibilità deriva solo dall’aver idealmente campionato xa(t) da -∞ a +∞: il fatto didisporre di infiniti campioni ci garantisce, in base al teorema del campionamento, che possiamorecuperare, tramite i campioni, tutto il contenuto informativo del segnale analogico di partenza.

Potevamo procedere anche in altro modo, lavorando nel dominio del tempo. Infatti, data lasequenza x(n) (di lunghezza infinita) rappresentativa del segnale xa(t), la sua funzione diautocorrelazione è definita come

∑+∞

−∞=

+=n

*XX )kn(x)n(x)k(r

e la trasformata di Fourier di questa funzione è

∑+∞

−∞=

π−=k

kf2jXX

2e)k(r)f(X

6 Non avrebbe senso fare il calcolo su tutto l’asse delle frequenze, dato che lo spettro si ripete periodicamente a passo fS.

Page 5: Elaborazione numerica dei segnali - users.libero.itusers.libero.it/sandry/download/DSPdownload/DSP_06a.pdf · una frequenza di campionamento f S che rispetti il teorema di campionamento

Stima spettrale (parte I) - Metodi non parametrici

Autore: Sandro Petrizzelli5

Quindi, partendo da x(n), ci calcoliamo rXX(k), poi trasformiamo per ricavare 2

)f(X ed infine

applichiamo la relazione 2

2S

2

a )f(Xf

1)f(X = per risalire a

2

a )f(X .

Quindi, mentre prima, lavorando nel dominio della frequenza (metodo diretto), si trattavasemplicemente di calcolare la trasformata di Fourier della sequenza x(n) e poi di calcolarne ilmodulo quadro, adesso, lavorando nel dominio del tempo (metodo indiretto), si tratta di calcolarela funzione di autocorrelazione di x(n), poi di trasformarla ed infine di dividerla per 2

Sf .

In entrambi i casi, comunque, si giunge a valutare proprio 2

a )f(X , il che è dovuto essenzialmente

alla possibilità di disporre di infiniti campioni di xa(t) (nonché anche all’aver evitato aliasing infrequenza).

Nella pratica, però, noi possiamo disporre solo di sequenze x(n) di durata finita. Vediamo allora, alivello essenzialmente qualitativo, come cambiano le cose.

Sia N la durata di x(n), il che significa che x(n) consta di N campioni numerati da 0 ad N-1. Ilfatto di disporre di un numero finito di campioni equivale a considerare la sequenza x(n), di duratainfinita, moltiplicata per una finestra rettangolare (sempre discreta) di durata N:

≤≤

==altrimenti 0

1-Nn0 )n(x)n(w)n(x)n(x~

Lo spettro della sequenza finestrata )n(x~ è la convoluzione degli spettri di x(n) e w(n):

∫−

−==2/f

2/f

S

S

da)af(W)a(X)f(W*)f(X)f(X~

Lo spettro di x(n) è la ripetizione periodica dello spettro di xa(t), mentre lo spettro di w(n) è laripetizione periodica dello spettro di un rettangolo esteso da t=0 a t=NT; si tratta cioè dellaperiodicizzazione di un sin(f)/f. Si può agevolmente ricavare che l’espressione analitica dello spettrodella finestra rettangolare campionata è

( )( )fsin

fNsine)f(W 2

1Nf2j

ππ

=−

π−

dove il fattore di fase deriva semplicemente dal fatto che il rettangolo non è centrato in 0.Per capire quale sia l’effetto della convoluzione di X(f) con W(f), facciamo riferimento ad un

esempio concreto.Per quanto riguarda X(f), consideriamo il caso semplice in cui tale spettro è costante, sul valore

unitario, per |f|<0.1fS e 0 altrove:

=altrimenti 0

f1.0f 1)f(X S

Abbiamo dunque considerato una frequenza di campionamento pari a 10 volte la massimafrequenza del segnale in esame, il che ci garantisce il rispetto del teorema del campionamento.

Andiamo a convolvere X(f) per una finestra rettangolare di lunghezza N=61. Indicato sempre conW(f) lo spettro della finestra rettangolare campionata, il suo modulo risulta essere del tipo seguente:

Page 6: Elaborazione numerica dei segnali - users.libero.itusers.libero.it/sandry/download/DSPdownload/DSP_06a.pdf · una frequenza di campionamento f S che rispetti il teorema di campionamento

Appunti di “Elaborazione numerica dei segnali” - Capitolo 6

Autore: Sandro Petrizzelli6

Modulo (in dB) dello spettro della finestra rettangolare campionata, diagrammato sia per N=M=31campioni sia per N=M=61 campioni. Sull’asse delle frequenze è riportata la frequenza normalizzata,

ossia il rapporto tra la frequenza effettiva e la frequenza di campionamento fS

Convolvendo, in base alla relazione )f(W*)f(X)f(X~ = , lo spettro appena riportato con X(f), si

ottiene uno spettro )f(X~

la cui densità di energia è del tipo riportato nella figura seguente:

Si nota immediatamente che l’energia ha subito una dispersione nella banda di frequenzacompresa tra 0.1fS e 0.5fS, dove invece lo spettro del nostro segnale è nullo. Questa dispersioneprende il nome di leakage ed è dovuta principalmente ai lobi secondari dello spettro W(f). Perquanto riguarda, invece, il lobo principale di W(f), esso contribuisce ad attenuare lo spettro X(f)all’interno della banda in cui è definito.

Così come si è visto a proposito del progetto dei filtri FIR, possiamo ridurre il leakage dovuto ailobi secondari facendo in modo che W(f) abbia lobi secondari meno pronunciati. Per ottenere questo,dobbiamo usare una finestra che non sia più rettangolare, cioè con una discontinuità brusca alla fine,ma presenti una transizione più dolce dalla banda passante alla banda arrestata. Riferendoci ancoraall’esempio appena illustrato, se usiamo una finestra di Blackman, ancora con N=61, l’esito della

convoluzione risulta essere uno spettro )f(X~

la cui densità di energia è del tipo seguente:

Page 7: Elaborazione numerica dei segnali - users.libero.itusers.libero.it/sandry/download/DSPdownload/DSP_06a.pdf · una frequenza di campionamento f S che rispetti il teorema di campionamento

Stima spettrale (parte I) - Metodi non parametrici

Autore: Sandro Petrizzelli7

Il leakage dovuto ai lobi laterali è sicuramente diminuito.Da questi discorsi, più che altro qualitativi, emerge dunque che lo spettro di energia della

sequenza )n(w)n(x)n(x~ = finestrata è solo una approssimazione dello spettro di energia di x(n).Analiticamente, la densità spettrale di energia di )n(x~ è

21N

0n

fn2j2

XX e)n(x~)f(X~

)f(S~ ∑

=

π−==

Questo spettro di energia può anche essere valutato numericamente tramite la DFT, ossiaconsiderando un dominio campionato delle frequenze: infatti, a partire da )n(x~ , possiamo calcolarnela DFT, ossia il vettore

∑−

=

π−=

1N

0n

nN

k2j

e)n(x~)k(X~

e successivamente possiamo calcolare il modulo quadro degli elementi di tale vettore, ottenendo una

versione campionata di )f(S~

XX :

== ∑

=

π−

N

kS~

e)n(x~)k(X~

XX

21N

0n

nN

k2j2

Il vettore

N

kS~

XX così ottenuto è una versione campionata distorta del reale spettro di energia

( )fSXX , ossia una versione distorta di

N

kSXX (=reale spettro di energia sottoposto a

campionamento).

STIMA DELLA DENSITÀ SPETTRALE DI POTENZA PER SEGNALI ALEATORI

TEMPO-CONTINUI

I segnali ad energia finita considerati nel precedente paragrafo hanno tutti una trasformata diFourier (proprio perché sono ad energia finita) e sono caratterizzati, in frequenza, dal loro spettro dienergia (o densità spettrale di energia). Al contrario, i segnali appartenenti alla famiglia dei

Page 8: Elaborazione numerica dei segnali - users.libero.itusers.libero.it/sandry/download/DSPdownload/DSP_06a.pdf · una frequenza di campionamento f S che rispetti il teorema di campionamento

Appunti di “Elaborazione numerica dei segnali” - Capitolo 6

Autore: Sandro Petrizzelli8

processi casuali stazionari non hanno energia finita e quindi non possiedono una trasformata diFourier. Questi segnali hanno una potenza media finita, per cui sono caratterizzabili solo tramite unadensità spettrale di potenza. Facciamo allora anche qui qualche richiamo.

Se x(t) è un generico processo casuale stazionario, la sua funzione di autocorrelazione èdefinita come

[ ])t(x)t(xE)( *XX τ+=τγ

Si tratta cioè della media (o aspettazione) statistica del prodotto di due processi casuali7.Trasformando secondo Fourier la funzione γXX(τ), si ottiene la densità spettrale di potenza

del processo considerato:

∫+∞

∞−

τπ− ττγ=Γ de)()f( f2jXXXX

Come già detto in precedenza, noi disponiamo di una sola realizzazione del processo stocastico inesame e ci stiamo ponendo il problema di risalire, da tale realizzazione, allo spettro di potenza delprocesso. Il fatto di avere un’unica realizzazione ci dice che non possiamo conoscere la verafunzione di autocorrelazione γXX(τ) del processo e quindi, come conseguenza, non possiamonemmeno calcolare )f(XXΓ , cioè la vera densità spettrale di potenza del processo.

Al contrario, basandoci proprio sulla singola realizzazione, possiamo calcolare la cosiddettafunzione di autocorrelazione media temporale (time-average autocorrelation function), che ècosì definita:

∫−

τ+=τ2/T

2/T

*

0XX

0

0

dt)t(x)t(xT

1)(R

dove T0 è il tempo di osservazione, cioè la durata finestra di osservazione della realizzazione.Nell’ipotesi che il processo, oltre che stazionario, sia ergodico sui momenti del primo ordine (la

media) e del secondo ordine (la funzione di autocorrelazione), sappiamo di poter scrivere che

∫−

∞→∞→τ+=τ=τγ

2/T

2/T

*

0T

XXT

XX

0

000

dt)t(x)t(xT

1lim)(Rlim)(

Questa relazione dice in pratica che, facendo tendere ad ∞ il tempo di osservazione, la funzioneRXX(τ) tende a coincidere con la reale funzione di autocorrelazione γXX(τ) del processo. In altreparole, per T0→∞, la funzione di autocorrelazione media temporale tendea coincide con la funzione di autocorrelazione media statistica8.Quindi, la funzione RXX(τ) è uno stimatore della funzione γXX(τ).

Possiamo ripetere un identico discorso nel dominio della frequenza. Infatti, data RXX(τ), possiamofarne la trasformata di Fourier:

∫−

τπ− ττ=2/T

2/T

f2jXXXX

0

0

de)(R)f(P

7 La stazionarietà del processo (che d’ora in poi sarà sempre sottintesa) garantisce che la funzione di autocorrelazione sia funzione

solo dell’ampiezza τ dell’intervallo (t,t+τ) prescelto e non degli estremi dell’intervallo stesso

8 E’ qui che interviene l’ipotesi di ergodicità del processo: essa garantisce che la media temporale tenda, con probabilità 1, allamedia statistica (o media di insieme).

Page 9: Elaborazione numerica dei segnali - users.libero.itusers.libero.it/sandry/download/DSPdownload/DSP_06a.pdf · una frequenza di campionamento f S che rispetti il teorema di campionamento

Stima spettrale (parte I) - Metodi non parametrici

Autore: Sandro Petrizzelli9

Sostituendo l’espressione di RXX(τ), otteniamo che

22/T

2/T

f2j

0

2/T

2/T

f2j2/T

2/T

*

0

2/T

2/T

f2j2/T

2/T

*

0XX

0

0

0

0

0

0

0

0

0

0

dte)t(xT

1

...dedt)t(x)t(xT

1dedt)t(x)t(x

T

1)f(P

∫ ∫∫ ∫

τπ−

τπ−

−−

τπ−

=

=ττ+=τ

τ+=

In base a quanto ottenuto, la densità spettrale di potenza del processo sarà la media di PXX(f),ottenuta per T0→∞:

[ ]

==Γ ∫

τπ−

∞→∞→

22/T

2/T

f2j

0T

XXT

XX

0

000

dte)t(xET

1lim)f(PElim)f(

Abbiamo dunque trovato, ancora una volta, due distinti modi di procedere:

• con il metodo diretto, cioè quello appena esposto, calcoliamo PXX(f), ne facciamo la media epoi calcoliamo il limite per T0→∞;

• con il metodo indiretto, invece, calcoliamo prima RXX(τ), ne facciamo il limite per T0→∞ inmodo da risalire a γXX(τ) ed infine trasformiamo per calcolare )(XX τΓ

Appare evidente, d’altra parte, che siamo ancora ai margini del problema che ci siamo posti, inquanto abbiamo detto all’inizio che intendiamo effettuare le nostre stime non su una realizzazioneanalogica, sia pure osservata per un tempo finito, ma su una sequenza x(n) ottenuta campionando lasuddetta realizzazione. Dobbiamo quindi rianalizzare il problema partendo da una sequenza x(n) didurata finita. E’ quello che faremo nei prossimi paragrafi.

USO DEL BANCO DI FILTRI

Volendo stimare lo spettro di potenza di una sequenza x(n) di lunghezza N, senza fare alcun altraipotesi sulla sequenza stessa (analisi non parametrica), dobbiamo necessariamente adottare unoschema di misura del tipo seguente:

x(n)n=0,1,...N-1

H1(f)

H2(f)

HK(f)

y1(n)

y2(n)

yK(n)

Page 10: Elaborazione numerica dei segnali - users.libero.itusers.libero.it/sandry/download/DSPdownload/DSP_06a.pdf · una frequenza di campionamento f S che rispetti il teorema di campionamento

Appunti di “Elaborazione numerica dei segnali” - Capitolo 6

Autore: Sandro Petrizzelli10

Abbiamo cioè un banco di K filtri passa-banda, ciascuno con una propria banda Bk centrata su unapropria frequenza centrale fk. Il segnale9 passa attraverso ognuno di questi filtri, dando origine a Ksequenze filtrate y1(n), y2(n),...,yK(n).

Ogni sequenza yk(n) rappresenta il contenuto spettrale della sequenza di partenza x(n) all’internodella banda Bk del corrispondente filtro Hk(f). Queste sequenze vanno a loro volta in ingresso a deimisuratori di potenza, i quali quindi forniscono una stima della potenza che x(n) ha allocata nellevarie bande:

Hk(f)

Bk

fk

yk(n)realizzazionedel processo

x(n) stima dellapotenza in Bk

Posta la questione in questo modo, si presentano una serie di problemi, che analizziamo neiseguenti paragrafi.

Risoluzione in frequenza del sistema di misura

Il primo problema cui va incontro lo schema di misura dell’ultima figura (cioè l’analizzatore dispettro numerico) riguarda la cosiddetta risoluzione in frequenza.

Supponiamo, infatti, che il segnale in ingresso sia costituito solamente da due sinusoidi, afrequenza f1 ed f2. Per risoluzione in frequenza (o risoluzione spettrale)di un qualsiasisistema di misura si intende il valore minimo della distanza ∆f=f2-f1 alla quale il sistema riesceancora a interpretare i due termini spettrali come distinti.

E’ ovvio che a noi interessa avere una risoluzione in frequenza il piùbassa possibile, cioè interessa avere una minima distanza ∆f alla quale due componentispettrali distinte possono essere effettivamente separate.

Il valore di ∆f, cioè appunto della risoluzione in frequenza, dipende dalla larghezza di banda deifiltri utilizzati; essendo tale larghezza di banda inversamente proporzionale alla durata della rispostaall’impulso, deduciamo che, per avere buona risoluzione in frequenza,dobbiamo avere filtri stretti, ossia con funzioni di rispostaall’impulso lunghe. Vedremo tra un attimo che questo requisito è in contrasto con altrirequisiti, per cui dovremo richiedere il solito compromesso tra opposte esigenze.

9 Val la pena ricordare che x(n), essendo la versione campionata di un certo x(t) analogico, ha uno spettro che è ancora continuo ed

è una ripetizione periodica dello spettro di x(t). Se è stato rispettato il teorema del campionamento, tale ripetizione periodica èesente da aliasing.

Page 11: Elaborazione numerica dei segnali - users.libero.itusers.libero.it/sandry/download/DSPdownload/DSP_06a.pdf · una frequenza di campionamento f S che rispetti il teorema di campionamento

Stima spettrale (parte I) - Metodi non parametrici

Autore: Sandro Petrizzelli11

Dispersione della stima spettrale

Il secondo problema che emerge nello schema di misura prima descritto è il seguente: all’uscitadai vari filtri noi otteniamo delle sequenze yk(n), i cui campioni sono delle variabili casualirappresentative dell’unica realizzazione, del processo in esame, che abbiamo a nostra disposizione;di queste sequenze noi vogliamo calcolare la potenza, ossia il valore quadratico medio (o varianza),che indicheremo con 2

kσ :

[ ])n(yVar k2k =σ

In effetti, però, avendo a che fare con una sola realizzazione del processo, non siamo in grado dicalcolare 2

kσ , ma ci dobbiamo limitare ad una sua stima (che quindi chiameremo stima spettrale).

Hk(f)

Bk

fk

yk(n)realizzazionedel processo

x(n) stima dellapotenza 2

Per effettuare tale stima, usiamo il classico medio usato in statistica:

• consideriamo una variabile aleatoria kS che abbia, come proprietà fondamentale, quella per cui

il suo valore medio è pari proprio a 2kσ ; in questo senso, kS è uno stimatore della quantità

2kσ ;

• questa variabile aleatoria kS diventa un numero, che sarà appunto la nostra stima di 2kσ , nel

momento in cui la calcoliamo in corrispondenza dei valori della realizzazione che abbiamoosservato.

Richiami di statisticaFacciamo alcuni richiami di statistica per comprendere meglio quanto appena detto.

Supponiamo di osservare un certo fenomeno, quale può essere il lancio di un dado o di unamoneta o l’estrazione di palline da un’urna. Supponiamo che vengano effettuate varierealizzazione di tale fenomeno (per esempio vari lanci del dado o varie estrazioni dall’urna). Peresaminare tali realizzazioni, associamo a ciascuna di essere una variabile aleatoria che ne indichiil risultato: per n realizzazioni avremo perciò le variabili aleatorie X1,X2,....,Xn. Ad esempio, nelcaso del lancio di un dado, la variabile aleatoria X1 corrisponde all’esito del primo lancio, lavariabile aleatoria X2 all’esito del secondo lancio e così via fino all’ n-simo lancio.

Dato che riguardano realizzazioni di uno stesso fenomeno, queste variabili aleatorie avrannoragionevolmente la stessa distribuzione. Se allora diciamo che una generica variabile aleatoria Xha la stessa distribuzione delle variabili X1,X2,....,Xn, in pratica stiamo dicendo che Xrappresenta il modello probabilistico che usiamo per la descrizione del fenomeno.

Un’altra condizione importante è che la le variabili X1,X2,....,Xn siano tra loro indipendenti: adesempio, nel caso dell’estrazione di una pallina da un’urna, l’indipendenza è verificata solo se,dopo ogni estrazione, riponiamo nell’urna la pallina estratta.

Page 12: Elaborazione numerica dei segnali - users.libero.itusers.libero.it/sandry/download/DSPdownload/DSP_06a.pdf · una frequenza di campionamento f S che rispetti il teorema di campionamento

Appunti di “Elaborazione numerica dei segnali” - Capitolo 6

Autore: Sandro Petrizzelli12

Nell’ipotesi fondamentale che le variabili X1,X2,...,Xn sianoindipendenti, noi diremo che esse costituiscono un campionealeatorio per la variabile aleatoria X.Una cosa è subito evidente: le variabili aleatorie X1,X2,....,Xn sono i modelli probabilistici chedescrivono le n distinte (ed indipendenti) realizzazioni del fenomeno; come tali, esse sono unqualcosa di teorico; al contrario, l’aver supposto di aver realmente osservato le n realizzazionidel fenomeno, significa che per ciascuna variabile è stato osservato un preciso valore10: indicatiallora con x1,x2,...,xn i valori osservati nelle n prove, diremo che essi costituiscono uncampione per X.

La differenza tra un campione aleatorio di X ed un campione di X è dunque evidente: ilcampione aleatorio è un insieme di variabili aleatorie, cioè sostanzialmente di modellimatematici, mentre il campione è un insieme di numeri, corrispondenti ai valori realmenteassunti dalle variabili aleatorie.

Può essere utile un esempio chiarificatore. Supponiamo che il fenomeno in esame sia il lancio diuna moneta: per ogni lancio sappiamo di poter utilizzare una variabile aleatoria di Bernoulli Bi

che assuma valore 1 in caso esca testa e valore 0 in caso esca croce. Supponiamo di effettuare 3lanci. Il modello probabilistico del fenomeno è dunque una variabile X con distribuzione diBernoulli e parametro p=0.5; un campione aleatorio è costituito da 3 variabili aleatorie(B1,B2,B3), le quali descrivono teoricamente l’esito delle 3 prove, mentre invece un campionepotrebbe essere costituito da (1,0,1) nel caso siano uscite testa, croce, testa, oppure (1,1,0) otutte le altre possibili combinazioni.

Il concetto fondamentale della statistica può essere così riassunto: si vuole sfruttare leinformazioni a disposizione su una certa variabile aleatoria X, ossia la conoscenza di uncampione aleatorio e di un eventuale campione, per stimare uno o più parametri legati alladistribuzione di probabilità di X.

A tale scopo, si introduce il fondamentale concetto di stimatore: qualitativamente, unostimatore è una variabile aleatoria che ci da una idea, più o menoapprossimata, di quanto valga il parametro che andiamo ricercando.Il problema che ci si pone è allora duplice: da un lato, individuare le caratteristiche che deveavere uno stimatore per essere un BUON stimatore del parametro in questione e, in secondoluogo, stabilire un criterio in base al quale scegliere, tra più stimatori disponibili, quellomigliore.

Diamo una definizione più rigorosa di stimatore per un parametro: sia X una variabile aleatoriacon una certa distribuzione di probabilità dipendente da un parametro θ incognito; sia (X1,..,Xn)un campione aleatorio di lunghezza n per X e sia (x1,...,xn) l’insieme dei corrispondenti valoriosservati. Se g(X1,...,Xn) è una funzione, dipendente dal campione aleatorio, usata per stimare θ,

noi diremo che g è uno stimatore del parametro θ. Porremo in generale )X,..,X,X(gˆn21=θ .

Il valore che θ (cioè la funzione g) assume in corrispondenza del campione (x1,...,xn) si diràinvece stima del parametro θ. In altre parole, mentre g(X1,...,Xn) è una variabile aleatoria,g(x1,...,xn) è il particolare valore che tale variabile aleatoria ha assunto per la realizzazioneconsiderata.

Uno stimatore θ si dice corretto (o non polarizzato) se la sua media è esattamente pari alparametro che lo stimatore deve stimare. In caso contrario, si parla di stimatore polarizzato.

10 Tale valore, ovviamente, sarà tra quelli appartenenti al rango della variabile aleatoria. Ad esempio, nel caso del lancio di una

moneta, il rango della corrispondente variabile aleatoria sarà (testa,croce), dato che sono gli unici risultati osservabili.

Page 13: Elaborazione numerica dei segnali - users.libero.itusers.libero.it/sandry/download/DSPdownload/DSP_06a.pdf · una frequenza di campionamento f S che rispetti il teorema di campionamento

Stima spettrale (parte I) - Metodi non parametrici

Autore: Sandro Petrizzelli13

Uno stimatore è detto consistente se e solo se è verificata la condizione

( ) 1ˆPlimn

≅ε<θ−θ∞→

ossia se, fissato un qualsiasi ε>0, la probabilità che θ coincida con θ (parametro da stimare)tende a diventare unitaria se la lunghezza n del campione osservato tende a diventare moltogrande.

Si tratta allora di capire quanto buona risulti la stima così ottenuta.Il parametro che ci fornisce l’attendibilità di una stima è la deviazione standard

kSσ dello

stimatore kS : infatti, questo parametro ci dice quanto, mediamente, i valori assunti kS si discostano

dal valore medio. Quanto minore è kS

σ , cioè quanto minore è il

discostamento di kS rispetto al suo valore medio 2kσ , tanto migliore

è la stima che stiamo facendo.Dobbiamo allora chiederci da cosa sia influenzata la stima11. Sicuramente, i parametri critici, da

questo punto di vista, sono la durata totale del tempo di osservazione (pari al prodotto del numero Ndi campioni per il passo di campionamento T) e le larghezze di banda Bk dei vari filtri12.

Vedremo che la dispersione (kS

σ ) della stima e la risoluzione spettrale (∆f) sono entrambe

limitate dal valore Ttot=NT del tempo di osservazione e, inoltre, al migliorare della risoluzione

spettrale peggiorerà (cioè crescerà) la dispersione di kS . Vediamo.

Come vedremo tra un attimo, lo stimatore che abbiamo indicato con kS avrà la seguenteespressione:

∑−

=

=1N

0n

2kk )n(y

N

1S

Il motivo per cui calcoliamo la media dei quadrati dei campioni all’uscita dei filtri è che solo cosìriusciamo ad ottenere delle misure ripulite dalle fluttuazioni statistiche. Per ottenere una varianzadella media che sia più piccola possibile, dobbiamo eseguire delle misure che siano statisticamenteindipendenti, cioè dobbiamo ottenere, all’uscita dei filtri, dei campioni yk(n) che siano tra loroindipendenti. Il problema, allora, è duplice: bisogna stabilire in quanto tempo mediare e quantemisure effettuare per ridurre la varianza al di sotto di un certo limite prefissato.

Per ottenere misure (campioni) statisticamente indipendenti all’uscita di un filtro, ci basta tenerconto della lunghezza della risposta all’impulso del filtro; infatti, dobbiamo ricordarci che filtrareequivale a fare una media mobile: il filtro, eseguendo una convoluzione tra ingresso e rispettoall’impulso, fornisce in pratica ciascun campione di uscita come media pesata di un certo numero dicampioni (successivi) dell’ingresso.

11 Si badi che, a voler essere rigorosi, le parole stima e stimatore non sono sinonimi, in quanto la prima è un numero, mentre il

secondo è una variabile aleatoria. Di fatto, però, nei nostri discorsi possiamo confondere le due cose, a patto comunque diricordare la differenza.

12 Tali larghezze di banda non necessariamente sono uguali tra loro

Page 14: Elaborazione numerica dei segnali - users.libero.itusers.libero.it/sandry/download/DSPdownload/DSP_06a.pdf · una frequenza di campionamento f S che rispetti il teorema di campionamento

Appunti di “Elaborazione numerica dei segnali” - Capitolo 6

Autore: Sandro Petrizzelli14

Ad esempio, nel caso di un filtro passa-banda ideale, di banda Bk e frequenza centrale fk, lacorrispondente risposta all’impulso è un sin(f)/f cui è aggiunto un termine di fase dovuto al fatto cheil rettangolo, in frequenza, non è centrato nell’origine:

Hk(f)

risposta all'impulso:( ) tf2j

k

kk

ketB

tBsin)t(h π

ππ

=

Bk

fk

A prescindere dal termine di fase, abbiamo una funzione ( )

k

k

tB

tBsin

ππ

di durata infinita, ma con i

lobi secondari che diventano sempre più trascurabili. Di conseguenza, ciò che determina l’intervallotemporale nel quale il filtro effettua la suddetta media pesata dell’ingresso è sostanzialmente

l’apertura del lobo principale. Dato che la funzione ( )

k

k

tB

tBsin

ππ

presenta il primo zero positivo in 1/Bk

ed il primo zero negativo in -1/Bk, deduciamo l’intervallo di pesatura è all’incirca 2/Bk.Di conseguenza, se vogliamo che il filtro fornisca due campioni praticamente indipendenti,

dobbiamo effettuare le due misure ad una distanza ∆Tk≥ 2/Bk.Quindi, il tempo totale necessario per misurare, all’uscita del filtro, un numero N di campioni

complessi indipendenti è

kktot B

2NTNT ⋅≥∆⋅=

Questo è dunque il tempo necessario per avere una buona stima del valore quadratico medio deicampioni yk(n) e quindi della densità spettrale di potenza nella banda Bk.

Questo tempo Ttot sarà dunque tanto maggiore quanto più piccola è Bk e quanto maggiore è N. Se,invece, Ttot è stato già fissato a priori, allora dovremo scegliere se fare grande N, al fine di ottenereuna misura precisa (in quanto è piccola la varianza) oppure fare piccolo Bk (per avere una migliorerisoluzione in frequenza).

In definitiva, per Ttot fissato, se stringiamo la banda dei filtri, abbiamo una migliore risoluzionespettrale ma otteniamo misure in cui la statistica indipendenza non è completa, nelle quali cioè lavarianza della media non è piccola; se invece allarghiamo la banda dei filtri, la varianza si riduce, maabbiamo una peggiore risoluzione in frequenza, cioè distinguiamo meno le singole componentispettrali.

Page 15: Elaborazione numerica dei segnali - users.libero.itusers.libero.it/sandry/download/DSPdownload/DSP_06a.pdf · una frequenza di campionamento f S che rispetti il teorema di campionamento

Stima spettrale (parte I) - Metodi non parametrici

Autore: Sandro Petrizzelli15

Stima della potenza di segnali correlati gaussiani

Può essere utile cominciare lo studio della stima del valore quadratico medio di un segnale x(n)mediante un semplice esempio.

Consideriamo un segnale x(n), di lunghezza N, composto da campioni reali che supponiamoessere, per semplicità, indipendenti e gaussiani13. In particolare, supponiamo che i campioni siano avalor medio nullo e con uguale valore quadratico medio σ2.

N.B. E’ importante fare una precisazione sulla simbologia: il segnale x(n) cui stiamo facendoriferimento adesso non è il segnale x(n) dei paragrafi precedenti, cioè la versione campionatadella realizzazione del processo, ma è il segnale all’uscita del generico filtro passa-banda delbanco.

Hk(f)

Bk

fk

x(n)realizzazionedel processo

stima dellapotenza in Bk

Di conseguenza, nei discorsi che seguono c’è sempre da sottintendere un pedice “k” che indicache il ragionamento è relativo alla generica banda Bk centrata su fk.

Il nostro obbiettivo, come più volte ripetuto, è quello di ottenere una stima di σ2. Il migliorestimatore che si possa considerare per tale parametro è notoriamente rappresentato dalla seguentevariabile casuale:

∑−

=

=1N

0n

2)n(xN

1S

Dobbiamo cioè prendere i campioni della sequenza, farne il quadrato e poi mediare.Dato che gli x(n) sono, per ipotesi, gaussiani, esiste un noto risultato di probabilità in base al

quale la variabile casuale S ha distribuzione χ2 con N gradi di libertà. Applicando, però, il noto

teorema del limite centrale, possiamo affermare che S è approssimativamente gaussiana.

Indichiamo con S

µ la media di S e con S

σ la deviazione standard di S .

Per calcolare S

σ (si tratta della dispersione della nostra stima), procediamo nel modo classico:

calcoliamo S

σ come radice della varianza 2

Sσ , la quale, a sua volta, si ottiene tramite la relazione

[ ] [ ]SESE 222

S−=σ .

La media di S è data da

[ ] [ ] [ ]( )

221N

0n

2

1N

0n

221N

0n

21N

0n

21N

0n

2

S

NN

1

N

1

)n(xEN

1)n(xE

N

1)n(xE

N

1)n(x

N

1ESE

σ=σ⋅=σ=

=+σ==

=

==µ

∑∑∑∑−

=

=

=

=

=

13 Parliamo di campioni gaussiani perché, in questa fase, non possiamo conoscere a priori il loro valore, per cui dobbiamo

necessariamente considerare ciascuno di essi come una variabile aleatoria, della quale possiamo perciò specificare ladistribuzione ed i relativi parametri.

Page 16: Elaborazione numerica dei segnali - users.libero.itusers.libero.it/sandry/download/DSPdownload/DSP_06a.pdf · una frequenza di campionamento f S che rispetti il teorema di campionamento

Appunti di “Elaborazione numerica dei segnali” - Capitolo 6

Autore: Sandro Petrizzelli16

Abbiamo trovato che la media di S coincide con il parametro (σ2) che S stesso deve stimare.

Questo ci dice che lo stimatore S è corretto14 ed è una condizione irrinunciabile, che deve esseresempre verificata in questo tipo di discorsi.

Abbiamo inoltre che

[ ] [ ]

=

⋅=⋅= ∑∑∑∑

=

=

=

=

1N

0n

21N

0m

22

1N

0m

21N

0n

22 )n(x)m(xEN

1)m(x

N

1)n(x

N

1ESSESE

Supponendo, per semplicità, che i campioni x(n) siano incorrelati e facendo gli opportunipassaggi, si ottiene che

[ ] [ ] [ ] [ ]( )22mn

242

1N

0n

1N

0m

222

2 )n(x)m(xE)NN()n(xNEN

1)n(x)m(xE

N

1SE ≠

=

=

−+== ∑∑

Ci possiamo inoltre ricordare che, per una variabile gaussiana X a valor medio nullo e varianzaσ2, sussiste la relazione E[X4]=3σ4, per cui

[ ] ( ) ( ) ( ) 4444242

4242

2

N

2N2

N

1NN2

N

1)NN(3N

N

1SE σ+σ=σ+=σ+σ=σ⋅−+σ⋅=

Possiamo dunque scrivere che la varianza di S vale

[ ] [ ]N

2

N

2SESE

4444222

S

σ=σ−

σ+σ=−=σ

Infine, possiamo concludere che la deviazione standard di S è

N

222

SSσ=σ=σ

E’ evidente che, all’aumentare di N (numero di campioni indipendenti),la deviazione standard (cioè la dispersione) della nostra stimadiminuisce, cioè la stima è più attendibile.

Una applicazione di questa relazione potrebbe essere la seguente: vogliamo calcolare il numero N

di campioni indipendenti che permette una stima S che, con una probabilità dell’ 85%, differisca daσ2 (cioè dal valore vero) di una quantità inferiore a ±0.05σ2.

In primo luogo, possiamo sfruttare il seguente risultato: una variabile casuale gaussiana, quale è

approssimativamente S , è contenuta, entro un intervallo S

kσ± attorno al suo valore medio, con una

probabilità che è pari a )2/k(erf . Affinché la quantità )2/k(erf assuma il valore 0.85, ènecessario che risulti k=1.43.

Quindi, con probabilità dell’85%, la variabile S è contenuta entro un intervallo S

43.1 σ± attorno

al suo valore medio. Noi vogliamo che questo intervallo coincida con quello indicato dal problema,ossia vogliamo che risulti

2

S05.043.1 σ±=σ±

14 o, ciò che è lo stesso, che lo stimatore è non polarizzato.

Page 17: Elaborazione numerica dei segnali - users.libero.itusers.libero.it/sandry/download/DSPdownload/DSP_06a.pdf · una frequenza di campionamento f S che rispetti il teorema di campionamento

Stima spettrale (parte I) - Metodi non parametrici

Autore: Sandro Petrizzelli17

Questa condizione, ponendo N

22

Sσ=σ , si ottiene per N≅1635.

Campioni correlati

Nell’esempio appena illustrato, abbiamo fatto, come ipotesi di partenza, quella per cui i campionix(n) fossero, ciascuno, variabili aleatorie con distribuzione gaussiana (a media nulla e varianza σ2),indipendenti e incorrelate tra loro. Nel caso in cui, invece, tali campioni siano correlati, si puòverificare, con calcoli un po' laboriosi, che, a parità del numero N di campioni, ladispersione della stima è maggiore. Per renderci conto di questo, senza fare calcolitroppo onerosi, facciamo riferimento ad un esempio semplice.

Siano X e Y due variabili aleatorie (ad esempio due campioni successivi delle nostrarealizzazione) gaussiane, entrambe a valor medio nullo e varianza σ2 (il che significa che

[ ] [ ] 222 YEXE σ== ). Supponiamo che tali variabili siano tra loro correlate, in base ad una relazionedel tipo seguente:

WXY +µ=

dove W è un’altra variabile aleatoria, incorrelata con X e quindi anche con Y, mentre µ è un genericonumero reale.

E’ evidente che

[ ] [ ] [ ] [ ] [ ] ( ) ( )[ ]( ) ( )[ ] ( ) [ ] [ ] ( ) ( )222222222222

222222222

1021WXE2XE21WXXE21

WXXE21XYE2XEYEXY2XYEWE

µ−σ=−σµ−µ+σ=µ−µµ−µ+σ=+µµ−µ+σ=

=+µµ−µ+σ=µ−µ+=µ−µ+=

dove evidentemente abbiamo si è tenuto conto che E[X⋅W]=0 data la incorrelazione tra le duevariabili.

Altri risultati importanti, che si possono ricavare e che sfrutteremo più avanti sono

[ ][ ][ ] ( ) 4222

2

44

21YXE

XYE

3XE

σµ+=

µσ=

σ=

Adesso vediamo di sfruttare questi risultati, considerando una sequenza x(m) di campionicorrelati. In particolare, supponiamo che

[ ]22

2

]E[x(n)

mn )m(x)n(xE

σ=

≠µσ=

Volendo stimare la varianza (o valore quadratico medio) dei campioni della sequenza, dobbiamo

ancora una volta considerare lo stimatore ∑−

=

=1N

0n

2)n(xN

1S .

Page 18: Elaborazione numerica dei segnali - users.libero.itusers.libero.it/sandry/download/DSPdownload/DSP_06a.pdf · una frequenza di campionamento f S che rispetti il teorema di campionamento

Appunti di “Elaborazione numerica dei segnali” - Capitolo 6

Autore: Sandro Petrizzelli18

La media di S è ancora σ2, mentre cambia la varianza, data la presenza della correlazione: si puòverificare che tale varianza risulta essere

( ) 2222

SS21

N

2µ+µ−σ=σ=σ

E’ evidente, allora che, al contrario di quanto trovato prima, questa varianza non decresce a 0 perN che tende all’infinito: essa ha infatti un valore minimo (per N→∞) che è

( ) 22 222

minSµσ=µσ=σ

Consideriamo adesso un altro caso e cioè quello in cui i campioni della sequenza x(n) sonocorrelati, ma questa correlazione diminuisce all’aumentare della differenza di indice; in pratica,supponiamo che due generici campioni x(m) ed x(m+k) della sequenza risultino tanto meno correlatiquanto maggiore è il valore di k, cioè tra la distanza temporale tra i due campioni. Per esempio,possiamo supporre che valga una legge del tipo

[ ] 1 )km(x)m(xE 2k ≤ρσρ=+

Essendo ρ una quantità non superiore ad 1, la correlazione tra i due campioni diminuisceesponenzialmente all’aumentare di k.

Andiamo allora a ripetere il calcolo della varianza di S sotto questa nuova condizione: dopoqualche passaggio, si trova che

2

22

S 1

1

N

2

ρ−ρ+

⋅σ≅σ

Anche in questo caso, notiamo che S

σ non decresce a 0 per N che tende ad infinito. Tuttavia,

confrontando questa espressione della varianza con quella ottenuta in assenza di correlazione tra i

campioni, si nota che la differenza è tutta nel termine 2

2

1

1

ρ−ρ+

. Possiamo allora esprimerci nel modo

seguente: avendo da stimare la varianza di una serie di N campioni non indipendenti, ma correlaticon un decadimento esponenziale (di ragione ρ) della correlazione, è come avere a disposizione unnumero equivalente di campioni indipendenti pari a

2

2

eq 1

1MM

ρ+ρ−

Questo numero è chiaramente inferiore ad M, il che significa che, per ottenere un prefissato valoredi

Sσ , dobbiamo considerare, in questo caso, più campioni rispetto al caso dei campioni non

correlati. Più campioni significa, evidentemente, un maggiore tempo di osservazione e si tratta di unrisultato assolutamente intuitivo: maggiore è la correlazione tra i campioni, più tempo dovremoaspettare per ottenere delle misure che siano significative per la nostra stima.

In altre parole, quando i campioni del segnale sono correlati, ladispersione della stima della loro varianza cresce, in quantodiminuisce il numero equivalente di campioni indipendenti checontribuiscono alla stima.

Page 19: Elaborazione numerica dei segnali - users.libero.itusers.libero.it/sandry/download/DSPdownload/DSP_06a.pdf · una frequenza di campionamento f S che rispetti il teorema di campionamento

Stima spettrale (parte I) - Metodi non parametrici

Autore: Sandro Petrizzelli19

Per esprimerci in modo molto pratico, è come dire che, volendo stimare la distribuzione delle etàdei cittadini in una nazione, non è opportuno raccogliere i dati solo nelle scuole elementari o nellecase di riposo.

METODO DEL PERIODOGRAMMA

Vediamo un metodo con cui si può realizzare un analizzatore di spettro numerico. Per semplicità,supponiamo che il processo del quale vogliamo stimare la densità spettrale sia a media nulla. Oltre aquesto, siamo sempre nelle ipotesi che il processo sia stazionario ed ergodico.

Per stimare la densità spettrale di potenza, abbiamo già osservato che possiamo seguire unmetodo indiretto, che consiste nello stimare prima la funzione di autocorrelazione15 delprocesso [ ])t(x)t(xE)( *

XX τ+=τγ e poi nel ricavare la densità spettrale di potenza come trasformata

di γXX(τ):

[ ]

)(Rlim)(dt)t(x)t(xT

1)(R :stimate funzioni

de)()f()t(x)t(xE)( :effettive funzioni

XXT

XX

2/T

2/T

*

0XX

f2jXXXX

Fourier*XX

0

0

0

τ=τγ→τ+=τ

ττγ=Γ →τ+=τγ

∞→−

+∞

∞−

τπ−

In pratica, quindi, possiamo stimare γXX(τ) tramite il limite, per un tempo di osservazione T0

tendente all’∞, della funzione di autocorrelazione media RXX(τ). Stimata γXX(τ), ci basta trasformareper ottenere una stima della densità spettrale di potenza ΓXX(f).

Questo vale, però, nell’ipotesi di ragionare con una realizzazione tempo-continua del processo edè il punto cui siamo arrivati nei precedenti paragrafi. Al contrario, noi abbiamo a che fare con unaosservazione campionata (oltre che limitata nel tempo) di una realizzazione del processo. In questocaso, possiamo pensare di approssimare γXX(τ) mediante una funzione tempo-discreta così definita:

1Nk0 )ki(x)i(xN

1)k(r

1kN

0i

*XX −≤≤+= ∑

−−

=

Il motivo per cui l’indice superiore della sommatoria vale N-k-1 è evidente: data lapresenza del termine x(i+k) all’interno della sommatoria, l’indice i+k non puòsuperare il valore N-1, dato che i campioni che abbiamo a disposizione sono solo N(numerati da 0 ad N-1), per cui il massimo valore di k non può che essere N-k-1.

Stiamo dunque approssimando la funzione tempo-continua γXX(τ) tramite la funzione tempo-discreta rXX(k).

Per semplicità, possiamo supporre che il processo considerato sia a valori reali, per cui anche ivalori campionati dell’unica realizzazione da noi osservata saranno reali e quindi l’operatore dicomplesso coniugato non avrà alcun esito pratico:

1Nk0 )ki(x)i(xN

1)k(r

1kN

0iXX −≤≤+= ∑

−−

=

15 Ricordiamo due note proprietà della funzione di autocorrelazione di un processo stazionario: la prima, come già ricordato in

precedenza, è che γXX(τ) dipende solo da τ e non da t; la seconda è che γXX(τ) è una funzione pari, il che significa cheγXX(τ)=γXX(-τ)

Page 20: Elaborazione numerica dei segnali - users.libero.itusers.libero.it/sandry/download/DSPdownload/DSP_06a.pdf · una frequenza di campionamento f S che rispetti il teorema di campionamento

Appunti di “Elaborazione numerica dei segnali” - Capitolo 6

Autore: Sandro Petrizzelli20

Questa sequenza rXX(k) è uno stimatore per γXX(k), cioè per la versione campionata di γXX(τ).Possiamo però renderci conto che si tratta di uno stimatore non sempre buono. Calcoliamo la suamedia, ossia la media dei valori del vettore rXX(k):

[ ] [ ] )k(N

|k|1)k(

N

|k|N...?...)ki(x)i(xE

N

1)ki(x)i(x

N

1E)k(rE XXXX

1kN

0i

1kN

0iXX γ

−=γ

==+=

+= ∑∑

−−

=

−−

=

dove ricordiamo che il valore assoluto posto sul ritardo k serve in quanto il calcolo della media deveincludere anche i campioni di rXX(k) per k<0.

Questo risultato ci dice che la media di rXX(k) non coincide con γXX(k), ma differisce da esso perun termine -(|k|/N)γXX(k). Si esprime questo dicendo che la stima presenta una polarizzazione: ilvalore atteso dello stimatore è diverso dal valore vero. Nel nostro caso, la polarizzazione (bias) èproprio il termine -(|k|/N)γXX(k).

D’altra parte, però, si osserva che tale polarizzazione decresce a zero per N→∞:

[ ] )k()k(rElim 0 )k(N

|k|lim che dato XXXX

NXX

Nγ=→=γ

∞→∞→

Si dice allora che lo stimatore rXX(k) è asintoticamente non polarizzato. In aggiunta aquesto, si può anche verificare che la sua varianza è molto piccola e, soprattutto, tende ad ∞ perN→∞.

Sulla base di questi risultati, possiamo comunque dire che rXX(k) è uno stimatoreconsistente per γXX(k).

In realtà, abbiamo un’altra possibilità, ossia possiamo pensare ad un nuovo stimatore per γXX(k):si tratta della funzione

1Nk0 )ki(x)i(xkN

1)k('r

1kN

0i

*XX −≤≤+

−= ∑

−−

=

ottenuta semplicemente da rXX(k) dividendo la sommatoria non per N, ma per N-k.Si può verificare che, con questo semplice accorgimento, la media di r’XX(k) coincide con γXX(k):

[ ] [ ] )k(...?...)ki(x)i(xEkN

1)ki(x)i(x

kN

1E)k('rE XX

1kN

0i

*1kN

0i

*XX γ==+

−=

+

−= ∑∑

−−

=

−−

=

Deduciamo dunque che r’XX(k) è uno stimatore non polarizzato per la funzione diautocorrelazione effettiva rXX(k) del nostro processo. Calcoliamo allora la varianza di questostimatore, per valutarne l’attendibilità: si dimostra (formula di Jenkins e Watts) che

[ ]( )

[ ]∑+∞

−∞=

+γ−γ+γ−

≅n

XX*

XX

2

XX2XX )mn()mn()n(kN

N)k('rVar

Questa quantità decresce a zero, per N→∞, nel momento in cui risulta convergente la sommatoria

∑+∞

−∞=

γn

2

XX )n( : deve cioè risultare

∞<γ∑+∞

−∞=n

2

XX )n(

Page 21: Elaborazione numerica dei segnali - users.libero.itusers.libero.it/sandry/download/DSPdownload/DSP_06a.pdf · una frequenza di campionamento f S che rispetti il teorema di campionamento

Stima spettrale (parte I) - Metodi non parametrici

Autore: Sandro Petrizzelli21

Sotto questa condizione, quindi, abbiamo un nuovo stimatore consistente per γXX(k).Rispetto alla funzione rXX(k) vista prima, la funzione r’XX(k) ha dunque il pregio di avere un

valore atteso coincidente con il valore vero (stimatore non polarizzato). A fronte di questo, però, sipuò verificare che, per alti valori di k, la funzione r’XX(k) risulta avere una grande varianza,specialmente quando k si avvicina ad N. Dato che questo problema non si presenta con rXX(k), neiprossimi discorsi faremo riferimento a rXX(k).

Se prendiamo rXX(k) come stimatore per la funzione di autocorrelazione del processo,automaticamente stiamo considerando, come stimatore per la densità spettrale di potenza ΓXX(f), latrasformata di Fourier di rXX(k):

∑−

−=

π−=1N

N1m

fmT2jXXXX e)m(r)f(P

dove la sommatoria è estesa anche per m<0 in quanto bisogna sempre includere anche i campioni dirXX(m) per m<0.

Che questo sia uno stimatore per ΓXX(f) si può vedere facilmente, sostituendo nella sommatorial’espressione di rXX(m):

221N

0n

fn2j

1N

0m

1kN

0i

fm2j*1N

0m

fm2j1kN

0i

*XX

)f(XN

1e)n(x

N

1

....?...e)mi(x)i(xN

1e)mi(x)i(x

N

1)f(P

==

==+=

+=

∑ ∑∑ ∑−

=

π−

=

−−

=

π−−

=

π−−−

=

dove ricordiamo che X(f) è la trasformata di Fourier della sequenza tempo-discreta x(n).La funzione PXX(f) così costruita è dunque uno stimatore della densità spettrale di potenza del

processo in esame ed è chiamata periodogramma. Il periodogramma si trasforma ovviamente dastimatore in stima nel momento in cui sostituiamo, nella sua espressione, i campioni x(n) dellarealizzazione che abbiamo effettivamente osservato.

E’ interessante osservare che anche il periodogramma, come la sua antitrasformata rXX(m), è unostimatore polarizzato. Possiamo infatti facilmente calcolare la media di PXX(f):

[ ] [ ] [ ]

∑∑∑−

−=

π−

−=

π−−

−=

π−−

−=

π−

γ

−=

===

=

1N

N1m

fmT2jXX

1N

N1m

fmT2jXX

1N

N1m

fmt2jXX

1N

N1m

fmT2jXXXX

e)m(N

m1

e)m(rEe)m(rEe)m(rE)f(PE

Come si vede, questo valore atteso differisce dal valore vero ΓXX(f), il che significa appunto chelo stimatore è polarizzato.

Ad ogni modo, così come accade per rXX(m), anche PXX(f) risulta asintoticamente non polarizzato:infatti, calcolandone il limite per N→∞, la polarizzazione si azzera e quindi il valor medio di PXX(f)tende a coincidere con ΓXX(f). Al contrario, la varianza di PXX(f) in generale non decresce a zero perN→∞, per cui, in generale, PXX(f) non è uno stimatore consistente delladensità spettrale di potenza ΓXX(f) del processo.

Possiamo dunque riepilogare quanto detto in questo paragrafo: la funzione rXX(m) è uno stimatoreconsistente della funzione di autocorrelazione γXX(m) del processo, mentre la sua trasformata diFourier PXX(f), che abbiamo chiamato periodogramma, non è uno stimatore consistente delladensità spettrale di potenza ΓXX(f) del processo.

Page 22: Elaborazione numerica dei segnali - users.libero.itusers.libero.it/sandry/download/DSPdownload/DSP_06a.pdf · una frequenza di campionamento f S che rispetti il teorema di campionamento

Appunti di “Elaborazione numerica dei segnali” - Capitolo 6

Autore: Sandro Petrizzelli22

Come vedremo tra poco, la polarizzazione che compare sia in rXX(m) sianel periodogramma PXX(f) è dovuta alla finestratura, ossia al fattodi usare sequenze temporali di durata finita.

IMPLEMENTAZIONE DEL PERIODOGRAMMA TRAMITE LA DFTDobbiamo adesso chiederci come implementare praticamente la stima fornita dal periodogramma.Possiamo partire dal fatto che il periodogramma PXX(f) è definito come la trasformata di Fourier

della seguente funzione:

1Nk0 )ki(x)i(xN

1)k(r

1kN

0i

*XX −≤≤+= ∑

−−

=

Per semplicità, supponiamo che x(n) sia reale, per cui l’operatore complesso coniugato si puòeliminare:

1Nk0 )ki(x)i(xN

1)k(r

1kN

0iXX −≤≤+= ∑

−−

=

Questa espressione può anche essere scritta in altro modo:

0ki

Nkiper 0k)con x(i )ki(x)i(x

N

1)k(r

1N

0iXX

<+∀>+∀

=++= ∑−

=

Posta in questa forma, quella sommatoria è rappresentativa di una particolare convoluzione tradue sequenze di lunghezza N; si tratta di convolvere x(n) per se stessa ribaltata: applicando infatti ladefinizione di convoluzione, abbiamo che

( ) ∑∑+∞

−∞=

+∞

−∞=

+=−−=−mm

)mn(x)m(x)m(nx)m(x)n(x*)n(x

Quindi per effettuare una correlazione tra x(n) e se stessa (cioèuna autocorrelazione), ci basta effettuare la convoluzione tra x(n)e se stessa ribaltata16.

Dovremmo perciò calcolare la suddetta convoluzione e poi trasformare secondo Fourier.In realtà, possiamo procedere più velocemente se ragioniamo direttamente nel dominio dellafrequenza e applichiamo il concetto della DFT.

Infatti, avendo detto che il periodogramma risulta essere

221N

0n

fn2jXX )f(X

N

1e)n(x

N

1)f(P == ∑

=

π−

è evidente che, campionando a passo 1/NT, otteniamo

221N

0n

nNT

k2j

XX NT

kX

N

1e)n(x

N

1

NT

kP

==

=

π−

16 Se la sequenza x(n) fosse complessa, per cui la correlazione necessiterebbe dell’applicazione dell’operatore di coniugio,

basterebbe convolvere x(n) con se stessa sia ribaltata sia coniugata.

Page 23: Elaborazione numerica dei segnali - users.libero.itusers.libero.it/sandry/download/DSPdownload/DSP_06a.pdf · una frequenza di campionamento f S che rispetti il teorema di campionamento

Stima spettrale (parte I) - Metodi non parametrici

Autore: Sandro Petrizzelli23

In base a questa relazione, se ci calcoliamo la DFT di x(n), ne facciamo il modulo quadro e poidividiamo per N, otteniamo proprio la versione campionata del periodogramma.

In conclusione, il tipo di procedura appena visto (metodo del periodogramma) per stimare ladensità spettrale di potenza di un processo consiste dei seguenti passi fondamentali:

• si considera la sequenza x(n) di partenza, rappresentativa di una realizzazione osservata su untempo NT e campionata a passo T;

• si calcola la DFT, basata su N campioni, di x(n);

• si calcola il modulo quadro del vettore della DFT e si divide per N.

Il risultato è la versione campionata del periodogramma.Naturalmente, viene spontaneo porsi una domanda: dato che stiamo considerando una versione

campionata del periodogramma, cioè sostanzialmente un campionamento in frequenza, possonoesserci dei problemi legati alla corrispondente periodicizzazione nei tempi? In effetti, il metodoappena descritto consente solo di stimare la densità spettrale di potenza del processo, mentre valeggermente modificato se si vuole usare il periodogramma per risalire alla funzione diautocorrelazione. Vediamo quali sono queste modifiche.

Abbiamo visto che2

XX NT

kX

N

1

NT

kP

=

Esplicitando in modo esteso il modulo quadro della DFT di x(n), abbiamo quanto segue:

=

∑∑

=

−−

=

1N

0n

nkN

1N

0n

nkN

*XX W)n(xW)n(x

N

1

NT

kP

Stiamo dunque eseguendo il prodotto della DFT di x(n) con la DFT di x(n) ribaltata nei tempi econiugata (se complessa). Nei tempi, stiamo perciò convolvendo circolarmente tali due sequenze. Mala funzione rXX(k), antitrasformata del periodogramma, è la convoluzione lineare di tali sequenze. Inaltre parole, finché usavamo il dominio continuo delle frequenze, cioè consideravamo PXX(f), nonavevamo problemi; dal momento in cui siamo passati a considerare il dominio campionato dellefrequenze, il prodotto in frequenza è diventato equivalente alla convoluzione circolare nel tempo enon più alla convoluzione lineare.

Per evitare, allora, ogni problema, dobbiamo forzare la convoluzione circolare a dare lo stessorisultato della convoluzione lineare e sappiamo che questo si ottiene semplicemente allungando x(n)con N zeri (zero padding) prima di calcolarne la DFT. Solo con questa ipotesi di partenza, siamo

certi che la IDFT del vettore

NT

kPXX potrà darci la funzione rXX(k) senza alcun errore.

Possiamo dunque concludere dicendo che, nel calcolare la DFT di x(n), dobbiamo scegliere tradue possibilità:

• la prima è quella di voler stimare solo lo spettro di potenza del processo: in questo caso, nondobbiamo eseguire alcuna IDFT e quindi possiamo lasciare invariata la x(n) e calcolarne laDFT su N campioni;

• la seconda è quella invece di voler stimare anche la funzione di autocorrelazione del processo,ottenuta come IDFT di PXX(k). In questo caso, prima di calcolare la DFT di x(n), dobbiamoallungare x(n) con N zeri in coda (zero padding), in modo da forzare la convoluzione

Page 24: Elaborazione numerica dei segnali - users.libero.itusers.libero.it/sandry/download/DSPdownload/DSP_06a.pdf · una frequenza di campionamento f S che rispetti il teorema di campionamento

Appunti di “Elaborazione numerica dei segnali” - Capitolo 6

Autore: Sandro Petrizzelli24

circolare a dare gli stessi risultati della convoluzione lineare, dopo di che dobbiamo calcolare laDFT basata su 2N campioni.

Osservazione: allungamento con zeri della x(n)

Abbiamo appena visto che, se decidiamo di usare la DFT per ottenere la versione campionata delperiodogramma, essendo x(n) lunga N ci basterebbe calcolare la DFT basata su N punti17:otterremmo il periodogramma alle frequenze fk=k/NT, per k=0,...,N-1. Nella pratica, però, un siffattocampionamento non fornisce una buona rappresentazione grafica dello spettro continuo checostituisce il periodogramma PXX(f). Per migliorare la rappresentazione grafica, è allora sufficienteaumentare il numero delle frequenze in cui si effettua il campionamento e sappiamo che questoequivale, ancora una volta, ad allungare x(n) con un certo numero di zeri posti in coda.

Ad esempio, se allunghiamo x(n) con M zeri, per cui calcoliamo la DFT di una sequenza lungaL=N+M, otteniamo

1-L0,1,....,k LT

kX

N

1e)n(x

N

1

LT

kP

221N

0n

nNT

k2j

XX =

==

=

π−

Il caso considerato prima è quello in cui gli zeri aggiunti sono N, per cui L=2N:

1-N0,1,....,2k NT2

kX

N

1e)n(x

N

1

NT2

kP

221N

0n

nNT

k2j

XX =

==

=

π−

In pratica, calcolare la DFT della sequenza x(n) allungata con N zeri equivale semplicemente adinterpolare, su un intervallo di ampiezza 1/2NT, la DFT della x(n) originale e tale interpolazioneviene effettuata con l’interpolatore ideale (cioè il sin(f)/f). Questo significa che i valori assunti dallaDFT nelle frequenze fk=k/NT sono identici a quelli assunti dalla DFT della sequenza x(n) originale,non bordata con zeri; ci stiamo limitando ad aggiungere i valori di PXX(f) alle frequenze k/2NT.

E’ opportuno ancora una volta ribadire che l’interpolazione deve essere fatta su X(k/NT) e non sulmodulo quadro; in altre parole, prima si fa l’interpolazione e poi il moduloquadro e non viceversa.

E’ anche importante sottolineare un’altra cosa: allungare la sequenza x(n) con Mzeri e calcolare la DFT basata su L>N punti non migliora larisoluzione spettrale della stima che stiamo effettuando. Si trattasemplicemente di un metodo semplice per incrementare la risoluzione con cui rappresentiamo ilperiodogramma. Al contrario, la risoluzione spettrale della nostra stima è univocamente determinatadal numero N di campioni da cui è composta x(n) e cioè sostanzialmente dalla durata della finestra diosservazione.

Ad ogni modo, dunque, l’allungamento di x(n) con N zeri in coda, pur non essendo semprenecessario (lo è solo quando si prevede di dover fare una IDFT del periodogramma), di fatto si rivelacomunque utile, in quanto fornisce una migliore rappresentazione grafica del periodogramma econsente di calcolare la IDFT del periodogramma stesso senza incorrere in errori.

17 Ricordiamo che una sequenza x(n) lunga N è del tutto equivalente alla sua DFT se e solo se quest’ultima è calcolata su almeno N

punti, cioè se si considerano almeno N campioni (in frequenza) nel periodo fondamentale [-fC/2,fC/2]; in caso contrario, se cioè siprendono meno di N campioni in frequenza, la corrispondente periodicizzazione nel tempo porta aliasing, il che significa che nonsarebbe possibile risalire ad x(n) partendo da X(k):

Page 25: Elaborazione numerica dei segnali - users.libero.itusers.libero.it/sandry/download/DSPdownload/DSP_06a.pdf · una frequenza di campionamento f S che rispetti il teorema di campionamento

Stima spettrale (parte I) - Metodi non parametrici

Autore: Sandro Petrizzelli25

Esempio

Un esempio molto significativo, relativo alla questione dell’allungamento con 0 della sequenzax(n) di partenza, è il seguente.

Consideriamo un segnale analogico x(t) costituito da due sole componenti spettrali

( ) ( )tf2costf2sin)t(x 21 ⋅⋅π+⋅⋅π=

Supponiamo che tali componenti spettrali siano distanziate di una quantità ∆f=f2-f1>0.Campioniamo questo segnale con una frequenza di campionamento opportuna:

( ) ( )nTf2cosnTf2sin)nT(x 21 ⋅⋅π+⋅⋅π=

Supponiamo, per esempio, che il rapporto tra f1 e la frequenza di campionamento sia pari a 0.135.Possiamo allora scrivere che

( ) ( )( )nf135.02cosn135.02sin)n(x ⋅∆+⋅π+⋅⋅π=

Supponiamo di aver preso N=16 campioni di questo segnale, per cui n varia tra 0 e 15.

Vogliamo adesso calcolare la densità spettrale 2

)f(XN

1)f(P = alle frequenze fk=k/L con

k=0,1,...,L-1 e vogliamo vedere cosa cambia prendendo L=8,16,32,128 e ∆f=0.06,0.01.A meno del fattore di scala 1/N, i risultati che si ottengono, al variare di L, sono rappresentati

nella figura seguente, nel caso di ∆f=0.06:

Si nota che la risoluzione (intesa come numero di punti riportati) con cui rappresentiamo P(k)aumenta, ovviamente, all’aumentare di L, mentre invece L non influisce sulla risoluzione spettrale

Page 26: Elaborazione numerica dei segnali - users.libero.itusers.libero.it/sandry/download/DSPdownload/DSP_06a.pdf · una frequenza di campionamento f S che rispetti il teorema di campionamento

Appunti di “Elaborazione numerica dei segnali” - Capitolo 6

Autore: Sandro Petrizzelli26

della stima. In questo caso, il valore ∆f=0.06 è sufficiente a risolvere le due componenti spettrali,cioè ad individuarle come distinte.

Adesso ripetiamo il calcolo prendendo invece ∆f=0.01. Si ottiene quanto segue:

In questo caso, le due componenti spettrali sono troppo vicine per essere risolte, a prescindere dalnumero L di campioni che abbiamo calcolato. Tale numero L contribuisce solo a migliorare larisoluzione della rappresentazione dello spettro.

Metodo del periodogramma, DFT e banco di filtri

Abbiamo visto che il metodo del periodogramma, per stimare la densità spettrale di potenzadi un processo, consiste dei seguenti passi fondamentali:

• data la sequenza x(n) di partenza, rappresentativa di una realizzazione osservata su un tempoNT e campionata a passo T, la allunghiamo con N zeri in coda;

• calcoliamo la DFT, basata su 2N campioni, di x(n);

• facciamo il modulo quadro del vettore della DFT e dividiamo per N.

Il risultato di questi passaggi è la versione campionata del periodogramma.E’ facile rendersi conto di come questo discorso si ricongiunga allo schema del banco di filtri da

cui siamo partiti all’inizio:

Page 27: Elaborazione numerica dei segnali - users.libero.itusers.libero.it/sandry/download/DSPdownload/DSP_06a.pdf · una frequenza di campionamento f S che rispetti il teorema di campionamento

Stima spettrale (parte I) - Metodi non parametrici

Autore: Sandro Petrizzelli27

x(n)n=0,1,...N-1

H1(f)

H2(f)

HK(f)

y1(n)

y2(n)

yK(n)

Abbiamo detto che, in questo schema, ogni filtro deve isolare il contenuto spettrale di x(n) in unaassegnata banda Bk , in modo che sia poi possibile stimare la potenza allocata in tale banda. Le Kuscite del banco di filtri costituiscono, in tal modo, una stima della densità spettrale di potenza delprocesso in osservazione. Ovviamente, l’intervallo di frequenza in osservazione è quellocorrispondente al periodo fondamentale [-fC/2,fC/2], date le note proprietà periodiche dello spettro dix(n).

Se i filtri passa-banda sono filtri a campionamento in frequenza, complessi, con frequenzecentrali

2

N,......,0,....,1

2

N-k

NT

kfk +==

allora il filtro k-simo ha la seguente risposta all’impulso:

1-N0,....,k We)n(h nkN

N

nk2j

k ===π

Funzione di trasferimento di un cosiddetto filtro a campionamento in frequenza: è un filtropassa-banda la cui funzione di trasferimento è un sin(f)/f centrato a frequenza fk e con banda 2/NT, il che

significa che il lobo principale è delimitato dai due zeri in posizioni fk+1/NT e fk-1/NT. Questo tipo difunzione di trasferimento deriva dall’interpolazione (con un filtro ideale nel dominio della frequenza) di

un segnale costituito da una sequenza di impulsi (cioè una sequenza di sinusoidi a frequenza fk=k/NT). Laconvoluzione con il sin(f)/f nel dominio delle frequenze corrisponde alla finestratura rettangolare nel

dominio dei tempi.

Page 28: Elaborazione numerica dei segnali - users.libero.itusers.libero.it/sandry/download/DSPdownload/DSP_06a.pdf · una frequenza di campionamento f S che rispetti il teorema di campionamento

Appunti di “Elaborazione numerica dei segnali” - Capitolo 6

Autore: Sandro Petrizzelli28

Questo filtro k-simo riceve in ingresso la sequenza x(n) di lunghezza N e da una uscita yk(n) che èla convoluzione di x(n) con hk(n):

∑∑∑−

=

−−

=

−+∞

−∞=

==−==1N

0m

nkN

mkN

1N

0m

k)mn(N

mkkk WW)m(xW)m(x)mn(h)m(x)n(h*)n(x)n(y

Per i noti problemi legati alla convoluzione tra due sequenze di durata finita18, l’uscita di questofiltro è valida solo nell’istante N-1, quando cioè il transitorio iniziale di caricamento dei dati nelfiltro è terminato, ma anche i dati sono terminati (visto che x(n) è lunga N). Deduciamo che l’unicocampione di cui disponiamo all’uscita del filtro k-simo è il campione

∑−

=

−==−1N

0m

mkNk W)m(x....(?)....)1N(y

A ben vedere, questo campione non è altro che il campione di posto k nel vettore X(k) checostituisce la DFT della sequenza x(n).

Quindi, determinare le uscite yk(N-1) del banco di N filtri passa-basso equivale al calcolo dellaDFT della sequenza di ingresso. Con questo ci riconduciamo evidentemente al metodo delperiodogramma: infatti, nello schema ridisegnato prima, in uscita dal k-simo filtro abbiamo un unicocampione, X(k); questo viene elevato al quadrato in modo da conoscerne la potenza |X(k)|2; talepotenza, divisa per N, è la stima P(k) della densità spettrale di potenza del processo in quella banda.Prendendo le N uscite e dividendole tutte per N, si ottiene quindi il periodogramma, cioè la stimadella densità spettrale di potenza del processo in tutte le frequenze fk=k/NT dell’intervallo da -fC/2 afC/2.

In effetti, si nota una differenza con il discorso iniziale da cui siamo partiti: infatti, quandoabbiamo introdotto la struttura a banco di filtri, abbiamo detto che, in uscita da ciascun filtro, ciaspettiamo una sequenza yk(n) di N campioni che dobbiamo quadrare e poi mediare per ottenere lastima della potenza in quella banda:

Hk(f)

Bk

fk

yk(n)realizzazionedel processo

x(n)stima della

potenza in Bk

∑−

=

=1N

0n

2kk )n(y

N

1S

Al contrario, abbiamo appena visto che, con la DFT, l’uscita di ciascun filtro è rappresentato daun unico campione: il quadrato di tale campione, diviso per N, è la stima della potenza nella bandaconsiderata.

18 Ci si riferisce ai problemi che sono stati già ampiamente esaminati nel capitolo sulla DFT: mentre hk(n) è una funzione

identicamente nulla nell’intervallo in cui non è specificata, x(n) non lo è, in quanto il numero limitato di campioni dipende solodall’aver effettuato l’osservazione (ed il campionamento) su un intervallo di durata finita; nel fare la convoluzione con hk(n),sono coinvolti i campioni x(-1), x(-2),..., che non sono a nostra disposizione; di conseguenza, quei campioni yk(n) dell’uscita checoinvolgono campioni non noti di x(n) sono necessariamente non significativi; se h(n) ha lunghezza N, come supposto, allora iprimi N-1 campioni dell’uscita non sono significativi, per cui resta un solo campione valido, che è quello relativo all’istante N-1:si tratta infatti dell’unico campione di uscita che coinvolge solo campioni noti di x(n). Teoricamente, sarebbero significativianche i campioni successivi di yk(n), da quello nell’istante N fino a quello calcolato nell’istante 2N-2 (in tutto, sono teoricamente2N-1 campioni di uscita). Tuttavia, l’ingresso è applicato fino all’istante N-1, per cui non ha senso calcolare le uscite oltre taleistante.

Page 29: Elaborazione numerica dei segnali - users.libero.itusers.libero.it/sandry/download/DSPdownload/DSP_06a.pdf · una frequenza di campionamento f S che rispetti il teorema di campionamento

Stima spettrale (parte I) - Metodi non parametrici

Autore: Sandro Petrizzelli29

In effetti, quello appena descritto non è l’unico modo di procedere: anche se lo vedremo megliopiù avanti, possiamo infatti già anticipare che non è sempre opportuno calcolare la DFT dell’interasequenza x(n) di partenza; al contrario, se si vuole ridurre la varianza della stima ottenuta colperiodogramma, è più opportuno spezzare x(n) in diverse sottosequenze, di ciascuna delle qualicalcolare la DFT. In tal modo, con il calcolo cioè di più DFT, ci troviamo, in uscita dai filtri, piùcampioni (tanti quante sono le DFT calcolate, cioè il numero di sottosequenze), che andremo aquadrare e mediare in base al discorso di prima.

Per concludere, osserviamo quanto segue: calcolando la DFT della sequenza di ingresso (o di unagenerica sottosequenza), è possibile dimostrare che i campioni X(k) (cioè le N uscite degli N filtripassa-banda) sono sostanzialmente indipendenti tra di loro; dato che tali campioni si riferiscono afrequenze distanti 1/NT una dall’altra, questo garantisce che la risoluzione spettrale siaeffettivamente 1/NT, cioè la massima possibile (in base a quanto detto nei paragrafi precedenti).

ALLUNGAMENTO DELLA FINESTRA DI OSSERVAZIONE

E’ noto che la stima del valor medio di una sequenza campioni (cioè di variabili casuali) è tantomigliore quanto maggiore è il numero di campioni: il motivo è che, all’aumentare del numero dicampioni, mentre rimane invariato il valore medio dello stimatore, diminuisce progressivamente lasua varianza. Si potrebbe allora pensare di considerare un numero maggiore di campioni anche per lastima della densità spettrale di potenza di un processo. In altre parole, si potrebbe pensare dimigliorare la dispersione della stima aumentando il numero N di campioni e quindi, di conseguenza,la durata Ttot=NT della finestra di osservazione del processo (o meglio della sua realizzazione)

Al contrario, l’aumento dei campioni non produce alcun miglioramento(cioè in questo caso una riduzione) della dispersione. Il motivo è che icampioni in più non vengono usati per mediare una stessa misura ripetuta più volte, ma vannosemplicemente ad allungare il vettore di cui si effettuerà la trasformata di Fourier.

L’effetto concreto dell’aumento del numero N di campioni riguarda solo due aspetti:

• in primo luogo, risulta aumentato, della stessa quantità, il numero dei campioni della stimadella densità spettrale di potenza, il che consente semplicemente una migliore rappresentazionegrafica di tale stima;

• in secondo luogo, come già sappiamo, risulta migliorata la risoluzione spettrale, ossia laminima distanza alla quale due componenti spettrali possono essere individuate come distinte:infatti, tale risoluzione spettrale è B=1/NT=1/Ttot ed è quindi inversamente proporzionale alladurata delle finestra di osservazione.

In ogni caso, la dispersione della stima non cambia all’aumentare di N: sostanzialmente, ladispersione della stima non diminuisce in quanto ciò che diminuisce è proprio la risoluzionespettrale.

Possiamo analizzare la cosa analiticamente.Stiamo supponendo di effettuare un campionamento della generica realizzazione del nostro

processo durante un intervallo Ttot=NT, dove N è il numero di campioni presi e T il passo dicampionamento. Tale sequenza campionata viene filtrata da un banco di filtri ed all’uscita di ciascunfiltro viene effettuata una media dei quadrati dei campioni, in modo da avere una stima della potenzache il processo ha nella corrispondente banda:

Page 30: Elaborazione numerica dei segnali - users.libero.itusers.libero.it/sandry/download/DSPdownload/DSP_06a.pdf · una frequenza di campionamento f S che rispetti il teorema di campionamento

Appunti di “Elaborazione numerica dei segnali” - Capitolo 6

Autore: Sandro Petrizzelli30

Hk(f)

Bk

fk

yk(n)realizzazionedel processo

x(n) stima dellapotenza in Bk

∑−

=

=1N

0n

2kk )n(y

N

1S

E’ evidente che, se vogliamo ottenere un valore B per la risoluzionespettrale della stima, dovremo usare un banco di filtri ciascuno conbanda passante pari proprio a B ; in questo modo, ciascun filtro consente una stimadello spettro di potenza del segnale attorno ad una frequenza qualsiasi e la risoluzione spettrale nonviene a dipendere dalla frequenza.

La stima che ciascun ramo del banco di filtri tira fuori è caratterizzata da una varianza che, nelcaso di campioni con distribuzione gaussiana, abbiamo visto essere

N

22

Skσ=σ

dove [ ] 2kSE σ= .

D’altra parte, abbiamo anche visto che dobbiamo prendere k

tot B

2NT ⋅≥ , ossia che 2N/Bk è il più

piccolo intervallo di tempo necessario ad avere una risoluzione spettrale Bk. Ponendo alloraN=TtotB/2, nell’ipotesi che i filtri passa-banda abbiano tutti la stessa banda, scriviamo che

2

tottot

2

tot

2

tot

2

S2

T

1T

4

BT

4

2/BT

2k

σ=σ=σ=σ=σ

In base a questa relazione, se aumentiamo Ttot, la dispersione kS

σ della stima rimane invariata in

quanto diminuisce contemporaneamente la risoluzione spettrale B.In altre parole, con il metodo del periodogramma noi raggiungiamo la

migliore risoluzione spettrale possibile, ma rinunciamo a migliorarela dispersione della nostra stima.

POLARIZZAZIONE DELLA STIMA SPETTRALE

Abbiamo considerato, nei paragrafi precedenti, il seguente stimatore per la funzione diautocorrelazione γXX(τ) del processo in esame:

0ki

Nkiper 0k)con x(i )ki(x)i(x

N

1)k(r

1N

0iXX

<+∀>+∀

=++= ∑−

=

Abbiamo osservato anche che questo stimatore è polarizzato, nel senso che la sua media risultadiversa da γXX(k), ossia da ciò che vogliamo stimare:

Page 31: Elaborazione numerica dei segnali - users.libero.itusers.libero.it/sandry/download/DSPdownload/DSP_06a.pdf · una frequenza di campionamento f S che rispetti il teorema di campionamento

Stima spettrale (parte I) - Metodi non parametrici

Autore: Sandro Petrizzelli31

[ ] )k(N

m1)k(rE XXXX γ

−=

Se considerassimo la trasformata di Fourier di rXX(k), ossia PXX(f) (cioè il periodogramma) e necalcolassimo la media, otterremmo lo stesso risultato, ossia che anche PXX(f) è uno stimatorepolarizzato (con la stessa polarizzazione di rXX(k)) per la densità spettrale di potenza del processo.

Possiamo chiederci a cosa sia dovuta questa polarizzazione. La risposta è che lapolarizzazione, sia nel tempo sia in frequenza, dipende dallafinestratura dei dati.

Per rendercene conto, consideriamo una sequenza z(n) di durata illimitata e effettuiamo unafinestratura, tramite una finestra generica w(n), su un intervallo di durata N:

1-Nn0 )n(w)n(z)n(x ≤≤⋅=

In pratica, z(n) è la realizzazione, di durata teoricamente infinita, del nostro processo e noi nestiamo osservando solo una parte, cioè x(n).

Consideriamo adesso la funzione di autocorrelazione di x(n) così come definita nei paragrafiprecedenti e calcoliamone la media:

[ ] [ ]∑∑∑−

=

=

=

+=

+=→

−≤≤

+=1N

0i

1N

0iXX

1N

0i

*XX )ki(x)i(xE

N

1)ki(x)i(x

N

1E)k(rE

1Nk0

)ki(x)i(xN

1)k(r

Sostituendo l’espressione di x(n), otteniamo che

[ ] [ ]∑−

=

++=1N

0iXX )ki(w)ki(z)i(w)i(zE

N

1)k(rE

I campioni w(n) della finestra sono deterministici (anche se per il momento non fissati), per cuipossiamo portarli fuori dalla media:

[ ] [ ] ∑∑∑−

=

=

=

+=+=++=1N

0iZZ

1N

0iZZ

1N

0iXX )ki(w)i(w)k(r

N

1)k(r)ki(w)i(w

N

1)ki(z)zikE)ki(w)i(w

N

1)k(rE

La sommatoria rimasta, moltiplicata per 1/N, non è altro che la funzione di autocorrelazione dellafinestra, per cui concludiamo che

[ ] )k(r)k(r)k(rE WWZZXX ⋅=

Questa espressione mostra dunque che il valore atteso di rXX(k) è pari alprodotto della funzione di autocorrelazione della sequenza originaleper quella della finestra. E’ dunque rWW(k), cioè la finestratura, che determina lapolarizzazione dei dati, cioè il fatto che [ ])k(rE XX differisca da rZZ(k).

Serviamoci ancora una volta di un esempio, considerando il caso semplice di una finestra w(n) ditipo rettangolare. Le figure seguenti mostrano quello che succede:

Page 32: Elaborazione numerica dei segnali - users.libero.itusers.libero.it/sandry/download/DSPdownload/DSP_06a.pdf · una frequenza di campionamento f S che rispetti il teorema di campionamento

Appunti di “Elaborazione numerica dei segnali” - Capitolo 6

Autore: Sandro Petrizzelli32

Esito della finestratura di una sequenza illimitata: nella figura (a) è rappresentata la sequenzaillimitata z(n); nella figura (b) è mostrata invece la finestra rettangolare w(n); l’esito dellamoltiplicazione tra le due sequenze è in figura (c); nella figura (d) è indicata la funzione di

autocorrelazione della sequenza z(n) considerata, mentre nella figura (e) è riportata la funzione diautocorrelazione di w(n), che notoriamente è di tipo triangolare; infine, nella figura (f) è riportata la

funzione di autocorrelazione di x(n)

Notiamo dunque che l’effetto della finestratura è quello per cui l’autocorrelazione vera rZZ(k) èmoltiplicata per una funzione triangolare rWW(k) che va a 0 dopo N campioni.

Una analisi assolutamente analoga si potrebbe fare nel dominio della frequenza considerando lafunzione PXX(f) e calcolandone la media:

[ ] [ ] [ ]

∑∑∑−

−=

π−

−=

π−−

−=

π−−

−=

π−

γ

−=

===

=

1N

N1m

fmT2jXX

1N

N1m

fmT2jXX

1N

N1m

fmT2jXX

1N

N1m

fmT2jXXXX

e)m(N

m1

e)m(rEe)m(rEe)m(rE)f(PE

La polarizzazione è identica a quella nel dominio del tempo e può essere ancora una voltainterpretata alla luce di una finestratura nel tempo, ossia di una convoluzione in frequenza tra ladensità spettrale vera del processo e lo spettro della finestra utilizzata. L’effetto di questaconvoluzione, come si è già visto in precedenza, è nella attenuazione della densità spettrale vera enei fenomeni di leakage, cioè di dispersione dell’energia del segnale al di fuori della banda in cui ilsegnale è definito. Per minimizzare questo tipo di inconvenienti, è opportuno scegliere finestre i cuilobi laterali, in frequenza, siano più bassi possibile. D’altra parte, questo tipo di finestre, purcontribuendo a ridurre la polarizzazione (dato appunto che i lobi laterali hanno ampiezza inferiorerispetto a quelli della trasformata della finestra rettangolare), riducono anche la risoluzione spettrale,dato che il lobo centrale della trasformata della finestra avrà larghezza maggiore che non quello

Page 33: Elaborazione numerica dei segnali - users.libero.itusers.libero.it/sandry/download/DSPdownload/DSP_06a.pdf · una frequenza di campionamento f S che rispetti il teorema di campionamento

Stima spettrale (parte I) - Metodi non parametrici

Autore: Sandro Petrizzelli33

corrispondente alla finestra rettangolare di uguale durata. Non solo, ma l’uso di finestre nonrettangolari, come si vedrà nei paragrafi successivi, porta a pesare i campioni di z(n) in modo diverso(alcuni hanno peso elevato ed altri peso praticamente nullo), il che contribuisce ad accrescere ladispersione della stima.

MODIFICA DEL PERIODOGRAMMA

Abbiamo in precedenza verificato che, aumentando la durata Ttot=NT della finestra diosservazione (tramite un aumento del numero N di campioni), la dispersione della stima ottenuta colperiodogramma rimane invariata, ma migliora (si riduce) la risoluzione spettrale. Allora, possiamopensare di procedere come segue: allunghiamo la finestra di osservazione, ad esempio, a 2NT, mapoi, anziché considerare una sequenza lunga 2N, consideriamo due sequenze lunghe N; da tali duesequenze otteniamo due periodogrammi, ossia due stime della densità spettrale di potenza, chesaranno distinti e indipendenti (dato che fanno riferimento a misure diverse e temporalmentedisgiunte dello stesso processo):

x(n) 43421x0(n) di lunghezza N

43421x1(n) di lunghezza N

X(k) 43421P0(k)

43421P1(k)

Così facendo, a pari risoluzione spettrale19, abbiamo due misure indipendenti. Se ne facciamo lamedia (il che significa che, per ognuna delle frequenze di campionamento, mediamo i valori fornitidai due periodogrammi), abbiamo evidentemente dimezzato la varianza della stima.

In generale, dunque, per ottenere una riduzione della fluttuazionestatistica della stima, dobbiamo effettuare delle misure ripetutestatisticamente indipendenti che poi vanno mediate. L’indipendenzagarantisce che questa operazione di media comporti una riduzione della varianza.

Una volta fissato il numero di campioni della sequenza di partenza, si tratta di decidere qualecompromesso si vuole ottenere tra la risoluzione spettrale e la dispersione della stima: dividendo lasequenza in più pezzi, la dispersione diminuisce, ma la risoluzione spettrale aumenta (quindipeggiora).

19 I due periodogrammi hanno la stessa risoluzione spettrale del periodogramma che si otterrebbe con un unica sequenza di

lunghezza N e questo perché la risoluzione spettrale, nel metodo del periodogramma, coincide praticamente con l’inverso delladurata della finestra di osservazione.

Page 34: Elaborazione numerica dei segnali - users.libero.itusers.libero.it/sandry/download/DSPdownload/DSP_06a.pdf · una frequenza di campionamento f S che rispetti il teorema di campionamento

Appunti di “Elaborazione numerica dei segnali” - Capitolo 6

Autore: Sandro Petrizzelli34

Media di periodogrammi di sottosequenze: metodo di Bartlett

Un altro modo per diminuire la dispersione della stima della densità spettrale di potenza è quellodi suddividere l’intervallo Ttot=NT di osservazione in L sottointervalli disgiunti (non aventi cioèistanti in comune) e cioè quindi di dividere la sequenza totale in L sottosequenze, ciascuna dilunghezza M=N/L:

x(n) 43421x0(n) di lunghezza M

43421 43421x1(n) di lunghezza M xL-1(n) di lunghezza M

X(k) 43421X0(k) di lunghezza M

43421 43421X1(n) di lunghezza M XL-1(n) di lunghezza M

..........

..........

Di ogni sottosequenza si calcola la DFT e successivamente il periodogramma (basta dividere laDFT per N). Si ottengono in tal modo L periodogrammi: si fa una media di tali periodogrammi (ossiauna media dei quadrati dei campioni relativi a ciascuna frequenza) e quello che viene fuori è la stimadella densità spettrale di potenza del processo.

Vediamo a livello analitico.Partendo da x(n), di lunghezza N, il primo passo è quello di ottenere L sottosequenze ciascuna di

lunghezza M:

1-M0,1,.....,n

1-L0,1,...,i )iMn(x)n(x i

==+=

Con ciascuna sequenza si costruisce il periodogramma:

1-L0,1,...,i )f(XM

1e)n(x

M

1)f(P

221M

0n

fn2ji

)i(XX === ∑

=

π−

La corrispondente versione campionata, ottenuta cioè tramite la DFT, è

1-L0,1,...,i NT

kX

M

1e)n(x

M

1

NT

kP

221M

0n

nNT

k2j

i)i(

XX =

==

=

π−

Noti gli L periodogrammi, se ne fa una media, il che significa che la stima finale risulta essere

∑∑∑−

=

=

=

=

=

=

1L

0i

21L

0i

21L

0i

)i(XXXX NT

kX

ML

1

NT

kX

M

1

L

1

NT

kP

L

1

NT

kP

E’ immediato verificare che la media del vettore PXX(k) appena costruito coincide con la mediadel generico periodogramma:

( )[ ] ( ) ( )[ ] ( )[ ] ( )[ ]kPEkPELL

1kPE

L

1kP

L

1EkPE )i(

XX)i(

XX

1L

0i

)i(XX

1L

0i

)i(XXXX =⋅⋅==

= ∑∑

=

=

Page 35: Elaborazione numerica dei segnali - users.libero.itusers.libero.it/sandry/download/DSPdownload/DSP_06a.pdf · una frequenza di campionamento f S che rispetti il teorema di campionamento

Stima spettrale (parte I) - Metodi non parametrici

Autore: Sandro Petrizzelli35

Per verificare che la stima sia migliorata, dobbiamo calcolare la varianza di PXX(k): si ha che

( )[ ] ( )[ ] ( )[ ] ( )[ ]kPEL

1kPVarL

L

1kPVar

L

1kPVar )i(

XX)i(

XX2

1L

0i

)i(XX2XX =⋅⋅== ∑

=

Abbiamo dunque conservato la media, ma abbiamo ridotto la varianza di un fattore L. La stima èdunque più precisa.

Metodo WOSA (finestre a coseno rialzato sovrapposte)

Nel metodo visto nel paragrafo precedente, noi scegliamo di suddividere x(n) in L sottosequenze,per ciascuna delle quali calcoliamo il periodogramma, dopo di che mediamo i vari periodogrammi edotteniamo la stima ricercata. In questo modo, si è visto che la varianza della stima diminuisce,mentre la media rimane inalterata. Nonostante questo apparente miglioramento, ci sono comunquedei problemi, legati al fatto che la polarizzazione aumenta.Sappiamo, infatti, che il periodogramma è uno stimatore

polarizzato per la densità spettrale di energia, nel senso che ilsuo valore atteso differisce dal valore vero. Nel momento in cuimediamo gli L periodogrammi brevi ottenuti troncando x(n) confinestre corte, la polarizzazione della stima peggiora, a causa deitroncamenti e delle periodicizzazioni più fitte. Questi troncamentideterminano eccessive infiltrazioni spettrali (spectral leaks), pari a quelle che si ottengono, peresempio, quando si tronca una sinusoide, nei tempi, mediante una finestra rettangolare: il prodottodella sinusoide per il rettangolo equivale, in frequenza, alla convoluzione tra un impulso (spettrodella sinusoide) ed un sin(f)/f (spettro del rettangolo) ed il risultato è quindi un sin(f)/f posizionato acavallo dell’impulso. Si passa cioè da uno spettro impulsivo (idealmente) ad uno spettro del tiposin(f)/f, che presenta cioè componenti spettrali anche lontane dalla frequenza della sinusoide dipartenza.

Per ridurre al minimo questo tipo di problema, è necessario usare, quando si formano lesottosequenze nel tempo, delle finestre abbastanza “dolci”, senza cioè transizioni nette. Tipicamente,è opportuno usare delle finestre a coseno rialzato.

Ovviamente, usando finestre più dolci di quella rettangolare paghiamo un prezzo non indifferente:dato l’andamento non uniforme delle finestre, abbiamo campioni che contribuiscono maggiormentealla stima (sono quelli in corrispondenza dei valori maggiori delle finestre) e campioni checontribuiscono solo in minima misura (quelli in corrispondenza dei valori minori delle finestre). Inaltre parole, pesiamo poco alcuni dei dati a nostra disposizione, ossia evitiamo di utilizzarli perridurre le fluttuazioni statistiche.

Viene allora da pensare di utilizzare, in ciascuna misura, anche i campioni che sono intervenutisolo in minima parte nella misura precedente. In altre parole, si dispongono le finestre non più suintervalli adiacenti e disgiunti, ma su intervalli in sovrapposizione parziale (window overlap),come indicato nella figura seguente:

Page 36: Elaborazione numerica dei segnali - users.libero.itusers.libero.it/sandry/download/DSPdownload/DSP_06a.pdf · una frequenza di campionamento f S che rispetti il teorema di campionamento

Appunti di “Elaborazione numerica dei segnali” - Capitolo 6

Autore: Sandro Petrizzelli36

Overlap di finestre (temporali) a coseno rialzato: ciascuna finestra include una parte dei campioniinclusi dalla precedente e, in particolare, quei campioni che, nella misura precedente, sono stati pesati di

meno.

In questo modo, si sfrutta al massimo tutta l’informazione presente nella sequenza dei dati, fermorestando il concetto che è possibile migliorare la risoluzione spettrale pagando con una minorediminuzione delle fluttuazioni statistiche.

Questo è il metodo WOSA (Window Overlap Spectral Analisys).

Smussamento dello spettro: metodo di Blackman-Tukey

Entrambi i metodi visti in precedenza si basano su una opportuna finestratura dei dati, ossia dellasequenza x(n) rappresentativa dell’unica realizzazione a nostra disposizione. C’è anche unaltro modo di procedere, che consiste nel finestrare non più x(n),ma la sua funzione di autocorrelazione. Vediamo di che si tratta.

Data la sequenza x(n), abbiamo visto che la funzione

0ki

Nkiper 0k)con x(i )ki(x)i(x

N

1)k(r

1N

0i

*XX

<+∀>+∀

=++= ∑−

=

può essere usata come stimatore (polarizzato) per la funzione di autocorrelazione γXX(τ) del processoin esame. Questa funzione fornisce una stima di γXX(τ) fino a ritardi τ≤N-1:

∑∑

=

=

=

=

−+=→−=

+=→=

==→=

1N

0i

*XX

1N

0i

*XX

1N

0i

21N

0i

*XX

)1Ni(x)i(xN

1)N(r1Nk

...

)1i(x)i(xN

1)1(r1k

)i(xN

1)i(x)i(x

N

1)0(r0k

Page 37: Elaborazione numerica dei segnali - users.libero.itusers.libero.it/sandry/download/DSPdownload/DSP_06a.pdf · una frequenza di campionamento f S che rispetti il teorema di campionamento

Stima spettrale (parte I) - Metodi non parametrici

Autore: Sandro Petrizzelli37

Si nota, però, che, quando si considera un valore piccolo di k, intervengono nella sommatoria varicampioni, mentre invece, all’aumentare di k, il numero dei campioni coinvolti diminuisce: adesempio, per k=N-1, abbiamo che

)1N(x)0(x

)2N2(x)1N(x...)1N(x)2(x)N(x)1(x)1N(x)0(x)1Ni(x)i(xN

1)1N(r

*

****1N

0i

*XX

−=

=−−+++++−=−+=− ∑−

=

Come si vede, essendo x(n) di lunghezza finita N, c’è un solo termine presente della sommatoria.Quindi, per ritardi grandi, l’autocorrelazione γXX(τ) viene stimata utilizzando ben pochi campioni

(un solo campione se τ=N-1) e questo comporta che la stima abbia grande dispersione. Diventa alloraopportuno modificare il periodogramma per ridurre la dispersione ma anche la risoluzione spettrale.Un possibile modo per ridurre sia la risoluzione spettrale sia la

dispersione è quello di finestrare la funzione di autocorrelazionerXX(k) prima di effettuarne la trasformata di Fourier. In pratica, quindi,anziché calcolare direttamente la trasformata di rXX(k), effettuiamo prima la finestratura di rXX(k),mediante una opportuna funzione finestra w(k), e poi trasformiamo, in modo che lo stimatoredella densità spettrale di potenza risulti essere

∑−

−=

π−=1N

N1k

kf2jXX

BTXX e)k(w)k(r)f(P

L’apice BT che abbiamo usato in quest’ultima simbologia evidenzia che questo metodo èdovuto a Blackman e Tukey.

Dobbiamo capire come è fatta la funzione finestra e perché risulta utile.In primo luogo, dovendo finestrare una funzione, rXX(k), di lunghezza 2N-1, la durata M della

finestra non sarà evidentemente superiore ad 2N-1. La sua durata sarà 2N-1.Un requisito opportuno di cui deve godere w(k) è quello di garantire che )f(PBT

XX risulti positivo intutte le frequenze20. Dato che la moltiplicazione tra rXX(k) e w(k) equivale alla convoluzione diPXX(f) (periodogramma) e W(f) (trasformata di w(k)), dobbiamo dunque garantire che taleconvoluzione dia luogo ad uno spettro )f(PBT

XX sempre positivo.

∫−

αα−α=2/f

2/f

XXBTXX

C

C

d)f(W)(P)f(P

La convoluzione tra il periodogramma e lo spettro della finestra determina chiaramente unosmussamento del periodogramma ed è per questo che si parla anche di metodo di smussamentodello spettro. Uno smussamento equivale ad una regolarizzazione delladensità spettrale stimata il che si traduce in una riduzione dellarisoluzione spettrale della stima stessa e questo è il grossovantaggio di questo metodo.

20 Questo problema non si presentava quando finestravamo direttamente i campioni della realizzazione: in quel caso, infatti, dopo la

finestratura ed il calcolo della DFT, effettuavamo sempre il modulo quadro e questo ci garantiva di ottenere solo valori positiviper la stima della densità di potenza.

Page 38: Elaborazione numerica dei segnali - users.libero.itusers.libero.it/sandry/download/DSPdownload/DSP_06a.pdf · una frequenza di campionamento f S che rispetti il teorema di campionamento

Appunti di “Elaborazione numerica dei segnali” - Capitolo 6

Autore: Sandro Petrizzelli38

Un altro requisito da richiedere è che )f(PBTXX risulti essere una funzione reale e questo accade se

w(k) è una funzione reale e soprattutto pari, cioè simmetrica rispetto all’origine.In generale, si usa una finestra triangolare (nei tempi) di durata

esattamente pari alla riduzione della risoluzione spettrale che sivuole ottenere. Cerchiamo allora di capire, qualitativamente, ciò che accade.

Una finestra triangolare nei tempi ha uno spettro che è del tipo sin2(f)/f2, per cui presenta un loboprincipale e dei lobi secondari molto meno pronunciati21. Se, ad esempio, la finestra triangolare va a0 per τmax=NT/10, come spesso viene consigliato, il lobo principale dello spettro della finestra ha ilprimo zero in ±10/NT; ciò significa che tale lobo combina insieme i contenuti spettrali su unintervallo di ampiezza 20/NT. In realtà, facendo i conti con maggiore precisione, si trova chel’intervallo i cui contenuti spettrali vengono combinati è largo circa 15/NT.

Quindi, l’effetto della convoluzione è una media pesata di contenuti spettrali adiacenti, il cheporta alla regolarizzazione della densità spettrale stimata e quindi alla diminuzione della dispersione.

RIEPILOGO E CONCLUSIONI

Lo scopo che ci siamo posti in questi paragrafi è la stima della densità spettrale di potenza (o, ciòche è lo stesso, della funzione di autocorrelazione) di un processo stazionario. Dato che possiamoaccedere ad una sola realizzazione del processo, per cui non possiamo calcolare la media statisticadel processo, dobbiamo necessariamente ipotizzare l’ergodicità oltre alla stazionarietà. Questaipotesi consente di sostituire alla media statistica la media temporale.

Per ottenere una buona stima, sarebbe necessario avere a disposizione una sequenza di campioni,della realizzazione del processo, la cui lunghezza tenda all’infinito, ma questo non è ovviamentepossibile. Avendo necessariamente a che fare con una sequenza di lunghezza finita, bisogna perciòrassegnarsi ad avere una stima che, per quanto buona possa essere, avrà sempre una propria varianzanon nulla.

Appurato questo, si possono evidentemente usare tutti i campioni a disposizione per calcolareun’unica stima. In questo caso, si considera un segmento temporale che è il più lungo possibile.

Il fatto di usare necessariamente una versione finestrata della realizzazione comporta che lospettro di potenza fornito dalla stima sia la convoluzione dello spettro vero con lo spettro di potenzadella finestra. Ad esempio, se la finestra scelta è rettangolare, la sua trasformata è del tipo sin(f)/f equindi lo spettro di potenza è del tipo sin2f/f2.

In generale, si procede richiedendo, a priori, che la varianza della stima ottenuta risulti inferiore,con una certa probabilità, ad un valore massimo tollerabile.

Dato che il numero dei campioni è fissato a priori, l’unico modo per migliorare (cioè ridurre) ladispersione è rinunciare alla migliore risoluzione spettrale. Si deve cioè procedere a segmentare lasequenza di N campioni in L sottosequenze di durata N/L. Così facendo, la risoluzione spettralepeggiora (in quanto ciascuna sottosequenza è ottenuta con una finestra che, essendo di durataminore, ha il lobo centrale del sin(f)/f più largo), ma si ottengono delle misure ripetute. Se questemisure sono statisticamente indipendenti, è possibile farne una media, avendo la garanzia che lavarianza di questa media è la media delle varianze. In tal modo, la stima è diventata più accurata diun fattore L.

21 Notiamo subito che la finestra triangolare ha il pregio di avere uno spettro sempre positivo, il che risponde al requisito, prima

citato, per cui la convoluzione di tale spettro con quello del processo debba risultare positivo. Questo è proprio il motivo per cuiBlackman e Tukey pensarono ad una finestra triangolare.

Page 39: Elaborazione numerica dei segnali - users.libero.itusers.libero.it/sandry/download/DSPdownload/DSP_06a.pdf · una frequenza di campionamento f S che rispetti il teorema di campionamento

Stima spettrale (parte I) - Metodi non parametrici

Autore: Sandro Petrizzelli39

A livello pratico, quindi, si procede nel modo seguente:

• data la sequenza x(n) di partenza, la si spezza in L sottosequenze, per ciascuna delle quali sicalcola la DFT;

• si prende poi il modulo quadro delle varie DFT e lo si divide per M (durata della genericasottosequenza): così facendo si ottengono L periodogrammi indipendenti tra loro;

• infine si mediano i periodogrammi, per ottenere la stima conclusiva.

Autore: SANDRO PETRIZZELLIe-mail: [email protected]

sito personale: http://users.iol.it/sandrysuccursale: http://digilander.iol.it/sandry1