Cancellazione Eco Acustico

Cancellazione eco acustico Roberto Patrizi

Professore: Aurelio Uncini

Roberto Patrizi Frosinone, 2 maggio 2009 e-mail: [email protected] Tesina per il corso di: Circuiti e algoritmi per l’elaborazione dei segnali Prof. Aurelio Uncini Composto in LATEX

INDICE

v

PREFAZIONE

1

BASI 1.1

1

ALGORITMO LMS

1

1.1.1 Modello SISO 1 1.1.2 Algoritmo aggiornamento coefficienti

1.2 1.3 1.4

2

PROBLEMA DI NON UNICITÀ METODI DI DECORRELAZIONE ALGORITMO IPNLMS 5

I M P L E M E N TA Z I O N E 2.1

FILTER DESIGN ROOMSIM

4 5

9

9

2.1.1 Comandi comuni

2.2

2

10

12

2.2.1 Descrizione del toolbox 12 2.2.2 Risposte impulsive 14

2.3 2.4 2.5

PREPARAZIONE DELLA SIMULAZIONE 2.3.1 Algoritmo LMS 18 IMPLEMENTAZIONE IPNLMS 19 PRESTAZIONI 21 2.5.1 Troncamento della risposta impulsiva 2.5.2 Inseguimento 23 2.5.3 Coerenza delle sorgenti 24

BIBLIOGRAFIA

16

22

27

i

ACRONIMI INFOCOM

Dipartimento di Scienza e Tecnica dell’Informazione e della Comunicazione Dipartimento della facoltà di Ingegneria Elettronica della Sapienza situato in via Eudossiana 18 (Colle Oppio), Roma.

FA

Filtro Adattativo Filtro i cui coefficienti che vengono aggiornati ad ogni campione elaborato.

SISO

Single Input Single Output Ingresso Singolo Uscita Singola. Sistemi monodimensionali in cui si ha un segnale di uscita ed un ingresso.

MISO

Multiple Input Single Output Ingresso Multiplo Uscita Singola. Sistema con più ingressi ed una sola uscita.

MIMO

Multiple Input Multiple Output Ingressi Multipli Uscite Multiple. Sistemi a più dimensioni in cui si hanno molteplici segnali di uscita e di ingresso.

LMS

Least Mean Square Media Quadratica Minima. Algoritmo autoadattativo basato sulla minimizzazione della funzione errore quadratica.

IPNLMS

Improved Proportional LMS Algoritmo LMS Proporzionale Migliorato. Algoritmo analogo al LMS adattato al caso multicanale, in cui la funzione costo viene aggiornata con un peso variabile.

ARMA

Auto Regressive Moving Average Media Mobile Autoregressiva. Filtro per generare segnali di tipo stocastico a partire da un rumore normale opportunamente modellato dal filtro stesso.

MSE

Mean Square Error Errore Quadratico Medio. Criterio per l’individuazione deterministica della media dei quadrati di un segnale di errore.

SSE

Sum of Squared Error Somma degli Errori Quadratici. Criterio per l’individuazione stocastica della media dei quadrati di un segnale di errore.

MCAEC

Multichannel Acoustic Eco Canceller Cancellatore Acustico d’Eco Multicanale. Algoritmo per la rimozione tramite filtro adattativo del segnale che ritorna dagli altoparlanti al microfono.

MATLAB ©

Matlab © Matlab © è un programma della MathWorks © Inc, 1994-2009.

iii

iv

Indice FFT

Fast Fourier Transform Trasformata di Fourier Veloce. Algoritmo efficiente per il calcolo numerico della trasformata di Fourier.

LTI

Lineari Tempo Invarianti Classe di sistemi descritti da un’equazione differenziale lineare a coefficienti costanti della forma (2.1), che godono delle proprietà di linearità (la funzione e tutte le sue differenze finite figurano solo con esponente unitario o nullo) e permanenza (tutti i coefficienti sono costanti). Nei testi in lingua inglese tale proprietà è definita come LSI, Linear Shift-Invariant.

FIR

Finite Impulse Response Risposta Impulsiva Finita. Filtri con risposta impulsiva di lunghezza finita. Sono semplici da implementare, la stabilità è immediata se i coefficienti del filtro sono limitati.

IIR

Infinite Impulse Response Risposta Impulsiva Infinita. Filtri con risposta impulsiva infinita, che presentano un meccanismo di feedback per cui l’uscita dipende dall’uscita stessa negli istanti precedenti.

ERLE

Echo Return Loss Enhancement Perdita del Ritorno d’Eco. Indice per la misurazione dell’attenuazione dell’eco acustico di un filtro adattativo. È definito come il rapporto fra la potenza istantanea del segnale desiderato d[n] e la potenza istantanea del segnale di errore residuo e[n] che si ottiene dalla cancellazione.

SNR

Signal to Noise Ratio Rapporto Segnale Rumore. Rapporto in decibel tra il segnale utile ed il disturbo.

P R E FA Z I O N E In questo scritto è sviluppato un Multichannel Acoustic Eco Canceller (MCAEC) con algoritmo Improved Proportional LMS (IPNLMS), cioé un Filtro Adattativo (FA) in grado di attenuare, ed al limite di cancellare del tutto, il percorso di ritorno tra altoparlanti e microfoni. É possibile eliminare dal segnale registrato dal microfono (segnale desiderato), il segnale emesso dagli altoparlanti, ricercando la risposta impulsiva dell’ambiente. Tramite la risposta impulsiva è possibile simulare il percorso tra l’ingresso (l’altoparlante) e l’uscita (microfono) del sistema. Sottraendo dall’uscita del microfono, il segnale ricostruito elettronicamente si ha la cancellazione dell’eco acustico. Vediamo in questo studio di approfondire il problema descritto, e di realizzare un sistema in grado di eliminarlo. Il lavoro è svolto passo passo, a cominciare dai casi più semplici con l’algoritmo Least Mean Square (LMS) fino ad arrivare all’implementazione definitiva, con molteplici analisi prestazionali. Gli algoritmi sono simulati in linguaggio Matlab; sebbene attualmente sia in fase di studio al dipartimento Dipartimento di Scienza e Tecnica dell’Informazione e della Comunicazione (InfoCom) della Sapienza un programma basato sul linguaggio C per l’elaborazione online in tempo reale dei segnali audio, allo stato attuale dei lavori non è ancora possibile effettuare dei test realistici. Una delle applicazioni più diffuse dei sistemi di cancellazione dell’eco, è la realizzazione di sistemi di teleconferenza. Attualmente i sistemi più sviluppati per teleconferenze sono a singolo canale: si ha un solo microfono ed un solo altoparlante. Il limite di un sistema a singolo canale è l’assenza della spazialità della ripresa sonora, senza la quale, non è possibile percepire la provenienza dei suoni. Nel caso di videoconferenza, l’unico modo per capire chi sta parlando è quello di osservare i movimenti delle labbra e la gestualità, che non è affatto naturale. La maggiore naturalezza permessa da un sistema multicanale con almeno 2 altoparlanti, consente invece la ricostruzione spaziale della scena sonora e l’individuazione della posizione della sorgente, rendendo la comunicazione più istintiva e quindi meno affaticante, ciò spiega la spinta verso la realizzazione di sistemi multicanale. Quello che si vuole ottenere è un sistema che sia utilizzabile con microfoni fissi (per la massima naturalezza di impiego e libertà di movimenti), che sia di semplice utilizzo, che non richieda alcuna fase di configurazione, si potrebbe dire plug & play, e che sia sufficientemente realistico da permettere la localizzazione delle sorgenti. In un sistema di questo tipo, l’accoppiamento acustico tra l’altoparlanti e microfoni causerebbe fastidiosi echi dei suoni emessi e lunghi riverberi, con il ritorno dei suoni tra gli ambienti messi in comunicazione dal sistema. Nei casi peggiori, la presenza di un percorso ciclico chiuso nella comunicazione full-duplex, per l’instaurarsi di una controreazione positiva, potrebbe portare il sistema all’instabilità, che si manifesta sotto forma di fastidiosi fischi che impedirebbero la comunicazione.

v

vi

