Laboratorio per l'Elaborazione MultiMediale Lezione 7 - La...

IntroduzioneICA

Altre Tecniche di Separazione

Laboratorio per l’ElaborazioneMultiMediale

Lezione 7 -La Separazione delle Sorgenti

Prof. Michele Scarpiniti

Dipartimento di Ingegneria dell’Informazione, Elettronica e Telecomunicazioni“Sapienza” Universita di Roma

http://ispac.diet.uniroma1.it/scarpiniti/index.htm

[email protected]

M. Scarpiniti Laboratorio per l’Elaborazione MultiMediale Lezione 7 - La Separazione delle Sorgenti 1 / 64

http://ispac.diet.uniroma1.it/scarpiniti/index.htm

IntroduzioneICA


1 IntroduzioneIntroduzioneTecniche di Pre-filtraggio

2 ICAIntroduzioneIl FastICAInfoMax

3 Altre Tecniche di SeparazioneBSS Sotto-determinata ed Estrazione di SorgentiSeparazione in Ambiente ConvolutivoSeparazione in Ambiente Non Lineare


IntroduzioneICA


IntroduzioneTecniche di Pre-filtraggio

Introduzione

Introduzione


IntroduzioneICA



Introduzione

L’obbiettivo di questa lezione e di illustrare alcune tecniche statistiche, notecome Analisi delle Componenti Indipendenti o ICA (Independent Compo-nent Analysis), utilizzate in algoritmi per la Separazione Cieca di Sorgen-ti o Blind Source Separation (BSS). Il termine “cieco” (blind) significaessenzialmente “senza informazioni a priori”.

In questa lezione verranno descritti ed esaminati gli algoritmi piu noti diseparazione (cieca) di sorgenti.


IntroduzioneICA



Il modello

Il modello piu semplice e quello lineare ed istantaneo: osserviamo N variabili aleato-rie x1, . . . , xN , che sono modellate come combinazione lineare delle N variabili aleatories1, . . . , sN :

x = As (1)

dove gli aij sono i coefficienti della matrice A.

Lo scopo degli algoritmi di Separazione Cieca di Sorgenti o BSS e di identificare unamatrice di de-mixing W, tale che le componenti del vettore di uscita

u = Wx (2)

siano quanto piu possibile statisticamente indipendenti .


IntroduzioneICA



Il modello lineare

Dando uno sguardo piu in profondita al modello lineare (ed istantaneo),possiamo esplicitamente porre:

x1 = a11s1 + a12s2 + . . .+ a1NsNx2 = a21s1 + a22s2 + . . .+ a2NsN

...xN = aN1s1 + aN2s2 + . . .+ aNNsN

Si desidera un modello di de-mixing del tutto equivalente.M. Scarpiniti Laboratorio per l’Elaborazione MultiMediale Lezione 7 - La Separazione delle Sorgenti 6 / 64

IntroduzioneICA



Eliminare la media

L’operazione piu semplice di pre-processing e l’eliminazione della media ocentramento dei dati: dato il vettore dei dati x e possibile sottrargli lamedia m = E {x}, tale che il nuovo vettore x sia a media nulla:

x = x−m

Tale pre-processamento e effettuato per semplificare molti algoritmi e nonperche la media non possa essere stimata. Dopo la stima dei parametriincogniti con i dati centrati, e possibile completare la stima sommando lamedia m precedentemente stimata.


IntroduzioneICA



Eliminare la media

In Matlab, la media puo essere rimossa semplicemente con il comando mean.Se x e un vettore di dati, scrivo semplicemente:

xn = x - mean(x);

Se invece x e una matrice, e i dati sono organizzati per righe:

m = mean(x,2);

xn = x - m*ones(1,length(x));


IntroduzioneICA



Normalizzazione

In generale i dati possono essere normalizzati in un intervallo [−xmax xmax],ma in pratica e molto utilizzata la normalizzazione che rende i dati a me-dia nulla e varianza unitaria. La normalizzazione, cioe, dopo il classicocentramento, e fatta rispetto alla deviazione standard σ dei dati:

xi =xi − µiσi

dove µi indica la media dell’i-esimo segnale xi , e σi e la deviazione standarddi xi .

Nel caso di semplice normalizzazione, per esempio dividendo gli elementi delvettore dei dati per il suo valore massimo, si ha che i nuovi dati risultanoappartenere all’intervallo [−1 1]. Inoltre, se esiste un grande divario tragli “individui” del vettore, la dinamica puo essere notevolmente ridottaapplicando una trasformazione logaritmica.


IntroduzioneICA



Normalizzazione

In Matlab la normalizzazione nell’intervallo [−1 1] e fatta attraverso:

m = max(x,[] ,2);

for i=1:N,

xn(i,:) = x(i,:)/m(i);

end

Nel caso di normalizzazione rispetto alla varianza (energia) ottengo:

m = mean(x,2);

ds = std(x,0,2);

for i=1:N,

xn(i,:) = (x(i,:)-m(i))/ds(i);

end


IntroduzioneICA



Sbiancamento

Un ulteriore strategia di pre-processing consiste nello sbiancamento dellevariabili osservate. Questo significa che prima di applicare un algoritmo diBSS e dopo la rimozione della media, trasformiamo linearmente il vettoredelle osservazioni x in un nuovo vettore x tale che sia bianco. In altre parole,la matrice di covarianza di x e la matrice identita: Rxx = E {xx} = I.

