% BICOHERE computes the bicoherence of two input time-series
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% DESCRIPTION: BICOHERE computes the bicoherence of two input time-series.
%
% CALL:        bs = bicohere(a1,a2,pl)
%
% INPUTS:      aN    - input analysis objects
%              a1,a2 - input analysis objects array
%              pl    - input parameter list
%
% OUTPUTS:     bs   - xyz data analysis object
%
% <a href="matlab:utils.helper.displayMethodInfo('bicohere', 'psd')">Parameters Description</a>
%
% VERSION:    $Id: psd.m,v 1.59 2010/12/17 17:45:08 ingo Exp $
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


function varargout = bicohere(varargin)
  
  % Check if this is a call for parameters
  if utils.helper.isinfocall(varargin{:})
    varargout{1} = getInfo(varargin{3});
    return
  end
  
  import utils.const.*
  utils.helper.msg(msg.PROC3, 'running %s/%s', mfilename('class'), mfilename);
  
  % Collect input variable names
  in_names = cell(size(varargin));
  for ii = 1:nargin,in_names{ii} = inputname(ii);end
  
  % Collect all AOs and plists
  [as, ao_invars] = utils.helper.collect_objects(varargin(:), 'ao', in_names);
  pl              = utils.helper.collect_objects(varargin(:), 'plist', in_names);
  
  if numel(as) ~= 2
    error('bicohere only works with 2 time-series at the moment.');
  end
  
  % Get data
  a  = as(1);
  b  = as(2);
  
  % same fs?
  if a.data.fs ~= b.data.fs
    error('### Two time-series have different sample rates.');
  end
  
  % Same length vectors?
  if a.len ~= b.len
    error('### Two time-series must have same length');
  end
  
  % Combine plists
  pl = combine(pl, getDefaultPlist); 
  usepl = utils.helper.process_spectral_options(pl, 'lin', a.len, a.fs);
  
  win     = find(usepl, 'Win');
  nfft    = find(usepl, 'Nfft');
  olap    = find(usepl, 'Olap')/100;
  Xolap   = round(olap*nfft);
  order   = find(usepl, 'order');
    
  fs = a.data.fs;
  x  = {a.data.y; b.data.y};
    
  [x,M,isreal_x,y,Ly,win,winName,winParam,noverlap,k,L,options] = ...
    ao.welchparse(x,'',win.win, Xolap, nfft, fs);
  
  select = 1:(nfft+1)/2;

  % loop over segments
  LminusOverlap = L-noverlap;
  xStart = 1:LminusOverlap:k*LminusOverlap;
  xEnd   = xStart+L-1;
  
  m = zeros(length(select), length(select));

  for i = 1:k
    if order < 0
      Xseg = x(xStart(i):xEnd(i));
      Yseg = x(xStart(i):xEnd(i));
    else
      [Xseg,coeffs] = ltpda_polyreg(x(xStart(i):xEnd(i)), order);
      [Yseg,coeffs] = ltpda_polyreg(y(xStart(i):xEnd(i)), order);
    end
     
    % window
    xw = Xseg.*win.';
    yw = Yseg.*win.';

    % FFT
    xx2s = fft(xw);    
    xx   = xx2s(select);
    yy2s = fft(yw);
    yy   = yy2s(select);
    
    scalex = abs(xx);
    scaley = abs(yy);
    sc = scalex * scaley';
    m = m + (xx * yy')./sc;
        
  end
    
  m = m./k;
  
  f = psdfreqvec('npts',nfft,'Fs',fs);
  f = f(select);
  do = xyzdata(f, f, m);
  do.setXunits('Hz');
  do.setYunits('Hz');
  do.setZunits(a.yunits*b.yunits);
  
  ma = ao(do);
  ma.setName(sprintf('%s, %s', a.name, b.name));
  ma.addHistory(getInfo('None'), usepl, ao_invars, [a.hist b.hist]);

  
  % Set output
  varargout{1} = ma;
end

%--------------------------------------------------------------------------
% Get Info Object
%--------------------------------------------------------------------------
function ii = getInfo(varargin)
  if nargin == 1 && strcmpi(varargin{1}, 'None')
    sets = {};
    pl   = [];
  else
    sets = {'Default'};
    pl   = getDefaultPlist;
  end
  % Build info object
  ii = minfo(mfilename, 'ao', 'ltpda', utils.const.categories.sigproc, '$Id: psd.m,v 1.59 2010/12/17 17:45:08 ingo Exp $', sets, pl);
end

%--------------------------------------------------------------------------
% Get Default Plist
%--------------------------------------------------------------------------
function plout = getDefaultPlist()
  persistent pl;  
  if exist('pl', 'var')==0 || isempty(pl)
    pl = buildplist();
  end
  plout = pl;  
end

function pl = buildplist()
  
  % General plist for Welch-based, linearly spaced spectral estimators
  pl = plist.WELCH_PLIST;
  
  
end
% END