Home > m > sigproc > frequency_domain > ltpda_lcohere.m

ltpda_lcohere

PURPOSE ^

LTPDA_LCOHERE implement weighted coherence estimation computed on a log frequency axis.

SYNOPSIS ^

function varargout = ltpda_lcohere(varargin)

DESCRIPTION ^

 LTPDA_LCOHERE implement weighted coherence estimation computed on a log frequency axis.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

 LTPDA_LCOHERE implement weighted cross-power estimation computed on a 
 log frequency axis.
 
 Usage:
   >> txy = ltpda_lcohere(as)
   >> txy = ltpda_lcohere(as, pl)
   >> [txy, sig2] = ltpda_lcohere(as, pl)  % returns variance estimates
 
   Inputs:
     as  - array of analysis objects
     pl  - parameter list (see below)
 
   Outputs:
     bs  - array of analysis objects, one for each input
 
   Parameter list:
     Kdes  - desired number of averages   [default: 100]
     Kmin  - minimum number of averages   [default: 10]
     Jdes  - number of spectral frequencies to compute [default: 1000]
     Win   - a specwin window object
             Only the design parameters of the window object are used; the
             window is recomputed for each DFT length inside the lpsd_core
             algorithm. [default: Kaiser with -200dB PSLL]
     Nolap - desired overlap [default: taken from window]
     Order - order of detrending.
             -1 - no detrending
              0 - subtract mean
              1 - subtract linear fit
              N - subtract fit of polynomial, order N  
 
 The following call returns a parameter list object that contains the
 default parameter values:
 
 >> pl = ltpda_lcohere('Params')
 
 
 M Hewitson 02-02-07
 
 $Id: ltpda_ltfe.m,v 1.5 2007/06/12 12:32:14 hewitson Exp $
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