La trasformazione di sbiancamento e sempre possibile. Un metodo popolareper ottenerla si basa sulla decomposizione autovalore-autovettore o eigen-value decomposition (EVD) della matrice di covarianza E

{xxT}

= EDET ,dove E e una matrice ortogonale che contiene gli autovettori per colon-na, e D e una matrice diagonale che contiene gli autovalori sulla diagonaleprincipale:

x = ED−1/2ETx


IntroduzioneICA



Sbiancamento

Lo sbiancamento trasforma la matrice di mixing A nella nuova matrice A:

x = ED−1/2ETx = ED−1/2ETAs = As

L’utilita dello sbiancamento risiede nel fatto che la nuova matrice A e ortogonale. Ciopuo essere visto da

Rxx = E {xx} = AE{

ssT}

AT = AAT

= I

E’ evidente che lo sbiancamento riduce il numero di parametri da stimare. Invece di doverstimare gli n2 parametri della matrice A, dobbiamo semplicemente stimare la matriceortogonale A che contiene solo n(n − 1)/2 parametri. In questo senso si dice che losbiancamento risolve mezzo problema.


IntroduzioneICA



Sbiancamento

In Matlab possiamo effettuare

R = cov(x.’);

[E,D] = eig(R);

xn = E*sqrt(inv(D))*E.’*x;

Riportiamo un esempio su un vettore x ottenuto dal vettore s = rand(2,1000),a cui e stata applicata la matrice A = [2, 3; −1, 4]


IntroduzioneICA



Altre Tecniche di Pre-processing

Se i dati sono segnali temporali e possibile utilizzare anche un filtraggiopassa-banda. Se filtriamo linearmente il segnale osservato xi (t) ottenia-mo un nuovo segnale x∗i (t), che puo essere espresso attraverso l’operatorelineare:

x∗ = Mx

Piu in generale, l’operazione di pre-filtraggio puo essere fatto in modo nonlineare:

x∗ = M(x)


IntroduzioneICA


IntroduzioneIl FastICAInfoMax

ICA

Independent Component Analysis


IntroduzioneICA



L’indipendenza statistica

Il problema che ci prefiggiamo di risolvere e di tipo blind, cioe senza in-formazioni a priori . L’unica cosa che e possibile affermare e che i segnaliprovengono da processi distinti e quindi sono statisticamente indipendenti.Allora un possibile modo di agire e quello di rintracciare nelle misture le com-ponenti indipendenti che possono essere buoni candidati per le stime dellesorgenti originali. Queste considerazioni giustificano l’adozione dell’Analisidelle componenti indipendenti o ICA, che e un insieme di tecniche stati-stiche utilizzate per “rintracciare” componenti indipendenti (o IC) in uninsieme di dati.


IntroduzioneICA



Limiti dell’ICA

Affinche l’ICA possa stimare efficientemente il modello di de-mixing delnostro problema, bisogna introdurre alcune condizioni e vincoli sulla naturadei segnali in gioco e del modello di mixing.

1 I segnali originali sono supposti statisticamente indipendenti.

2 I segnali originali, eccetto al massimo uno, non devono averedistribuzione gaussiana. Infatti la statistica di ordine superiore (al secondo)e sempre nulla per la distribuzione gaussiana, mentre tale informazione sullastatistica di ordine superiore e essenziale per l’indipendenza statistica.

3 Per semplicita assumiamo che la matrice di mixing sia quadrata. Cioe,il numero di sensori utilizzati e pari al numero di sorgenti. Assumiamo ancheche la matrice di mixing sia non singolare.


IntroduzioneICA



Ambiguita dell’ICA

E’ facile verificare che il modello ICA descritto dall’eq. (1), e soggetto alle seguentiambiguita:

1 E’ impossibile determinare le varianze (energie) delle componentiindipendenti. Il motivo e che entrambi s ed A sono incogniti, quindi ogni scalareche moltiplica una componenti si puo sempre essere cancellata dividendo larelativa colonna ai di A per lo stesso scalare, per esempio αi :

x =∑i

(1

αiai

)(αi si )

Si noti che esiste anche l’ambiguita del segno: e possibile moltiplicare lacomponente indipendente per −1 senza perturbare il modello.

2 E’ impossibile determinare l’ordine delle componenti indipendenti. Il motivo eche nuovamente entrambi s e A sono incogniti, quindi possiamo liberamentescambiare l’ordine dei termini nella somma in (1),e chiamare qualsiasi componentecome prima componente. Formalmente, e possibile sostituire al modello ICA unamatrice di permutazione P e la sua inversa x = AP−1Ps. Gli elementi di Ps sonole variabili indipendenti sj , ma in un altro ordine.


IntroduzioneICA



Perche le variabili Gaussiane sono proibite

Cerchiamo di capire il perche le variabili gaussiane sono proibite nel modello ICA.Assumiamo che la distribuzione congiunta di due sorgenti s1 e s2 sia gaussiana,cioe

p (s1, s2) =1

2πexp

(− s2

1 + s22

2

)=

1

2πexp

(−‖s‖

2

2

)ed assumiamo che al matrice di mixing A sia ortogonale, cioe A−1 = AT e det A =1. Allora la distribuzione congiunta delle osservazioni x1 e x2 e data da

p (x1, x2) =1