PREFAZIONE Il ruolo del Multichannel Acoustic Eco Canceller, è di aprire l’anello, bloccando i suoni che tornano dopo essere passati nell’altro ambiente di comunicazione. Attualmente sono diffusi sistemi Single Input Single Output (SISO), con un solo microfono ed un solo altoparlante, ottimizzati proprio per la cancellazione dell’eco. Un semplice filtro LMS è in grado di soddisfare le specifiche. In un sistema Multiple Input Multiple Output (MIMO) debbono essere impiegati più altoparlanti e più microfoni, almeno 2 altoparlanti e 2 microfoni per ciascuno degli ambienti messi in contatto. Sebbene apparentemente siano poco evidenti le complicazioni che incorrono nei sistemi multicanale, queste rendono inefficaci i metodi applicabili al caso SISO. Può accadere che canali riprodotti siano tra loro altamente correlati, anzi, è esattamente ciò che avviene ad esempio per i segnali di 2 microfoni in posizioni diverse nella stessa stanza con un unica sorgente, per cui il sistema di equazioni associato al MCAEC diviene singolare o quantomeno malposto. Assicurare una rapida convergenze in queste condizioni richiede la decorrelazione delle sorgenti, alterando i segnali riprodotti in qualche modo, ad esempio con delle distorsioni non lineari, senza distruggere però l’informazione che questi portano! Dalla qualità dei filtri di decorrelazione dipende la fedeltà della percezione sonora: meno percepibile è l’effetto, migliore è il filtro di decorrelazione (a parità di decorrelazione ovviamente). In questo testo vedremo quindi come affrontare e risolvere i problemi elencati. Nel primo capitolo saranno descritti gli algoritmi e le problematiche annesse, in modo analitico. La bontà dei metodi di decorrelazione è data in primo luogo dalla percezione soggettiva dei segnali che produce, per cui si tratta di un’elemento empirico piuttosto che matematico del problema, la cui progettazione si avvicina concettualmente all’arte. Nel secondo capitolo verranno studiate le tecniche per la simulazione dei sistemi, che consentono un’analisi delle prestazioni degli algoritmi. Gli algoritmi saranno implementati in Matlab, per cui si introdurranno i concetti fondamentali del filter design toolbox.

1

BASI

In letteratura, sono molteplici gli algoritmi adattativi studiati, la cui efficacia dipende dall’ambito di impiego. Lo scopo di un filtro adattativo è la stima di una risposta impulsiva reale h, ottenuta elaborando i campioni del segnale x man mano che questi si rendono disponibili al sistema. In base all’algoritmo scelto per questo compito, convergenza, capacità di inseguimento, complessità, robustezza, possono essere molto differenti! Un criterio di classificazione, è in base al tipo di segnale su cui operano: esistono algoritmi deterministici ed algoritmi stocastici. L’algoritmo LMS appartiene a questa seconda categoria. La differenza principale tra le due classi di algoritmi è data dalla funzione costo, che nel caso deterministico sarà del tipo Sum of Squared Error (SSE), mentre nel caso stocastico l’errore si definisce come il valore atteso dei quadrati della sequenza, cioè l’errore sarà del tipo Mean Square Error (MSE).

1.1

ALGORITMO LMS

Nel 1960 il prof. Bernard Widrow con il dottorando Marcian Edward Ted Hoff diedero vita all’algoritmo del gradiente stocastico o Least Mean Square (LMS), uno degli algoritmi adattativi più diffusi, probabilmente per la sua semplicità di implementazione e robustezza. Il principio di base con il quale viene costruito il filtro, è la minimizzazione dell’errore quadratico MSE: durante il percorso tra l’altoparlante ed il microfono, il segnale sonoro viene riprodotto dall’altoparlante, che introduce le prime distorsioni1 del segnale, poi sarà riflesso dalle pareti dell’ambiente, che introducono ulteriori distorsioni, l’eco, infine giunge al microfono che a sua volta sarà fonte di ulteriore distorsione, sia per l’inserzione fisica del microfono stesso che altera il campo sonoro, sia per l’equalizzazione introdotta; in definitiva, si può considerare una risposta impulsiva dell’ambiente, che porta dal segnale di ingresso, al segnale desiderato2 in uscita dal microfono. Se il filtro adattativo riproduce esattamente la risposta impulsiva del sistema fisico, posso creare un segnale elettronico identico al segnale desiderato, la loro differenza permette di togliere dall’uscita del microfono la parte di segnale che proviene dagli altoparlanti.

1.1.1

Modello SISO

In un modello Single Input Single Output (SISO), ho un solo microfono ed un solo altoparlante. Sia x [n] il segnale che pilota l’altoparlante, d[n] 1 Qui distorsioni è inteso come una qualsiasi alterazione del segnale, anche una equalizzazione, riverbero o cammini multipli. 2 In uscita dal filtro vorrei avere lo stesso segnale che ho in uscita dal microfono, quindi nella realizzazione del filtro il segnale di riferimento del microfono è indicato come desiderato.

1

2

BASI

x₁

nl

e

x₂

h₁

w₁ min{Ĵ }

nl

y₁ -

y₂

d

w₂

S

h₂

Figura 1: Schema completo di un MCAEC. Per ciascun microfono occorre ripetere gli elementi principali del filtro presentati su sfondo colorato. La minimizzazione della funzione costo dipende dai P segnali in ingresso al sistema, ed agisce sui filtri adattati w. Se tali filtri rispecchiano le risposte impulsive h dell’ambiente, i segnali y sono identici ai segnali catturati dal microfono provenienti dagli altoparlanti, per cui il segnale di errore ne sarà privo! Nell’ambiente di sinistra quindi non ci sarà l’eco dei segnali x, ma soltanto i suoni provenienti dalla sorgente S.

il segnale in uscita dal microfono. Stiamo considerando segnali campionati, un convertitore Digitale/Analogico a monte dell’altoparlante, ed un convertitore Analogico/Digitale a valle del microfono completano il sistema, in generale tali componenti non saranno ulteriormente considerati poiché inessenziali ai fini della trattazione. Il segnale desiderato si ottiene dalla relazione: d [ n ] = h [ n ] ⊗ x [ n ] + b [ n ],

(1.1)

in cui ovviamente h[n] è la risposta impulsiva del sistema, che comprende tutti gli elementi (lineari) che intervengono nel percorso tra x [n] e d[n], mentre b[n] raccoglie tutti gli ulteriori segnali che non provengono dall’altoparlante, che possono essere intesi sia come rumore, che come segnale, ad esempio una voce da trasmettere.

1.1.2

Algoritmo aggiornamento coefficienti

Gli algoritmi adattativi sono basati sulla minimizzazione di una funzione costo J (w) che dipende indirettamente dal filtro adattato, tipicamente indicato con w. L’uscita y del filtro w quando in ingresso ho il segnale x, dipende ovviamente dal filtro stesso. Se consideriamo l’uscita desiderata o misurata d d [ n ] = x [ n ] ⊗ h [ n ], e la confrontiamo con l’uscita y stimata dal filtro adattativo w y [ n ] = x [ n ] ⊗ w [ n ], è possibile definire l’errore del filtro, banalmente, come la differenza tra i due segnali e [ n ] ≡ d [ n ] − y [ n ] = d [ n ] − x [ n ] ⊗ w [ n ],

(1.2)

1.1 ALGORITMO LMS in cui si è esplicitato solo il segnale stimato, poiché in generale d[n] è un segnale noto del sistema, ad esempio il segnale proveniente da un microfono. L’algoritmo LMS ha il compito di stimare la risposta impulsiva h, e di ricrearne una sua replica che comunemente è indicata con w. L’algoritmo confronta il segnale desiderato d = x ⊗ h con il segnale ottenuto dalla simulazione y = x ⊗ w, se questi segnali non sono uguali ho un segnale di errore e[n] non nullo, dato dalla loro differenza. Gli algoritmi adattativi si possono classificare in base alla funzione costo, che può essere di tipo deterministico, definita come Sum of Squared Error (SSE): J (w) ≡

∑|e[n]|2 ,

o di tipo stocastico se la funzione costo è definita come Mean Square Error (MSE)  J (w) ≡ ∑ E |e[n]|2 , (1.3) come nel caso dell’algoritmo LMS. Minimizzare la funzione costo, vuol dire minimizzare la distanza3 tra il segnale desiderato ed il segnale stimato. É proprio dalla minimizzazione della funzione costo che si ricava l’algoritmo di aggiornamento dei coefficienti del filtro. L’algoritmo LMS considera, analogamente alle tecniche ls, l’errore quadratico istantaneo e2 [n] al posto della sua aspettazione. Per ciascuna iterazione calcolo l’errore in un instante n fissato che chiameremo k, per cui la (1.3) si riduce ad una semplice Jˆ(w) = e2 [k]