0001 function varargout = ltpda_lcohere(varargin)
0002 % LTPDA_LCOHERE implement weighted coherence estimation computed on a log frequency axis.
0003 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0004 %
0005 % LTPDA_LCOHERE implement weighted cross-power estimation computed on a
0006 % log frequency axis.
0007 %
0008 % Usage:
0009 %   >> txy = ltpda_lcohere(as)
0010 %   >> txy = ltpda_lcohere(as, pl)
0011 %   >> [txy, sig2] = ltpda_lcohere(as, pl)  % returns variance estimates
0012 %
0013 %   Inputs:
0014 %     as  - array of analysis objects
0015 %     pl  - parameter list (see below)
0016 %
0017 %   Outputs:
0018 %     bs  - array of analysis objects, one for each input
0019 %
0020 %   Parameter list:
0021 %     Kdes  - desired number of averages   [default: 100]
0022 %     Kmin  - minimum number of averages   [default: 10]
0023 %     Jdes  - number of spectral frequencies to compute [default: 1000]
0024 %     Win   - a specwin window object
0025 %             Only the design parameters of the window object are used; the
0026 %             window is recomputed for each DFT length inside the lpsd_core
0027 %             algorithm. [default: Kaiser with -200dB PSLL]
0028 %     Nolap - desired overlap [default: taken from window]
0029 %     Order - order of detrending.
0030 %             -1 - no detrending
0031 %              0 - subtract mean
0032 %              1 - subtract linear fit
0033 %              N - subtract fit of polynomial, order N
0034 %
0035 % The following call returns a parameter list object that contains the
0036 % default parameter values:
0037 %
0038 % >> pl = ltpda_lcohere('Params')
0039 %
0040 %
0041 % M Hewitson 02-02-07
0042 %
0043 % $Id: ltpda_ltfe.m,v 1.5 2007/06/12 12:32:14 hewitson Exp $
0044 %
0045 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0046 
0047 % Check if this is a call for parameters
0048 if nargin == 1
0049   in = char(varargin{1});
0050   if strcmp(in, 'Params')
0051     varargout{1} = getDefaultPL();
0052     return
0053   end
0054 end
0055 
0056 % capture input variable names
0057 invars = {};
0058 as     = [];
0059 ps     = [];
0060 for j=1:nargin
0061   if isa(varargin{j}, 'ao')
0062     as = [as varargin{j}];
0063     invars = [invars cellstr(inputname(j))];
0064   end
0065   if isa(varargin{j}, 'plist')
0066     ps = [ps varargin{j}];
0067   end
0068 end
0069 
0070 ALGONAME = mfilename;
0071 VERSION  = '$Id: ltpda_ltfe.m,v 1.5 2007/06/12 12:32:14 hewitson Exp $';
0072 
0073 % initialise output array
0074 bs = []; 
0075 
0076 % check plist
0077 if isempty(ps)
0078   pl = getDefaultPL();
0079 else
0080   pl = combine(ps, getDefaultPL);
0081 end
0082 
0083 % Go through each input AO and build a matrix of signals
0084 na    = length(as);
0085 
0086 %----------------- Resample all AOs
0087 
0088 disp('*** Resampling AOs...');
0089 fsmax = findFsMax(as);
0090 fspl  = plist(param('fsout', fsmax));
0091 bs = [];
0092 for i=1:na
0093   a = as(i);
0094   % check this is a time-series object
0095   if ~isa(a.data, 'tsdata')
0096     error('### ltpda_lcohere requires tsdata (time-series) inputs.');
0097   end  
0098   % Check Fs
0099   if a.data.fs ~= fsmax
0100     warning(sprintf('!!! Resampling AO %s/%s to %f Hz', a.name, a.data.name, fsmax));
0101     a = resample(a, fspl);
0102   end
0103   bs = [bs a];
0104 end
0105 as = bs;
0106 clear bs;
0107 disp('*** Done.');
0108 
0109 %----------------- Truncate all vectors
0110 
0111 % Get shortest vector
0112 disp('*** Truncating all vectors...');
0113 lmin = findShortestVector(as);
0114 nsecs = lmin / fsmax;
0115 bs = [];
0116 for i=1:na
0117   a = as(i);
0118   if len(a) ~= lmin
0119     warning(sprintf('!!! Truncating AO %s/%s to %d secs', a.name, a.data.name, nsecs));
0120     bs = [bs  select(a, 1:lmin)];
0121   else
0122     bs = [bs a];
0123   end
0124 end
0125 as = bs;
0126 clear bs;
0127 disp('*** Done.');
0128 
0129 
0130 %----------------- Build signal Matrix
0131 N     = len(as(1)); % length of first signal
0132 iS    = zeros(N,na);
0133 for j=1:na
0134   a = as(j);
0135   iS(:,j) = a.data.x;
0136 end
0137 
0138 %----------------- check input parameters
0139 
0140 % Averaging
0141 Kdes = find(pl, 'Kdes');   % desired averages
0142 Kmin = find(pl, 'Kmin');   % minimum averages
0143 Jdes = find(pl, 'Jdes');   % num desired spectral frequencies
0144 Win   = find(pl, 'Win');   % Window object
0145 Nolap = find(pl, 'Nolap');    % desired overlap
0146 if Nolap == -1
0147   % use window function rov
0148   Nolap = Win.rov/100;
0149   disp(sprintf('! Using recommended overlap of %d', Nolap))
0150   pl = set(pl, 'Nolap', Nolap);
0151 end
0152 Order = find(pl, 'Order');    % desired overlap
0153 
0154 % Debug
0155 debug = find(pl, 'DEBUG');
0156 
0157 %----------------- Get frequency vector
0158 [f,r,m,L] = ltpda_compute_f(fsmax, lmin, Kdes, Kmin, Jdes, Nolap, debug);
0159   
0160 %----------------- compute CPSD Estimates
0161 
0162 if nargout == 2
0163   disp('** Computing Coherence matrix with variance estimates')
0164   tic
0165   [Txy, ENBW, Sig] = mlcohere(iS, f, r, m, L, fsmax, Win, Order);
0166   toc
0167 else
0168   disp('** Computing Coherence matrix')
0169   tic
0170   [Txy, ENBW]  = mlcohere(iS, f, r, m, L, fsmax, Win, Order);
0171   toc
0172 end
0173   
0174 %----------------- Build output Matrix of AOs
0175 disp('** Building output AO matrix')
0176 bs  = cell(na);
0177 sig = cell(na);
0178 for j=1:na  % Loop over input signal
0179   bi = as(j);
0180   for k=1:na  % Loop over output signal
0181     bo = as(k);   
0182     % create new output fsdata
0183     fs = fsdata(f, squeeze(Txy(j,k,:)), fsmax);
0184     fs = set(fs, 'name', sprintf('LCOHERE(%s->%s)', bi.data.name, bo.data.name));
0185     fs = set(fs, 'yunits', [bo.data.yunits '/' bi.data.yunits]);
0186     fs = set(fs, 'enbw', squeeze(ENBW(j,k,:)));
0187     % create new output history
0188     h = history(ALGONAME, VERSION, pl, [bi.hist bo.hist]);
0189     if na == length(invars)
0190       h = set(h, 'invars', invars([j k]));
0191       nk = invars{k};
0192       nj = invars{j};
0193     else
0194       nk = bo.name;
0195       nj = bi.name;
0196     end
0197     % make output analysis object
0198     b = ao(fs, h);    
0199     % set name
0200     b = set(b, 'name', sprintf('LCOHERE(%s->%s)', nj, nk));    
0201     % add to outputs
0202     bs(j,k) = {b};
0203     
0204     % Make variance outputs
0205     if nargout == 2
0206       % create new output fsdata
0207       fs = fsdata(f, squeeze(Sig(j,k,:)), fsmax);
0208       fs = set(fs, 'name', ['\sigma^2' sprintf('(%s->%s)', bi.data.name, bo.data.name)]);
0209       fs = set(fs, 'yunits', [bo.data.yunits '/' bi.data.yunits]);
0210       fs = set(fs, 'enbw', squeeze(ENBW(j,k,:)));
0211       % create new output history
0212       h = history(ALGONAME, VERSION, pl, [bi.hist bo.hist]);
0213       if na == length(invars)
0214         h = set(h, 'invars', invars{[j k]});
0215       end
0216       % make output analysis object
0217       b = ao(fs, h);
0218       % set name
0219       b = set(b, 'name', ['\sigma^2' sprintf('(%s->%s)', nj, nk)]);
0220       % add to outputs
0221       sig(j,k) = {b};
0222     
0223     end
0224   end % End loop over output signal
0225 end % End loop over input signal
0226 
0227 % now we have a cell matrix of AOs but
0228 % we want a normal matrix
0229 b  = [bs{:}];
0230 bs = reshape(b, na, na);
0231 varargout{1} = bs;
0232 if nargout == 2
0233   b  = [sig{:}];
0234   s  = reshape(b, na, na);
0235   varargout{2} = s;
0236 end
0237 %--------------------------------------------------------------------------
0238 % Returns the length of the shortest vector in samples
0239 function lmin = findShortestVector(as)
0240 
0241 lmin = 1e20;
0242 for j=1:length(as)
0243   if len(as(j)) < lmin
0244     lmin = len(as(j));
0245   end
0246 end
0247 
0248 
0249 %--------------------------------------------------------------------------
0250 % Returns the max Fs of a set of AOs
0251 function fs = findFsMax(as)
0252 
0253 fs = 0;
0254 for j=1:length(as)
0255   a = as(j);
0256   if a.data.fs > fs
0257     fs = a.data.fs;
0258   end
0259 end
0260 
0261 %--------------------------------------------------------------------------
0262 % Get default params
0263 function plo = getDefaultPL()
0264 
0265 disp('* creating default plist...');
0266 plo = plist();
0267 plo = append(plo, param('Kdes', 100));
0268 plo = append(plo, param('Jdes', 1000));
0269 plo = append(plo, param('Kmin', 2));
0270 plo = append(plo, param('Win',  specwin('Kaiser', 1, 200)));
0271 plo = append(plo, param('Nolap', -1));
0272 plo = append(plo, param('Order', 0));
0273 disp('* done.');
0274 
0275 
0276 %--------------------------------------------------------------------------
0277 % Compute all TF combinations
0278 function varargout = mlcohere(varargin)
0279 
0280 % MLCOHERE
0281 %
0282 % Txy = mlcohere(X,f,r,m,L,fs,win)
0283 %
0284 % M Hewitson 19-05-07
0285 %
0286 %
0287 
0288 % Get inputs
0289 X     = varargin{1};
0290 f     = varargin{2};
0291 r     = varargin{3};
0292 m     = varargin{4};
0293 L     = varargin{5};
0294 fs    = varargin{6};
0295 win   = varargin{7};
0296 order = varargin{8};
0297 
0298 % do we need to compute variance?
0299 if nargout == 2
0300   computeVariance = 1;
0301 else
0302   computeVariance = 0;
0303 end
0304 
0305 twopi    = 2.0*pi;
0306 twopiofs = 2.0*pi/(1.0*fs);
0307 
0308 si   = size(X);
0309 nx   = si(1);
0310 nc   = si(2);
0311 nf   = length(f);
0312 Txy  = zeros(nc,nc,nf);
0313 Sig  = zeros(nc,nc,nf);
0314 ENBW = zeros(nc,nc,nf);
0315 
0316 
0317 switch order
0318   case -1
0319     detrender = '';
0320   case 0
0321     detrender = 'xi = xi - mean(xi);xo = xo - mean(xo);';
0322   case 1
0323     detrender = 'xi = detrend(xi);xo = detrend(xo);';
0324   otherwise
0325     detrender = 'xi = polydetrend(xi, order); xo = polydetrend(xo, order);';
0326 end
0327 
0328 disp_each = round(nf/100)*10;
0329 
0330 for fi=int32(1:nf)
0331   
0332   % compute DFT exponent and window
0333   l = L(fi);
0334   switch win.name
0335     case 'Kaiser'       
0336       win = specwin(win.name, l, win.psll);
0337     otherwise
0338       win = specwin(win.name, l);
0339   end    
0340       
0341   p     = 1i * twopi * m(fi)/l.*[0:l-1];
0342   C     = win.win .* exp(p);
0343   Xolap = (1-win.rov/100);
0344   
0345   if mod(fi,disp_each) == 0 || fi == 1 || fi == nf
0346     disp(sprintf('++ computing frequency %03d of %03d: %f Hz', fi, nf, f(fi)));
0347   end
0348   
0349   %-------------- Count segments
0350   nsegs = 0;
0351   start = 1;
0352   lm1   = l-1;
0353   while start + l - 1 <= nx
0354     start = floor(start + l*Xolap);
0355     nsegs = nsegs + 1;
0356   end
0357 
0358   % now allocate indices for segments
0359 %   idxs  = repmat(struct('idx', zeros(1,l)), 1, nsegs);
0360   idxs  = zeros(nsegs,l, 'uint32');
0361   nsegs = 0;
0362   start = 1;
0363   while start + l - 1 <= nx
0364     idxs(1+nsegs,:) = start:start+lm1;
0365     start = floor(start + l*Xolap);
0366     nsegs = nsegs + 1;
0367   end
0368   Sxy = zeros(nsegs,1);
0369   i2  = zeros(nsegs,1);
0370   xi  = zeros(l,1);
0371   xo  = zeros(l,1);
0372   
0373   %------------------ Process the channels
0374   
0375   % Loop over input channels
0376   for j=1:nc        
0377     % loop over output channels
0378     for k=1:nc      
0379       % compute lower triangle of matrix
0380       if j<k
0381         
0382         % do segments
0383         Sxx    = 0;
0384         Syy    = 0;      
0385         Sxy    = 0;
0386         i2(:)  = 0;
0387         
0388         % Loop over segments
0389         for s=1:nsegs           
0390           % get segments
0391           xi(:) = X(idxs(s,:),j);
0392           xo(:) = X(idxs(s,:),k);
0393 
0394           % detrend segments
0395           eval(detrender);
0396 
0397           % compute DFT
0398           ti      = C*xi;
0399           to      = C*xo;
0400           cto     = conj(to);
0401           Sxy     = Sxy + ti.*cto;
0402           Sxx     = Sxx + ti.*conj(ti);
0403           Syy     = Syy + to.*cto;
0404 
0405         end % End segment loop
0406         disp(['-processed ' num2str(nsegs) ' segments']);
0407         
0408         Sxy = Sxy / nsegs;
0409         Sxx = Sxx / nsegs;
0410         Syy = Syy / nsegs;
0411         
0412         S1      = win.ws;
0413         S12     = S1*S1;
0414         S2      = win.ws2;
0415         ENBW(j,k,fi) = fs*S2/S12;
0416         Txy(j,k,fi)  = (abs(Sxy)^2)/(Sxx*Syy);
0417 
0418         %----------- compute variance
0419         if computeVariance
0420           A     = 0;
0421           scale = fs*S2;
0422           for s=1:nsegs
0423             S  = (Sxy(s)/scale-Txy(j,k,fi))^2;
0424             A  = A + S * i2(s);
0425           end
0426           Sig(j,k,fi) = A/sum(i2);
0427         end % End compute variance
0428         
0429       end % End if j<k
0430     end % End loop over output channels
0431   end % End loop over input channels
0432 end % End loop over frequencies
0433 
0434 % Fill the holes
0435 aoOnes = ones(1,nf);
0436 for j=1:nc
0437   for k=1:nc
0438     if j == k % the answer is 1
0439       Txy(j,k,:) = aoOnes;
0440     end
0441     if j>k
0442       Txy(j,k,:) = Txy(k,j,:);
0443       if computeVariance
0444         Sig(j,k,:) = Sig(k,j,:);
0445       end
0446     end
0447   end
0448 end
0449 
0450 
0451 varargout{1} = Txy;
0452 if nargout == 2
0453   varargout{2} = ENBW;
0454 elseif nargout == 3
0455   varargout{3} = Sig;
0456   varargout{2} = ENBW;
0457 end
0458 
0459 
0460 % END

Generated on Fri 02-Nov-2007 19:39:27 by m2html © 2003