Home > classes > @ao > lxspec.m

lxspec

PURPOSE ^

LXSPEC performs log-scale cross-spectral analysis of various forms.

SYNOPSIS ^

function varargout = lxspec(varargin)

DESCRIPTION ^

 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
                       'cohere'   - estimate cross-coherence
             mi     - minfo object for calling method
             invars - invars variable from the calling higher level script

 OUTPUTS:    b  - Na x Na matrix of output AOs

 VERSION:    $Id: lxspec.m,v 1.10 2008/09/03 17:05:51 ingo Exp $

 HISTORY:    16-05-2008 M Hewitson
                Creation


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 % LXSPEC performs log-scale cross-spectral analysis of various forms.
0002 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0003 %
0004 % DESCRIPTION: LXSPEC performs log-scale cross-spectral analysis of various forms.
0005 %              The function is a helper function for various higher level
0006 %              functions. It is meant to be called from other functions
0007 %              (e.g., ltfe).
0008 %
0009 % CALL:       b = lxspec(a, pl, method, iALGO, iVER, invars);
0010 %
0011 % INPUTS:     a      - vector of input AOs
0012 %             pl     - input parameter list
0013 %             method - one of
0014 %                       'cpsd'     - compute cross-spectral density
0015 %                       'tfe'      - estimate transfer function between inputs
0016 %                       'cohere'   - estimate cross-coherence
0017 %             mi     - minfo object for calling method
0018 %             invars - invars variable from the calling higher level script
0019 %
0020 % OUTPUTS:    b  - Na x Na matrix of output AOs
0021 %
0022 % VERSION:    $Id: lxspec.m,v 1.10 2008/09/03 17:05:51 ingo Exp $
0023 %
0024 % HISTORY:    16-05-2008 M Hewitson
0025 %                Creation
0026 %
0027 %
0028 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0029 
0030 function varargout = lxspec(varargin)
0031 
0032   import utils.const.*
0033   
0034   % unpack inputs
0035   as       = varargin{1};
0036   pl       = varargin{2};
0037   method   = varargin{3};
0038   mi       = varargin{4};
0039   invars   = varargin{5};
0040 
0041   VERSION  = '$Id: lxspec.m,v 1.10 2008/09/03 17:05:51 ingo Exp $';
0042   % Set the method version string in the minfo object
0043   mi.setMversion([VERSION '-->' mi.mversion]);
0044 
0045   na     = numel(as);
0046   copies = zeros(size(as));
0047 
0048   %----------------- Resample all AOs
0049   fsmax = ao.findFsMax(as);
0050   fspl  = plist(param('fsout', fsmax));
0051   for i=1:na
0052     % check this is a time-series object
0053     if ~isa(as(i).data, 'tsdata')
0054       warning('### xspec requires tsdata (time-series) inputs. Skipping AO %s', invars{i});
0055       as(i) = [];
0056     else
0057       % Check Fs
0058       if as(i).data.fs ~= fsmax
0059         utils.helper.msg(msg.PROC2, 'resampling AO %s to %f Hz', as(i).name, fsmax);
0060         % Make a deep copy so we don't
0061         % affect the original input data
0062         as(i) = copy(as(i), 1);
0063         copies(i) = 1;
0064         as(i) = resample(as(i), fspl);
0065       end
0066     end
0067   end
0068 
0069   %----------------- Truncate all vectors
0070   % Get shortest vector
0071   lmin = ao.findShortestVector(as);
0072   nsecs = lmin / fsmax;
0073   for i=1:na
0074     if len(as(i)) ~= lmin
0075       utils.helper.msg(msg.PROC1, 'truncating AO %s to %d secs', as(i).name, nsecs);
0076       % do we already have a copy?
0077       if ~copies(i)
0078         % Make a deep copy so we don't
0079         % affect the original input data
0080         as(i) = copy(as(i), 1);
0081         copies(i) = 1;
0082       end
0083       as(i) = select(as(i), 1:lmin);
0084     end
0085   end
0086 
0087   %----------------- Build signal Matrix
0088   N     = len(as(1)); % length of first signal
0089   iS    = zeros(na, N);
0090   for j=1:na
0091     iS(j,:) = as(j).data.getY;
0092   end
0093 
0094   %----------------- check input parameters
0095 
0096   % Averaging
0097   Kdes = find(pl, 'Kdes');    % desired averages
0098   Jdes = find(pl, 'Jdes');    % num desired spectral frequencies
0099   Win   = find(pl, 'Win');    % Window object
0100   Nolap = find(pl, 'Olap');  % desired overlap
0101   if Nolap == -1
0102     % use window function rov
0103     Nolap = Win.rov/100;
0104     utils.helper.msg(msg.PROC1, 'using recommended overlap of %d', Nolap);
0105     pl = pset(pl, 'Olap', Nolap);
0106   end
0107   Order = find(pl, 'Order');    % desired overlap
0108   Lmin = find(pl, 'Lmin');  % Minimum segment length
0109 
0110   %----------------- Get frequency vector
0111   [f, r, m, L, K] = ao.ltf_plan(lmin, fsmax, Nolap, 1, Lmin, Jdes, Kdes);
0112 
0113   %----------------- compute TF Estimates
0114   Txy = ao.mltfe(iS, f, r, m, L, fsmax, Win, Order, Nolap*100, Lmin, method);
0115 
0116   %----------------- Build output Matrix of AOs
0117   utils.helper.msg(msg.PROC1, 'building output AO matrix');
0118   bs(na,na) = ao;
0119   for j=1:na  % Loop over input signal
0120     for k=1:na  % Loop over output signal
0121       % Keep the data shape of the first AO
0122       y = squeeze(Txy(j,k,:));
0123       if size(as(1).data.y,1) == 1
0124         f = f.';
0125         y = y.';
0126       end
0127       % create new output fsdata
0128       fsd = fsdata(f, y, fsmax);
0129       fsd.setYunits(as(k).data.yunits / as(j).data.yunits);
0130       % make output analysis object
0131       bs(j,k) = ao(fsd);
0132       % add history
0133       bs(j,k).addHistory(mi, pl, [invars(j) invars(k)], [as(j).hist, as(k).hist]);
0134       % set name
0135       bs(j,k).setName(sprintf('%s(%s->%s)', method, invars{j}, invars{k}), 'internal');
0136     end % End loop over output signal
0137   end % End loop over input signal
0138 
0139   % Set output
0140   varargout{1} = bs;
0141 end
0142 % END

Generated on Mon 08-Sep-2008 13:18:47 by m2html © 2003