(1.4)

stimata, come indicato dal simbolo di accento circonflesso (ˆ) sulla funzione J. Ogni iterazione i coefficienti del filtro w[n] vengono aggiornati, può essere quindi conveniente indicare il filtro come un vettore usando una lettera in neretto, ed aggiungere un pedice che indica l’iterazione piuttosto che il coefficiente n-esimo. Indichiamo quindi il filtro all’iterazione k-esima con wk , in modo da poter scrivere l’espressione di adattamento come:
1 wk = wk−1 + µ −∇ Jˆ(w) , 2

(1.5)

che a questo punto dovrebbe risultare ovvia, in cui µ è un coefficiente costante detto passo di aggiornamento. Il gradiente della funzione costo si calcola sostituendo alla funzione costo la sua definizione con la (1.2), per cui  ∂ d[n] − x [n] ⊗ w[n] ∂e[n] ˆ = = −2e[n] x [n], (1.6) ∇ J (w) = 2e[n] ∂w ∂w poi applico l’approssimazione, o meglio, lo stimatore della (1.4), calcolando la (1.6) all’iterazione k-esima:

∇ Jˆ(w)|n=k = −2e[k] x [n]. 3 In generale posso minimizzare una qualsiasi norma p cioè || x [n]|| ≡ ∑n norma 2 coincide con la distanza euclidea.

p p

x p [n]. La

3

4

BASI Sostituendo infine nella (1.5), si ottiene l’equazione di aggiornamento dei coefficienti dell’algoritmo LMS: wk = wk−1 + µe[k] xk , in cui il vettore xk contiene gli ultimi punti del segnale di ingresso, in numero pari alla lunghezza del filtro adattativo w.

1.2

PROBLEMA DI NON UNICITÀ

L’algoritmo LMS in assenza di sorgenti nella stanza adatta i coefficienti del filtro affinché il segnale registrato dal microfono sia identico al segnale ricostruito elettronicamente dal filtro. L’operazione descritta è lineare nel dominio della frequenza, in cui l’operazione di filtraggio, che nel tempo corrisponde ad una convoluzione, si trasforma in un semplice prodotto. La risposta impulsiva di una sala, così come il filtri in grado di modellarne i canali acustici, possono essere considerati come sistemi Lineari Tempo Invarianti (LTI), semplificando notevolmente la trattazione e la realizzazione di filtri adattativi. Il segnale che posso ottenere dal filtro adattativo è una qualunque trasformazione lineare dell’ingresso. Una sorgente nella stanza (diversa dall’altoparlante), in generale, non sarà affatto correlata con il segnale emesso dall’altoparlante, il filtro adattativo quindi può simulare l’effetto del canale acustico sul segnale emesso dall’altoparlante, ma non può ricavare il segnale della sorgente indipendente. Quello che si può ottenere è la cancellazione del solo segnale riprodotto, quello che rimane è il solo segnale della sorgente indipendente. Se il sistema è multicanale, può accadere che i segnali in ingresso siano correlati, ad esempio generati da un’unica sorgente xs , filtrata secondo le P risposte della stanza che chiameremo g p . Questi P segnali correlati poi, attraverso un qualsiasi canale di comunicazione, raggiungono il secondo ambiente del Multichannel Acoustic Eco Canceller (MCAEC). I segnali da cancellare possono essere espressi dalla seguente: x p = g p ⊗ xs ,

p = 1, 2, . . . , P,

convoluzione tra la risposta impulsiva g e l’unica sorgente xs . Studiando la matrice di autocorrelazione R xx , che sarà data da: R xx = E{ x [k] ⊗ x [k]} = GE{ xs ⊗ xs } G T = GR xs xs G T , in cui G è la matrice delle risposte impulsive e R xs xs è la matrice di autocorrelazione della sorgente, si osserva che il rango della matrice di autocorrelazione non può essere completo, ma si ha al massimo:  Rg{ R xx } = min Rg{ G }, Rg{ R xs xs } Ker{ R xx } = PQ − Rg{ R xx }, in cui Ker e Rg indicano, rispettivamente, il nucleo ed il rango. Il sistema considerato quindi, avrebbe dimensione inferiore al numero delle equazioni, per cui non esiste una soluzione unica del problema, ed il MCAEC può non convergere alla risposta impulsiva dell’ambiente.

1.3 METODI DI DECORRELAZIONE

5

In generale, più i segnali x p sono correlati, peggiore è il condizionamento della matrice di autocorrelazione, quindi la convergenza dell’algoritmo sarà sempre più lenta o non convergerà affatto. È possibile aggirare questo ostacolo introducendo una distorsione non lineare che non sia percettibile, ma sufficente a ridurre la correlazione dei segnali, è quello che vedremo nella prossima sezione.

1.3

METODI DI DECORRELAZIONE

Esistono molteplici approcci per incrementare la decorrelazione dei segnali coinvolti, principalmente basati sull’introduzione di qualche tipo di non-linearità o distorsione. Un rettificatore può essere utilizzato a tale scopo in modo tale che la trasformazione si possa esprimere come:

x 0p (n) = x p (n) + α

x p (n) + | x p (n)| . 2

(1.7)

Risultati sperimentali mostrano che la percezione sonora è debolmente affetta da tali distorsioni che vengono assimilate a suoni di seconda armonica. Inoltre si hanno degli effetti di mascheramento nel parlato. Altri approcci sfruttano una rettificazione di semi-onda graduale che risulta più appropriata per la musica, sebbene presenti maggior complessità computazionale. Un altro modo per ridurre la correlazione è quello di inserire rumore opportunamente sagomato in modo da disturbare il meno possibile. Infine un ulteriore tecnica consiste nell’utilizzare filtri tempo varianti in modo da variare casualmente la posizione dei poli del filtro.

1.4

percezione sonora delle distorsioni

ALGORITMO IPNLMS

Quando le caratteristiche del filtro da costruire sono note si possono migliorare gli algoritmi esistenti. Sappiamo infatti che l’eco ha un inviluppo esponenziale: la maggior parte dell’energia si trova nella prima parte che contiene le prime riflessioni o early reflections , il resto della risposta impulsiva dato dalla radiazione diffusa o late reflections spesso contiene meno energia della prima parte, sebbene sia molto più lunga! In generale per tutti quei filtri sparse , cioè con energia concentrata in una porzione piccola della risposta impulsiva, è possibile scalare il fattore di aggiornamento dei coefficienti in proporzione ai coefficienti del filtro. L’algoritmo che si ottiene è l’LMS Proporzionale e Normalizzato. Come accennato nell’introduzione, in base all’algoritmo scelto si possono avere performance molto diverse. Considerando l’algoritmo LMS, le sue qualità sono la semplicità e la robustezza, il suo tallone d’Achille è la velocità di convergenza, ed anche l’inseguimento di conseguenza non è molto rapido. Si osserva (cfr. figura 2 nella pagina seguente) che la rapidità di convergenza dipende dal passo di aggiornamento dei coefficienti µ. La figura è stata ottenuta con l’algoritmo LMS, nel § 2.3.1 ne è riportata l’implementazione in Matlab. In ingresso è presente rumore uniforme bianco. Il passo più piccolo si ha per µ = 1/20 del massimo valore che consente la convergenza (1.8), si ottiene la curva in blu. L’errore in dB,

peculiarità della risposta impulsiva

importanza del passo di aggiornamento sulla convergenza dell’algoritmo

6

0

BASI

Disallineamento

Errore

0 −50

−40

−100

−60

−150

−80 0

−200 0

10

ERLE

200

−20

5

250 150 100 50 0

5

10

−50 0

5

10

Figura 2: Dipendenza dell’algoritmo LMS dal passo di aggiornamento µ, pari ad 1/20 (—), 1/10 (—) ed 1/5 (—) del massimo valore µ che consente la convergenza (1.8).

vantaggi per µ crescente . . .

. . . e svantaggi

