PRÀCTICA 3

SISTEMES ADAPTATIUS

1.- FILTRAT DE WIENER I SISTEMES ADAPTATIUS

FILTRAT DE WIENER

El filtrat de Wiener és un cas particular dels filtrats adaptatius. El filtrat de Wiener intenta anar modificant els pesos d'un filtre FIR de tal forma que la potència del senyal de sortida e[n] sigui mínima. Així doncs, el problema queda limitat a poder calcular de forma continuada el valor dels coeficients del filtre FIR tals que aconsegueixen reduir al màxim la potència a la sortida. Aquest càlcul es fa a partir dels senyals S[n] i r[n].

Tenim un filtre FIR de n coeficients que el representem com:

Per tant tenim que la senyal a la sortida del filtre FIR serà:

On SN[n] és una finestra de N mostres de la senyal d'entrada S[n] i la podem representar com:

Per altre banda tenim que:

Si calculem la potència de la senyal d'error e[n] tenim el què s'anomena la funció de cost:

Si volem minimitzar l'error quadràtic cal minimitzar la funció de cost i per fer-ho podem calcular-ne la derivada en totes les n dimensions, és a dir, el gradient.

El gradient de la funció de cost té la forma d'una paràbola de n dimensions. És a dir, tenim un problema de minimització amb tantes dimensions com pesos té el filtre. En el cas d'utilitzar un filtre de dos pesos tenim que el problema es pot esquematitzar com:

En el cas general de n pesos tenim que en calcular el mínim igualant el gradient a zero i obtenim:

IDENTIFICACIÓ DE SISTEMES

El filtrat de Wiener es pot utilitzar per identificar sistemes incògnits. Per tal de poder identificar sistemes incògnits ens cal muntar una estructura del següent tipus:

Figura 1. Esquema utilitzat per la identificació de sistemes

Tenim que el desenvolupament obtingut per al filtre de Wiener és un desenvolupament que aconsegueix obtenir una potència a la sortida mínima. Amb la topologia de la figura anterior podem veure que si aconseguim tenir a la sortida una potència nul·la tindrem que el filtre de Wiener és idèntic al sistema incògnita. En el cas que el sistema incògnita sigui un filtre FIR podrem arribar a fer una identificació perfecte si l'ordre del filtre de Wiener és superior al del filtre incògnita. En el cas que el sistema incògnita sigui un filtre IIR no el podrem arribar a identificar perfectament mai encara que tinguem molts coeficients.

A continuació vaig a demostrar que es pot identificar un sistema introduint soroll a l'entrada del sistema incògnita i utilitzant un filtre de Wiener d'ordre superior al sistema incògnita.

Per a realitzar aquesta demostració he utilitzat el següent programa en Matlab:

Ordre=40; %Ordre del filtre identificat

FirInc=randn(1,50); %Coeficients d'un FIR incògnita

x=randn(1,500); %Entrem com a senyal soroll

r=filter(FirInc,1,x);

rx=xcorr(x,x,'biased'); %Calculem Rxx real

rxx=rx(length(rx)/2:length(rx)/2 + Ordre-1);%Ens quedem amb part dreta de autocorrelació

rxr=xcorr(x,r,'biased'); %Calculem correlació creuada

rxr=rxr(length(rxr)/2:length(rxr)/2 + Ordre-1);

Rxx=toeplitz(rxx); %Matriu rx(0) diagonal principal rx(1) segones diagonals

H=inv(Rxx)*rxr'; %Calculem els coeficients que minimitzen la sortida

figure(1); %Mostrem el filtre incògnit i el filtre identificat

subplot(2,1,1);

plot(FirInc);

hold on

plot(H,'g');

hold off

XLABEL('Coeficients dels filtres');

YLABEL('Valors dels coeficients');

subplot(2,1,2); %Mostrem l'error comès en la identificació dels coeficients

if length(H)<length(FirInc)

