Home > classes > @ao > xspec.m

xspec

PURPOSE ^

XSPEC performs cross-spectral analysis of various forms.

SYNOPSIS ^

function varargout = xspec(varargin)

DESCRIPTION ^

 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., ltpda_tfe).

 CALL:        b = xspec(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 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  - Na x Na matrix of output AOs

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

 HISTORY:     30-05-2007 M Hewitson
                 Creation

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

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 % XSPEC performs cross-spectral analysis of various forms.
0002 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0003 %
0004 % DESCRIPTION: XSPEC performs 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., ltpda_tfe).
0008 %
0009 % CALL:        b = xspec(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 %                       'mscohere' - estimate cross-coherence
0017 %              iALGO  - ALGONAME from the calling higher level script
0018 %              iVER   - VERSION from the calling higher level script
0019 %              invars - invars variable from the calling higher level script
0020 %
0021 % OUTPUTS:     b  - Na x Na matrix of output AOs
0022 %
0023 % VERSION:     $Id: xspec.m,v 1.15 2008/09/03 17:05:51 ingo Exp $
0024 %
0025 % HISTORY:     30-05-2007 M Hewitson
0026 %                 Creation
0027 %
0028 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0029 
0030 function varargout = xspec(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: xspec.m,v 1.15 2008/09/03 17:05:51 ingo Exp $';
0042 
0043   % Set the method version string in the minfo object
0044   mi.setMversion([VERSION '-->' mi.mversion]);
0045 
0046   na = numel(as);
0047   if na < 2
0048     error('### XSPEC needs at least two time-series AOs.');
0049   end
0050 
0051   %----------------- Resample all AOs
0052   copies = zeros(size(as));
0053 
0054   fsmax = ao.findFsMax(as);
0055   fspl  = plist(param('fsout', fsmax));
0056   for i=1:na
0057     % check this is a time-series object
0058     if ~isa(as(i).data, 'tsdata')
0059       error('### ao.xspec requires tsdata (time-series) inputs.');
0060     end
0061     % Check Fs
0062     if as(i).data.fs ~= fsmax
0063       utils.helper.msg(msg.PROC1, 'resampling AO %s to %f Hz', as(i).name, fsmax);
0064       % Make a deep copy so we don't
0065       % affect the original input data
0066       as(i) = copy(as(i), 1);
0067       copies(i) = 1;
0068       resample(as(i), fspl);
0069     end
0070   end
0071 
0072   %----------------- Truncate all vectors
0073 
0074   % Get shortest vector
0075   lmin = ao.findShortestVector(as);
0076   nsecs = lmin / fsmax;
0077   for i=1:na
0078     if len(as(i)) ~= lmin
0079       utils.helper.msg(msg.PROC2, 'truncating AO %s to %d secs', as(i).name, nsecs);
0080       % do we already have a copy?
0081       if ~copies(i)
0082         % Make a deep copy so we don't
0083         % affect the original input data
0084         as(i) = copy(as(i), 1);
0085         copies(i) = 1;
0086       end
0087       as(i).select(1:lmin);
0088     end
0089   end
0090 
0091   %----------------- check input parameters
0092 
0093   % points in FFT
0094   Nfft  = find(pl, 'Nfft');
0095   if ~isempty(Nfft)
0096     if Nfft < 0
0097       Nfft = abs(Nfft) * lmin;
0098       utils.helper.msg(msg.PROC2, 'using default Nfft of %g', Nfft);
0099       pl.pset('Nfft', Nfft);
0100     end
0101   end
0102 
0103   % Window object
0104   Win   = find(pl, 'Win');
0105   if ~isempty(Win)
0106     if length(Win.win)~=Nfft
0107       switch Win.type
0108         case 'Kaiser'
0109           Win = specwin(Win.type, Nfft, Win.psll);
0110         otherwise
0111           Win = specwin(Win.type, Nfft);
0112       end
0113       utils.helper.msg(msg.PROC2, 'reset window to %s(%d)', strrep(Win.type, '_', '\_'), length(Win.win));
0114       pl.pset('Win', Win);
0115     end
0116   end
0117 
0118   % desired overlap
0119   Olap = find(pl, 'Olap');
0120   if ~isempty(Olap)
0121     if Olap < 0
0122       % use window function rov
0123       Olap = Win.rov;
0124       utils.helper.msg(msg.PROC2, 'using recommended overlap of %2.1f%%', Olap);
0125       pl.pset('Olap', Olap);
0126     end
0127   end
0128 
0129   % desired detrending order
0130   order = find(pl, 'Order');
0131   if isempty(order)
0132     order = 0;
0133     utils.helper.msg(msg.PROC1, 'using default detrending order 0 (mean)');
0134     pl.append('Order', order);
0135   end
0136 
0137   % Loop over input AOs
0138   bs(na,na) = ao;
0139   for j=1:na
0140     for k=1:na
0141       if k>=j
0142         utils.helper.msg(msg.PROC1, 'computing (%d,%d) %s -> %s', j, k, invars{j}, invars{k});
0143 
0144         % -------- Make Xspec estimate
0145 
0146         % Compute xspec using welch
0147         [txy, f, info] = ao.welch(as(j), as(k), method, pl);
0148 
0149         % Keep the data shape of the first input AO
0150         if size(as(1).data.y,1) == 1
0151           txy = txy.';
0152           f   = f.';
0153         end
0154 
0155         % create new output fsdata
0156         fs = fsdata(f, txy, fsmax);
0157         fs.setEnbw(info.enbw);
0158         fs.setXunits('Hz');
0159         fs.setYunits(info.units);
0160         fs.setNavs(info.navs);
0161 
0162         % make output analysis object
0163         bs(j,k) = ao(fs);
0164 
0165         %----------- Add history
0166 
0167         % we need to get the input histories in the same order as the inputs
0168         % to this function call, not in the order of the input to tfestimate;
0169         % otherwise the resulting matrix on a 'create from history' will be
0170         % mirrored.
0171         if j<k
0172           bs(j,k).addHistory(mi, pl, [invars(j) invars(k)], [as(j).hist, as(k).hist]);
0173         else
0174           bs(j,k).addHistory(mi, pl, [invars(j) invars(k)], [as(k).hist, as(j).hist]);
0175         end
0176 
0177         % set name
0178         bs(j,k).setName(sprintf('%s(%s->%s)', method, invars{j}, invars{k}), 'internal');
0179 
0180       end % End upper triangle
0181     end % End second loop over AOs
0182   end % End first loop over AOs
0183 
0184   % Fill remainder
0185   for j=1:na
0186     for k=1:na
0187       if k<j
0188         % Do we need a copy here? Perhaps
0189         % copying the pointer is enough...
0190         bs(j,k) = bs(k,j);
0191       end
0192     end
0193   end
0194 
0195   % Set output
0196   varargout{1} = bs;
0197 end
0198

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