decresce lentamente. così come il disallineamento (2.3). Al contrario l’Echo Return Loss Enhancement (ERLE)4 cresce sia con il tempo sia con il passo di aggiornamento µ, proprio ad indicare il tendere dell’uscita del filtro adattato al segnale desiderato. Ovviamente, mostrati gli aspetti positivi di un passo di aggiornamento più ampio, ci sono ottime ragioni per operare con un coefficiente µ molto piccolo, contrariamente a quanto appena detto, e sono:

• La convergenza dell’algoritmo LMS è dimostrata per processi gaussiani, stazionari, a media nulla e statisticamente indipendenti5 ! Si ottiene la condizione di convergenza: µ<

2  ; L E x2 [n]

(1.8)

• un coefficiente più elevato comporta un errore a regime più elevato, in particolare l’errore MSE a regime sarà: Jr =

µ 2 2 Lσ σ , 2 x b

(1.9)

l’errore dipende quindi dalle potenze dei segnali coinvolti nell’algoritmo, x2 e b2 . soluzione bilanciata: passo di aggiornamento µ variabile

passo di aggiornamento proporzionale ai coefficienti del filtro

Il primo accorgimento che porta ad un miglioramento della velocità di convergenza e riduce l’errore a regime, è di variare il passo di aggiornamento µ per ciascuna iterazione. Nelle prime iterazioni un passo di aggiornamento più rapido velocizza la convergenza, successivamente un passo più stretto permette di affinare la soluzione. L’algoritmo che si ottiene in questo modo è detto LMS Normalizzato. Perché allora non cambiare il coefficiente di aggiornamento, non solo per ciascuna iterazione, ma anche per ciascun elemento del filtro adattativo all’istante corrente? L’idea, vincente, ha originato l’algoritmo Proporzionale: per elementi di ampiezza elevata è possibile adattare con un passo elevato, per gli elementi più piccoli del filtro sarà più efficace un coefficiente di aggiornamento ad essi proporzionato. L’algoritmo proporzionale converge bene quando la risposta è sparsa o diffusa, cioè quando c’è molta disparità tra le ampiezze del filtro, come avviene, ad esempio, tra la parte iniziale e finale di una risposta impulsiva di una stanza. 4 La definizione (2.4) dell’Echo Return Loss Enhancement (ERLE) si trova al § 2.5 nella pagina 21. 5 Le ipotesi elencate semplificano la dimostrazione, ma non sono del tutto vincolanti.

1.4 ALGORITMO IPNLMS Tabella 1: Specifica algoritmo IPNLMS

Inizializzazione: Parametri:

Errore: Aggiornamento:

h = 0| L × P 0<α≤1 −1 ≤ β ≤ −1 1−β 0 < δ  σs2 2L ε>0 e(k ) = d − ∑ p w Tp (k − 1) x p

|w p (k)| 1−β + (1 + β ) 2PL 2kw p (k )k1 + ε α µ(k) = ∑ p ∑n x2p (k) g p (k) + δ w p ( k ) = w p ( k − 1) + µ ( k ) g p ( k ) x p e ( k ) g p (k) =

Quando la risposta non è sparsa, la proporzionalità del passo di aggiornamento deteriora le caratteristiche del filtro adattato, ottenendo un comportamento peggiore rispetto ad un filtro LMS normalizzato. L’approccio che permette di superare queste limitazioni si ha usando dei coefficienti che siano in parte pesati, ma che conservano hanno una parte fissa. L’algoritmo ottenuto con questa miglioria (Improved) prende il nome di Improved Proportional LMS (IPNLMS), ed attualmente è l’algoritmo derivato dalla famiglia di filtri LMS che converge più rapidamente e segue più accuratamente le variazioni della risposta impulsiva, comunque sia fatta: omogenea o sparsa. I parametri elencati, nell’ordine sono: α Passo base. È il passo di aggiornamento dell’algoritmo LMS, punto di partenza per trovare il passo di aggiornamento dell’IPNLMS. β Parametro di stabilizzazione. Controlla la proporzionalità dell’algoritmo, rendendo l’algoritmo più o meno proporzionale. Per β = 1 si annulla la parte costante del vettore di proporzionalità g, rimane l’algoritmo LMS Normalizzato e Proporzionale, ma non Improved. Tipicamente β ∈ [−0.5, 0]. δ Parametro di regolarizzazione. Corregge la parte costante del passo di aggiornamento Normalizzato µ. e Soglia. È una costante molto piccola che serve ad evitare divisioni per zero nell’implementazione dell’algoritmo.

7

2

I M P L E M E N TA Z I O N E

In questa seconda parte andiamo ad analizzare le soluzioni adottate per l’implementazione dell’algoritmo, in Matlab. Un toolbox per il filtraggio adattativo è già presente in Matlab, e fa parte dal pacchetto filterdesign. Per prima cosa daremo quindi uno sguardo al toolbox, poi analizzaremo il toolbox roomsim per la simulazione delle risposte impulsive, infine andremo ad analizzare più in dettaglio l’implementazione dei nuovi algoritmi.

2.1

FILTER DESIGN

Come già accennato nell’introduzione, in Matlab (Versione 7.5) è presente un toolbox (in versione 4.2) chiamato filterdesign, per l’implementazione di filtri deterministici ed adattativi. L’algoritmo LMS con la sua variante normalizzata sono pronti per l’uso, mancano però le varianti più complesse fino alla IPNLMS, che perciò dovranno essere implementate manualmente. Il toolbox esistente è basato sul metodo di programmazione ad oggetti, che permette di costruire i filtri adattativi come oggetti, ossia come classe di dati strutturati che contengono tutte le informazioni relative al filtro. Per ciascuna classe esiste uno specifico metodo che implementa in modo opportuno, dipendente dalla classe, alcune elaborazioni per l’oggetto considerato. Chiariamo con un esempio. Supponiamo di essere interessati al filtraggio di un segnale. Per fare questa operazione esiste il comando filter. Nel caso in cui il filtro sia statico, il metodo della convoluzione, come noto, risolve il problema. Lo stesso comando filter può operare su una classe diversa, ad esempio se il filtro è di tipo adattativo con algoritmo LMS, le operazioni richieste sono quelle elencate nel listato 2.2 nella pagina 18. Nel gergo della programmazione ad oggetti si dice che si ha un overloading dell’operatore, cioè si hanno diversi metodi selezionati dal tipo di classe con la quale si sta operando. In Matlab, il meccanismo di selezione del metodo appropriato, in base alla classe, avviene in modo automatico, rispettando la semplice regola della struttura delle directory: gli operatori per una data classe devono trovarsi all’interno di una cartella con il nome formato dal carattere @ seguito dal nome della classe (@nomeclasse). È possibile anche annidare più classi, ad esempio per i filtri adattativi esiste la classe genitore adaptfilt i cui metodi si trovano nella cartella /@adaptfilt,1 al cui interno sono annidate le classi specifiche a ciascun algoritmo. Quindi i file che implementano il filtro adattativo con algoritmo LMS vanno collocati nella cartella /@adaptfilt/@LMS. Per espandere il toolbox aggiungendo anche l’algoritmo IPNLMS i relativi 1 Tale cartella sarà posizionata nella sotto cartella /toolbox/filterdesign che solitamente si trova nel percorso di installazione di Matlab.

9

IMPLEMENTAZIONE

10

file di implementazione andranno posti in una cartella /@adaptfilt/ @IPNLMS.

Tuttavia, in questo testo non si andrà ad espandere il toolbox di Matlab come descritto, ma ci è sembrato opportuno dare dei cenni di base per capire da subito come partire con il piede giusto per cominciare a lavorare e per capire il lavoro svolto.

2.1.1 Comandi comuni Sempre a titolo informativo, una breve descrizione dei comandi più comuni permette di capire da subito gli algoritmi proposti, maggiori informazioni ovviamente sono disponibili sulla guida di Matlab adaptfilt.lms Permette di generare dei filtri adattativi LMS o di altri tipi, specificandone il nome (scelto tra quelli disponibili per default) dopo il punto. La sintassi è: 1

ha=adaptfilt.lms(L,stepsize,leakage,coeffs,states)

in cui i parametri sono: L Lunghezza del filtro adattativo. È l’unico parametro che deve essere impostato STEPSIZE è il passo di adattamento µ, se omesso la funzione calcola un valore opportuno. LEAKAGE Ë un fattore di perdita compreso in [0, 1], agisce moltiplicando il filtro adattato corrente in ciascuna iterazione. Per default è 1, cioè non si ha leakage. COEFFS Condizioni iniziali dei coefficienti del filtro. Per default è 0. STATES Condizioni iniziali stati del filtro. Nulli per default. Eseguendo il comando specificando la sola lunghezza L pari a 1024, nel prompt di Matlab sono mostrate alcune delle informazioni caratteristiche per l’oggetto adaptfilt.lms appena creato: 1 2 3 4 5

