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 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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