%Matlab program for use in conjunction with Homework no. 1
%on adaptive noise cancellation.  Filtered noise is added 
%to a sinewave burst.  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]);
h =[0.0034    0.0125   -0.0030   -0.0188   -0.0096    0.0270    0.0351 ...
-0.0231   -0.0812   -0.0153    0.1925    0.3864    0.3864    0.1925 ...
-0.0153   -0.0812   -0.0231    0.0351    0.0270   -0.0096   -0.0188 ...
-0.0030    0.0125    0.0034];
N=512;
A=1;  %amplitude of desired sinewave
dsize=N;
omega0=pi/8;
ordinate=0:dsize-1;
datar=cos(omega0*ordinate);
%rdd0=datar*datar'/dsize;
%noisepwr=24*rdd0;
noisepwr=24*(A^2/2);
noise=randn(1,dsize+2*M); remote=noisepwr*noise;
%Option:  Simulate correlated noise at remote sensor,
%as opposed to white noise, by running white noise
%through a filter.  Program uses a filter with a triangular impulse response.
Nmax=5;
corrfilt=[1:Nmax Nmax-1:-1:1]/Nmax;
filtnoise=conv(randn(1,dsize+2*M),corrfilt);
corrnoise=filtnoise(1,Nmax+1:dsize+2*M-Nmax);
remote=noisepwr*corrnoise;
%
for n=1:dsize;
x(1,n)=datar(1,n)+h*remote((n+M-1):-1:n).';
end
plot(ordinate,datar,'k','Linewidth',3);
hold on
plot(ordinate,x,'r','Linewidth',3);
legend('Desired Signal','Noise Corrupted Signal');
title('Desired Signal and Desired Signal Plus Noise');
hold off
pause
%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) ; 0.5; zeros((M/2-1),1)];
hrls=[zeros(M/2,1) ; 0.5; zeros((M/2-1),1)];
%hrls=zeros(M,1);
rww0=remote*remote.'/max(size(remote));
mu=1/(2*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).';
errhestlms(1,n)=(hlms-h')'*(hlms-h');
%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;
errhestrls(1,n)=(hrls-h')'*(hrls-h');
end
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('Adaptively Estimated vs True Sinewave');
pause
plot(datar,'k','Linewidth',3);
hold on
plot(ylms,'r','Linewidth',3);
plot(yrls,'m','Linewidth',3);
hold off
legend('Uncorrupted Sinewave','LMS Estimate','RLS Estimate');
title('Estimated versus True Impulse Response');
pause
plot(errhestlms,'r','Linewidth',3);
hold on
plot(errhestrls,'m','Linewidth',3);
hold off
legend('LMS Error','RLS Error');
title('Filter Estimation Error vs Time');