2πexp

(−∥∥ATx

∥∥2

2

)∣∣det AT∣∣ =

1

2πexp

(−‖x‖

2

2

)

Si vede che la matrice ortogonale non cambia la pdf: le distribuzioni delle sorgentie delle misture sono identiche. Il fenomeno per cui una matrice ortogonale non puoessere stimata da variabili gaussiane e legata alla proprieta che variabili gaussianeincorrelate sono anche indipendenti.


IntroduzioneICA



Perche le variabili Gaussiane sono proibite

Graficamente, il fenomeno puo essere osservato graficando la distribuzione dellemisture, che e la stessa della distribuzione delle due sorgenti indipendenti. La figuraseguente mostra che la pdf ha una simmetria di rotazione. Di conseguenza noncontiene alcuna informazione sulle direzioni delle colonne della matrice di mixingA: per questo motivo la matrice A non puo essere stimata.

Quindi, nel caso di variabili con distribuzione gaussiana, la matrice A nonpuo essere identificata.

Con variabili gaussiane, l’unica cosa che possiamo fare e lo sbiancamento. E’per questi motivi che e possibile accettare al massimo una sola componentecon distribuzione gaussiana.


IntroduzioneICA



Diverse soluzioni per l’ICA

Il problema ICA puo essere risolto con approcci differenti.

1 Come notato in precedenza, il concetto chiave dell’ICA e la nongaussianita delle sorgenti. Un primo modo di risolvere il problema e distimare le direzioni che massimizzano la non gaussianita. Una misuradi non gaussianita e dato dal cumulante del quarto ordine chiamatocurtosi o dagli altri cumulanti di ordine superiore. La curtosi e moltosensibile agli outliers, quindi e possibile utilizzare alternativamente,come misura di non gaussianita, la negentropy. Questi approcci sonodefiniti Higher Order Statistics o HOS.

2 Una metodologia alternativa utilizza misure di indipendenza statisticabasate sulla teoria dell’informazione (Mutua Informazione, DivergenzaKullback-Leibler, ecc.). Questi approcci sono definiti InformationTheoretic.


IntroduzioneICA



L’algoritmo FastICA

Uno degli algoritmi piu performanti per la separazione di sorgenti e il FastI-CA, denominato cosı perche ha una velocita di convergenza molto elevata.Esso appartiene alla prima classe di algoritmi, cioe massimizza la non gaus-sianita dei dati. Utilizza uno schema di ottimizzazione basato sul punto fissodella funzione costo. Una prima versione utilizzava come funzione costo lacurtosi

kurt(y) = E{y4}− 3

Una versione piu efficiente si ottiene utilizzando come funzione costo lanegentropy

J(y) = H(ygauss)− H(y)

o meglio una sua approssimazione attraverso una funzione non lineare G (·)

J(y) ∝ [E {G (y)} − E {G (ν)}]2

dove ν e una variabile gaussiana a media nulla e varianza unitaria.M. Scarpiniti Laboratorio per l’Elaborazione MultiMediale Lezione 7 - La Separazione delle Sorgenti 22 / 64

IntroduzioneICA




Lo scopo dell’algoritmo FastICA e quello di trovare le direzioni, cioe il vet-tore w, tale che la proiezione wTx massimizzi la non gaussianita misurataattraverso l’approssimazione della negentropy J(wTx). Applicando il me-todo dei moltiplicatori di Lagrange, l’ottimo di E

{G (wTx)

}soggetto al

vincolo E{

(wTx)2}

= ‖w‖2 = 1 e calcolato nei punti in cui si annulla laderivata della funzione Lagrangiana:

F = E{

xg(wTx)}

+ λw

dove g(·) e la derivata della funzione G (·). Valutando il gradiente:

∂F

∂w= E

{xxTg ′(wTx)

}+ λI = 0

Dopo semplici calcoli, si ottiene l’algoritmo:

w = E{

xg(wTx)}− E

{g ′(wTx)

}w (3)


IntroduzioneICA




Se non si e effettuato lo sbiancamento, conviene “ortogonalizzare” i vettoriw. Dopo aver stimato p− 1 componenti indipendenti, il vettore wp e orto-gonalizzato alla Gram-Schmidt, utilizzando i p− 1 vettori precedentementestimati:

wp = wp −∑p−1

j=1

(wT

p wj

)wj ,

wp =wp

‖wp‖ .

Comunque dopo la stima della p-esima componente, e necessario toglieredall’insieme delle misture x, la componente indipendente yp appena stimata.Tale procedura, denominata deflation, puo essere ottenuta semplicementenel seguente modo:

x′ = x−wy.


IntroduzioneICA




L’algoritmo FastICA puo essere schematizzato come segue:

1 Rimuovere la media.

2 Sbiancare i dati, ottenendo x.

3 Scegliere una condizione iniziale (random) w a norma unitaria.

4 Valutare w = E{

xg(wTx)}− E

{g ′(wTx)

}w, per una data scelta di g .

5 Normalizzare il vettore w = w/ ‖w‖.6 Tornare al punto 4 fino a convergenza.

7 Applicare la deflation per eliminare la stima corrente dall’insieme dellemisture e tornare al punto 3 per stimare una nuova componente.

Hyvarinen et al. hanno messo a disposizione un pacchetto Matlab che implemental’algoritmo FastICA. E’ possibile scaricare questo software dal link: http://www.cis.hut.fi/projects/ica/fastica/