Algorithm: FilterLength: StepSize: Leakage: PersistentMemory:

’Direct-Form FIR LMS Adaptive Filter’ 1024 0.1 1 false

tra le strutture di dati non mostrate, l’elemento più importante probabilmente è coefficients, che contiene tutti i punti del filtro calcolato, accessibili con il comando ha.coefficients. filter L’operazione di filtraggio avviene tramite il comando filter, che può accettare semplici vettori o oggetti più complessi come il filtro adattativo creato precedentemente. Se ho un semplice vettore posso scrivere 1

y = filter(b,a,x);

2.1 FILTER DESIGN per calcolare il filtro Infinite Impulse Response (IIR) con l’equazione N

∑

M

ak y[n − k ] =

k =0

∑ bk x [ n − k ] .

(2.1)

k =0

Se ho un filtro adattativo invece, posso scrivere 1

[y,e] = filter(ha,x,d);

per calcolare sia l’uscita del filtro adattativo per l’ingresso x in riferimento al segnale desiderato d, che l’errore e tra il segnale generato dal filtro y ed il segnale desiderato d. xcorr Calcola la crosscorrelazione tra dei vettori x, y, o l’autocorrelazione se un solo vettore viene passato alla funzione. La sintassi è 1 2

c = xcorr(x,y); % Correlazione c = xcorr(x); % Autocorrelazione

mfilt Crea un filtro multirate, permettendo di scegliere tra vari metodi, tra cui cascata di filtri, interpolatori e decimatori. In particolare un semplice decimatore Finite Impulse Response (FIR) si ottiene con la sintassi 1

hd = mfilt.firdecim(m)

con m fattore di decimazione. Anche in questo caso, dopo aver creato il filtro, la decimazione vera e propria si ottiene filtrando il segnale da decimare con l’istruzione y = filter(hd,x). norm Una semplice istruzione per calcolare la norma di un vettore x modulo è

p, 1

normap = norm(x,p);

File audio Per operare con i file audio sono disponibili i comandi 1

x = wavread(file);

che assegna al vettore x il contenuto del file indicato dal vettore di tipo stringa file, eventualmente passato come testo tra singoli apici e non come variabile. Il comando inverso, per salvare un vettore elaborato come file .wav è 1

wavwrite(x,path)

Per ascoltare un vettore è necessario per prima cosa creare un oggetto player 1

plx = audioplayer(x,fs)

11

IMPLEMENTAZIONE

12

Tabella 2: Principali cartelle del toolbox roomsim

Directory

Contenuto

\cipic \Excel_setups \text_setups \sensor\types \Impulse_response

Risposte impulsive del capo File di impostazione in formato excel File di impostazione in formato testo Risposte direzionali dei microfoni Risposte impulsive calcolate da roomsim

che necessita anche della frequenza di campionamento oltre che del segnale, poi è possibile ascoltare l’oggetto player creato con l’istruzione 1

play(plx)

2.2

ROOMSIM

Il toolbox roomsim2 è un insieme di programmi Matlab per calcolare la risposta impulsiva di un ambiente. Non appena viene installato il toolbox, occorre aggiornare le directory Matlab, successivamente sarà possibile avviare l’interfaccia grafica mostrata in figura 3a nella pagina successiva, digitando al prompt di Matlab: 1

roomsim

purché sia impostata come directory di lavoro la directory radice di roomsim.

2.2.1 Descrizione del toolbox Le risposte impulsive si possono calcolare selezionando il pulsante Set -up and Run the Room simulation. Nel menù successivo è possibile specificare il modo in cui sono passati a roomsim i parametri di input. Si è scelto di impostare la simulazione tramite fogli excel, per cui occorre selezioniare il pulsante Read from an Excel spreadsheet, come mostrato in figura 3b a fronte. Prima di proseguire, selezionando l’opportuno file excel, occorre spendere alcune parole sui file che vengono copiati da roomsim. Alcune delle directory generate durante l’installazione di roomsim sono elencate nella tabella 2. In particolare la cartella contenente i file excel permette di impostare velocemente la simulazione, in base ad una serie di parametri. Nella stessa cartella è disponibile il file absorption_table. xls che mostra i valori tipici dei coefficienti di assorbimento in base al materiale delle pareti. Non appena si seleziona il file di setup ha inizio la simulazione, che dopo pochi istanti mostra i grafici: 1. Reverberation time vs frequency 2. Absorption coefficients vs frequency 2 É possibile scaricare roomsim sul sito Matlab all’indirizzo http://www.mathworks.com/ matlabcentral/fileexchange/5184.

2.2 ROOMSIM

(a)

(b)

Figura 3: Menù principali del programma roomsim

(a)

(b)

Figura 4: Menù per l’impostazione delle unità di ascisse ed ordinate del grafico della risposta impulsiva roomsim

3. Room receiver & source geometry Per quanto riguarda la risposta impulsiva, mostrata nel 4° grafico, ci sono ben 2 successivi menù (fig. 4) che permettono la personalizzazione delle ascisse ed ordinate. Il 5° grafico rappresenta la Fast Fourier Transform (FFT) della risposta impulsiva, ma prima un’ulteriore finestra di dialogo chiede di selezionare la lunghezza dell’algoritmo FFT. La risposta impulsiva è calcolata con il metodo delle immagini, proposto da Allen a Berkley. La soluzione dell’equazione delle onde viene ricavata in modo approssimato con i principi dell’ottica geometrica: una parete piatta e rigida riflette il suono, che si propaga in linea retta. Ogni volta che il suono viene riflesso, fino ad arrivare al microfono, è possibile tracciare una linea retta fino ad una sorgente immagine posta esternamente alla stanza. Il problema della ricerca degli infiniti percorsi dalla sorgente al microfono, può essere rovesciato, considerando un solo percorso rettilineo ma moltiplicando il numero di sorgenti. Le sorgenti si ottengono in modo molto semplice specchiando la stanza su ciascuna parete, e reiterando il procedimento fino ad allontanarsi oltre una certa distanza. Il 6° grafico di roomsim mostra proprio le sorgenti immagini.3 Può 3 Sia il grafico delle sorgenti immagini che il successivo, possono essere omessi tramite

13

14

IMPLEMENTAZIONE

0.15

h11, ξ=0.53117

h12, ξ=0.50478

h , ξ=0.51977

h , ξ=0.57511

0.1 0.05 0 -0.05 21

0.15

22

0.1 0.05 0 -0.05

0

0.1

0.2

0.3

0

0.1

0.2

0.3

Figura 5: I grafici delle 4 risposte impulsive calcolate

essere difficile da decifrare, specialmente a prima vista, dato l’elevato numero di punti coinvolti nel calcolo, ma fornisce a colpo d’occhio un’idea dell’approssimazione applicata nel calcolo della risposta impulsiva. Prima del grafico, come sempre, alcune finestre di dialogo permettono di impostare alcuni parametri di visualizzazione. Occorre tenere presente che il sistema delle immagini approssima la risposta impulsiva, ma soffre di alcune limitazioni. In particolare in ambito audio, le lunghezze d’onda sono spesso spesso elevate, soprattutto alle basse frequenze com’è ovvio. Quando le lunghezze d’onda sono più elevate delle dimensioni della parete subentrano fenomeni di diffusione sonora o scattering, che vengono ignorati dall’algoritmo, rendendo i risultati meno accurati. Al termine delle elaborazioni, roomsim chiede di salvare i dati calcolati come unico file in formato .mat infine torna al menù iniziale. Nella directory impulse response invece è possibile trovare il file della sola risposta impulsiva, più comodo da utilizzare nelle simulazioni. La risposta impulsiva di un ambiente in generale presenta delle caratteristiche ricorrenti, come un inviluppo esponenziale, picchi distinti nella parte iniziale ed una coda molto lunga acusticamente più omogenea. La lunghezza della coda necessita un’approssimazione in fase di calcolo, difficilmente si potranno gestire tutti i campioni della risposta reale, sarà necessario troncare la risposta stimata quando l’ampiezza scende al di sotto di un certo valore, che per roomsim è pari a −45 dB. L’approssimazione commessa dovrebbe essere accettabile, cioè tale da non introdurre errori elevati!

