--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/m-toolbox/classes/+utils/@math/loglikelihood.m	Wed Nov 23 19:22:13 2011 +0100
@@ -0,0 +1,99 @@
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%
+% 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
+