http://www.cis.hut.fi/projects/ica/fastica/

http://www.cis.hut.fi/projects/ica/fastica/

IntroduzioneICA



Matlab: l’algoritmo FastICA

Il pacchetto FastICA proposto da Hyvarinen si utilizza nel seguente modo:

y = fastica(x);

[y,A,W] = fastica(x);

[y,A,W] = fastica(x,’varargin ’);

in cui x e il vettore delle misture, y e la stima delle sorgenti separate, Ala stima della matrice di mixing e W la stima della matrice di de-mixing(cioe tale che W · A = P · Λ). E’ possibile specificare diverse opzioni, conil relativo valore, tramite la stringa ’varargin’.

Le opzioni piu comuni sono:

’approach’: puo assumere i due valori:1 ’symm’: stima tutte le componenti in una volta sola;2 ’defl’: stima le componenti una dopo l’altra (opzione di default).

’numOfIC’: specifica il numero di componenti indipendenti dastimare. Di default e pari alla dimensionalita dei dati.


IntroduzioneICA




Ulteriori opzioni sono:’g’: scelta della non linearita utilizzata; e possibile scegliere tra:

1 ’pow3’: g(u) = u3 (e l’opzione di default);2 ’tanh’: g(u) = tanh(a1u);3 ’gauss’: g(u) = u exp(−a2

u2

2 );4 ’skew’: g(u) = u2.

’a1’: valore del parametro a1.

’a2’: valore del parametro a2.

’epsilon’: criterio di stop. Il valore di default e ε = 0.0001.

’maxNumIterations’: numero massimo di iterazioni primadell’arresto. Il valore di default e Nit = 1000.

’verbose’: vale ’on’ (default) o ’off’ e mostra le informazioni sulprompt dei comandi.

’displayMode’: vale ’on’ o ’off’ (default) e disegna la stima deisegnali recuperati.


IntroduzioneICA




Ulteriori opzioni comuni sono:

’firstEig’: specifica l’indice dell’autovalore piu ampio. Di default e 1.

’lastEig’: specifica l’indice dell’autovalore piu piccolo. Di default e parialla dimensionalita dei dati.

’interactivePCA’: vale ’on’ o ’off’ (default) e calcola la PCA(riduzione della dimensionalita e sbiancamento) prima di effettuare l’ICA.

’only’: puo assumere tre valori:1 ’pca’: calcola solo la PCA e restituisce le matrici degli autovettori E

e degli autovalori D;2 ’white’: calcola solo lo sbiancamento e restituisce le matrici di

whitening e di de-whitening ;3 ’all’: applica l’ICA dopo lo sbiancamento. E’ l’opzione di default.

Ad esempio:

y = fastica(x,’approach ’,’symm’,’g’,’tanh’);

y = fastica(x,’lastEig ’,10,’numOfIC ’ ,3);

y = fastica(x,’verbose ’,’off’,’displayMode ’,’on’);


IntroduzioneICA




Riportiamo un esempio di separazione di tre sorgenti: la radio, la poesia e il rumore.

[s1 ,Fs] = audioread(’noise.wav’);

s2 = audioread(’radio.wav’).’;

s3 = audioread(’poesia.wav’).’;

s = [s1.’;s2;s3];

A = [1.5 -0.4 0.5; 0.5 0.4 -0.1; 0.2 -0.6 1.3];

x = A*s;

y = fastica(x,’displayMode ’,’on’);

for i=1:3,

soundsc(y(i,:),Fs);

end


IntroduzioneICA



BSS: valutare le prestazioni

Esistono diversi modi per valutare la “bonta” di un algoritmo di separazione. I piu immediati sonoi seguenti due:

1 Rapporto Segnale Interferenza o SIR:

SIRj = 10 log

E {(|y |σ(j),j

)2}

/ E

∑k 6=j

(|y |σ(j),k

)2

indica quanto segnale utile e presente sul j-esimo canale rispetto alle altre componenticonsiderate come segnali interferenti. La separazione e tanto piu perfetta, quanto piu eelevato del valore del SIR per la j-esima stima.

2 Performance Index (PI) o Indice di Vicinanza al Modello Inverso:

PI =N∑i=1

N∑

k=1|qik |2

maxp

[∣∣qip∣∣2] − 1

+N∑

k=1

N∑i=1|qik |2

maxp

[∣∣qpk ∣∣2] − 1

indica quanto il prodotto Q = WA e simile al prodotto di una matrice di permutazione Ped una matrice diagonale Λ, ovvero Q = PΛ. La separazione e tanto piu perfetta, quantoil PI tende a zero.


IntroduzioneICA




Un’implementazione pratica del SIR e la seguente:

function SIR = evalSIR(s,Q)

N = size(s,1);

SIR = zeros(1,N);

P = permutation(Q);

for i=1:N,

[m,q] = max(P(i,:));

sn = zeros(size(s));

sn(i,:) = s(i,:);

un = Q*sn;

sig = 0;

noise = 10^ -17;

for k=1:N,

if k == q

sig = abs(sum(un(k,:) .^2));

else

noise = noise + abs(sum(un(k,:) .^2));

end

end

SIR(i) = sig./noise;

end


IntroduzioneICA




