% LINEDETECT find spectral lines in the ao/fsdata objects.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% DESCRIPTION: LINEDETECT find spectral lines in the ao/fsdata objects.
%
% CALL:        b = linedetect(a, pl)
%
% <a href="matlab:utils.helper.displayMethodInfo('ao', 'linedetect')">Parameters Description</a>
%
% VERSION:    $Id: linedetect.m,v 1.15 2011/04/08 08:56:16 hewitson Exp $
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function varargout = linedetect(varargin)

  % Check if this is a call for parameters
  if utils.helper.isinfocall(varargin{:})
    varargout{1} = getInfo(varargin{3});
    return
  end

  if nargout == 0
    error('### cat cannot be used as a modifier. Please give an output variable.');
  end

  % Collect input variable names
  in_names = cell(size(varargin));
  for ii = 1:nargin,in_names{ii} = inputname(ii);end

  % Collect all AOs
  [as, ao_invars] = utils.helper.collect_objects(varargin(:), 'ao', in_names);
  [pli, pl_invars] = utils.helper.collect_objects(varargin(:), 'plist', in_names);

  % Decide on a deep copy or a modify
  bs = copy(as, nargout);

  Na = numel(bs);
  if isempty(bs)
    error('### Please input at least one AO.');
  end

  % Combine plists
  if ~isempty(pli)
    pl = parse(pli, getDefaultPlist());
  else
    pl = getDefaultPlist();
  end

  % Get parameters from plist
  N       = find(pl, 'N');
  fsearch = find(pl, 'fsearch');
  thresh  = find(pl, 'thresh');

  % Loop over input AOs
  for jj = 1:Na
    if isa(bs(jj).data, 'fsdata')

      % Make noise-floor estimate
      nf = smoother(bs(jj), pl);

      % Make ratio
      r = bs(jj)./nf;

      % find lines
      lines = findLines(bs(jj).data.getY, r.data.getX, r.data.getY, thresh, N, fsearch);

      f = [lines(:).f];
      y = [lines(:).a];

      % Keep the data shpare of the input AO
      if size(bs(jj).data.y, 2) == 1
        f = f.';
        y = y.';
      end

      % Make output data: copy the fsdata object so to inherit all the feautures
      fs = copy(bs(jj).data, 1);
      
      % Make output data: set the values
      fs.setX(f);
      fs.setY(y);

    else
      error('### I can only find lines in frequency-series AOs.');
    end

    %------- Make output AO

    % make output analysis object
    bs(jj).data = fs;

    bs(jj).name = sprintf('lines(%s)', ao_invars{1});
    bs(jj).addHistory(getInfo('None'), pl, ao_invars(jj), bs(jj).hist);
  end

  % Set output
  if nargout == numel(bs)
    % List of outputs
    for ii = 1:numel(bs)
      varargout{ii} = bs(ii);
    end
  else
    % Single output
    varargout{1} = bs;
  end
end

%--------------------------------------------------------------------------
% find spectral lines
function lines = findLines(ay, f, nay, thresh, N, fsearch)

  % look for spectral lines
  l       = 0;
  pmax    = 0;
  pmaxi   = 0;
  line    = [];
  idx     = find( f>=fsearch(1) & f<=fsearch(2) );
  for jj = 1:length(idx)
    v = nay(idx(jj));
    if v > thresh
      if v > pmax
        pmax  = v;
        pmaxi = idx(jj);
      end
    else
      if pmax > 0
        % find index when we have pmax
        fidx = pmaxi; %(find(nay(1:idx(jj))==pmax));
        l = l+1;
        line(l).idx = fidx;
        line(l).f   = f(fidx);
        line(l).a   = ay(fidx);
      end
      pmax = 0;
    end
  end

  % Select largest peaks
  lines = [];
  if ~isempty(line)
    [bl, lidx] = sort([line.a], 'descend');
    lidxs = lidx(1:min([N length(lidx)]));
    lines = line(lidxs);
    disp(sprintf('   + found %d lines.', length([lines.f])));
  else
    disp('   + found 0 lines.');
  end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                               Local Functions                               %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION:    getInfo
%
% DESCRIPTION: Get Info Object
%
% HISTORY:     11-07-07 M Hewitson
%                Creation.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

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: linedetect.m,v 1.15 2011/04/08 08:56:16 hewitson Exp $', sets, pl);
  ii.setModifier(false);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION:    getDefaultPlist
%
% DESCRIPTION: Get Default Plist
%
% HISTORY:     11-07-07 M Hewitson
%                Creation.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function plout = getDefaultPlist()
  persistent pl;  
  if exist('pl', 'var')==0 || isempty(pl)
    pl = buildplist();
  end
  plout = pl;  
end

function pl = buildplist()

  pl = plist();
  
  % N
  p = param({'N', 'The maximum number of lines to return.'}, {1, {10}, paramValue.OPTIONAL});
  pl.append(p);
  
  % fsearch
  p = param({'fsearch', 'The frequency search interval.'}, {1, {[0 1e10]}, paramValue.OPTIONAL});
  pl.append(p);
  
  % thresh
  p = param({'thresh', 'A threshold to test normalised amplitude against. (A sort-of SNR threshold.)'}, {1, {2}, paramValue.OPTIONAL});
  pl.append(p);
  
  % BW
  p = param({'bw', ['The bandwidth of the running median filter used to<br>'...
                    'estimate the noise-floor.']}, {1, {20}, paramValue.OPTIONAL});
  pl.append(p);                
                  
  % HC
  p = param({'hc', 'The cutoff used to reject outliers (0-1).'}, {1, {0.8}, paramValue.OPTIONAL});
  pl.append(p);
  
  
end

% PARAMETERS:  N        - max number of lines to return  [default: 10]
%              fsearch  - freqeuncy search interval  [default: all]
%              thresh   - a threshold to test normalised amplitude spectrum against
%                         [default: 2]
%              bw       - bandwidth over which to compute median [default: 20 samples]
%              hc       - percent of outliers to exclude from median estimation (0-1)
%                         [default: 0.8]