--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/m-toolbox/classes/@ao/psdconf.m	Wed Nov 23 19:22:13 2011 +0100
@@ -0,0 +1,861 @@
+% PSDCONF Calculates confidence levels and variance for psd
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%
+% DESCRIPTION: PSDCONF Input a spectrum estimated with psd or lpsd (Welch's
+% Overlapped Segment Averaging Method) and calculates confidence levels and
+% variance for the spectrum.
+%
+% CALL:         [lcl,ucl] = psdconf(a,pl)
+%               [lcl,ucl,var] = psdconf(a,pl)
+%
+% INPUTS:
+%               a  -  input analysis objects containing power spectral
+%                     densities calculated with psd or lpsd.
+%               pl  - input parameter list
+%
+% OUTPUTS:                
+%               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.
+%
+% 
+% NOTE1: PSDCONF checks the navs field to distinguish between psd and lpsd
+% power spectra. If a.data.navs is NaN then it assumes to dealing with lpsd
+% power spectrum at input.
+%
+% NOTE2: Copied directly from MATLAB chi2conf function (Signal Processing
+% Toolbox) and extended to do degrees of fredom calculation, variance
+% calculation, to input AOs and to input plist.
+%   Copyright 2007 The MathWorks, Inc.
+%   Revision: 1.6.4.4   Date: 2008/05/31 23:27:28 
+%
+%
+% <a href="matlab:utils.helper.displayMethodInfo('ao', 'psdconf')">Parameters Description</a>
+%
+% VERSION:    $Id: psdconf.m,v 1.14 2011/05/18 07:53:57 mauro Exp $
+%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+function varargout = psdconf(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, PSDCONF can process only one AO per time !!!')
+  end
+  
+  %%% check that fsdata is input
+  if ~isa(as.data, 'fsdata')
+    error('!!! Non-fsdata input, PSDCONF can process only fsdata !!!')
+  end
+  
+  %%% avoid input modification
+  if nargout == 0
+    error('!!! PSDCONF cannot be used as a modifier. Please give an output variable !!!');
+  end
+  
+  %%% Parse plists
+  pl = parse(pl, getDefaultPlist());
+  
+  %%% Find parameters
+  conf    = find(pl, 'conf');
+  conf = conf/100; % go from percentage to fractional
+  
+  %%% getting metho
+  if  isnan(as.data.navs)
+    mtd = 'lpsd';
+  else
+    mtd = 'psd';
+  end
+  
+  %%% switching over methods
+  switch mtd
+    case 'psd'
+      %%% confidence levels for spectra calculated with psd
+      % get number of averages
+      navs = as.data.navs;
+      % get window object
+      w = as.hist.plistUsed.find('WIN');
+      % percentage of overlap
+      olap = as.hist.plistUsed.find('OLAP')./100;
+      % number of bins in each fft
+      nfft = as.hist.plistUsed.find('NFFT');
+      
+      % Normalize window data in order to be square integrable to 1
+      win = w.win ./ sqrt(w.ws2);
+      
+      % Calculates total number of data in the original time-series
+      Ntot = ceil(navs*(nfft-olap*nfft)+olap*nfft);
+      
+      % defining the shift factor
+      n = floor((Ntot-nfft)./(navs-1));
+      
+      % Calculating correction factor
+      sm = 0;
+      % Padding to zero the window
+      % willing to work with columns
+      [ii,kk] = size(win);
+      if ii<kk
+        win = win.';
+      end
+      win = [win; zeros(navs*n,1)];
+      for mm = 1:navs-1
+        pc = 0;
+        for tt = 1:nfft
+          pc = pc + win(tt)*win(tt + mm*n);
+        end
+        sm = sm + (1 - mm/navs)*(abs(pc)^2);
+      end
+      % Calculating degrees of freedom
+      dof = (2*navs)/(1+2*sm);
+      dof = round(dof);
+      
+      % Calculating Confidence Levels factors
+      alfa = 1 - conf;
+      c = 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
+      % get window used
+      uwin = as.hist.plistUsed.find('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
+      nx = L(1);
+      
+      % windows overlap
+      olap = as.hist.plistUsed.find('OLAP');
+      
+      dofs = ones(nf,1);
+      cl = ones(nf,2);
+      for jj=1:nf
+        l = L(jj);
+        % compute window
+        switch uwin.type
+          case 'Kaiser'
+            w = specwin(uwin.type, l, uwin.psll);
+          otherwise
+            w = specwin(uwin.type, l);
+        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);
+        
+        % Compute steps between segments
+        if navg == 1
+          shift = 1;
+        else
+          shift = (nData - segLen) / (navg - 1);
+        end
+        if shift < 1 || isnan(shift)
+          shift = 1;
+        end
+        n = floor(shift);
+        
+        % Calculating correction factor
+        sm = 0;
+        % Padding to zero the window
+        % willing to work with columns
+        [ii,kk] = size(owin);
+        if ii<kk
+          owin = owin.';
+        end
+        win = [owin; zeros(navg*n,1)];
+        for mm = 1:navg-1
+          pc = 0;
+          for tt = 1:segLen
+            pc = pc + win(tt)*win(tt + mm*n);
+          end
+          sm = sm + (1 - mm/navg)*(abs(pc)^2);
+        end
+        % Calculating degrees of freedom
+        dof = (2*navg)/(1+2*sm);
+        dof = round(dof);
+
+        % Calculating Confidence Levels factors
+        alfa = 1 - conf;
+        c = chi2inv([1-alfa/2 alfa/2],dof);
+        c = dof./c;
+        
+        % storing c and dof
+        dofs(jj) = dof;
+        cl(jj,1) = c(1);
+        cl(jj,2) = c(2);
+      end % for jj=1:nf
+      
+      % 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
+      
+  end %switch mtd
+  
+  % Output data
+  
+  % defining units
+  inputunit = get(as.data,'yunits');
+  varunit = unit(inputunit.^2);
+  varunit.simplify;
+  levunit = inputunit;
+  levunit.simplify;
+  
+  
+  % variance
+  dvar = fsdata();
+  dvar.setFs(as.data.fs);
+  dvar.setT0(as.data.t0);
+  dvar.setEnbw(as.data.enbw);
+  dvar.setNavs(as.data.navs);
+  dvar.setXunits(copy(as.data.xunits,1));
+  dvar.setYunits(varunit);
+  dvar.setX(as.data.x);
+  dvar.setY(expvar);
+  ovar = ao(dvar);
+  % Set output AO name
+  ovar.name = sprintf('var(%s)', ao_invars{:});
+  % Add history
+  ovar.addHistory(getInfo('None'), pl, [ao_invars(:)], [as.hist]);
+  
+  % lower confidence level
+  dlwb = fsdata();
+  dlwb.setFs(as.data.fs);
+  dlwb.setT0(as.data.t0);
+  dlwb.setEnbw(as.data.enbw);
+  dlwb.setNavs(as.data.navs);
+  dlwb.setXunits(copy(as.data.xunits,1));
+  dlwb.setYunits(levunit);
+  dlwb.setX(as.data.x);
+  dlwb.setY(lwb);
+  olwb = ao(dlwb);
+  % Set output AO name
+  clev = [num2str(conf*100) '%'];
+  olwb.name = sprintf('%s_low_conf_level(%s)', clev, ao_invars{:});
+  % Add history
+  olwb.addHistory(getInfo('None'), pl, [ao_invars(:)], [as.hist]);
+  
+  % upper confidence level
+  dupb = fsdata();
+  dupb.setFs(as.data.fs);
+  dupb.setT0(as.data.t0);
+  dupb.setEnbw(as.data.enbw);
+  dupb.setNavs(as.data.navs);
+  dupb.setXunits(copy(as.data.xunits,1));
+  dupb.setYunits(levunit);
+  dupb.setX(as.data.x);
+  dupb.setY(upb);
+  oupb = ao(dupb);
+  % Set output AO name
+  oupb.name = sprintf('%s_up_conf_level(%s)', clev, ao_invars{:});
+  % Add history
+  oupb.addHistory(getInfo('None'), pl, [ao_invars(:)], [as.hist]);
+  
+  % output
+  if nargout == 2
+    varargout{1} = olwb; % lower conf level
+    varargout{2} = oupb; % upper conf level
+  elseif nargout == 3
+    varargout{1} = olwb; % lower conf level
+    varargout{2} = oupb; % upper conf level
+    varargout{3} = ovar; % expected variance
+  end
+  
+end
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%                               Local Functions                               %
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+%----- Matlab Functions ---------------------------------------------------
+
+%--------------------------------------------------------------------------
+function x = chi2inv(p,v)
+  %CHI2INV Inverse of the chi-square cumulative distribution function (cdf).
+  %   X = CHI2INV(P,V)  returns the inverse of the chi-square cdf with V
+  %   degrees of freedom at the values in P. The chi-square cdf with V
+  %   degrees of freedom, is the gamma cdf with parameters V/2 and 2.
+  %
+  %   The size of X is the common size of P and V. A scalar input
+  %   functions as a constant matrix of the same size as the other input.
+  
+  %   References:
+  %      [1]  M. Abramowitz and I. A. Stegun, "Handbook of Mathematical
+  %      Functions", Government Printing Office, 1964, 26.4.
+  %      [2] E. Kreyszig, "Introductory Mathematical Statistics",
+  %      John Wiley, 1970, section 10.2 (page 144)
+  
+  
+  if nargin < 2,
+    error(generatemsgid('Nargchk'),'Requires two input arguments.');
+  end
+  
+  [errorcode p v] = distchck(2,p,v);
+  
+  if errorcode > 0
+    error(generatemsgid('InvalidDimensions'),'Requires non-scalar arguments to match in size.');
+  end
+  
+  % Call the gamma inverse function.
+  x = gaminv(p,v/2,2);
+  
+  % Return NaN if the degrees of freedom is not a positive integer.
+  k = find(v < 0  |  round(v) ~= v);
+  if any(k)
+    tmp  = NaN;
+    x(k) = tmp(ones(size(k)));
+  end
+  
+end
+
+%--------------------------------------------------------------------------
+
+function x = gaminv(p,a,b)
+  %GAMINV Inverse of the gamma cumulative distribution function (cdf).
+  %   X = GAMINV(P,A,B)  returns the inverse of the gamma cdf with
+  %   parameters A and B, at the probabilities in P.
+  %
+  %   The size of X is the common size of the input arguments. A scalar input
+  %   functions as a constant matrix of the same size as the other inputs.
+  %
+  %   GAMINV uses Newton's method to converge to the solution.
+  
+  %   References:
+  %      [1]  M. Abramowitz and I. A. Stegun, "Handbook of Mathematical
+  %      Functions", Government Printing Office, 1964, 6.5.
+  
+  %   B.A. Jones 1-12-93
+  %   Was: Revision: 1.2, Date: 1996/07/25 16:23:36
+  
+  if nargin<3,
+    b=1;
+  end
+  
+  [errorcode p a b] = distchck(3,p,a,b);
+  
+  if errorcode > 0
+    error(generatemsgid('InvalidDimensions'),'The arguments must be the same size or be scalars.');
+  end
+  
+  %   Initialize X to zero.
+  x = zeros(size(p));
+  
+  k = find(p<0 | p>1 | a <= 0 | b <= 0);
+  if any(k),
+    tmp = NaN;
+    x(k) = tmp(ones(size(k)));
+  end
+  
+  % The inverse cdf of 0 is 0, and the inverse cdf of 1 is 1.
+  k0 = find(p == 0 & a > 0 & b > 0);
+  if any(k0),
+    x(k0) = zeros(size(k0));
+  end
+  
+  k1 = find(p == 1 & a > 0 & b > 0);
+  if any(k1),
+    tmp = Inf;
+    x(k1) = tmp(ones(size(k1)));
+  end
+  
+  % Newton's Method
+  % Permit no more than count_limit iterations.
+  count_limit = 100;
+  count = 0;
+  
+  k = find(p > 0  &  p < 1 & a > 0 & b > 0);
+  pk = p(k);
+  
+  % Supply a starting guess for the iteration.
+  %   Use a method of moments fit to the lognormal distribution.
+  mn = a(k) .* b(k);
+  v = mn .* b(k);
+  temp = log(v + mn .^ 2);
+  mu = 2 * log(mn) - 0.5 * temp;
+  sigma = -2 * log(mn) + temp;
+  xk = exp(norminv(pk,mu,sigma));
+  
+  h = ones(size(pk));
+  
+  % Break out of the iteration loop for three reasons:
+  %  1) the last update is very small (compared to x)
+  %  2) the last update is very small (compared to sqrt(eps))
+  %  3) There are more than 100 iterations. This should NEVER happen.
+  
+  while(any(abs(h) > sqrt(eps)*abs(xk))  &  max(abs(h)) > sqrt(eps)    ...
+      & count < count_limit),
+    
+    count = count + 1;
+    h = (gamcdf(xk,a(k),b(k)) - pk) ./ gampdf(xk,a(k),b(k));
+    xnew = xk - h;
+    % Make sure that the current guess stays greater than zero.
+    % When Newton's Method suggests steps that lead to negative guesses
+    % take a step 9/10ths of the way to zero:
+    ksmall = find(xnew < 0);
+    if any(ksmall),
+      xnew(ksmall) = xk(ksmall) / 10;
+      h = xk-xnew;
+    end
+    xk = xnew;
+  end
+  
+  
+  % Store the converged value in the correct place
+  x(k) = xk;
+  
+  if count == count_limit,
+    fprintf('\nWarning: GAMINV did not converge.\n');
+    str = 'The last step was:  ';
+    outstr = sprintf([str,'%13.8f'],h);
+    fprintf(outstr);
+  end
+end
+
+%--------------------------------------------------------------------------
+function x = norminv(p,mu,sigma)
+  %NORMINV Inverse of the normal cumulative distribution function (cdf).
+  %   X = NORMINV(P,MU,SIGMA) finds the inverse of the normal cdf with
+  %   mean, MU, and standard deviation, SIGMA.
+  %
+  %   The size of X is the common size of the input arguments. A scalar input
+  %   functions as a constant matrix of the same size as the other inputs.
+  %
+  %   Default values for MU and SIGMA are 0 and 1 respectively.
+  
+  %   References:
+  %      [1]  M. Abramowitz and I. A. Stegun, "Handbook of Mathematical
+  %      Functions", Government Printing Office, 1964, 7.1.1 and 26.2.2
+  
+  %   Was: Revision: 1.2, Date: 1996/07/25 16:23:36
+  
+  if nargin < 3,
+    sigma = 1;
+  end
+  
+  if nargin < 2;
+    mu = 0;
+  end
+  
+  [errorcode p mu sigma] = distchck(3,p,mu,sigma);
+  
+  if errorcode > 0
+    error(generatemsgid('InvalidDimensions'),'Requires non-scalar arguments to match in size.');
+  end
+  
+  % Allocate space for x.
+  x = zeros(size(p));
+  
+  % Return NaN if the arguments are outside their respective limits.
+  k = find(sigma <= 0 | p < 0 | p > 1);
+  if any(k)
+    tmp  = NaN;
+    x(k) = tmp(ones(size(k)));
+  end
+  
+  % Put in the correct values when P is either 0 or 1.
+  k = find(p == 0);
+  if any(k)
+    tmp  = Inf;
+    x(k) = -tmp(ones(size(k)));
+  end
+  
+  k = find(p == 1);
+  if any(k)
+    tmp  = Inf;
+    x(k) = tmp(ones(size(k)));
+  end
+  
+  % Compute the inverse function for the intermediate values.
+  k = find(p > 0  &  p < 1 & sigma > 0);
+  if any(k),
+    x(k) = sqrt(2) * sigma(k) .* erfinv(2 * p(k) - 1) + mu(k);
+  end
+end
+
+%--------------------------------------------------------------------------
+function p = gamcdf(x,a,b)
+  %GAMCDF Gamma cumulative distribution function.
+  %   P = GAMCDF(X,A,B) returns the gamma cumulative distribution
+  %   function with parameters A and B at the values in X.
+  %
+  %   The size of P is the common size of the input arguments. A scalar input
+  %   functions as a constant matrix of the same size as the other inputs.
+  %
+  %   Some references refer to the gamma distribution with a single
+  %   parameter. This corresponds to the default of B = 1.
+  %
+  %   GAMMAINC does computational work.
+  
+  %   References:
+  %      [1]  L. Devroye, "Non-Uniform Random Variate Generation",
+  %      Springer-Verlag, 1986. p. 401.
+  %      [2]  M. Abramowitz and I. A. Stegun, "Handbook of Mathematical
+  %      Functions", Government Printing Office, 1964, 26.1.32.
+  
+  %   Was: Revision: 1.2, Date: 1996/07/25 16:23:36
+  
+  if nargin < 3,
+    b = 1;
+  end
+  
+  if nargin < 2,
+    error(generatemsgid('Nargchk'),'Requires at least two input arguments.');
+  end
+  
+  [errorcode x a b] = distchck(3,x,a,b);
+  
+  if errorcode > 0
+    error(generatemsgid('InvalidDimensions'),'Requires non-scalar arguments to match in size.');
+  end
+  
+  %   Return NaN if the arguments are outside their respective limits.
+  k1 = find(a <= 0 | b <= 0);
+  if any(k1)
+    tmp   = NaN;
+    p(k1) = tmp(ones(size(k1)));
+  end
+  
+  % Initialize P to zero.
+  p = zeros(size(x));
+  
+  k = find(x > 0 & ~(a <= 0 | b <= 0));
+  if any(k),
+    p(k) = gammainc(x(k) ./ b(k),a(k));
+  end
+  
+  % Make sure that round-off errors never make P greater than 1.
+  k = find(p > 1);
+  if any(k)
+    p(k) = ones(size(k));
+  end
+end
+
+%--------------------------------------------------------------------------
+function y = gampdf(x,a,b)
+  %GAMPDF Gamma probability density function.
+  %   Y = GAMPDF(X,A,B) returns the gamma probability density function
+  %   with parameters A and B, at the values in X.
+  %
+  %   The size of Y is the common size of the input arguments. A scalar input
+  %   functions as a constant matrix of the same size as the other inputs.
+  %
+  %   Some references refer to the gamma distribution with a single
+  %   parameter. This corresponds to the default of B = 1.
+  
+  %   References:
+  %      [1]  L. Devroye, "Non-Uniform Random Variate Generation",
+  %      Springer-Verlag, 1986, pages 401-402.
+  
+  %   Was: Revision: 1.2, Date: 1996/07/25 16:23:36
+  
+  if nargin < 3,
+    b = 1;
+  end
+  
+  if nargin < 2,
+    error(generatemsgid('Nargchk'),'Requires at least two input arguments');
+  end
+  
+  [errorcode x a b] = distchck(3,x,a,b);
+  
+  if errorcode > 0
+    error(generatemsgid('InvalidDimensions'),'Requires non-scalar arguments to match in size.');
+  end
+  
+  % Initialize Y to zero.
+  y = zeros(size(x));
+  
+  %   Return NaN if the arguments are outside their respective limits.
+  k1 = find(a <= 0 | b <= 0);
+  if any(k1)
+    tmp = NaN;
+    y(k1) = tmp(ones(size(k1)));
+  end
+  
+  k=find(x > 0 & ~(a <= 0 | b <= 0));
+  if any(k)
+    y(k) = (a(k) - 1) .* log(x(k)) - (x(k) ./ b(k)) - gammaln(a(k)) - a(k) .* log(b(k));
+    y(k) = exp(y(k));
+  end
+  k1 = find(x == 0 & a < 1);
+  if any(k1)
+    tmp = Inf;
+    y(k1) = tmp(ones(size(k1)));
+  end
+  k2 = find(x == 0 & a == 1);
+  if any(k2)
+    y(k2) = (1./b(k2));
+  end
+end
+
+function [errorcode,out1,out2,out3,out4] = distchck(nparms,arg1,arg2,arg3,arg4)
+  %DISTCHCK Checks the argument list for the probability functions.
+  
+  %   B.A. Jones  1-22-93
+  %   Was: Revision: 1.2, Date: 1996/07/25 16:23:36
+  
+  errorcode = 0;
+  
+  if nparms == 1
+    out1 = arg1;
+    return;
+  end
+  
+  if nparms == 2
+    [r1 c1] = size(arg1);
+    [r2 c2] = size(arg2);
+    scalararg1 = (prod(size(arg1)) == 1);
+    scalararg2 = (prod(size(arg2)) == 1);
+    if ~scalararg1 & ~scalararg2
+      if r1 ~= r2 | c1 ~= c2
+        errorcode = 1;
+        return;
+      end
+    end
+    if scalararg1
+      out1 = arg1(ones(r2,1),ones(c2,1));
+    else
+      out1 = arg1;
+    end
+    if scalararg2
+      out2 = arg2(ones(r1,1),ones(c1,1));
+    else
+      out2 = arg2;
+    end
+  end
+  
+  if nparms == 3
+    [r1 c1] = size(arg1);
+    [r2 c2] = size(arg2);
+    [r3 c3] = size(arg3);
+    scalararg1 = (prod(size(arg1)) == 1);
+    scalararg2 = (prod(size(arg2)) == 1);
+    scalararg3 = (prod(size(arg3)) == 1);
+    
+    if ~scalararg1 & ~scalararg2
+      if r1 ~= r2 | c1 ~= c2
+        errorcode = 1;
+        return;
+      end
+    end
+    
+    if ~scalararg1 & ~scalararg3
+      if r1 ~= r3 | c1 ~= c3
+        errorcode = 1;
+        return;
+      end
+    end
+    
+    if ~scalararg3 & ~scalararg2
+      if r3 ~= r2 | c3 ~= c2
+        errorcode = 1;
+        return;
+      end
+    end
+    
+    if ~scalararg1
+      out1 = arg1;
+    end
+    if ~scalararg2
+      out2 = arg2;
+    end
+    if ~scalararg3
+      out3 = arg3;
+    end
+    rows = max([r1 r2 r3]);
+    columns = max([c1 c2 c3]);
+    
+    if scalararg1
+      out1 = arg1(ones(rows,1),ones(columns,1));
+    end
+    if scalararg2
+      out2 = arg2(ones(rows,1),ones(columns,1));
+    end
+    if scalararg3
+      out3 = arg3(ones(rows,1),ones(columns,1));
+    end
+    out4 =[];
+    
+  end
+  
+  if nparms == 4
+    [r1 c1] = size(arg1);
+    [r2 c2] = size(arg2);
+    [r3 c3] = size(arg3);
+    [r4 c4] = size(arg4);
+    scalararg1 = (prod(size(arg1)) == 1);
+    scalararg2 = (prod(size(arg2)) == 1);
+    scalararg3 = (prod(size(arg3)) == 1);
+    scalararg4 = (prod(size(arg4)) == 1);
+    
+    if ~scalararg1 & ~scalararg2
+      if r1 ~= r2 | c1 ~= c2
+        errorcode = 1;
+        return;
+      end
+    end
+    
+    if ~scalararg1 & ~scalararg3
+      if r1 ~= r3 | c1 ~= c3
+        errorcode = 1;
+        return;
+      end
+    end
+    
+    if ~scalararg1 & ~scalararg4
+      if r1 ~= r4 | c1 ~= c4
+        errorcode = 1;
+        return;
+      end
+    end
+    
+    if ~scalararg3 & ~scalararg2
+      if r3 ~= r2 | c3 ~= c2
+        errorcode = 1;
+        return;
+      end
+    end
+    
+    if ~scalararg4 & ~scalararg2
+      if r4 ~= r2 | c4 ~= c2
+        errorcode = 1;
+        return;
+      end
+    end
+    
+    if ~scalararg3 & ~scalararg4
+      if r3 ~= r4 | c3 ~= c4
+        errorcode = 1;
+        return;
+      end
+    end
+    
+    
+    if ~scalararg1
+      out1 = arg1;
+    end
+    if ~scalararg2
+      out2 = arg2;
+    end
+    if ~scalararg3
+      out3 = arg3;
+    end
+    if ~scalararg4
+      out4 = arg4;
+    end
+    
+    rows = max([r1 r2 r3 r4]);
+    columns = max([c1 c2 c3 c4]);
+    if scalararg1
+      out1 = arg1(ones(rows,1),ones(columns,1));
+    end
+    if scalararg2
+      out2 = arg2(ones(rows,1),ones(columns,1));
+    end
+    if scalararg3
+      out3 = arg3(ones(rows,1),ones(columns,1));
+    end
+    if scalararg4
+      out4 = arg4(ones(rows,1),ones(columns,1));
+    end
+  end
+end
+
+%----- 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: psdconf.m,v 1.14 2011/05/18 07:53:57 mauro 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 = plist({'conf', 'Required percentage confidence level. [0-100]'}, paramValue.DOUBLE_VALUE(95));
+end
+% END
+
+% PARAMETERS: 
+%             'conf'  - Required percentage confidence level. It is a
+%                       number between 0 and 100 [Default: 95,
+%                       corresponding to 95% confidence].