--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/m-toolbox/classes/@ao/lxspec.m	Wed Nov 23 19:22:13 2011 +0100
@@ -0,0 +1,208 @@
+% LXSPEC performs log-scale cross-spectral analysis of various forms.
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%
+% DESCRIPTION: LXSPEC performs log-scale 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., ltfe).
+%
+% CALL:       b = lxspec(a, pl, method, iALGO, iVER, invars);
+%
+% INPUTS:     a      - vector of 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
+%             mi     - minfo object for calling method
+%             invars - invars variable from the calling higher level script
+%
+% OUTPUTS:    b  - output AO
+%
+% VERSION:    $Id: lxspec.m,v 1.44 2011/04/05 08:12:40 mauro Exp $
+%
+% HISTORY:    16-05-2008 M Hewitson
+%                Creation
+%
+%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+function varargout = lxspec(varargin)
+  
+  import utils.const.*
+  
+  % unpack inputs
+  as       = varargin{1};
+  pl       = varargin{2};
+  method   = varargin{3};
+  mi       = varargin{4};
+  invars   = varargin{5};
+  
+  VERSION  = '$Id: lxspec.m,v 1.44 2011/04/05 08:12:40 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
+  % Check if there are some AOs left
+  if numel(tsao) ~= 2
+    error('### LXSPEC 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
+  
+  copies = zeros(size(tsao));
+  
+  %----------------- Resample all AOs
+  fsmax = ao.findFsMax(tsao);
+  fspl  = plist(param('fsout', fsmax));
+  for ll = 1:numel(tsao)
+    % Check Fs
+    if tsao(ll).data.fs ~= fsmax
+      utils.helper.msg(msg.PROC2, 'resampling AO %s to %f Hz', tsao(ll).name, fsmax);
+      % Make a deep copy so we don't
+      % affect the original input data
+      tsao(ll) = copy(tsao(ll), 1);
+      copies(ll) = 1;
+      tsao(ll).resample(fspl);
+    end
+  end
+  
+  %----------------- Truncate all vectors
+  % Get shortest vector
+  lmin = ao.findShortestVector(tsao);
+  nsecs = lmin / fsmax;
+  for ll = 1:numel(tsao)
+    if len(tsao(ll)) ~= lmin
+      utils.helper.msg(msg.PROC1, '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
+  
+  %----------------- Build signal Matrix
+  N     = len(tsao(1)); % length of first signal
+  iS    = zeros(numel(tsao), N);
+  for jj = 1:numel(tsao)
+    iS(jj,:) = tsao(jj).data.getY;
+  end
+  
+  %----------------- check input parameters
+  pl = utils.helper.process_spectral_options(pl, 'log');
+  
+  % Desired number of averages
+  Kdes = find(pl, 'Kdes');
+  % num desired spectral frequencies
+  Jdes = find(pl, 'Jdes');
+  % Minimum segment length
+  Lmin = find(pl, 'Lmin');
+  % Window function
+  Win = find(pl, 'Win');
+  % Overlap
+  Nolap = find(pl, 'Olap')/100;
+  % Order of detrending
+  Order = find(pl, 'Order');
+  
+  %----------------- Get frequency vector
+  [f, r, m, L, K] = ao.ltf_plan(lmin, fsmax, Nolap, 1, Lmin, Jdes, Kdes);
+  
+  %----------------- compute TF Estimates
+  [Txy dev]= ao.mltfe(iS, f, r, m, L,K,fsmax, Win, Order, Nolap*100, Lmin, method);
+  
+  % Keep the data shape of the first AO
+  if size(tsao(1).data.y, 1) == 1
+    f   = f.';
+    Txy = Txy.';
+    dev = dev.';
+  end
+  
+  %----------------- Build output Matrix of AOs
+  
+  % create new output fsdata
+  fsd = fsdata(f, Txy, fsmax);
+  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 {'cohere','mscohere'}
+      fsd.setYunits(unit());
+    otherwise
+      error(['### Unknown method:' method]);
+  end
+  
+  % 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 standard deviation to dy field
+  bs.data.dy = dev;
+  % simplify the units
+  if strcmp(method, 'cpsd')
+    bs.simplifyYunits(...
+      plist('prefixes',false,'exceptions','Hz'),...
+      'internal');
+  end
+  
+  % set name
+  bs.name = sprintf('L%s(%s->%s)', upper(method), invars{1}, invars{2});
+  % set procinfo
+  bs.procinfo = combine(bs.procinfo,plist('r', r, 'm', m, 'l', L, 'k', K));
+  % Propagate 'plotinfo'
+  if ~isempty(tsao(1).plotinfo) || ~isempty(tsao(2).plotinfo)
+    bs.plotinfo = combine(tsao(1).plotinfo, tsao(2).plotinfo);
+  end
+  % Add history
+  bs.addHistory(mi, pl, [invars(:)], inhists);
+  
+  
+  % Set output
+  varargout{1} = bs;
+end
+% END