line source
+ − % XSPEC performs cross-spectral analysis of various forms.
+ − %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ − %
+ − % DESCRIPTION: XSPEC performs cross-spectral analysis of various forms.
+ − % The function is a helper function for various higher level
+ − % functions. It is meant to be called from other functions
+ − % (e.g., tfe).
+ − %
+ − % CALL: b = xspec(a, pl, method, iALGO, iVER, invars);
+ − %
+ − % INPUTS: a - vector of 2 input AOs
+ − % pl - input parameter list
+ − % method - one of
+ − % 'cpsd' - compute cross-spectral density
+ − % 'tfe' - estimate transfer function between inputs
+ − % 'mscohere' - estimate magnitude-squared cross-coherence
+ − % 'cohere' - estimate complex cross-coherence
+ − % iALGO - ALGONAME from the calling higher level script
+ − % iVER - VERSION from the calling higher level script
+ − % invars - invars variable from the calling higher level script
+ − %
+ − % OUTPUTS: b - output AO
+ − %
+ − % VERSION: $Id: xspec.m,v 1.61 2011/05/22 22:09:44 mauro Exp $
+ − %
+ − % HISTORY: 30-05-2007 M Hewitson
+ − % Creation
+ − %
+ − %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ −
+ − function varargout = xspec(varargin)
+ −
+ − import utils.const.*
+ −
+ − % unpack inputs
+ − as = varargin{1};
+ − pl = varargin{2};
+ − method = varargin{3};
+ − mi = varargin{4};
+ − invars = varargin{5};
+ −
+ − VERSION = '$Id: xspec.m,v 1.61 2011/05/22 22:09:44 mauro Exp $';
+ −
+ − % Set the method version string in the minfo object
+ − mi.setMversion([VERSION '-->' mi.mversion]);
+ −
+ − %----------------- Select all AOs with time-series data
+ − tsao = [];
+ − for ll = 1:numel(as)
+ − if isa(as(ll).data, 'tsdata')
+ − tsao = [tsao as(ll)];
+ − else
+ − warning('### xspec requires tsdata (time-series) inputs. Skipping AO %s. \nREMARK: The output doesn''t contain this AO', invars{ll});
+ − end
+ − end
+ − if numel(tsao) ~= 2
+ − error('### XSPEC needs two time-series AOs.');
+ − end
+ −
+ − %----------------- Gather the input history objects
+ − inhists = [tsao(:).hist];
+ −
+ − %----------------- Check the time range
+ − time_range = mfind(pl, 'split', 'times');
+ − for ll = 1:numel(tsao)
+ − if ~isempty(time_range)
+ − switch class(time_range)
+ − case 'double'
+ − tsao(ll) = split(tsao(ll), plist(...
+ − 'times', time_range));
+ − case 'timespan'
+ − tsao(ll) = split(tsao(ll), plist(...
+ − 'timespan', time_range));
+ − case 'time'
+ − tsao(ll) = split(tsao(ll), plist(...
+ − 'start_time', time_range(1), ...
+ − 'end_time', time_range(2)));
+ − case 'cell'
+ − tsao(ll) = split(tsao(ll), plist(...
+ − 'start_time', time_range{1}, ...
+ − 'end_time', time_range{2}));
+ − otherwise
+ − end
+ − end
+ − if tsao(ll).len <= 0
+ − error('### The object is empty! Please revise your settings ...');
+ − end
+ − end
+ −
+ − %----------------- Resample all AOs
+ − copies = zeros(size(tsao));
+ −
+ − fsmin = ao.findFsMin(tsao);
+ − fspl = plist(param('fsout', fsmin));
+ − for ll = 1:numel(tsao)
+ − % Check Fs
+ − if tsao(ll).data.fs ~= fsmin
+ − utils.helper.msg(msg.PROC1, 'resampling AO %s to %f Hz', tsao(ll).name, fsmin);
+ − % Make a deep copy so we don't
+ − % affect the original input data
+ − tsao(ll) = copy(tsao(ll), 1);
+ − copies(ll) = 1;
+ − resample(tsao(ll), fspl);
+ − end
+ − end
+ −
+ − %----------------- Truncate all vectors
+ −
+ − % Get shortest vector
+ − lmin = ao.findShortestVector(tsao);
+ − nsecs = lmin / fsmin;
+ − for ll = 1:numel(tsao)
+ − if len(tsao(ll)) ~= lmin
+ − utils.helper.msg(msg.PROC2, 'truncating AO %s to %d secs', tsao(ll).name, nsecs);
+ − % do we already have a copy?
+ − if ~copies(ll)
+ − % Make a deep copy so we don't
+ − % affect the original input data
+ − tsao(ll) = copy(tsao(ll), 1);
+ − copies(ll) = 1;
+ − end
+ − tsao(ll).select(1:lmin);
+ − end
+ − end
+ −
+ − %----------------- check input parameters
+ −
+ − usepl = utils.helper.process_spectral_options(pl, 'lin', lmin, fsmin);
+ −
+ − % Loop over input AOs
+ − utils.helper.msg(msg.PROC1, 'computing %s(%s -> %s)', method, invars{1}, invars{2});
+ −
+ − % -------- Make Xspec estimate
+ −
+ − % Compute xspec using welch and always scale to PSD
+ − [txy, f, info, dev] = welch(tsao(1), tsao(2), method, usepl.pset('Scale', 'PSD'));
+ −
+ − % Keep the data shape of the first input AO
+ − if size(tsao(1).data.y,1) == 1
+ − txy = txy.';
+ − f = f.';
+ − dev = dev.';
+ − end
+ −
+ − % create new output fsdata
+ − fsd = fsdata(f, txy, fsmin);
+ − fsd.setEnbw(info.enbw);
+ − fsd.setXunits('Hz');
+ −
+ − switch lower(method)
+ − case 'tfe'
+ − fsd.setYunits(tsao(2).data.yunits / tsao(1).data.yunits);
+ − case 'cpsd'
+ − fsd.setYunits(tsao(2).data.yunits * tsao(1).data.yunits / unit('Hz'));
+ − case {'mscohere','cohere'}
+ − fsd.setYunits(unit());
+ − otherwise
+ − error(['### Unknown method:' method]);
+ − end
+ − fsd.setNavs(info.navs);
+ −
+ − % set object t0
+ − if eq(tsao(1).t0, tsao(2).t0)
+ − fsd.setT0(tsao(1).t0);
+ − end
+ −
+ − % make output analysis object
+ − bs = ao(fsd);
+ − % add variance
+ − bs.data.dy = dev;
+ −
+ − % simplify the units in the case of cpsd calculation
+ − if strcmp(method, 'cpsd')
+ − bs.simplifyYunits(...
+ − plist('prefixes', false, 'exceptions', 'Hz'));
+ − end
+ −
+ − %----------- Add history
+ −
+ − % set name
+ − bs.name = sprintf('%s(%s->%s)', upper(method), invars{1}, invars{2});
+ −
+ − % we need to get the input histories in the same order tsao the inputs
+ − % to this function call, not in the order of the input to xspec;
+ − % otherwise the resulting matrix on a 'create from history' will be
+ − % mirrored.
+ − bs.addHistory(mi, usepl, [invars(:)], inhists);
+ −
+ − % Propagate 'plotinfo'
+ − if ~isempty(tsao(1).plotinfo) || ~isempty(tsao(2).plotinfo)
+ − bs.plotinfo = combine(tsao(1).plotinfo, tsao(2).plotinfo);
+ − end
+ −
+ − % Set output
+ − varargout{1} = bs;
+ − end
+ −