2.2.2 Risposte impulsive Se consideriamo 2 sorgenti e 2 microfoni, le possibili risposte impulsive sono 4. Consideriamo quindi una stanza con dimensioni di 5.10 × 4.25 × 3.4 m. La maggior parte dei parametri sono comuni a tutte opportune impostazioni del file excel di setup, tuttavia sono fornite in questa parte alcune informazioni di interesse.

2.2 ROOMSIM Tabella 3: Parametri generali di impostazione del toolbox roomsim

Valore 32 000.00 40.00 20.00 −1 8192 H_11

Nome

Descrizione

Fs humidity TEMP order H_length H_filename

Frequenza di campionamento Percentuale di umidità relativa Temperatura in °C Limite massimo di riflessioni (automatico) Lunghezza risposta impulsiva (4 × 2048) Nome del file della risposta impulsiva Schema stanza

y 4 3

R₂ 2

h₁₂

R₁

h₂₁

h₂₂

1 0

h₁₁

S₁

S₂ 0

1

2

3

4

5

x

Figura 6: Pianta dell’ambiente simulato, per il calcolo delle risposte impulsive. Con S sono indicate le sorgenti, con R i sensori

le simulazioni, In particolare sono stati applicati i parametri elencati nella tabella 3. La frequenza di lavoro adottata, adeguata all’elaborazione vocale che ha una banda più stretta della musica, è di 8 KHz, ma roomsim richiede un numero maggiore di punti per poter calcolare la risposta impulsiva. La simulazione avviene quindi a 32 KHz, frequenza quadrupla rispetto alla frequenza operativa, per decimare poi i campioni di un fattore 4. I coefficienti di riflessione sono gli stessi del file setup_6_surfaces fornito con roomsim. Le simulazioni sono state effettuate con un solo sensore di tipo omnidirezionale, soltanto la posizione delle sorgenti e dei sensori cambiano. Con 4 simulazioni sono state ottenute le 4 risposte impulsive h pq di figura 5 a fronte tra la sorgente p-esima ed il sensore q-esimo, come riportato nella figura 6. I dati dei sensori e delle sorgenti sono indicati nella tabella 4 nella pagina seguente. Risposta impulsiva sparsa L’ambiente simulato è molto riverberante, le risposte impulsive infatti non sono molto sparse, al contrario, si presentano dense. Compreso il concetto intuitivo di distribuzione sparsa, proviamo darne una definizione formale che sia in grado di associare un numero al segnale che quantifichi il fenomeno. Tale definizione dovrebbe essere:

15

16

IMPLEMENTAZIONE Tabella 4: Specifica dei parametri delle 4 simulazioni

1 x y z

Sensori 2

1 2 1.2

1 2.25 1.2

Microfoni 1 2 4 2 1.2

3 1 1.2

• invariante rispetto al prodotto per una costante non nulla; • indipendente dall’ordine dei coefficienti della risposta impulsiva e dipendente dalla loro ampiezza relativa; • limitata, l’indice di densità o diffusione dovrebbe essere un numero compreso in un intervallo finito. Per rispettare le linee di principio elencate, è stata proposta la seguente definizione dell’indice di densità,   L k h k1 √ √ ξ (h) ≡ 1− , L− L L k h k2

(2.2)

che gode proprio delle proprietà descritte, in particolare in formule 0 ≤ ξh ≤ 1,

∀ a 6= 0,

ξ ( ah) = ξ (h).

L’estremo superiore si ha per la funzione di dirac che è la distribuzione sparsa per eccellenza, quindi ξ (δ) = 1, al contrario la densità di una distribuzione uniforme di infiniti elementi identici è indicata con ξ (u) = 0. La resa acustica delle risposte impulsive presentate in figura 5 nella pagina 14 è piuttosto riverberante, indice di densità elevata di una coda che decade lentamente. L’indice di densità per queste risposte è indicato sopra ciascuna di esse, ed è poco al di sopra di 0.5, con lievi differenze. La risposta più densa è la h12 , mentre la risposta tipicamente utilizzata quando non diversamente specificato è la h11 . Un’ulteriore risposta impulsiva è stata calcolata per ottenere un risultato con un indice ξ ancora maggiore. Per questo, è stata fatta una simulazione tra la sorgente 2 ed il ricevitore 1 cambiando solamente i coefficienti di assorbimento dell’ambiente. Per le pareti sono stati usati i coefficienti di assorbimento medi delle mattonelle, pavimento e soffitto sono stati considerati in legno. I coefficienti sono elencati nel file absorption_table, rispettivamente con il nome di acoustic_tile_average e platform_floor_wooden. Tali valori sono riportati nella tabella 5 a fronte, mentre il grafico della risposta impulsiva è mostrato in figura 7.

2.3

PREPARAZIONE DELLA SIMULAZIONE

Per preparare la simulazione occorre il segnale proveniente da uno dei microfoni. Una prima simulazione con un semplice algoritmo LMS

2.3 PREPARAZIONE DELLA SIMULAZIONE Risposta impulsiva,

0.1

ξ=0.632

0.08 0.06 0.04 0.02 0 -0.02 -0.04

0

0.05

0.1

0.15

0.2

0.25

0.3

0.35

Figura 7: Risposta impulsiva con un coefficiente di densità più elevato. Tabella 5: Coefficienti di assorbimento usati per la simulazione della risposta impulsiva mostrata in figura 7

frequenza (Hz)

125

250

500

1000

2000

4000

piastrelle legno

0.10 0.40

0.30 0.30

0.80 0.20

0.85 0.20

0.75 0.15

0.65 0.10

può servire per capire se la simulazione è stata impostata correttamente. Sono stati presi 2 file .wav a 8 KHz, il primo di una voce maschile (malespeech.wav), il secondo una voce femminile (femalespeech.wav). La voce femminile coincide con la sorgente 1 e si suppone sia originata da un altoparlante, la voce maschile è posta nella sorgente 2, entrambe in riferimento allo schema in figura 6 nella pagina 15. Ciò è ottenuto in pratica associando le opportune risposte impulsive! In questo modo è stato ottenuto un segnale che simula il suono catturato da un ipotetico microfono posto nella posizione R1, salvato nel file micrec.wav. 1 2 3 4 5 6 7 8 9 10

% Carico le risposte impulsive (e la frequenza di campionamento Fs) load([path,’\Matlab\ImpulseResponse\H_11_S1’]); h11=data; load([path,’\Matlab\ImpulseResponse\H_21_S1’]); h21=data; clear data; % Decimo le risposte impulsive di 4 volte, da 32KHz a 8KHz hm = mfilt.firdecim(4); % Filtro di decimazione di un fattore 4, h11 = filter(hm,h11); % Risposta decimata h21 = filter(hm,h21); Fs=Fs /4; clear hm; % Aggiorno Fs e cancello hm

11 12 13 14 15 16 17 18

x = wavread([path,’\Matlab\Sound\femalespeech’]); % Segnale speaker b = wavread([path,’\Matlab\Sound\malespeech’]); % Segnale oratore x(10*Fs +1:length(x))=[]; % Tronco i segnali dopo 10 secondi b(10*Fs +1:length(b))=[]; d= conv(x,h11); % Segnale x con eco b= conv(b,h21); d= d+b; d=d./max(d); % Segnale microfono normalizato

17

IMPLEMENTAZIONE

18

19 20

%wavwrite(d,[path,’\Matlab\Sound\micrec’])

Listato 2.1: Comandi Matlab per generare il segnale registrato del microfono.

Il codice è mostrato nel listato 2.1, ci sono poche cose da aggiungere:

• il comando load carica 2 variabili dette data e Fs, rispettivamente la risposta impulsiva del filtro e la frequenza di campionamento; • il segnale d catturato dal microfono è dato dalla somma dei 2 ingressi filtrati ciascuno secondo il proprio percorso. Il segnale x dell’altoparlante non viene sovrascritto al contrario del segnale b che corrisponde alla voce maschile. Entrambi sono troncati dopo 10 s, in questo modo è possibile sommarli per calcolare l’uscita del microfono, che è stata normalizzata tra −1 e 1.

