--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/m-toolbox/classes/@ao/confint.m	Wed Nov 23 19:22:13 2011 +0100
@@ -0,0 +1,427 @@
+% CONFINT Calculates confidence levels and variance for psd, lpsd, cohere, lcohere and curvefit parameters
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%
+% DESCRIPTION: CONFINT Input psd, mscohere (magnitude square coherence) 
+% and return confidence levels and variance for them.
+% Spectra are assumed to be calculated with the WOSA method (Welch's
+% Overlapped Segment Averaging Method)
+%
+% CALL:         out = confint(a,pl)
+%               
+%
+% INPUTS:
+%               a  -  input analysis objects containing power spectral
+%                     densities or magintude squared coherence.
+%               pl  - input parameter list
+%
+% OUTPUTS:         
+%               out - a collection object containing:
+%                 lcl - lower confidence level
+%                 ucl - upper confidence level
+%                 var - expected spectrum variance
+%
+%
+%              If the last input argument is a parameter list (plist).
+%              The following parameters are recognised.
+%
+% <a href="matlab:utils.helper.displayMethodInfo('ao', 'confint')">Parameters Description</a>
+% 
+% 
+% 
+%
+%
+% VERSION:    $Id: confint.m,v 1.20 2011/04/29 13:54:36 luigi Exp $
+%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+function varargout = confint(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, CONFINT can process only one AO per time !!!')
+  end
+  
+  %%% avoid input modification
+  if nargout == 0
+    error('!!! CONFINT 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'));
+  conf = find(pl, 'conf');
+  dof  = find(pl, 'dof');
+  conf = conf/100; % go from percentage to fractional
+  Ntot = find(pl,'DataLength');
+  
+  %%% check that fsdata is input
+  if ~isa(as.data, 'fsdata')
+    error('!!! Non-fsdata input, CONFINT can process only fsdata !!!')
+  end
+
+  
+  % looking to dof
+  if isempty(dof)
+    calcdof = true;
+  else
+    if isa(dof, 'ao')
+      dof = dof.data.y;
+      calcdof = false;
+    else
+      calcdof = false;
+    end
+  end
+    
+  
+  %%% switching over methods
+  switch mtd
+    case 'psd'
+      %%% confidence levels for spectra calculated with psd
+      
+      % calculating dof
+      if calcdof
+        dofs = getdof(as,plist('method',mtd,'DataLength',Ntot));
+        dof = dofs.y;
+      end % if calcdof
+      dof = round(dof);
+      if length(dof)~=1
+        error('!!! CONFINT for ao/psd method, dof must be a single number')
+      end
+      
+      % Calculating Confidence Levels factors
+      alfa = 1 - conf;
+      c = utils.math.Chi2inv([1-alfa/2 alfa/2],dof);
+      c = dof./c;
+      
+      % calculating variance
+      expvar = ((as.data.y).^2).*2./dof;
+      
+      % calculating confidence levels
+      lwb = as.data.y.*c(1);
+      upb = as.data.y.*c(2);
+      
+    case 'lpsd'
+      %%% confidence levels for spectra calculated with lpsd
+      
+      % calculating dof
+      if calcdof
+        dofs = getdof(as,plist('method',mtd,'DataLength',Ntot));
+        dofs = dofs.y;
+
+        % extract number of frequencies bins
+        nf = length(as.x);
+        
+        cl = ones(nf,2);
+        for jj=1:nf
+
+          % Calculating Confidence Levels factors
+          alfa = 1 - conf;
+          c = utils.math.Chi2inv([1-alfa/2 alfa/2],dofs(jj));
+          c = dofs(jj)./c;
+
+          % storing c
+          cl(jj,1) = c(1);
+          cl(jj,2) = c(2);
+        end % for jj=1:nf
+      else % if calcdof
+        if length(dof)~=length(as.x)
+          error('!!! CONFINT for ao/lpsd method, dof must be a vector of the same length of the frequencies vector')
+        end
+        dofs = round(dof);
+        cl = ones(length(as.x),2);
+        for jj = 1:length(as.x)
+          % Calculating Confidence Levels factors
+          alfa = 1 - conf;
+          c = utils.math.Chi2inv([1-alfa/2 alfa/2],dofs(jj));
+          c = dofs(jj)./c;
+
+          % storing c
+          cl(jj,1) = c(1);
+          cl(jj,2) = c(2);
+        end
+      end % if calcdof
+      % willing to work with columns
+      dy = as.data.y;
+      [ii,kk] = size(dy);
+      if ii<kk
+        dy = dy.';
+        rsp = true;
+      else
+        rsp = false;
+      end
+      % calculating variance
+      expvar = ((dy).^2).*2./dofs;
+      
+      % calculating confidence levels
+      lwb = dy.*cl(:,1);
+      upb = dy.*cl(:,2);
+      
+      % reshaping if necessary
+      if rsp
+        expvar = expvar.';
+        lwb = lwb.';
+        upb = upb.';
+      end
+      
+    case 'mscohere'
+      %%% confidence levels for mscohere calculated with ao/cohere
+      
+      % calculating dof
+      if calcdof
+        dofs = getdof(as,plist('method',mtd,'DataLength',Ntot));
+        dof = dofs.y;
+      end % if calcdof
+      dof = round(dof);
+      if length(dof)~=1
+        error('!!! CONFINT for ao/cohere method, dof must be a single number')
+      end
+      
+      % Defining Y variable
+      Y = atanh(sqrt(as.data.y));
+      
+      % Calculating Confidence Levels factor
+      alfa = 1 - conf;
+      c = -sqrt(2).*erfcinv(2*(1-alfa/2))./sqrt(dof);
+      Ylwb = Y - c;
+      Yupb = Y + c;
+      
+      % calculating confidence levels
+      lwb = tanh(Ylwb).^2;
+      upb = tanh(Yupb).^2;
+      
+      % calculating variance
+      expvar = ((1-(as.data.y).^2).^2).*((as.data.y).^2).*4./dof;
+      
+    case 'mslcohere'
+      %%% confidence levels for spectra calculated with lpsd
+      
+      % calculating dof
+      if calcdof
+        dofs = getdof(as,plist('method',mtd,'DataLength',Ntot));
+        dofs = dofs.y;
+
+        % extract number of frequencies bins
+        nf = length(as.x);
+       
+        % willing to work with columns
+        dy = as.data.y;
+        [ii,kk] = size(dy);
+        if ii<kk
+          dy = dy.';
+          rsp = true;
+        else
+          rsp = false;
+        end
+  
+        % Defining Y variable
+        Y = atanh(sqrt(dy));
+
+        cl = ones(nf,2);
+        for jj=1:nf
+
+          % Calculating Confidence Levels factors
+          alfa = 1 - conf;
+          c = -sqrt(2).*erfcinv(2*(1-alfa/2))./sqrt(dofs(jj));
+
+          % storing c and dof
+          cl(jj,1) = Y(jj) - c;
+          cl(jj,2) = Y(jj) + c;
+        end % for jj=1:nf
+        
+      else % if calcdof
+        if length(dof)~=length(as.x)
+          error('!!! CONFINT for ao/lcohere method, dof must be a vector of the same length of the frequencies vector')
+        end
+        dofs = round(dof);
+        
+        % willing to work with columns
+        dy = as.data.y;
+        [ii,kk] = size(dy);
+        if ii<kk
+          dy = dy.';
+          rsp = true;
+        else
+          rsp = false;
+        end
+        
+        % Defining Y variable
+        Y = atanh(sqrt(dy));
+        
+        cl = ones(length(as.x),2);
+        for jj = 1:length(as.x)
+          % Calculating Confidence Levels factors
+          alfa = 1 - conf;
+          c = -sqrt(2).*erfcinv(2*(1-alfa/2))./sqrt(dofs(jj));
+
+          % storing c
+          cl(jj,1) = Y(jj) - c;
+          cl(jj,2) = Y(jj) + c;
+        end
+      end % if calcdof
+      
+      % calculating variance
+      expvar = ((1-(dy).^2).^2).*((dy).^2).*4./dofs;
+      
+      % get not well defined coherence estimations
+      idd = dofs<=2;
+      
+      % calculating confidence levels
+      lwb = tanh(cl(:,1));
+      % remove negative elements
+      idx = lwb < 0;
+      lwb(idx) = 0;
+      % set lower bound to zero in points where coharence is not well defined 
+      lwb(idd) = 0;
+      upb = tanh(cl(:,2));
+      % set upper bound to one in points where coharence is not well
+      % defined
+      upb(idd) = 1;
+      lwb = lwb.^2;
+      upb = upb.^2;
+      
+      % reshaping if necessary
+      if rsp
+        expvar = expvar.';
+        lwb = lwb.';
+        upb = upb.';
+      end
+      
+      
+      
+  end %switch mtd
+  
+  % Output data
+ 
+    
+  % defining units
+  inputunit = get(as.data,'yunits');
+  varunit = unit(inputunit.^2);
+  varunit.simplify;
+  levunit = inputunit;
+  levunit.simplify;
+  
+  
+  % variance
+  plvar = plist('xvals', as.data.x, 'yvals', expvar, 'type', 'fsdata');
+  ovar = ao(plvar);
+  
+  ovar.setFs(as.data.fs);
+  ovar.setT0(as.data.t0);
+  ovar.data.setEnbw(as.data.enbw);
+  ovar.data.setNavs(as.data.navs);
+  ovar.setXunits(as.data.xunits);
+  ovar.setYunits(varunit);
+  % Set output AO name
+  ovar.name = sprintf('var(%s)', ao_invars{:});
+  
+  % lower confidence level
+  pllwb = plist('xvals', as.data.x, 'yvals', lwb, 'type', 'fsdata');
+  olwb = ao(pllwb);
+
+  olwb.setFs(as.data.fs);
+  olwb.setT0(as.data.t0);
+  olwb.data.setEnbw(as.data.enbw);
+  olwb.data.setNavs(as.data.navs);
+  olwb.setXunits(copy(as.data.xunits,1));
+  olwb.setYunits(levunit);
+  % Set output AO name
+  clev = [num2str(conf*100) '%'];
+  olwb.name = sprintf('%s_low_conf_level(%s)', clev, ao_invars{:});
+  
+  % upper confidence level
+  plupb = plist('xvals', as.data.x, 'yvals', upb, 'type', 'fsdata');
+  oupb = ao(plupb);
+
+  oupb.setFs(as.data.fs);
+  oupb.setT0(as.data.t0);
+  oupb.data.setEnbw(as.data.enbw);
+  oupb.data.setNavs(as.data.navs);
+  oupb.setXunits(copy(as.data.xunits,1));
+  oupb.setYunits(levunit);
+  % Set output AO name
+  oupb.name = sprintf('%s_up_conf_level(%s)', clev, ao_invars{:});
+  
+  
+  outobj = collection(olwb,oupb,ovar);
+  outobj.setName(sprintf('%s conf levels for %s', clev, ao_invars{:}));
+  outobj.addHistory(getInfo('None'), pl, [ao_invars(:)], [as.hist]);
+  
+  varargout{1} = outobj;
+    
+  
+end
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%                               Local Functions                               %
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+
+%----- LTPDA FUNCTIONS ----------------------------------------------------
+%--------------------------------------------------------------------------
+% 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: confint.m,v 1.20 2011/04/29 13:54:36 luigi Exp $', sets, pl);
+  ii.setModifier(false);
+  ii.setOutmin(2);
+  ii.setOutmax(3);
+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()
+
+ 
+  pl = ao.getInfo('getdof', 'Default').plists;
+  
+  % Conf
+  p = param({'Conf', ['Required percentage confidence level.<br>' ...
+    'It is a number between 0 and 100.']}, ...
+    {1, {95}, paramValue.OPTIONAL});
+  pl.pset(p);
+  
+  % DOF
+  p = param({'dof', ['Degrees of freedom of the estimator. If it is<br>'...
+    'left empty they are calculated.']}, paramValue.EMPTY_DOUBLE);
+  pl.pset(p);
+end
+% END
+