%Demo done during Session 4 illustrating
%adaptive noise cancellation.
%Filtered noise is added to a short speech segment.
%An unfiltered version of the noise is picked up
%at a "remote" sensor.
clear all
clf
set(0,'defaultaxesfontsize',22);
M=24;
h=remez(M-1,[0 .333 .5 1],[1 1 0 0]);
data=getspeech('00f1s2t0.wav');
Fs=12500;
plot(data)
end1=input('end1?');
end2=input('end2?');
datar=data(end1:end2);
datar=datar/max(abs(datar));
[d,dsize]=size(datar);
rdd0=datar*datar'/dsize;
noisepwr=24*rdd0;
remote=noisepwr*randn(1,dsize+2*M);
for n=1:dsize;
x(1,n)=datar(1,n)+h*remote((n+M-1):-1:n).';
end
input('Original Utterance');
soundsc(datar,Fs)
input('Utterance with noise');
soundsc(x,Fs)
%estimation of optimal filter coefficient vector
%based on block time-averaged estimates of Rxx and rdx
Rww=zeros(M,M);
rxw=zeros(M,1);
w=1; alpha=1/w;
%initial filter coefficient vectors for LMS and RLS
hlms=[zeros(M/2,1) ; 1;zeros((M/2-1),1)];
hrls=[zeros(M/2,1) ; 0.1;zeros((M/2-1),1)];
%hrls=zeros(M,1);
rww0=remote*remote.'/dsize;
%mu=2/(M*rxx0);
mu=2/(6*M*rww0);
Rinv=0.1*eye(M,M);
%iterating over successive data blocks of length M
for n=1:dsize;
%LMS adaptation
errlms=x(1,n)-remote((n+M-1):-1:n)*hlms;
ylms(1,n)=errlms;
hlmsold=hlms;
hlms=hlmsold+mu*errlms*remote((n+M-1):-1:n).';
%RLS adaptation
hrlsold=hrls;
errrls=x(1,n)-remote((n+M-1):-1:n)*hrlsold;
Rinvold=Rinv;
f=Rinvold*remote((n+M-1):-1:n).';
murls=remote((n+M-1):-1:n)*f;
Kvec=f/(w+murls);
Rinv=alpha*(Rinvold-Kvec*f.');
hrls=hrlsold+errrls*Kvec;
yrls(1,n)=x(1,n)-remote((n+M-1):-1:n)*hrls;
if n==1
lh=plot(h,'k','Linewidth',3);
hold on
lhlms=plot(hlms,'r','Linewidth',3,'erasemode','xor');
lhrls=plot(hrls,'m','Linewidth',3,'erasemode','xor');
hold off
legend('True h[n]','LMS Estimate','RLS Estimate');
title('Estimated versus True Impulse Response');
th = xlabel(sprintf('n=%g',n));
set(th,'erasemode','xor')
drawnow
elseif n< 500
    %set(lh,'ydata',h)
    set(lhlms,'ydata',hlms)
    set(lhrls,'ydata',hrls)
	set(th,'string',sprintf('n=%g',n))
end
end
input('Utterance after LMS adaptive filtering')
soundsc(ylms(300:dsize),Fs)
input('Utterance after RLS adaptive filtering')
soundsc(yrls(300:dsize),Fs)
plot(h,'k','Linewidth',3);
hold on
plot(hlms,'r','Linewidth',3);
plot(hrls,'m','Linewidth',3);
hold off
legend('True h[n]','LMS Estimate','RLS Estimate');
title('Estimated versus True Impulse Response');