La precedente funzione puo essere riscritta in modo piu furbo. Infatti leprestazioni devono dipendere solo dal modello e non dal particolare insiemedi segnali utilizzati. Questo significa che l’indice SIR deve potersi calcolarea partire dalla matrice Q = W · A nel seguente modo:

SIRj =maxp|qjp|2

N∑k=1

|qjk |2 −maxp|qjp|2

.

Il codice necessario ad implementare tale formula e il seguente:

function SIR = evalSIR2(Q)

SIR = max(Q.^2) ./( diag(Q.’*Q).’ - max(Q.^2));


IntroduzioneICA




Un’implementazione pratica del PI e la seguente:

function PI = evalPI(Q)

N = size(Q,1);

QP = abs(Q.^2);

PI1 = sum(QP ,1)./max(QP ,[],1);

PI2 = sum(QP ,2)./max(QP ,[],2);

PI = (sum(PI1) + sum(PI2) - 2*N)/(N*(N-1));


IntroduzioneICA




Nel caso dell’esempio precedente e possibile valutare le prestazioni ottenutedall’algoritmo FastICA.

[y,B,W] = fastica(x);

Q = W*A;

SIR = 10* log10(evalSIR(Q));

PI = evalPI(Q);

Si ottengono i seguenti valori:

SIR = 29.4422 28.1170 40.0449 dB

PI = 0.0056

Entrambi gli indici indicano quindi che e stata ottenuta un’ottima separa-zione, come era stato messo in evidenza dall’ascolto.


IntroduzioneICA




Se non si dispone della funzione fastica, e possibile implementarla facilmente. Si esegueil codice seguente comprensivo di centramento e sbiancamento:

x = x - mean(x,2)*ones(1,length(x));

R = cov(x.’);

[E,D] = eig(R);

x = E*sqrt(inv(D))*E.’*x;

for k=1:N,

w = ones(N,1)/sqrt(N);

for i=1: numIteration ,

w = mean(x*(x.’*w).^3 ,2) - mean (3*(x.’*w).^2)*w;

w = w/norm(w);

end

W(:,k) = E*sqrt(D)*w;

u = w.’*x;

y(k,:) = u;

x = x - w*u;

end


IntroduzioneICA



L’algoritmo InfoMax

Un altro algoritmo molto utilizzato e l’InfoMax, proposto da Bell & Sejnowski nel1995. L’InfoMax cerca di massimizzare la mutua informazione I (y; x), tra il vettoredi ingresso x ed una sua trasformazione non lineare ed invertibile y, ottenuta come

y = h(u) = h(Wx)

dove W e la N × N matrice di de-mixing, e h(u) = [h1(u1), . . . , hN(uN)]T e uninsieme di N funzioni non lineari.

Si noti che la rete utilizzata dal-l’algoritmo di Bell & Sejnow-ski e una rete neurale a singo-lo strato. Per questo motivo lefunzioni non lineari hi (ui ) sonosolitamente chiamate funzionidi attivazione o AF.


IntroduzioneICA




L’intero sistema utilizzato per la separazione puo essere schematizzato nellaseguente figura:

Una funzione costo ottimizza l’indipendenza statistica dell’uscita y della rete, cam-biando i parametri liberi della rete (pesi della matrice, o parametri della nonlinearita).


IntroduzioneICA




Analizziamo come funziona l’algoritmo InfoMax. Poiche il mapping tra ingresso ed uscitae deterministico, massimizzare I (y, x) equivale a massimizzare l’entropia congiunta H(y).Infatti:

I (y, x) = H (y)− H ( y| x)

Differenziando l’equazione precedente rispetto ai parametri w, coinvolti nel mapping trax e y, si ottiene:

∂

∂wI (y, x) =

∂

∂wH (y) (4)

dato che H(y|x) non dipende da w. Per questo motivo l’algoritmo InfoMax e equivalentealla massimizzazione dell’entropia o ME. Ricordando che

py (y) =px (x)

|det J|

dopo qualche passaggio e possibile valutare l’entropia congiunta H(y) delle uscite come:

H (y) = H (x) + ln det W +N∑i=1

ln∣∣h′i ∣∣. (5)


IntroduzioneICA




Applicando la regola della discesa del gradiente all’eq. (5), si ottiene la regola diapprendimento relativa alla matrice W di de-mixing:

∆W =∂H (y)

∂W= W−T + ΨxT (6)

dove W−T =(W−1

)T, Ψ = [Ψ1, . . . ,ΨN ]T e Ψk = h′′k (uk)/h′k (uk).

Un problema sussiste nella scelta delle non linearita hi (ui ). In ogni caso, si puoscrivere che

H(y) = −I (u) + E

{N∑i=1

ln|h′i |

pui (ui )

}(7)

quindi se si scegliesse |h′i | = pui (ui ) (∀i) allora massimizzare l’entropia congiuntaH(y) e equivalente alla minimizzazione della mutua informazione e quindi allarisoluzione del problema della ricerca delle componenti indipendenti.


IntroduzioneICA




Come prima approssimazione e possibile utilizzare una funzione di attiva-zione di tipo squashing (a saturazione) tra 0 e 1.

Anche la funzione tanh(·) va bene in molte applicazioni, soprattutto persegnali super-gaussiani, come lo speech. In questo caso, il vettore Ψ puoessere scritto esplicitamente e la regola di apprendimento diviene:

∆W =∂H (y)

