--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/m-toolbox/classes/@ao/xspec.m	Wed Nov 23 19:22:13 2011 +0100
@@ -0,0 +1,197 @@
+% 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
+