function [h,J]=mnsed_opt(x_1,x_2,N,r_1,r_2)
%MNSED_OPT	h=mnsed_opt(x_1,x_2,N,r_1,r_2)
%	omputes an optized FIR filter of size NxN and returns the filter 
%	coefficients h.r_1 and r_2 are the 2D autocorrelation matrices
%	of x_1 and x_2, and if supplied, the x_1 and x_2 inputs may be [].
%	The filter is optimal with respect to the Mean-Normaised
%	Square Error Difference (MNSED).
%
%	Example use
%	x_1=hipsread('brodatz/D92.hips',0:127,0:127);
%	x_2=hipsread('brodatz/D4.hips',0:127,0:127);
%	h=mnsed_opt(x_1,x_2,6);

% Compute the correlation matrices
if nargin<4,
  r_1=corrcoeff(x_1,x_1,N);
  r_2=corrcoeff(x_2,x_2,N);
end
R_1=r2R(r_1,N);
R_2=r2R(r_2,N);

% Solve the optimization equation
[V,D]=eig(inv(R_2)*R_1);

% Select the eigenvector yielding the maximum object function
J_max=-inf;
for i=1:N*N,
  h=V(:,i);
  J=(h'*R_1*h - h'*R_2*h)^2/((h'*R_1*h)*(h'*R_2*h));
  if J>J_max,
    h_max=h;
    J_max=J;
    i_max=i;
  end
end

J=J_max;
h=reshape(h_max,N,N)';