∂W= W−T − 2yxT (8)


IntroduzioneICA



L’algoritmo InfoMax : il Gradiente Naturale

La regola della discesa del gradiente stocastico lavora bene se i dati appartengonoad uno spazio Euclideo, altrimenti si ha il problema dei minimi locali , nonche unavelocita di convergenza molto bassa. Per evitare questo tipo di problemi, si applical’algoritmo di ottimizzazione nello spazio di Riemann. Amari ha dimostrato che lanuova regola (chiamata Gradiente Naturale) e semplicemente quella della discesadel gradiente moltiplicata a destra per il prodotto WTW:

∆W = η(

W−T + ΨxT)

WTW = η(

I + ΨxTWT)

W =

= η(

I + ΨuT)

W

ovvero, per hi (ui ) = tanh(ui ), si ottiene

∆W = η(I− 2yuT

)W

Il gradiente naturale e molto piu efficiente, in quanto:

1 evita il problema dei minimi locali;

2 evita l’inversione della matrice W.


IntroduzioneICA



Matlab: l’algoritmo InfoMax

L’implementazione in Matlab dell’algoritmo InfoMax non e problematica. Infatti il coredell’algoritmo, cioe l’eq. (8) puo essere scritta come segue.

N = size(x,1);

W = rand(N);


u = W*x;

y = tanh(u);

DW = inv(W.’) - 2*y*x.’;

W = W + eta*DW;

end

u = W*x;

oppure nel caso di gradiente naturale:


u = W*x;

y = tanh(u);

DW = eye(N) - 2*y*u.’;

W = W + eta*DW;

end

u = W*x;


IntroduzioneICA


BSS Sotto-determinata ed Estrazione di SorgentiSeparazione in Ambiente ConvolutivoSeparazione in Ambiente Non Lineare




IntroduzioneICA




Gli ambienti che abbiamo visto in precedenza sono molto semplici, in quantoviene considerato un modello di mixing lineare, istantaneo e con un numerodi sensori pari al numero di sorgenti (matrice di mixing quadrata). In realta,queste ipotesi possono essere messe in discussione:

1 BSS sotto-determinata: nel caso in cui ho piu sorgenti di sensori, equindi la matrice A di mixing e rettangolare;

2 BSS convolutiva: ciascuna mistura non e semplicemente unacombinazione lineare delle sorgenti, ma una loro versione filtrata:

xi (n) =N∑j=1

K−1∑k=0

aij (k)sj (n − k) .

3 BSS non lineare: le misture sono una versione non lineare del mixingprecedentemente visto:

x = f (As)


IntroduzioneICA



BSS sotto-determinata

Se il numero M di sensori e minore del numero N di sorgenti, il problema (notocome under-determined o over-complete BSS) e piu complicato. Un approcciomolto semplice e quello che utilizza un algoritmo geometrico.

Questo approccio e particolarmente performante se le sorgenti sono sparse.La sparsita e una misura di quanti valori del segnale sono nulli. Quindimaggiore e la sparsita, piu il segnale e super-gaussiano.

Se i segnali sono abbastanza sparsi , lo scatterplot delle misture evidenziaalcune direzioni “caratteristiche”. Da queste direzioni e possibile stimare lamatrice di mixing A. Riportiamo un esempio di 2 sensori ed N sorgenti:

A =

[cosα1 cosα2 · · · cosαN

sinα1 sinα2 · · · sinαN

]


IntroduzioneICA




Per aumentare la sparsita delle misture e possibile applicare delle trasformazioni tempo-frequenza, quali:

1 la Short-Time Fourier Transform (STFT):

STFT {x (t)} ≡ X (ω, τ) =

∞∫−∞

x (t)w (t − τ) e−jωtdt

2 una Trasformata Wavelet, ad esempio:

Xω (a, b) =1√a

∞∫−∞

x (t)Ψ∗(t − b

a

)dt

dove Ψ∗(t) e una funzione continua sia nel tempo che in frequenza ed e chiamatawavelet madre.


IntroduzioneICA




La figura seguente illustra la densita di probabilita di un segnale vocale,dopo l’applicazione della STFT

L’algoritmo di BSS nel caso sotto-determinato, puo essere schematizzatoin due passi:

1 Stima delle direzioni αi attraverso un algoritmo di clustering (K-means, adesempio): passo BMMR o Blind Mixing Model Recovery;

2 Ricostruzione delle sorgenti, attraverso la risoluzione del sistemasotto-determinato (non quadrato): passo BSR o Blind Source Recovery.


IntroduzioneICA



Estrazione di Sorgenti (BSE)

Se il numero di misture e molto elevato, ma siamo interessati ad un numero moltolimitato di sorgenti utili , allora e piu conveniente adoperare algoritmi di EstrazioneCieca di sorgenti o Blind Source Extraction (BSE): questi algoritmi estraggono unasorgente per iterazione.

Una singola unita di elaborazione (un neurone) e utilizzata per estrarre unasola sorgente indipendente con le proprieta statistiche specificate. L’unitasuccessiva implementa una tecnica di deflation utilizzata per rimuovere ilsegnale estratto dal resto delle misture.


IntroduzioneICA




La regola di apprendimento per l’estrazione della prima componente puo essere derivata,dopo lo sbiancamento, minimizzando la seguente funzione costo:

J1(w1) = −1

4|k(y1)| = −β

