%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Compute log-likelihood
%
% M Nofrarias 15-06-09
%
% $Id: loglikelihood.m,v 1.2 2011/03/24 19:59:50 ingo Exp $
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function loglk = loglikelihood(xn,in,out,noise,model,params)
  
  % parameters
  fs = noise(1).fs;
  N = length(noise(1).y);
  
  if numel(in) == 1 % 1 input / 1 output case
    f = in(1).x;
    % evaluate models
    eval = model(1).setParams(params,xn);
    h11 = double(eval);
    
    % spectra to variance
    C11 = (N*fs/2)*noise(1).y;
    
    % compute elements of inverse cross-spectrum matrix
    InvS11 = 1./C11;
    
    % compute log-likelihood terms first, all at once does not cancel the
    % imag part when multiplying x.*conj(x)
    v1v1 = conj(out(1).y - h11.*in(1).y).*(out(1).y - h11.*in(1).y);
    
    loglk = sum(InvS11.*v1v1);
    
  elseif numel(in) == 2
    
    if(numel(out) == 1) % 2 input / 1 output case
      f = in(1).x;
      
      % evaluate models
      h11 = model(1).setParams(params,xn).double;
      h12 = model(2).setParams(params,xn).double;
      
      % spectra to variance
      C11 = (N*fs/2)*noise(1).y;
      
      % compute elements of inverse cross-spectrum matrix
      InvS11 = 1./C11;
      
      % compute log-likelihood terms first, all at once does not cancel the
      % imag part when multiplying x.*conj(x)
      in1 = in(1).data.getY;
      f = in(1).data.getX;
      in2 = in(2).data.getY;
      out1 = out(1).data.getY;
      
      v1v1 = conj(out1 - h11.*in1 - h12.*in2).*(out1 - h11.*in1 - h12.*in2);
      
      loglk = sum(InvS11.*v1v1);
      
    elseif(numel(out) == 2) % 2 input / 2 output case
      f = in(1).x;      
      
      % evaluate models
      h11 = model(1).setParams(params,xn).double;
      h12 = model(2).setParams(params,xn).double;
      h21 = model(3).setParams(params,xn).double;
      h22 = model(4).setParams(params,xn).double;
      
      % spectra to variance
      C11 = (N*fs/2)*noise(1).y;
      C22 = (N*fs/2)*noise(2).y;
      C12 = (N*fs/2)*noise(3).y;
      C21 = (N*fs/2)*noise(4).y;
      
      % compute elements of inverse cross-spectrum matrix
      InvS11 = (C22./(C11.*C22 - C12.*C21));
      InvS22 = (C11./(C11.*C22 - C12.*C21));
      InvS12 = (C21./(C11.*C22 - C12.*C21));
      InvS21 = (C12./(C11.*C22 - C12.*C21));
      
      % compute log-likelihood terms first, all at once does not cancel the
      % imag part when multiplying x.*conj(x)
      in1 = in(1).data.getY;
      in2 = in(2).data.getY;
      out1 = out(1).data.getY;
      out2 = out(2).data.getY;
      
      v1v1 = conj(out1 - h11.*in1 - h12.*in2).*(out1 - h11.*in1 - h12.*in2);
      v2v2 = conj(out2 - h21.*in1 - h22.*in2).*(out2 - h21.*in1 - h22.*in2);
      v1v2 = conj(out1 - h11.*in1 - h12.*in2).*(out2 - h21.*in1 - h22.*in2);
      v2v1 = conj(out2 - h21.*in1 - h22.*in2).*(out1 - h11.*in1 - h12.*in2);
      
      loglk = sum(InvS11.*v1v1 + InvS22.*v2v2 - InvS12.*v1v2 - InvS21.*v2v1);
    end
  end
  
end