Chiarimenti sulla funzione xcorr del Matlab

[anton89]

Salve,
ho bisogno di capire bene come funziona la xcorr del Matlab.
Ho implementato questo codice:

a = [1, 0, 1, 1, 0,1,0];
b = [1,1];
[c, lags] = xcorr(a,b);

Facendo il plot(lags,c), la funzione mi fa vedere un grafico che ha un picco laddove le due sequenze si assomigliano di più. Infatti, come potete notare, il vettore [a(3),a(4)] è esattamente uguale a b, quindi il grafico mi dice che ho un picco in posizione 3.

Ora, il punto è l'ampiezza del picco. Voglio che mi restituisca un valore che mi quantifichi quanto le due sequenze siano correlate. Cioè, se modifico il codice in questo modo:

a = [1, 0, 1, 1, 0,1,0];
b = [0.8,0.8];
[c, lags] = xcorr(a,b);

voglio che l'ampiezza del picco sia proporzionale al grado di somiglianza delle due sequenze. E' ovvio che è un problema di normalizzazione, ma non riesco a capire come si faccia.

Mi aiutate?

[Russell]

premetto che ho fatto il test usando Octave, e non con Matlab
ma girando i tuoi esempi ho nel primo caso che il picco vale 2,.... ma nel secondo caso vale 1.6
ovvero sarebbe quello che ti aspetti in pratica
Con Matlab hai una normalizzazione automatica? ... strano

[anton89]

E' questo il punto, il Matalb non mi fa la normalizzazione automatica. Vorrei capire come si normalizza. Tutto qui

[Russell]

ah... ma quindi vuoi che sia normalizzato il massimo in entrambi i casi?

beh, basta che fai:

c=c/max(abs(c(:)));

prima di fare il plot

(solo che quando il segnale è un po' rumoroso rischi di amplificare schifezze... comunque prova)

[anton89]

No, non è quello che voglio. Io voglio calcolare l'ampiezza assoluta della cross-correlazione, cioè voglio un coefficiente proporzionale al grado di correlazione fra le due sequenze.
Ad esempio, io mi aspetto che in posizione 3 la correlazione abbia ampiezza 1 nel primo esempio e abbia una ampiezza minore di 1 nel secondo esempio; infatti, nel secondo esempio b = [0.8, 0.8] che è quasi 1, ma non è 1. Quindi saranno correlate ma un po' di meno rispetto al primo esempio.

Spero di essermi spiegato meglio. Se normalizzo come dici tu, non ottengo questa informazione. Anzi, la perdo.

[Russell]

allora dividi c per l'energia di a e b
quindi qualcosa tipo

ea=sqrt(sum(a(:).^2)/length(a));
eb=sqrt(sum(b(:).^2)/length(b));
c=c/(ea*eb);

[gill90]

La correlazione tra due segnali e è una funzione del tipo

e, se i segnali sono discretizzati (), diventa


Mi pare di aver capito che poi vuoi anche cercare quanto i segnali sono somiglianti nelle ampiezze nel punto di massima somiglianza.
Supponiamo che il secondo segnale siano una porzione (che va da a ) traslata del primo segnale (tuo primo esempio), per cui non c'è nessuno scalamento nelle ampiezze tra di loro. In questo modo la correlazione nel punto massimo sarà direttamente proporzionale all'energia del secondo segnale:
.
Ma questo corrisponde anche all'energia del primo segnale troncato anch'esso in , questo perché i due segnali sono uguali nelle ampiezze (abbiamo supposto ci fosse solo una traslazione).

Se ora hai segnali scalati nelle ampiezze oltreché traslati, è evidente che, nel punto di massima correlazione, non avrai più che le due energie (quella del secondo segnale e quella del primo troncato) sono uguali, proprio perché i due segnali sono ora diversi tra loro.
Supponi di fare una cosa del genere: prendi ora un terzo segnale molto simile al secondo (), che valga:
- 1 dove è diverso da 0;
- 0 dove è uguale a 0.
In pratica è un segnale "finestra" del secondo, che identifica solo dove ci sono dei "buchi" in .
Se ora calcoli la correlazione tra il quadrato del primo () e questo terzo segnale finestra (), ottieni proprio l'energia del segnale troncato, esattamente come fatto nel caso precedente.
Questa appena trovata è una misura della correlazione se i due segnali fossero perfettamente uguali nelle ampiezze, per cui valuti il massimo per vedere a quanto corrisponde tale valore e lo utilizzi per normalizzare la tua correlazione "normale" tra e .

In pratica:
- crei il segnale "finestra":

Codice: Seleziona tutto

w = y;
w(find(w)) = 1;

- calcoli l'energia "troncata" di x:

Codice: Seleziona tutto

ct = xcorr(x.^2,w);

- ne trovi il massimo:

Codice: Seleziona tutto

Cmax = max(ct);

- calcoli la correlazione solita tra e :

Codice: Seleziona tutto

c = xcorr(x,y)

- la normalizzi per il valore trovato precedentemente:

Codice: Seleziona tutto

Cn = c/Cmax

In questo modo ti trovi una misura della similitudine dei due segnali: laddove la correlazione normalizzata è più vicina a 1, allora in tal punto i due segnali saranno tanto più simili anche nelle ampiezze.

[anton89]

Ciao Gill90,
io non so come ringraziarti, mi salvi la vita

Senti, un'ultima cosa ti volevo chiedere. Se volessi utilizzare lo stesso "trucchetto" della funzione finestra per implementare la stima della cross-correlazione attraverso le fft, come potrei fare?

In altre parole, dati due segnali x e y, calcolo la fft(x), calcolo la fft(y) e poi faccio il prodotto dell'uno con il complesso e coniugato dell'altro; infine faccio la ifft del risultato.

Come si potrebbe normalizzare in questo caso?

[gill90]

Perdonami ma non capisco. perché vorresti passare da un dominio all'altro? La cross-correlazione ha senso nel dominio di appartenenza, viene fuori un casino se trasformi troppe volte questo segnale... A cosa ti serve più nel dettaglio?

[anton89]

Ciao,
allora il punto è questo: se guardi come è implementata la xcorr del Matlab, ti accorgi che è implementata esattamente come ti ho descritto io, con le fft dei segnali.
Capisco il tuo punto di vista e ora ti spiegherò il punto. Non ho alcun controllo di quello che fa la xcorr, non mi piace come funziona, perciò la sto re-implementando a modo io, usando sempre la stessa logica della fft. E mi sono accorto che la funzione che ho implementato io è molto più efficiente in termini di risultati rispetto alla xcorr del Matlab (che ripeto, implementa la stessa logica). Ragion per cui quella funzione non la voglio utilizzare più e utilizzare la mia.

Ultimo passo per completare la realizzazione della funzione è "aggiustare" le ampiezze, cioè fare in modo che l'ampiezza della cross-correlazione sia quella che voglio ottenere, cioè compresa tra -1 e 1 in dipendenza del grado di somiglianza dei due segnali.

Sai, non riesco a spiegarmi perché la xcorr del Matlab funzioni male ma ora che mi hai fatto notare che è un casino passare da un dominio all'altro, credo di capire perché xD In ogni caso, poiché la mia funzione funziona (scusa il gioco di parole), la lascio funzionare con questa logica. Mi serve solo correggere le ampiezze.
Se mi vuoi dare qualche dritta, te ne sarò molto grato