% MLTFE compute log-frequency space TF
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% DESCRIPTION: MLTFE compute log-frequency space TF
%
% CALL:        Txy = mltfe(X,f,r,m,L,fs,win,order,olap,Lmin,method,variance)
%
% VERSION:     $Id: mltfe.m,v 1.25 2009/11/27 15:18:45 miquel Exp $
%
% HISTORY:     19-05-2007 M Hewitson
%                 Creation
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function varargout = mltfe(varargin)
  
  import utils.const.*
  
  % Get inputs
  X      = varargin{1};
  f      = varargin{2};
  r      = varargin{3};
  m      = varargin{4};
  L      = varargin{5};
  K      = varargin{6};
  fs     = varargin{7};
  win    = varargin{8};
  order  = varargin{9};
  olap   = varargin{10};
  Lmin   = varargin{11};
  method = varargin{12};
  
  % --- Prepare some variables
  si         = size(X);
  nc         = si(1);
  nf         = length(f);
  Txy        = zeros(nf,1);
  dev        = zeros(nf,1);
  disp_each  = round(nf/100)*10;
  
  winType = win.type;
  winPsll = win.psll;
  
  % ----- Loop over Frequency
  for fi=1:nf
    [Txy(fi) dev(fi)]= computeTF(fs, L(fi), K(fi), m, winType, winPsll, X, olap, order, nc, f(fi), fi, nf, disp_each, method);
  end
  
  % Set output
  varargout{1} = Txy;
  varargout{2} = dev;
end

%--------------------------------------------------------------------------
% Function to run over channels
function  [Txy,dev]= computeTF(fs, l, K, m, winType, winPsll, X, olap, order, nc, ffi, fi, nf, disp_each, method)
  
  switch lower(winType)
    case 'kaiser'
      lwin = specwin(winType, l, winPsll);
    otherwise
      lwin = specwin(winType, l);
  end
  
  % Compute DFT coefficients
  twopi = 2.0*pi;
  p     = 1i * twopi * m(fi)/l.*[0:l-1];
  C     = lwin.win .* exp(p);
  
  if mod(fi,disp_each) == 0 || fi == 1 || fi == nf
    utils.helper.msg(utils.const.msg.PROC1, 'computing frequency %04d of %04d: %f Hz', fi, nf, ffi);
  end
  % Loop over input channels
  [Txy,dev] = in2out(l, X, K, olap, order, nc, method, C, lwin, fs);
  
end

%--------------------------------------------------------------------------
% Compute 1 input to multiple outputs
function [Txy,dev]= in2out(l, X, K, olap, order, nc, method, C, lwin, fs)
  
  % if no errors are required the function returns zero but errors are not
  % stored in the final ao
  dev = 0;
  
   switch lower(method)
    case 'tfe'
      % Core cross-DFT part in C-mex file
      % We need cross-spectrum and Power spectrum
      [XY, XX, YY, M2, nsegs] = ltpda_dft(X(1,:), X(2,:), l, C, olap, order);
      Txy  = conj(XY)/(XX);
      if nsegs == 1
        dev = Inf;
      else
        dev = sqrt((nsegs/(nsegs-1)^2)*(YY./XX).*(1 - (abs(XY).^2)./(XX.*YY)));
%       dP =  sqrt((k/(k-1)^2)*(Pyy./Pxx).*(1 - (abs(Pxy).^2)./(Pxx.*Pyy)));
      end
      
    case 'cpsd'
      % Core cross-DFT part in C-mex file
      [XY, XX, YY, M2, nsegs] = ltpda_dft(X(1,:), X(2,:), l, C, olap, order);
      S1      = lwin.ws;
      S2      = lwin.ws2;
      Txy  = 2.0*XY/fs/S2;
      Var = 4.0*M2/fs^2/S2^2/nsegs;
       if nsegs == 1
        dev = Inf;
      else
        dev = sqrt(Var);
       end
      
    case 'mscohere'
      % Core cross-DFT part in C-mex file
      % We need cross-spectrum and Power spectrum
      [XY, XX, YY, M2, nsegs] = ltpda_dft(X(1,:), X(2,:), l, C, olap, order);
      Txy  = (abs(XY).^2)./(XX.*YY); % Magnitude-squared coherence
      if nsegs == 1
        dev = Inf;
      else
        dev = sqrt((2*Txy/nsegs).*(1-Txy).^2);
      end
      
    case 'cohere'
      % Core cross-DFT part in C-mex file
      % We need cross-spectrum and Power spectrum
      [XY, XX, YY, M2, nsegs] = ltpda_dft(X(1,:), X(2,:), l, C, olap, order);
      Txy  = XY./sqrt(XX.*YY);  % Complex coherence
      if nsegs == 1
        dev = Inf;
      else
        dev = sqrt((2*abs(Txy)/nsegs).*(1-abs(Txy)).^2);
      end
  end
  
end