Home > m > sigproc > frequency_domain > ltpda_xspec.m

ltpda_xspec

PURPOSE ^

LTPDA_XSPEC performs cross-spectral analysis of various forms.

SYNOPSIS ^

function varargout = ltpda_xspec(varargin)

DESCRIPTION ^

 LTPDA_XSPEC performs cross-spectral analysis of various forms.

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

 DESCRIPTION: LTPDA_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 = ltpda_xspec(a, pl, method, iALGO, iVER, invars);

 INPUTS:     a      - vector of input AOs
             pl     - input parameter list
             method - one of 'TF', 'CPSD', 'COHERE'
             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: ltpda_xspec.m,v 1.13 2008/03/30 23:15:26 mauro Exp $

 HISTORY:    30-05-2007 M Hewitson
                Creation

 The following call returns a parameter list object that contains the
 default parameter values:

 >> pl = ltpda_xspec('Params')

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

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

0001 function varargout = ltpda_xspec(varargin)
0002 % LTPDA_XSPEC performs cross-spectral analysis of various forms.
0003 %
0004 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0005 %
0006 % DESCRIPTION: LTPDA_XSPEC performs cross-spectral analysis of various forms.
0007 %              The function is a helper function for various higher level
0008 %              functions. It is meant to be called from other functions
0009 %              (e.g., ltpda_tfe).
0010 %
0011 % CALL:       b = ltpda_xspec(a, pl, method, iALGO, iVER, invars);
0012 %
0013 % INPUTS:     a      - vector of input AOs
0014 %             pl     - input parameter list
0015 %             method - one of 'TF', 'CPSD', 'COHERE'
0016 %             iALGO  - ALGONAME from the calling higher level script
0017 %             iVER   - VERSION from the calling higher level script
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: ltpda_xspec.m,v 1.13 2008/03/30 23:15:26 mauro Exp $
0023 %
0024 % HISTORY:    30-05-2007 M Hewitson
0025 %                Creation
0026 %
0027 % The following call returns a parameter list object that contains the
0028 % default parameter values:
0029 %
0030 % >> pl = ltpda_xspec('Params')
0031 %
0032 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0033 
0034 % unpack inputs
0035 as     = varargin{1};
0036 pl     = varargin{2};
0037 method = varargin{3};
0038 iALGO  = varargin{4};
0039 iVER   = varargin{5};
0040 invars = varargin{6};
0041 
0042 ALGONAME = iALGO;
0043 VERSION  = [iVER '/' '$Id: ltpda_xspec.m,v 1.13 2008/03/30 23:15:26 mauro Exp $'];
0044 
0045 
0046 % initialise output array
0047 bs = [];
0048 
0049 na = length(as);
0050 if na < 2
0051   error('### LTPDA_XSPEC needs at least two AOs to make a transfer function.');
0052 end
0053 
0054 
0055 %----------------- Resample all AOs
0056 
0057 disp('*** Resampling AOs...');
0058 fsmax = findFsMax(as);
0059 fspl  = plist(param('fsout', fsmax));
0060 bs = [];
0061 for i=1:na
0062   a = as(i);
0063   % check this is a time-series object
0064   if ~isa(a.data, 'tsdata')
0065     error('### ltpda_xspec requires tsdata (time-series) inputs.');
0066   end
0067   % Check Fs
0068   if a.data.fs ~= fsmax
0069     warning(sprintf('!!! Resampling AO %s/%s to %f Hz', a.name, a.data.name, fsmax));
0070     a = resample(a, fspl);
0071   end
0072   bs = [bs a];
0073 end
0074 as = bs;
0075 clear bs;
0076 disp('*** Done.');
0077 
0078 
0079 %----------------- Truncate all vectors
0080 
0081 % Get shortest vector
0082 disp('*** Truncating all vectors...');
0083 lmin = findShortestVector(as);
0084 nsecs = lmin / fsmax;
0085 bs = [];
0086 for i=1:na
0087   a = as(i);
0088   if len(a) ~= lmin
0089     warning(sprintf('!!! Truncating AO %s/%s to %d secs', a.name, a.data.name, nsecs));
0090     bs = [bs  select(a, 1:lmin)];
0091   else
0092     bs = [bs a];
0093   end
0094 end
0095 as = bs;
0096 clear bs;
0097 disp('*** Done.');
0098 
0099 %----------------- check input parameters
0100 
0101 % points in FFT
0102 Nfft  = find(pl, 'Nfft');
0103 if ~isempty(Nfft)
0104   if Nfft < 0
0105     Nfft = abs(Nfft) * lmin;
0106     disp(sprintf('! Using default Nfft of %g', Nfft))
0107     pl = pset(pl, 'Nfft', Nfft);
0108   end
0109 end
0110 
0111 % Window object
0112 Win   = find(pl, 'Win');
0113 if ~isempty(Win)
0114   if length(Win.win)~=Nfft
0115     switch Win.name
0116       case 'Kaiser'
0117         Win = specwin(Win.name, Nfft, Win.psll);
0118       otherwise
0119         Win = specwin(Win.name, Nfft);
0120     end
0121     disp(sprintf('! Reset window to %s(%d)', strrep(Win.name, '_', '\_'), length(Win.win)))
0122     pl = pset(pl, 'Win', Win);
0123   end
0124 end
0125 
0126 % desired overlap
0127 Olap = find(pl, 'Olap');
0128 if ~isempty(Olap)
0129   if Olap < 0
0130     % use window function rov
0131     Olap = Win.rov;
0132     disp(sprintf('! Using recommended overlap of %2.1f%%', Olap))
0133     pl = pset(pl, 'Olap', Olap);
0134   end
0135 end
0136 
0137 % desired detrending order
0138 order = find(pl, 'Order');
0139 if isempty(order)
0140   order = 0;
0141   warning('! using default detrending order 0 (mean)');
0142   pl = append(pl, 'Order', order);
0143 end
0144 
0145 % Loop over input AOs
0146 bs = cell(na);
0147 for j=1:na
0148   aj = as(j);
0149   for k=1:na
0150     ak = as(k);
0151     % -------- Make Xspec estimate
0152 
0153     % Do detrending
0154     x = aj.data.y;
0155     y = ak.data.y;
0156     switch order
0157       case -1
0158         % do nothing
0159       case 0
0160         x = x - mean(x);
0161         y = y - mean(y);
0162       case 1
0163         x = detrend(x);
0164         y = detrend(y);
0165       otherwise
0166         y = polydetrend(y, order);
0167         y = polydetrend(y, order);
0168     end
0169 
0170     switch method
0171       case 'TF'
0172         % Compute TF using tfestimate
0173         [txy,f] = tfestimate(x,y,Win.win,round(Olap/100*Nfft),Nfft,fsmax);
0174       case 'COHERE'
0175         % Compute coherence using mscohere
0176         [txy,f] = mscohere(x,y,Win.win,round(Olap/100*Nfft),Nfft,fsmax);
0177       case 'CPSD'
0178         % Compute Cross-spectral density using CPSD
0179         [txy,f] = cpsd(x,y,Win.win,round(Olap/100*Nfft),Nfft,fsmax);
0180       otherwise
0181         error('### Unknown method.')
0182     end
0183 
0184     % create new output fsdata
0185     fs = fsdata(f, txy, fsmax);
0186     fs = set(fs, 'name', sprintf('%s(%s->%s)', method, aj.data.name, ak.data.name));
0187     fs = set(fs, 'enbw', Win.nenbw*fsmax/Nfft);
0188 
0189     %----------- create new output history
0190 
0191     % we need to get the input histories in the same order as the inputs
0192     % to this function call, not in the order of the input to tfestimate;
0193     % otherwise the resulting matrix on a 'create from history' will be
0194     % mirrored.
0195     if j<k
0196       h = history(ALGONAME, VERSION, pl, [aj.hist ak.hist]);
0197     else
0198       h = history(ALGONAME, VERSION, pl, [ak.hist aj.hist]);
0199     end
0200     if na == length(invars)
0201       h = set(h, 'invars', [invars(j) invars(k)]);
0202       nk = invars{k};
0203       nj = invars{j};
0204     else
0205       nk = ak.name;
0206       nj = aj.name;
0207     end
0208 
0209     % make output analysis object
0210     b = ao(fs, h);
0211 
0212     % set name
0213     b = setnh(b, 'name', sprintf('%s(%s->%s)', method, nj, nk));
0214 
0215     % add to outputs
0216     bs(j,k) = {b};
0217   end % End second loop over AOs
0218 end % End first loop over AOs
0219 
0220 % now we have a cell matrix of AOs but
0221 % we want a normal matrix
0222 b  = [bs{:}];
0223 varargout{1} = reshape(b, na, na);
0224 
0225 
0226 %--------------------------------------------------------------------------
0227 % Returns the length of the shortest vector in samples
0228 function lmin = findShortestVector(as)
0229 
0230 lmin = 1e20;
0231 for j=1:length(as)
0232   if len(as(j)) < lmin
0233     lmin = len(as(j));
0234   end
0235 end
0236 
0237 
0238 %--------------------------------------------------------------------------
0239 % Returns the max Fs of a set of AOs
0240 function fs = findFsMax(as)
0241 
0242 fs = 0;
0243 for j=1:length(as)
0244   a = as(j);
0245   if a.data.fs > fs
0246     fs = a.data.fs;
0247   end
0248 end
0249 
0250 
0251

Generated on Mon 31-Mar-2008 13:54:54 by m2html © 2003