4k(y1)

dove k(y1) e la curtosi della prima uscita y1 e β = ±1 determina il segno della curtosi delsegnale estratto. Applicando la regola della discesa del gradiente, si ottiene:

∆w1 = −µ1∂J1w1

∂w1= µ1β

m4(y1)

m32(y1)

[m2(y1)

m4(y1)E{y 3

1 x1

}− E {y1x1}

]dove µ1 > 0 e il learning rate, m2(y1) = E

{|y |2}

e m4(y1) = E{|y |4}

. Il terminem4(y1)/m3

2(y1) e sempre positivo e puo essere assorbito nel learning rate. Ponendo

ϕ(y1) = βm4(y1)

m32(y1)

[m2(y1)

m4(y1)y 3

1 − y1

]si ottiene la seguente regola di apprendimento:

∆w1 = µ1ϕ(y1)x1 (9)


IntroduzioneICA




Dopo l’estrazione della prima sorgente y1 ≈ s1, possiamo applicare una proceduradi deflation che rimuove il segnale appena estratto dalle misture. Questa procedurapuo essere applicata ricorsivamente fino all’estrazione di tutti i segnali di interesse:

xj+1 = xj − wjyj j = 1, 2, . . .

dove il filtro wj puo essere ottimamente stimato, minimizzando la funzione costo

Jj (wj) = E{

xTj+1xj+1

}= E

{xTj xj

}− 2wT

j E {xjyj}+ wTj wjE

{y2j

}rispetto a wj , che porta alla semplice regola di apprendimento:

wj =E {xjyj}

E{y2j

} =E{

xjxTj}

wj

E{y2j

}Si noti che se i dati sono stati sbiancati, allora E

{y2j

}= 1 e E

{xjxTj

}= I. Di

conseguenza wj = wj , cioe i filtri di deflation coincidono con i filtri di de-mixing.


IntroduzioneICA



Matlab: l’algoritmo BSE

Implementiamo in Matlab l’algoritmo di estrazione di sorgenti o BSE:

% Centramento e Sbiancamento

x = x - mean(x,2)*ones(1,length(x));

R = cov(x.’);

[E,D] = eig(R);

x = E*sqrt(inv(D))*E.’*x;

%-----------------------------------------------------------

% Estrazione

y = zeros(size(x));

for i=1:M,

w = ones(N,1)/sqrt(N);

for t=1: numIteration ,

u = w.’*x;

k = kurtosis(u) -3;

b = sign(k);

m4 = moment(u,4);

m2 = moment(u,2);


IntroduzioneICA



Matlab: l’algoritmo BSE

phi = b*(m4./(m2.^3))*((u.^3)*m2 ’/m4 ’ - u);

Dw = eta*x*phi.’;

w = w + Dw;

w = w/norm(w);

end

%-----------------------------------------------------------

% Deflation

W(:,i) = w;

u = w.’*x;

w = (cov(x.’)*w)/var(u);

x = x - w*u;

y(i,:) = u;

fprintf(’Estrazione sorgente %d su %d di %d misture \n’,i,M

,N);

end


IntroduzioneICA



Separazione in Ambiente Convolutivo

Un esempio di problema di separazione in ambiente convolutivo e il coctail-party problem, illustrato nella seguente figura:

x (t) = A (t) ∗ s (t)STFT⇔ X (f , t) = A (f ) S (f , t)

Un problema BSS convolutivo puo essere approcciato risolvendo N problemiistantanei, uno per ogni frequency bin.


IntroduzioneICA




Esempio di spettrogramma di due sorgenti e due misture.


IntroduzioneICA




La procedura da seguire, quindi, e la seguente:1 si finestra il segnale;2 si calcola la STFT su N bin;3 si applica un algoritmo di separazione per ogni bin;4 si trasforma nel tempo.

Poiche il risultato della trasformazione in frequenza e di natura complessa(C), l’algoritmo da utilizzare deve essere in grado di lavorare con i numericomplessi, il che puo essere effettuato con semplici modifiche.

Inoltre la separazione e il risultato di N problemi diversi, ognuno dei qualipuo essere soggetto alle ambiguita di scaling e permutazione. Si capisce chequeste ambiguita diventano molto problematiche per ambienti convolutivi.Il risultato della separazione puo quindi essere distorto.

Ulteriori problemi nel dominio complesso sono dovuti alla scelta della fun-zione di attivazione complessa.


IntroduzioneICA



Separazione in Ambiente Non Lineare

Sfortunatamente il modello di mixing lineare ed istantaneo e troppo pocorealistico in molte applicazioni. Un modello piu reale e ottenuto introdu-cendo delle distorsioni non lineari nei segnali. In generale il caso non linearenon garantisce una soluzione unica. La soluzione diviene piu semplice in uncaso particolare, chiamato separazione Post Nonlinear (PNL), ben nota inletteratura. In questo caso la soluzione esiste ed e unica.

La procedura di mixing e realizzata indue passi : prima e presente un mixinglineare, quindi ogni mistura e ottenu-ta come una versione non lineare di-storta della mistura lineare. Il modelloPNL e abbastanza realistico in moltesituazioni operative.


IntroduzioneICA



BSS Non Lineare

E’ stato detto che se sul modello non lineare non diamo vincoli particolari, come lastruttura PNL, allora non e garantita l’esistenza e/o l’unicita della soluzione. Per meglioillustrare il concetto, forniamo un esempio.

