line source
+ − % 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
+ −