%Demo done during Session 5 illustrating
%adaptive cancellation of narrowband interferer.
%A tone is added to a short speech segment.
clear all
clf
set(0,'defaultaxesfontsize',22);
M=16;
omega0=pi/4;
%fp=fopen('0ef1s1t0.wav','r');
%[data,length]=fread(fp,[1,1024],'char');
%[data,length]=fread(fp,[1,inf],'short');
%fclose(fp);
data = getspeech('0wf1s1t0.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;
tonepwr=24*rdd0;
x=datar+tonepwr*cos(omega0*(0:dsize-1));
input('Original Utterance');
soundsc(datar,Fs)
input('Utterance with tone');
soundsc(x,Fs)
D=200;
%estimation of optimal filter coefficient vector
%based on block time-averaged estimates of Rxx and rdx
[A,xsize]=size(x);
d=x(D-1:xsize);
Rxx=zeros(M,M);
rdx=zeros(M,1);
w=1; alpha=1/w;
for n=1:dsize-M-D;
Rxx=Rxx+x((n+M-1):-1:n).'*x((n+M-1):-1:n);
rdx=rdx+d(1,n)*x((n+M-1):-1:n).';
end
%initial filter coefficient vectors for LMS and RLS
hest=inv(Rxx)*rdx;
hlms=[zeros(M/2,1) ; 1;zeros((M/2-1),1)];
hrls=[zeros(M/2,1) ; 0.5;zeros((M/2-1),1)];
%hrls=zeros(M,1);
rxx0=x*x.'/xsize;
mu=2/(6*M*rxx0);
Rinv=0.1*eye(M,M);
%iterating over successive data blocks of length M
for n=1:dsize-D;
y(1,n)=d(1,n)-x((n+M-1):-1:n)*hest;
%LMS adaptation
errlms=d(1,n)-x((n+M-1):-1:n)*hlms;
ylms(1,n)=errlms;
hlmsold=hlms;
hlms=hlmsold+mu*errlms*x((n+M-1):-1:n).';
%RLS adaptation
hrlsold=hrls;
errrls=d(1,n)-x((n+M-1):-1:n)*hrlsold;
Rinvold=Rinv;
f=Rinvold*x((n+M-1):-1:n).';
murls=x((n+M-1):-1:n)*f;
Kvec=f/(w+murls);
Rinv=alpha*(Rinvold-Kvec*f.');
hrls=hrlsold+errrls*Kvec;
yrls(1,n)=d(1,n)-x((n+M-1):-1:n)*hrls;
end
input('Utterance with noise removed via block')
input('estimate of optimal adaptive filter');
soundsc(y,Fs)
input('Utterance after LMS adaptive filtering')
soundsc(ylms,Fs)
input('Utterance after RLS adaptive filtering')
soundsc(yrls,Fs)
hnotch=[1 zeros(1,D-1) -hest.'];
freqomega=linspace(-pi,pi,512);
H0=20*log10(abs(fftshift(fft(hnotch,512))));
H0mag=H0-max(H0);
plot(freqomega,H0mag,'k','LineWidth',2);
axis([-pi pi -50 10]);
title('DTFT of h[n;omega0]'),...
xlabel('Omega'), ylabel('Magnitude'), grid
hold on
plot([omega0 omega0],[-50 10],'b--','LineWidth',1);
plot([-omega0 -omega0],[-50 10],'b--','LineWidth',1);
hold off