% 2-dimensional MMSE demonstration
% ECE4703 5-Dec-2006
% Updated on 6-Dec-2006 to show
%  - exact gradient descent
%  - approximate gradient descent (LMS)

format compact
unknownsystem = [-1,1]	% generate random unknown system
btest = -3:0.1:3;		% range of filter coefficients

x = randn(5000,1);		% white noise input signal
x = filter([0.5 1 0.5 ],1,x);  % make input non-white
d = conv(unknownsystem,x);	% output of unknown system

% main loop to compute MSE
i1 = 0;
for b0=btest,
  i2 = 0;
  i1 = i1+1;
  for b1=btest,
    i2 = i2+1;
    b = [b0;b1];		% form filter
    y = conv(b,x);		% compute filter output
    e = d-y;			% error between unknown system and our filter
    mse(i2,i1) = mean(e.^2);	% mean squared error (MSE)
  end
end

% MMSE filter computation
r00 = mean(x.*x);
r11 = mean(x.*x);
r01 = mean(x(2:end).*x(1:end-1));
r10 = mean(x(1:end-1).*x(2:end));
p0 = mean(d(1:1000).*x(1:1000));
p1 = mean(d(2:1000).*x(1:999));
R = [r00 r01 ; r10 r11];
p = [p0 ; p1];
bmmse = inv(R)*p

% plots
contourf(btest,btest,mse,20);
hold on
plot(unknownsystem(1),unknownsystem(2),'rp');
xlabel('b0');
ylabel('b1');
title('Mean Squared Error contour');
axis image

% ----------------------------------------------------------
% Exact gradient descent
% ----------------------------------------------------------
iterations = 100;
b = [-2;-3];			% initial filter guess
mu = 0.1;			% step size
plot(b(1),b(2),'w+');
bnew = b;
%pause
for i=1:iterations,
        bold = bnew;
	bnew = bold-mu*(2*R*bold-2*p);
	plot([bold(1) bnew(1)],[bold(2) bnew(2)],'w-+');
	%pause
end
plot(unknownsystem(1),unknownsystem(2),'rp');

% ----------------------------------------------------------
% Approximate gradient descent (LMS)
% ----------------------------------------------------------
b = [-2;-3];			% initial filter guess
mu = 0.1;			% step size
plot(b(1),b(2),'yx');
bnew = b;
%pause
for i=2:iterations+1,
        bold = bnew;
	yn = bold(1)*x(i)+bold(2)*x(i-1);
	en = d(i)-yn;
	bnew = bold+mu*en*[x(i);x(i-1)];
	plot([bold(1) bnew(1)],[bold(2) bnew(2)],'y-x');
	%pause
end
plot(unknownsystem(1),unknownsystem(2),'rp');

hold off