plot(H'-FirInc(1:length(H)));

else

plot(FirInc' - H(1:length(FirInc)));

end;

axis([0 length(FirInc) min(FirInc) max(FirInc)]);

ylabel('Error en els coeficients');

XLABEL('Coeficients dels filtres');

Els resultats obtinguts són els esperats, és a dir, si utilitzem un ordre de filtre inferior a l'ordre del sistema a identificar tenim que cometem error, mentre que si utilitzem un ordre de filtre superior a l'ordre del sistema a identificat tenim que la identificació és quasi exacta. Si fem el mateix però en comptes de intentar identificar un sistema que es comporti com un filtre FIR ho fem amb un sistema que es comporti com un filtre IIR tenim que, independentment de l'ordre del sistema a identificar i de l'ordre del filtre que nosaltres utilitzem per identificar el sistema, la identificació no és mai exacta i sempre cometem una mica d'error.

Figura 2. Identificació d'un sistema incògnita d'ordre superior al filtre utilitzat en la identificació

Figura 3. Identificació d'un sistema incògnita d'ordre inferior al filtre utilitzat en la identificació

LIMITACIONS DEL FILTRAT DE WIENER

Les limitacions a l'hora d'utilitzar el filtrat de Wiener són bàsicament dos: en primer lloc la gran potència de càlcul que és necessària per tal de realitzar el càlcul d'una matriu inversa i per altre banda el fet d'haver de calcular una autocorrelació que està definida entre menys infinit i més infinit.

Cal tenir en compte que si els senyals d'entrada no són estacionaris (cas més probable) tenim que les autocorrelacions aniran variant en el temps i per tant tindrem que haurem de calcular hN* repetidament. Per aplicacions que necessiten el temps real no podem utilitzar dons el filtrat de Wiener i per tant s'utilitza un altre mètode que permet calcular els coeficients òptims del filtre (hN* ) mostra a mostra sense haver de calcular les autocorrelacions. Aquests mecanismes s'anomenen sistemes adaptatius.

Càlcul d'una matriu inversa

El fet d'haver d'invertir una matriu implica haver de disposar d'una potència de càlcul molt elevada que tant sols es veu reduïda si la matriu és del tipus toeplitz. Aquest inconvenient que pot semblar no gaire important és un factor de limitació important a l'hora de la seva utilització ja que impossibilita la seva utilització en totes aquelles aplicacions que necessiten el càlcul en temps real i per tant el filtrat de Wiener queda reduït a la utilització en llocs on el temps de càlcul no és molt rellevant.

Càlcul d'una autocorrelació

Com ja he dit anteriorment, per tal de calcular una autocorrelació verdadera cal realitzar una integral de menys a més infinit i per tant és impossible realitzar-la físicament. La solució està en fer el càlcul aproximat de l'autocorrelació mitjançant una finestra. El problema d'aquest mètode és que la solució no és una autocorrelació estrictament dita i, per tant, la matriu resultant no és toeplitz i per tant el càlcul de la matriu inversa es complica molt.

Per tal de fer una estimació de l'autocorrelació tenim dos mètodes:

* Càlcul de l'autocorrelació utilitzant una finestra constant

Aquest mètode es basa en utilitzar una finestra de mida constant del senyal a correlar. Utilitzem tota la senyal per realitzar l'autocorrelació però només contemplem les mostres que es troben dins la finestra. Aquest mètode ens dona una bona aproximació a l'autocorrelació però el problema és que no es tracta d'una autocorrelació pròpiament dita i per tant el resultat no és un vector parell i per tant no es pot utilitzar com una matriu toeplitz. El mecanisme es basa en:

Figura 4. Estimació d'una autocorrelació utilitzant finestra constant

* Càlcul de l'autocorrelació considerant tant sols el senyal enfinestrat

Aquest mètode té la avantatge que si que calculem una autocorrelació verdadera i per tant el resultat es un vector parell que pot utilitzar en una matriu toeplitz i ens facilita el càlcul de la matriu inversa. El problema és que estem calculant l'autocorrelació d'un senyal que no és el senyar real, sinó que és un senyal que coincideix amb el senyal real durant la finestra i la resta val zeros. El mètode que s'utilitza per calcular l'autocorrelació d'aquesta manera és:

Figura 5. Càlcul de l'autocorrelació d'una senyal prèviament enfinestrada

SISTEMES ADAPTATIUS

Un sistema adaptatius són aquells que adapten els seus coeficients mostra a mostra o bloc a bloc en funció del senyal d'entrada. L'esquema bàsic de tot sistema adaptatiu és:

Figura 6. Esquema d'un sistema adaptat

Un filtre adaptatiu és un filtre digital que té la capacitat d'adaptar les seves característiques en el temps.

Usualment desitgem separar senyal i soroll que ocupen bandes de freqüències diferents, en aquests casos un filtre lineal amb coeficients fixes ens solucionarà el problema. A vegades però, és necessari que les característiques del filtre siguin variables adaptant-se a la variació de les característiques del senyal en cada instant de temps. Aquests filtres tenen la característica que els coeficients es van modificant al llarg del temps i, per tant, també la seva resposta freqüencial es veu alterada al llarg del temps.

Un dels sistemes adaptatius millors és l'anomenat LMS (Least Mean Square) que es basa en anar ajustant els coeficients mostra a mostra de tal forma que es minimitzi l'error quadràtic mig.

La idea és senzilla: en comptes de calcular on es troba el mínim en la funció de cost es tracta de realitzar una aproximació recursiva fins a trovar-lo. Es tracta doncs de realitzar una espècie de aproximació per càlcul numèric (semblant a Boltzano) fins a trobar la hN* però sense calcular-la. El funcionament de l'algorisme es basa en anar modificant els pesos del filtre en sentit descendent del gradient i per tant anirem en direcció al mínim.

Això permet evitar els càlculs de les matrius d'autocorrelació i per tant estalviar molts càlculs de tal forma que permeten implementar aquest algoritme en temps real.

Els pesos del filtre es modifiquen mostra a mostra de la següent forma:

Fixar tots els pesos a un valor inicial com per exemple 0

Calcular la sortida del filtre com:

Calcular l'error estimat:

Modificar els següents pesos del filtre

(5) Tornar al punt (2)

A la pràctica tenim que mitjançant aquesta metodologia mai aconseguim el valor obtingut per Wiener, però la solució obtinguda fluctua entorn la solució òptima (la calculada mitjançant Wiener). Per altre banda, el sistema serà estable mentre el pas d'aprenentatge () estigui comprès entre 0 i el màxim valor de la matriu de covariança de la senyal d'entrada.

La funció que he utilitzat que calcula els coeficients òptims mitjançant el mètode de LMS és:

function [Hnou,err]=LMS(S,R,E,Ordre,h)

if nargin==4

h=zeros(1,Ordre); %Definim els coeficients si no els hem passat

VecSen=zeros(1,Ordre); %Definim els finestra de senyal

end;

%Implementació de l'algorisme LMS

for i=1:length(S) %Realitzem càlcul mostra a mostra per totes les mostres

if (i>Ordre) %Només podem començar quan i>Ordre per no accedir a vector(<0)

YN=0;

for j=1:Ordre

YN=YN+h(j)*S(i-j+1);

end;

err(i)=R(i)-YN; %Calculem error

for j=1:Ordre %Calculem nous coeficients

h(j)=h(j)+E*err(i)*S(i-j+1);

end;

end;

end;

Hnou=h; %Retornem valors calculats

Per tal de fer la identificació del sistema a part he fet:

Ordre=20; %Definim ordre del sistema identificació

E=0.005; %Valor pas d'aprenentatge

Senyal=randn(1,1000); %Entrem soroll

FirInc=randn(1,20); %Filtre incògnita

r=filter(FirInc,1,Senyal); %Obtenim la referència r

[h,err]=LMS(Senyal,r,E,Ordre); %Realitzem identificació

figure(1); %Mostrem els resultats

subplot(3,1,1);

plot(FirInc);

ylabel('Coeficients filtre incògnita');

subplot(3,1,2);

plot(h);

ylabel('Coeficients filtre identificat');

xlabel('Coeficients');

subplot(3,1,3);

plot(err);

xlabel('Mostres');

ylabel('Error');

El que he pogut observar és que si utilitzo un valor petit per al pas d'aprenentatge () tenim que el sistema tarda més en convergir però obtenim un error residual molt petit, mentre que si utilitzo un pas d'aprenentatge () gran el sistema és molt més ràpid en convergir però el valor residual d'error és més gran.

A continuació tenim una identificació d'un sistema FIR d'ordre 20 (utilitzant soroll blanc a l'entrada) mitjançant l'algorisme LMS i amb diferents valors per el pas d'aprenentatge:

Figura 7. Identificació d'un sistema FIR mitjançant l'algorisme LMS amb =0,005

En aquesta identificació obtenim un error residual de 1,094E-4

Figura 8. Identificació d'un sistema FIR mitjançant l'algorisme LMS amb =0,05

En aquest cas he obtingut un error residual molt major que en el cas anterior 0,77. Observem per altre banda que en aquest cas només calen 200 mostres per tal d'haver adaptat el sistema, mentre que en la gràfica anterior calien unes 700 mostres per adaptar el sistema.

En el cas d'introduir senyal de veu obtenim resultats similars, però amb la diferència que cal utilitzar un valor per al pas d'aprenentatge molt menor per tal que el sistema sigui estable.

A continuació tenim les gràfiques de com evoluciona la potència de la senyal d'error en funció del pas d'aprenentatge. La primera gràfica és quan apliquem a l'entrada soroll, mentre que la segona és quan apliquem senyal de veu a l'entrada.

Figura 9. Evolució de la potència d'error en funció de  amb soroll a l'entrada

Figura 10. Evolució de la potència d'error en funció de  amb senyal de veu a l'entrada

Observem que les dues gràfiques són idèntiques i amb la única diferència que per tal de trobar els coeficients amb senyal de veu a l'entrada cal un valor del pas d'aprenentatge molt inferior que en el cas d'introduir soroll. Per altre banda, tenim que tenen aquesta forma ja que quan menor és  tenim que el sistema tarda més en convergir i per tant tenim més tros amb un error important, mentre que quan  és gran tenim que el sistema convergeix ràpidament i per tant la potència d'error és molt menor. No arribem mai a potència d'error zero ja que tot i que el sistema convergeixi molt ràpidament ens queda un error residual que és degut a que els pesos del filtre calculat van fluctuant entorn als valors òptims.

Per realitzar aquestes gràfiques he realitzat els següents passos:

Ordre=20; %Definim ordre

Senyal=wavread('test.wav'); %Apliquem a l'entrada senyal de veu

Senyal=Senyal(length(Senyal)-10000:length(Senyal));

Senyal=Senyal-mean(Senyal);

FirInc=randn(1,20); %Definim el filtre a identificar

r=filter(FirInc,1,Senyal); %Calculem r

i=0;

for E=0.000005:0.000001:0.00005 %Realitzem escombrat en pas d'apren.

[h,err]=LMS(Senyal,r,E,Ordre); %Calculem h i err

i=i+1;

Pot(i)=err*err'; %Calculem potència d'error

pua(i)=E;

end;

figure(1);

plot(pua,Pot); %Mostrem resultat

xlabel(`Pas de aprenentatge');

ylabel(`Potència de error');

2.- CANCEL·LACIÓ CEGA D'ECOS

En aquest apartat he de realitzar els algorismes de generació i cancel·lació d'eco no reverberant. En primer lloc m'agradaria explicar els dos tipus d'eco que existeixen:

Eco reverberant

L'eco reverberant es caracteritza per tenir diverses repeticions de la senyal original retardades i amb potència cada cop inferior. Podem modelar l'eco reverberant com:

Com podem observar tenim que la senyal amb eco està composta per la senyal sense eco amb diverses rèpliques retardades i cada cop amb menys amplitud. Aquest eco apareix en tot recinte tancat ja que les reflexions de la veu en les parets del recinte provoquen que ens arribin amb retard rèpliques de la senyal.

Eco no reverberant

L'eco no reverberant es caracteritza per tenir tant sols una repetició retardada del senyal. Per tant podem modelar l'eco no reverberant com:

Per tant tenim que la funció de transferència d'un generador d'eco no reverberant és:

GENERACIÓ D'ECO NO REVERBERANT

Si volem realitzar una funció generadora d'eco no reverberant ens cal tant sols implementar aquest filtre i aplicar-lo a alguna senyal d'entrada. A continuació tenim el codi que he utilitzat per fer una funció que generi eco no reverberant:

function [y]=EcoNoRev(Senyal,RetEco,AmpEco)

% function [y]=EcoNoRev(Senyal,RetEco,AmpEco)

% Genera eco no revereberant

% RetEco => Retard de l'eco

% AmpEco => Amplitud de l'eco

Tmp=[1,zeros(1,RetEco-2),AmpEco]; % Defineix el filtre

y=filter(Tmp,1,Senyal); % Aplica el filtre

Si apliquem soroll blanc a l'entrada del filtre i fem l'autocorrelació a la sortida tenim:

Figura 11. Autocorrelació de soroll blanc amb eco no reverberant

Com podem observar en l'autocorrelació tenim que hi ha una repetició dels senyal d'entrada (Soroll) al cal de 100 mostres i que té una potència inferior a la potència del senyal pròpiament dit. Aquesta gràfica l'he obtingut introduint el següent codi:

Retard=100;

Amplitud=0.5;

Senyal=randn(1,1000); % Senyal d'entrada

SenyalEco=EcoNoRev(Senyal,Retard,Amplitud); % Realitza el filtrat per afegir eco

z=xcorr(SenyalEco,SenyalEco); % Calcula autocorrelació

pua=-length(z)/2+1:length(z)/2; % Realitza canvi d'eix per zero al centre

figure(1);

plot(pua,z); % Mostra els resultats

TITLE('Rxx DEL SENYAL AMB ECO NO REVERBERANT');

Com podem observar en la gràfica de la figura 11, tenim que la relació de potència entre el senyal i l'eco és d'aproximadament 0,25 i això és degut a que la relació d'amplitud entre el senyal i l'eco és de 0,5 i per tant la relació en potència és de 0,52.

CANCEL·LACIÓ D'ECO NO REVERBERANT

Per tal de generar una funció de cancel·lació d'eco no reverberant a ceges he utilitzat l'algorisme LMS però seguint l'estructura del cancel·lador que és la següent:

Figura 12. Esquema del cancel·lador d'eco no reverberant

El programa que he utilitzat que implementa aquest esquema és el següent:

function [SNet,hn]=Cancel(SBrut,RetMin,E,Ordre,h)

% [SNet,hn]=Cancel(SBrut,E,Ordre)

% Cancel·la eco no reverberant

if nargin==4

h=zeros(1,Ordre); %Definim els coeficients si no els hem passat

end;

err=zeros(1,length(SBrut));

%Implementació de l'algorisme LMS

for i=1:length(SBrut) %Realitzem calcul per totes les mostres

if (i>max(Ordre,RetMin))

S(i)=err(i-RetMin);

YN=0;

for j=1:Ordre

YN=YN+h(j)*S(i-j+1);

end;

err(i)=SBrut(i)-YN;

for j=1:Ordre

h(j)=h(j)+E*err(i)*S(i-j+1);

end;

end;

end;

SNet=err; %Retornem senyal amb eco cancel·lat

hn=h; %Retornem coeficients del filtre estimat

COMPORTAMENT DEL CANCEL·LADOR AMB SOROLL A L'ENTRADA

Tenim el nostre cancel·lador d'ecos al que li apliquem soroll blanc al que prèviament li hem afegit un eco no reverberant. Al afegir un eco no reverberant al soroll tenim que en l'autocorrelació d'aquest soroll apareixeran dos pics que corresponen un a la potència de soroll i l'altre a l'eco. Un cop hem aplicat aquest soroll al cancel·lador, aquest intentarà minimitzar la potència de senyal a la sortida i com que hi apliquem un retard mínim no podrà minimitzar-se mitjançant la cancel·lació amb el propi senyal per tant no tindrà altre solució que cancel·lar l'eco. Si fem aquest procés obtenim les següents gràfiques:

Figura 13. Autocorrelacions del senyal a l'entrada i a la sortida del cancel·lador

Figura 14. Filtre FIR identificat en el procés de cancel·lació

Com podem observar tenim que aquest filtre s'assembla al filtre que estem buscant que hauria de ser:

En aquest cas hem introduït un eco d'amplitud 0,5 (A=0,5) i de retard 30 (RetEco=30). Per altre banda he utilitzat un retard mínim de 10 mostres i per tant podem observar que en el coeficient 20 tenim un valor de 0,5 aproximadament. El coeficient 20 és el coeficient que fa que es cancel·li el retard situat a una distància de 30 mostres si el retard mínim és de 10 mostres. Com es pot observar tenim que la resta de coeficients que haurien de ser zero estan entorn a zero però no són exactament zero.

El programa utilitzat per realitzar aquestes gràfiques ha estat:

Ts=1/8000;

Retard=30; %Posició del eco

Amplitud=0.5;

%Lectura de la senyal d'entrada

Senyal=randn(1,1000);

SenyalEco=EcoNoRev(Senyal,Retard,Amplitud); %Apliquem eco no reverberant

Ordre=50; %Ordre del filtre a identificar

RetMin=10; %Retard mínim

E=0.005; %Pas d'aprenentatge

[SNet,hn]=Cancel(SenyalEco,RetMin,E,Ordre); %Cancel·lem eco

COMPORTAMENT DEL CANCEL·LADOR AMB SENYAL DE VEU AMB ECO

A continuació estudiem el comportament del cancel·lador mitjançant el senyal de veu a l'entrada. Aquí tenim un exemple dels resultats obtinguts amb aquest esquema:

Figura 15. Autocorrelacions del senyal amb eco i del senyal amb eco cancel·lat

Com es pot observar en les autocorrelacions hem aconseguit una cancel·lació bastant bona tot i que no és total. A continuació tenim una gràfica amb la senyal temporal amb eco i un cop se li ha aplicat el cancel·lador d'ecos.

Figura 16. Senyal amb eco no reverberant i senyal resultant de la cancel·lació

A continuació tenim els coeficients del filtre FIR que hem estimat amb l'algorisme:

Figura 17. Coeficients del filtre FIR estimat

Aquesta gràfica ha estat realitzada amb un retard mínim de 1400 i el retard de l'eco era de 1600 mostres (200ms). Per altre banda hi havia definit un ordre pel filtre a identificar de 400. Tenim que idealment hauríem de tenir uns coeficients que corresponguessin a la inversa del filtre de generació d'eco no reverberant i per tant hauriem de tenir:

Per una banda tenim que el retard total que obtenim és la suma del retard mínim, que en aquest cas és de 1400, més el retard estimat en els coeficients, que en aquest cas és de 200, i per tant tenim un retard total de 1600 que coincideix amb el retard de l'eco. Per altre banda tenim que l'amplitud del pic A hauria de ser igual a l'amplitud amb la que hem generat l'eco (en aquest cas era de 0,5) i per altre banda hem estimat una amplitud de 0,45 i per tant la cancel.lació de l'eco no pot ser perfecte. Per altre banda, tenim que tots els coeficients a excepció del situat al pes 1600 (Retard mínim + Retard estimat) haurien de ésser zero i en aquest cas tenim que són valors molt petits i que oscil·len entre el zero però no són zero.

Així doncs, tenim que idealment hauríem d'esperar obtenir els següents coeficients:

Figura 18. Coeficients del filtre FIR que esperava obtenir

Com podem observar tenim que els valors que esperàvem obtenir i els valors obtinguts tampoc són molt semblants però no iguals. És per això que no cancel·lem adequadament el 100% de l'eco que teniem en el senyal original.

Per obtenir aquestes gràfiques he utilitzat el següent codi de programa:

Ts=1/8000;

Retard=1600; %200ms a 8KHz

Amplitud=0.5;

%Lectura de la senyal d'entrada

Senyal=wavread('test.wav');

Senyal=Senyal-mean(Senyal);

SenyalEco=EcoNoRev(Senyal,Retard,Amplitud); %Generem eco no reverberant

Ordre=400;

RetMin=1400;

E=0.000001;

[SNet,hn]=Cancel(SenyalEco,RetMin,E,Ordre); %Cancel·lem l'eco no reverberant

figure(1); %Representem les autocorrelacions

subplot(2,1,1);

z=xcorr(SenyalEco,SenyalEco);

Dummy=-length(z)/2+1:length(z)/2;

plot(Dummy,z);

ylabel('xcorr del senyal amb eco');

subplot(2,1,2);

z=xcorr(SNet,SNet);

Dummy=-length(z)/2+1:length(z)/2;

plot(Dummy,z);

ylabel('xcorr del senyal netejat');

figure(2); %Representem els senyals en temps

t=Ts:Ts:length(SenyalEco)*Ts;

plot(t,SenyalEco,'g');

hold on

plot(t,SNet);

hold off;

xlabel('Temps (segons)');

figure(3); %Representem els coeficients del filtre obtingut

plot(hn);

TITLE('PESOS DEL FILTRE IDENTIFICAT');

XLABEL('Posició dels pesos');

sound(SenyalEco)

sound(SNet)

ESTUDIAR EL COMPORTAMENT AL VARIAR L'ESTIMACIÓ DE LA POTÈNCIA

En aquest cas ens cal estimar la potència del senyal d'entrada a partir d'una finestra de N mostres. Amb aquesta potència i una constant de proporcionalitat calcularé el valor del pas d'aprenentatge que s'ha d'utilitzar en cada mostra. L'avantatge d'aquest mètode és que es va adaptant el valor del pas d'aprenentatge a la potència del senyal de tal forma que en els llocs on la potència de senyal sigui gran tindrem un valor pel pas d'aprenentatge petit i per tant el sistema continuarà essent estable, mentre que en els punt on la potència del senyal d'entrada sigui petita tindrem que el valor del pas d'aprenentatge podrà ser superior mantenint la estabilitat del sistema i per tant tindrem un sistema millor.

He modificat el programa cancel·lador d'ecos per tal d'afegir l'estimació de la potència i el càlcul del valor d'aprenentatge mostra a mostra en funció de la potència estimada.

function [SNet,hn]=Cancel(SBrut,RetMin,K,Ordre,h)

% [SNet,hn]= Cancel(SBrut,RetMin,K,Ordre,h)

% Cancel·la eco no reverberant amb estimació de potència

Pot=0; %Inicialitza potència estimada

LongVent=100;

Beta=(LongVent-1)/LongVent;

............

err(i)=SBrut(i)-YN;

if i<LongVent %Estimació de potència

Pot=((i-1)*Pot+SBrut(i-RetMin)*SBrut(i-RetMin))/i;

else

Pot=Beta*Pot + (1 - Beta)*SBrut(i-RetMin)*SBrut(i-RetMin);

end;

if Pot>0 %Calcula pas d'aprenentatge en funció de Pot

E=K/(Ordre*Pot);

else

E=0;

end;

for j=1:Ordre

h(j)=h(j)+E*err(i)*S(i-j+1);

end;

..............

A continuació tenim una gràfica que mostra l'evolució dels coeficients del filtre FIR identificat en funció del valor del paràmetre K utilitzat en el valor del pas d'aprenentatge.

Figura 19. Coeficients del filtre FIR per diferents valors de K

Tenim que per valors petits de K ens cal una durada molt més elevada de mostres de la senyal d'entrada per tal d'obtenir una estimació acceptable dels coeficients. Per altre banda, al augmentar el valor de K tenim que cada cop obtenim una estimació més adequada del coeficient a on hi ha el retard, però per altre banda, tenim més dispersió entre els coeficients que haurien de ser zero.

El programa que he utilitzat per obtenir la gràfica anterior és:

Retard=30; %Característiques de l'eco introduït

Amplitud=0.5;

%Lectura de la senyal d'entrada

Senyal=randn(1,1500);

SenyalEco=EcoNoRev(Senyal,Retard,Amplitud); %Generació de l'eco

Ordre=50;

RetMin=10;

figure(1);

hold on;

for K=1/1000:1/100:1/10

[SNet,hn]=Cancel(SenyalEco,RetMin,K,Ordre); %Cancel·lació de l'eco

plot(hn); %Representació dels coeficients obtinguts

end;

TITLE('PESOS DEL FILTRE IDENTIFICAT');

XLABEL('Posició dels pesos');

hold off;

NETEJAR TANT COM SIGUI POSSIBLE LA SENYAL DE VEU

Ara tenim que no cal tenir temps real i per tant podem analitzar tota la senyal de veu primer i en segon lloc fer el processat. He decidit inventar una estratègia pròpia per tal de treure el eco. L'estratègia es basa en identificar automàticament l'eco a partir de l'autocorrelació i extreure'n els seus paràmetres (retard i amplitud) i un cop els tenim aplicar el filtre cancel·lador d'eco no reverberant suposant que coneixem aquests paràmetres. L'eliminació que s'aconsegueix és òptima de tal forma que es neteja totalment el senyal.

El programa utilitzat és:

function [SNet]=ideal(SBrut,RetMin)

z=xcorr(SBrut,SBrut); %Fa autocorrelació

[yMax,xMax]=max(z); %Detecta punt màxim

IncThr=abs( max(z)-min(z) )/100; %Resolució de detecció de posició eco

%Anem baixant threshold de comparació fins a trovar eco

xEco=[];

Threshold=yMax-IncThr;

while xEco==[]

Tmp=find(z>Threshold);

pua=find( (abs(Tmp-(length(z)/2) )>RetMin) & (Tmp>length(z)/2) );

xEco=Tmp(pua);

Threshold=Threshold-IncThr;

end;

%Calculem la posició del eco

RetEco=xEco-(length(z)/2)

AmpEco=sqrt(z(xEco)/yMax)

%Apliquem filtre IIR amb els valors estimats per eliminar eco

SNet=filter(1,[1 zeros(1,RetEco-2) AmpEco], SBrut);

Funcionament de l'algorisme. En primer lloc fem l'autocorrelació del senyal amb eco. En segon lloc detectem l'amplitud de l'autocorrelació a l'origen. A continuació calculem l'increment per al llindar de comparació com una fracció de l'amplitud màxima de l'autocorrelació. En segon lloc anem mirant tots els punts de l'autocorrelació que sobrepassen un cert nivell de comparació (Threshold) i anem baixant el nivell fins que ens trobem amb l'eco. Un cop tenim la posició de l'eco ja podem conèixer-ne els seus dos paràmetres: l'amplitud i el retard. Finalment apliquem el filtre cancel·lador d'ecos amb els valors d'amplitud i retard estimats i eliminem l'eco.

Figura 20. Autocorrelació del senyal amb eco i un cop eliminat l'eco

Com es pot observar en les autocorrelacions he obtingut un senyal totalment net d'eco i per tant les estimacions són correctes. Si observem els senyals en el temps tenim:

Figura 21. Senyal amb eco i senyal net d'eco

Per tal d'obtenir les dues gràfiques anteriors he utilitzat el següent programa:

Retard=1600; %200ms a 8KHz

Amplitud=0.8;

RetMin=20;

%Lectura de la senyal d'entrada

Senyal=wavread('test.wav');

Senyal=Senyal-mean(Senyal);

SenyalEco=EcoNoRev(Senyal,Retard,Amplitud); %Afegim eco al senyal

%Aplicació de la cancel·lació ideal amb estimació dels paràmetres de l'eco

[SNet]=ideal(SenyalEco,RetMin);

sound(SenyalEco);

sound(SNet);

%Representació de les autocorrelacions

figure(1);

subplot(2,1,1);

z=xcorr(SenyalEco,SenyalEco);

Dummy=(-length(z)/2+1:length(z)/2);

plot(Dummy,z);

ylabel('xcorr del senyal original');

subplot(2,1,2);

z=xcorr(SNet,SNet);

plot(Dummy,z);

Dummy=(-length(z)/2+1:length(z)/2);

plot(Dummy,z);

ylabel('xcorr del senyal net');

%Representació temporal del senyal net i el senyal amb eco

figure(2);

plot(SenyalEco);

hold on

plot(SNet,'g');

hold off;

L'algorisme que he dissenyat és molt ràpid i funciona molt correctament. Té l'únic inconvenient que necessita disposar de tota la senyal prèviament abans de poder-la netejar i per tant no es pot implementar en temps real. Apart d'aquest inconvenient tenim que és un mecanisme molt ràpid i que permet netejar totalment qualsevol senyal de veu amb eco no reverberant.

Una altre opció és utilitzar l'algorisme LMS i aplicar-lo iteradament de tal forma que cada cop els coeficients del filtre van tendint més al filtre cancel·lador ideal i per tant la cancel·lació és millor. Aquest mètode trobo que és massa lent i que no aconsegueix resultats tant òptims com els que he obtingut amb el meu mètode.

Processat del senyal

1

24

200

400

0

0

0,5

Senyal amb eco cancel·lat

err

Senyal amb eco

r

Y

eK[n]

XK[n]

YK[n]

FILTRE WIENER

r[n]

SISTEMA ?

eK[n]

Y[n]

S[n]

FILTRE WIENER

Jmin

hN*

h(1)

h(0)

J

e[n]

r[n]

S[n]

FILTRE FIR

y[n]

S

FIR

Z-D

ECO