% 2-dimensional MMSE demonstration
% ECE4703 5-Dec-2006

format compact
unknownsystem = 2*rand(2,1)-1	% generate random unknown system
btest = -3:0.1:3;		% range of filter coefficients

x = randn(1000,1);		% white noise input signal
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),'wp');
xlabel('b0');
ylabel('b1');
title('Mean Squared Error contour');
axis image
hold off