2.3.1 Algoritmo LMS Se l’algoritmo Least Mean Square (LMS) funziona, dovrebbe essere in grado di di ricostruire la risposta impulsiva h11 confrontando il segnale x ed il segnale ricevuto dal microfono. Il segnale y mi aspetto che vada a coincidere con il segnale ricevuto dal microfono se x fosse l’unica sorgente nella stanza. L’errore, vale a dire il segnale ricevuto dal microfono meno la parte relativa ad x, dovrebbe essere il solo segnale proveniente dalla sorgente 2, cioè la sola voce maschile. 1 2 3 4 5 6 7 8

for n=1:niter, X(2:L) = X(1:L-1); % parte di x traslata ed in ordine inverso X(1) = x(n); % aggiungo il campione corrente y(n) = W*X; % calcola ed assegna l’uscita corrente e(n) = d(n) - y(n); % calcola l’errore corrente W = W + mu*e(n)*X’; % aggiorna i coefficienti del filtro misa(n) = 20* log10(norm( (h(1:L)-W’),2 ) ./ norm(h(1:L),2)); end

Listato 2.2: Implementazione in Matlab dell’algoritmo LMS

Dall’algoritmo 2.2 si ottiene proprio, ben distinguibile, la voce maschile, con un lieve rumore di fondo, appena percepibile, dell’altra voce. Il principale problema è la bassa velocità di convergenza, per cui per ottenere il file malespeech_SISO_LMS-2nd_iter sono state necessarie 2 iterazioni, la prima permette di calcolare i coefficienti dell’algoritmo, la seconda permette di ottenere un output pulito. L’algoritmo necessita di essere inizializzato, vanno preallocati i vettori coinvolti nel ciclo principale, e devono essere calcolati i parametri necessari all’esecuzione. L’inizializzazione avviene tramite il codice: 1

% INIZIALIZZAZIONE dei coefficienti del filtro LMS

2 3 4 5 6 7 8 9

L=2048; Sx=size(x); niter = max(Sx); y = zeros(niter,1); e = y; X = zeros(L,P); clear Sx;

% Lunghezza FA

W = X’;

% inizializzo il filtro adattativo W

% % % %

numero di iterazioni creo un vettore di uscita di elem nulli creo il vettore errore di elem nulli L elem di x usati in ogni iterazione

10 11

2.4 IMPLEMENTAZIONE IPNLMS

12 13

mu=2/(L * max(var(x)) )/10; misa=zeros(L,P);

Listato 2.3: Implementazione in Matlab dell’algoritmo LMS

Per l’elaborazione principale dell’algoritmo, viene usato un vettore X(n) che contiene gli ultimi L punti di x, in ordine decrescente. Il 1° elemento infatti, contiene il campione corrente, col crescere dell’indice n si hanno i capioni precedenti. In questo modo, il prodotto riga × colonna (linea 16) corrisponde al calcolo della convoluzione per il solo punto corrente.

2.4

IMPLEMENTAZIONE IPNLMS

L’algoritmo IPNLMS è stato implementato considerando il caso più realistico di P ingressi ed una sola uscita, lo stesso algoritmo riassunto nella tabella 1 nella pagina 7. I nomi di alcuni parametri sono stati cambiati per adattarli a Matlab, il codice risultante è nel listato 2.4. 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18

for n=1:niter, X(2:L,:) = X(1:L-1,:);% parte di x traslata ed in ordine inverso X(1,:) = x(n,:); % aggiungo il campione corrente X2 = X.^2; % valore quadratico dei segnali for m=1:P, g(m,:)=g0+((1+b)*abs(W(m,:)))/(2*norm(W(m,:),1)+t);%vett scal y(n)=y(n)+W(m,:)*X(:,m); % uscita mu = alp / ( (g(m,:)*X2(:,m)) +delta);% cost d’aggiornamento end e(n) = d(n) - y(n); % errore for m=1:P, W(m,:)=W(m,:)+mu*e(n)* g(m,:).*(X(:,m))’; % aggiorna i coeff misa(n,m) = 20* log10(... norm( (h(1:L,m)- (W(m,:))’ ) ,2 ) ./ norm(h(1:L,m),2) ... ); end end beep

Listato 2.4: Implementazione in Matlab dell’algoritmo IPNLMS

Gli ultimi L campioni del vettore x sono memorizzati nella matrice X, in cui ciascuna colonna corrisponde ad un segnale in ingresso differente. Il prodotto tra il segnale ed il relativo filtro adattato viene effettuato con un ciclo interno, per ciascun ingresso. Anche in questo caso l’algoritmo deve essere inizializzato, i parametri da impostare sono molti di più, il codice che svolge questo compito è il seguente: 1 2 3 4 5 6

niter = max(Sx); clear Sx; y = zeros(niter,1); X = zeros(L,P); e = y; W = X’;

% numero di iterazioni % % % %

creo un vettore di uscita di elem nulli L elem di x usati in ogni iterazione creo il vettore errore di elem nulli inizializzo il filtro adattativo W

7 8 9 10

% Preallocazione variabili mu = e; misa = zeros(L,P);

19

IMPLEMENTAZIONE

20

ERLE

60 50 40

(dB)

30 20 10 0 −10 0

1

2 3 tempo(sec)

4

5

Figura 8: Prestazioni dell’IPNLMS al variare di β e della risposta impulsiva. In particolare beta = −0.5 per le curve in blu (—) e verde (—), le altre due curve si ottengono per β = 0.5. Le curve blu (—) e rossa (—) approssimano la risposta impulsiva h12 , per la quale ξ (h12 ) = 0.502, le curve verde (—) e magenta (—) sono calcolate con hs , ξ (hs ) = 0.632

11

g = zeros(P,L);

12 13 14 15 16 17 18

% Parametri specifici IPNLMS t = 0.001; %t = 5 * eps; % soglia (treshold), evita divisioni per 0 alp = 0.3; % b = -0.5; % parametro di stabilizzazione delta0 = 10 * mean((var(x)))^2; delta = delta0*(1-b)/(2*L); % parametro regolarizzazione

Listato 2.5: Implementazione in Matlab dell’algoritmo IPNLMS

Le lettere greche sono sostituite dalle equivalenti lettere romane, tranne δ che è chiamato delta, poiché d rappresenta il segnale desiderato. Confronto con gli algoritmi più semplici Si è discusso nel § 1.4 nella pagina 5 su come l’algoritmo IPNLMS costituisca un miglioramento rispetto alle altre varianti dell’LMS. Consideriamo quindi rumore bianco uniforme in ingresso del filtro adattativo, ed andiamo a vedere l’andamento temporale dell’ERLE4 nei primi 5 secondi, nel caso in cui la risposta impulsiva da approssimare sia più o meno sparsa. Consideriamo in primo luogo la risposta impulsiva h12 mostrata in figura 5 nella pagina 14, piuttosto densa poiché l’indice di dispersione ξ (h12 ) ≈ 0.502. I parametri tipicamente utilizzati sono:

• α = 0.3; • β = −0.5; 4 Per la definizione dell’ERLE si veda il § 2.5 nella pagina successiva

2.5 PRESTAZIONI MISA

10 0 −10

(dB)

−20 −30 −40 −50 −60 0

1

2 3 tempo(sec)

4

5

Figura 9: Disallineamento. Stesso test della figura 8 ma con misura del disallineamento invece dell’ERLE

• e = 0.001; • δ = 10σx2 . Per β = −0.5, si ha un buon compromesso tra proporzionalità e passo costante, mentre valori di β più elevati rendono l’algoritmo sempre più proporzionale. Per β = 0.5 l’IPNLMS si può considerare esclusivamente proporzionale, perdendo di fatto la migliorata capacità di convergenza per risposte non sparse. In definitiva per β = 0.5 e per risposte dense l’algoritmo dovrebbe convergere più lentamente. In figura 8 sono riportati i valori dell’ERLE ottenuti in simulazione, i parametri per ciascuna curva sono i seguenti: — Parametri standard. L’algoritmo converge rapidamente anche per risposte dense; — β = 0.5. Approssimazione dell’algoritmo PNLMS non migliorato. L”algoritmo non si comporta bene con le risposte dense, infatti è il più lento tra gli algoritmi presentati; — Risposta sparsa con parametri standard. Gli algoritmi convergono più rapidamente sulle risposte sparse, in particolare l’IPNLMS si dimostra particolarmente efficiente, ho la risposta migliore; — Risposta sparsa con β = 0.5. Anche in questo caso le prestazioni sono peggiori dell’IPNLMS, ma migliori del caso precedente. Come si osserva le prestazioni migliori si hanno per risposte sparse, in generale l’algoritmo IPNLMS converge più rapidamente.

