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