% GETDOF Calculates degrees of freedom for psd, lpsd, cohere and lcohere
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% DESCRIPTION: GETDOF Input psd or mscohere (magnitude square coherence)
% estimated with the WOSA (Welch's Overlapped Segment Averaging Method) and
% return degrees of freedom of the estimator.
%
% CALL:         dof = getdof(a,pl)
%
% INPUTS:
%               a  -  input analysis objects containing power spectral
%                     densities or magnitude squared coherence.
%               pl  - input parameter list
%
% OUTPUTS:
%               dof - cdata AO with degrees of freedom for the
%                     corresponding estimator. If the estimators are lpsd
%                     or lcohere then dof number of elements is the same of
%                     the spectral estimator
%
%
%              If the last input argument is a parameter list (plist).
%              The following parameters are recognised.
%
% 
%
% <a href="matlab:utils.helper.displayMethodInfo('ao', 'getdof')">Parameters Description</a>
%
% VERSION:    $Id: getdof.m,v 1.19 2011/05/23 20:36:47 mauro Exp $
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function varargout = getdof(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
  [as, ao_invars] = utils.helper.collect_objects(varargin(:), 'ao', in_names);
  pl              = utils.helper.collect_objects(varargin(:), 'plist', in_names);
  
  %%% avoid multiple AO at input
  if numel(as) > 1
    error('!!! Too many input AOs, GETDOF can process only one AO per time !!!')
  end
  
  %%% check that fsdata is input
  if ~isa(as.data, 'fsdata')
    error('!!! Non-fsdata input, GETDOF can process only fsdata !!!')
  end
  
  %%% avoid input modification
  if nargout == 0
    error('!!! GETDOF cannot be used as a modifier. Please give an output variable !!!');
  end
  
  %%% Parse plists
  pl = parse(pl, getDefaultPlist());
  
  %%% Find parameters
  mtd  = lower(find(pl, 'method'));
  if ~ischar(mtd)
    error('!!! Method must be a string !!!')
  end
  mtd  = lower(mtd);
  
  Ntot = find(pl,'DataLength');
  
  %%% switching over methods
  switch mtd
    case 'psd'
      % get hist
      hst = as.hist;
      % get nodes
      [n,a, nodes] = getNodes(hst);
      % get plists from nodes
      pls = [nodes(:).pl];
      if numel(pls) > 1
        plss = pls(1);
        for ii = 2:numel(pls)-1
          plss = parse(plss, pls(ii));
        end
      end
      % get number of averages
      navs = as.data.navs;
      % get window object
      w = find(plss, 'WIN');
      % percentage of overlap
      olap = find(plss, 'OLAP')./100;
      % number of bins in each fft
      nfft = find(plss, 'NFFT');
      psll = find(plss, 'psll');
      
      % Normalize window data in order to be square integrable to 1
      if ischar(w)
        switch lower(w)
          case 'kaiser'
            Win = specwin(w, nfft, psll);
          otherwise
            Win = specwin(w, nfft);
        end
      else
        Win = w;
      end
      win = Win.win ./ sqrt(Win.ws2);
      
      % Calculates total number of data in the original time-series
      if isempty(Ntot)
        Ntot = ceil(navs*(nfft-olap*nfft)+olap*nfft);
      end
      
      if navs == 1
        dofs = round(2*navs);
      else
        [R,n] = utils.math.overlapCorr(win,Ntot,navs);
        dof = 2*navs/(2*R*navs+1);
        dofs = round(dof);
      end
      
    case 'lpsd'
      % get hist
      hst = as.hist;
      % get nodes
      [n,a, nodes] = getNodes(hst);
      % get plists from nodes
      pls = [nodes(:).pl];
      if numel(pls) > 1
        plss = pls(1);
        for ii = 2:numel(pls)-1
          plss = parse(plss, pls(ii));
        end
      end
      % get window used
      uwin = find(plss, 'WIN');
      
      % extract number of frequencies bins
      nf = length(as.x);
      
      % dft length for each bin
      if ~isempty(as.procinfo)
        L = as.procinfo.find('L');
      else
        error('### The AO doesn''t have any procinfo with the key ''L''');
      end
      
      % set original data length as the length of the first window
      if isempty(Ntot)
        nx = L(1);
      else
        nx = Ntot;
      end
      
      % windows overlap
      olap = find(plss, 'OLAP')./100;
      psll = find(plss, 'psll');
      
      dofs = ones(nf,1);
      for jj = 1:nf
        l = L(jj);
        % compute window
        if ischar(uwin)
          switch lower(uwin)
            case 'kaiser'
              w = specwin(uwin, l, psll);
            otherwise
              w = specwin(uwin, l);
          end
        else
          w = uwin;
        end
        % Normalize window data in order to be square integrable to 1
        owin = w.win;
        owin = owin./sqrt(w.ws2);
        
        % Compute the number of averages we want here
        segLen = l;
        nData  = nx;
        ovfact = 1 / (1 - olap);
        davg   = (((nData - segLen)) * ovfact) / segLen + 1;
        navg   = round(davg);

        if navg == 1
          dof = 2*navg;
        else
          [R,n] = utils.math.overlapCorr(owin,nx,navg);
          dof = 2*navg/(2*R*navg+1);
        end
        
        % storing c and dof
        dofs(jj) = dof;
        
      end % for jj=1:nf
      
    case 'mscohere'
      % get hist
      hst = as.hist;
      % get nodes
      [n,a, nodes] = getNodes(hst);
      % get plists from nodes
      pls = [nodes(:).pl];
      if numel(pls) > 1
        plss = pls(1);
        for ii = 2:numel(pls)-1
          plss = parse(plss, pls(ii));
        end
      end
      % get number of averages
      navs = as.data.navs;
      % get window object
      w = find(plss,'WIN');
      % percentage of overlap
      olap = find(plss,'OLAP')./100;
      % number of bins in each fft
      nfft = find(plss,'NFFT');
      psll = find(plss, 'psll');
      
      % Normalize window data in order to be square integrable to 1
      if ischar(w)
        switch lower(w)
          case 'kaiser'
            Win = specwin(w, nfft, psll);
          otherwise
            Win = specwin(w, nfft);
        end
      else
        Win = w;
      end
      win = Win.win ./ sqrt(Win.ws2);
      
      % Calculates total number of data in the original time-series
      if isempty(Ntot)
        Ntot = ceil(navs*(nfft-olap*nfft)+olap*nfft);
      end
      
      if navs == 1
        dofs = round(2*navs);
      else
        [R,n] = utils.math.overlapCorr(win,Ntot,navs);
        dof = 2*navs/(2*R*navs+1);
        dofs = round(dof);
      end
      
    case 'mslcohere'
      % get hist
      hst = as.hist;
      % get nodes
      [n,a, nodes] = getNodes(hst);
      % get plists from nodes
      pls = [nodes(:).pl];
      if numel(pls) > 1
        plss = pls(1);
        for ii = 2:numel(pls)-1
          plss = parse(plss, pls(ii));
        end
      end
      % get window used
      uwin = find(plss,'WIN');
      
      % extract number of frequencies bins
      nf = length(as.x);
      
      % dft length for each bin
      if ~isempty(as.procinfo)
        L = as.procinfo.find('L');
      else
        error('### The AO doesn''t have any procinfo with the key ''L''');
      end
      
      % set original data length as the length of the first window
      if isempty(Ntot)
        nx = L(1);
      else
        nx = Ntot;
      end
      
      % windows overlap
      olap = find(plss, 'OLAP')./100;
      psll = find(plss, 'psll');
      
      dofs = ones(nf, 1);
      for jj = 1:nf
        l = L(jj);
        % compute window
        if ischar(uwin)
          switch lower(uwin)
            case 'kaiser'
              w = specwin(uwin, l, psll);
            otherwise
              w = specwin(uwin, l);
          end
        else
          w = uwin;
        end
        % Normalize window data in order to be square integrable to 1
        owin = w.win ./ sqrt(w.ws2);
        
        % Compute the number of averages we want here
        segLen = l;
        nData  = nx;
        ovfact = 1 / (1 - olap);
        davg   = (((nData - segLen)) * ovfact) / segLen + 1;
        navg   = round(davg);
        
        if navg == 1
          dof = 2*navg;
        else
          [R,n] = utils.math.overlapCorr(owin,nx,navg);
          dof = 2*navg/(2*R*navg+1);
        end
        
        % storing c and dof
        dofs(jj) = dof;
        
      end % for jj=1:nf
      
  end %switch mtd
  
  % Output data
  
  % dof
  ddof = cdata();
  ddof.setY(dofs);
  odof = ao(ddof);
  % Set output AO name
  odof.name = sprintf('dof(%s)', ao_invars{:});
  % Add history
  odof.addHistory(getInfo('None'), pl, [ao_invars(:)], [as.hist]);
  
  % output
  if nargout == 1
    varargout{1} = odof;
  else
    error('!!! getdof can have only one output')
  end
  
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: getdof.m,v 1.19 2011/05/23 20:36:47 mauro Exp $', sets, pl);
  ii.setModifier(false);
  ii.setOutmin(1);
  ii.setOutmax(1);
end

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

function pl = buildplist()

  pl = plist();
  
  p = param({'method', ['Set the desired method. Supported values are<ul>'...
    '<li>''psd'' power spectrum calculated with ao/psd, whatever the scale</li>'...
    '<li>''lpsd'' power spectrum calculated with ao/lpsd, whatever the scale</li>'...
    '<li>''mscohere'' magnitude square coherence spectrum calculated with ao/cohere</li>'...
    '<li>''mslcohere'' magnitude square coherence spectrum calculated with ao/lcohere</li>']}, ...
    {1, {'psd', 'lpsd', 'mscohere', 'mslcohere'}, paramValue.OPTIONAL});
  pl.append(p);
  
  p = param({'DataLength',['Data length of the time series.'...
    'It is better to input for more stable calculation.'...
    'Leave it empty if you do not know its value.']},...
      paramValue.EMPTY_DOUBLE);
  pl.append(p);
                       
                       
end



% END