2.5

PRESTAZIONI

Per finire vediamo come si comportano gli algoritmi presentati, quali sono cioè le prestazioni.

21

IMPLEMENTAZIONE 50

ERLE

10

40

0

30

−10

20

−20

(dB)

(dB)

22

10

−30

0

−40

−10 0

2 tempo(sec)

4

−50 0

MISA

2 tempo(sec)

4

Figura 10: Effetto del troncamento della risposta impulsiva sulle capacità di adattamento del filtro.

Una possibile misura della bontà dell’algoritmo si ottiene semplicemente dalla funzione errore: esprimendone il valore in dB è possibile apprezzarne la convergenza anche quando scende al di sotto di un millesimo del valore iniziale. Gli algoritmi adattativi, solitamente devono trovare in qualche modo una risposta impulsiva di un sistema, cercando di ricavarne l’andamento, con la migliore accuratezza possibile, nel più breve tempo possibile. Una stima della bontà di un filtro adattativo si può quindi ottenere dal confronto della risposta impulsiva stimata con la risposta impulsiva reale, solitamente normalizzato rispetto alla risposta impulsiva ed espresso in dB. Quello che si ottiene è il disalineamento che in formule si scrive ∆h ≡ 20 log10

k h − W k2 k h k2

(2.3)

Un parametro specifico per la cancellazione dell’eco, cioè per i Multichannel Acoustic Eco Canceller (MCAEC) si può definire considerando il segnale desiderato e l’errore, in particolare dal rapporto delle loro potenze in dB, η ≡ 10 log10

Pd Ed2 [n] = 10 log10 2 , Pe Ee [n]

(2.4)

si ottiene l’Echo Return Loss Enhancement (ERLE). Si può considerare l’ERLE come l’analogo del rapporto segnale rumore (Signal to Noise Ratio, SNR). Rispetto al disallineamento, l’ERLE si occupa dell’uscita del filtro, piuttosto che dell’ingresso, fornice quindi una stima più accurata della resa del filtro.

2.5.1 Troncamento della risposta impulsiva Le risposte impulsive tipicamente presentano lunghe code, della durata spesso superiore a qualche secondo. Tuttavia i primi coefficienti sono più significativi, poiché hanno un’ampiezza maggiore, mentre gli ultimi sono sempre più prossimi allo zero. Per semplificare l’elaborazione e ridurre il numero di operazioni necessarie ad ogni ciclo, la lunghezza del filtro adattato, spesso è inferiore alla reale lunghezza del filtro da ricavare. Vediamo come influisce

2.5 PRESTAZIONI Segnale desiderato

0.8 0.6 0.4 0.2 0 −0.2 −0.4 −0.6 −0.8

0

1

2

3

4 5 6 tempo(sec)

7

8

9

10

Figura 11: Andamento del segnale desiderato per via della variazione della risposta impulsiva dell’ambiente. Ciascun colore indica una risposta impulsiva.

questa approssimazione sulle capacità di adattamento del filtro, osservando il grafico in figura 10 nella pagina precedente, che mostra proprio l’andamento dell’ERLE e del disallineamento al variare della lunghezza del filtro adattato w. La prima esecuzione dell’algoritmo mostra l’adattamento nel caso in cui la lunghezza del filtro adattato coincide con quella della risposta impulsiva considerata, cioè per L = 2048. Successivamente è stata dimezzata la lunghezza del filtro adattato, portato prima a L = 1024 (curva rossa —), poi a L = 512 (curva verde —). Più il filtro adattato è breve, peggiore è l’approssimazione della risposta impulsiva, e maggiore è l’errore in uscita dal filtro. In questa simulazione è stata usata la risposta impulsiva h21 .

2.5.2 Inseguimento Un’altra proprietà delle risposte impulsive, è la variabilità. In generale i parametri di una sala non saranno mai costanti, basta aprire una finestra o una porta per variare considerevolmente l’acustica dell’ambiente, su cui agiscono anche molti altri fattori, come la presenza di persone in movimento, flussi d’aria e variazioni di temperatura. La simulazione della variazione della risposta impulsiva della stanza è stata sostituita con la variazione del percorso, sono state quindi considerate le risposte impulsive h12 , hs ed h11 , con coefficienti di densità rispettivamente di 0.505, 0.632 e 0.531. Il segnale desiderato è stato ottenuto a partire da un segnale tipo rumore bianco uniforme, variando la risposta impulsiva della sala ogni 3 secondi. Il segnale desiderato così generato per le simulazioni è mostrato in figura 11. Per ciascun segmento di segnale l’algoritmo è in grado di ricavare la risposta impulsiva, ma ciascuna variazione della risposta impulsiva comporta un nuovo periodo di inseguimento durante il quale l’algoritmo non è in grado di cancellare l’eco.

23

IMPLEMENTAZIONE ERLE

35 30 25 20 (dB)

24

15 10 5 0 −5 0

1

2

3

4 5 tempo(sec)

6

7

8

9

Figura 12: Inseguimento della risposta impulsiva che cambia ogni 3 secondi. Nell’ordine le risposte impulsive sono h12 , hs e h11 , rispettivamente con indice di densità ξ pari a 0.505, 0.632 e 0.531.

La curva mostrata in figura 12 è da ritenersi pessimistica, poiché sebbene in realtà la risposta impulsiva possa variare molto rapidamente, la variazione istantanea applicata nelle simulazioni è solo frutto di un’astrazione matematica.

2.5.3 Coerenza delle sorgenti Come ultima simulazione, si consideri un sistema i cui segnali in ingresso mostrano un diverso grado di correlazione. Tali segnali possono essere generati a partire da un rumore bianco uniforme ω, con un’operazione non lineare, il valore assoluto ad esempio si presta bene allo scopo. Due segnali non correlati possono essere dunque ricavati dal modulo, ω + |ω | 2 ω − |ω | x2 = 2

x1 =

come somma o come differenza. Tali segnali avranno un’ampiezza positiva pari all’ampiezza di ω, mentre il limite inferiore sarà l’asse delle ascisse. Per avere una correlazione variabile è possibile sommare lo stesso segnale di rumore bianco di partenza ω, un opportuno peso κ permette di variare la correlazione. Eliminando il valor medio κ/2 si ottiene ω + |ω | κ − 2 2 ω − |ω | κ x2 = ω + κ − , 2 2 x1 = ω + κ

2.5 PRESTAZIONI ERLE

50 40

(dB)

30 20 10 0 −10 0

1

2 3 tempo(sec)

4

5

Figura 13: Echo Return Loss Enhancement (ERLE) per un sistema MISO a 2 ingressi con correlazione variabile, per vari livelli di correlazione. L’equazione degli ingressi è la (2.5), i valori del parametro di correlazione κ sono rispettivamente κ = 5 (—), κ = 2 (—), κ = 1 (—), κ = 0(—) Disallineamento 0

(dB)

−2 −4 −6 −8 −10 0

1

2

3

4 0 tempo(sec)

1

2

3

4

Figura 14: Misura del disallineamento al variare della correlazione degli ingressi per un algoritmo IPNLMS MISO

e normalizzando infine i segnali si ottiene la seguente   ω + |ω | κ 1 x1 = ω + κ − κ 2 2 2 +1   ω − |ω | κ 1 x2 = ω + κ + , κ 2 2 2 +1

(2.5) (2.6)

espressione di due segnali normalizzati, con correlazione variabile in base al parametro κ, di ampiezza compresa tra −1 e 1.

25

BIBLIOGRAFIA [1] Benesty, J., Gänsler, T., Huang, Y. A. e Rupp, M. (2004), «Adaptive algorithms for MIMO acoustic echo cancellation», in «Audio signal processing for next generation multimedia communication systems», p. 120–147, Kluwer academic. [2] Huang, Y., Benesty, J. e Chen, J. (2006), Acoustic MIMO signal processing, Springer. [3] Uncini, A. (2006), Audio digitale, McGraw-Hill, Milano. [4] Uncini, A. (2008), «Filtraggio adattativo e circuiti intelligenti aggiornato», Appunti del corso di circuiti e algoritmi per l’elaborazione dei segnali.

27