Example

Date due variabili random indipendenti s1 con distribuzione uniforme [0, 2π) e s2 con

distribuzione di Rayleigh, cioe tale per cui ps2 (s2) = s2

σ22e−

s22/2 con varianza σ2

2 = 1.

Date le due trasformazioni non lineari y1 = s2 cos s1 e y2 = s2 sin s1, si puo verificare che y1 e y2

sono ancora indipendenti , ma con distribuzione Gaussiana, cosı non e possibile separarle comeconseguenza del Teorema di Darmois-Skitovich. Infatti lo Jacobiano J della trasformazione e:

det (J) = det

(∂y1∂s1

∂y1∂s2

∂y2∂s1

∂y2∂s2

)= det

(−s2 sin s1 cos s1

s2 cos s1 sin s1

)= −s2

quindi la pdf di y = [y1, y2] e espressa come

py1,y2 (y1, y2) =ps1,s2 (s1,s2)

|det(J)| = 12π

exp

(− y2

1 +y22

2

)=

=

(1√2π

exp

(− y2

12

))(1√2π

exp

(− y2

22

))≡ py1 (y1) · py2 (y2)


IntroduzioneICA



BSS Non Lineare: il modello mirror

Per risolvere il problema BSS nell’ambiente non lineare, di tipo PNL, si adopera ilseguente modello:

Questo modello e conosciuto come modello mirror , perche il modello di de-mixinge la versione speculare di quello di mixing. In formule, abbiamo:

v = Asx = F (v) , xi = fi (vi )r = G (x) , ri = gi (xi )u = Wr


IntroduzioneICA



BSS Non Lineare: il Sistema

L’intero sistema per la separazione PNL e il seguente:

L’architettura di learning deve adattare i pesi wij , le funzioni di attivazionihi (ui ) e le funzioni compensatrici gi (xi ) delle non linearita distorcenti fi (vi ).


IntroduzioneICA



La Funzione di Attivazione Flessibile: il Neurone Spline

Per ottenere una funzione adattativa da utilizzare per la funzione di attiva-zione o la non linearita compensatrice e possibile utilizzare il neurone spline.

Il neurone spline e descritto dalla seguente figura:

Nel secondo blocco i termini Mi sono le colonne della matrice M, chedipende dalla particolare base utilizzata.


IntroduzioneICA



La Funzione di Attivazione Flessibile: il Neurone Spline

Una formulazione matriciale del valore yi dell’uscita interpolata e data da:

yi = hi (u) = T ·M ·Qi

dove

T =[u3 u2 u 1

],

Qi =[

Qi Qi+1 Qi+2 Qi+3]T

e

M = 12

−1 3 −3 12 −5 4 −1−1 0 1 00 2 0 0

La matrice M determina la caratteristica dellacurva interpolante: CR-Spline o B-Spline

Per garantire una caratteristica di monotonia di tutta la funzione, bisognaimporre il seguente vincolo: Q1 < Q2 < . . . < QN .


IntroduzioneICA



Diverse basi spline

Le basi spline piu diffuse sono espresse dalle seguenti:

1 B-Spline:

M =1

6

−1 3 −3 13 −6 3 0−3 0 3 01 4 1 0

2 Bezier Spline:

M =1

6

−1 3 −3 13 −6 3 0−3 3 0 01 0 0 0

3 Catmull-Rom Spline:

M =1

2

−1 3 −3 12 −5 4 −1−1 0 1 00 2 0 0

4 Hermite Spline:

M =

2 −2 1 1−3 3 −2 −10 0 1 01 0 0 0


IntroduzioneICA



Diverse basi spline

Dato un ingresso u i valori dei parametri locali i e ν, possono essere calcolati da

z ′ = u∆u + N−2

2

z =

1 se z ′ < 1

z ′ se 1 ≤ z ′ ≤ N − 3

N − 3 se z ′ > N − 3

i = bzcν = z − i

dove ∆u e la distanza tra due punti di controllo della spline. In Matlab

z1 = u/Delta + (N/2-1);

if z1(k) <= 0

z(k) = 0;

elseif z1(k) >= (N-2)

z(k) = N-3;

else

z(k) = z1(k);

end

i = floor(z);

ni = z - i;


IntroduzioneICA



Bibliografia

A. Hyvarinen, J. Karhunen, E. Oja.

Independent Component Analysis.Wiley, 2001.

A. Hyvarinen and E. Oja.

Independent Component Analysis: algorithms and applications.Neural Networks, Vol. 13, pp. 411-430, 2000.

A. Hyvarinen2000.

Fast and Robust Fixed-Point Algorithms for Independent Component Analysis.IEEE Transaction on Neural Networks, Vol. 10, N. 3, pp. 626-634, May 1999.

P. Comon and C. Jutten.

Handbook of Blind Source Separation.Academic Press, 2010.

T.A. Davis.

MATLAB Primer.8-th Edition, CRC Press, 2010.

D.M. Smith.

Engineering Computation with MATLAB.2-nd Edition, Addison-Wesley, 2010.

A. Knight.

Basics of MATLAB and Beyond.CRC Press, 1999.


Laboratorio per l'Elaborazione MultiMediale Lezione 7 - La...

Documents

Transcript of Laboratorio per l'Elaborazione MultiMediale Lezione 7 - La...