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. 
 
 The function is a helper function for various higher level functions. It
 is meant to be called from other functions (e.g., ltpda_tfe).
 
 Usage: 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
 
 M Hewitson 30-05-07

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

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

Generated on Mon 02-Jul-2007 12:19:41 by m2html © 2003