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: fs/4]
     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: fs/4]
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   invars = [invars cellstr(inputname(j))];
0062   if isa(varargin{j}, 'ao')
0063     as = [as varargin{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(na, N);
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 if Jdes == -1
0145   Jdes = round(fsmax/4);
0146   disp(sprintf('! Using default Jdes of %d', Jdes))
0147   pl = set(pl, 'Jdes', Jdes);
0148 end
0149 Win   = find(pl, 'Win');   % Window object
0150 Nolap = find(pl, 'Nolap');    % desired overlap
0151 if Nolap == -1
0152   % use window function rov
0153   Nolap = Win.rov/100;
0154   disp(sprintf('! Using recommended overlap of %d', Nolap))
0155   pl = set(pl, 'Nolap', Nolap);
0156 end
0157 Order = find(pl, 'Order');    % desired overlap
0158 if isempty(Order)
0159   Order = 0;
0160   disp(sprintf('! Using default detrending order of %d', Order))
0161   pl = append(pl, param('Order', Order));
0162 end
0163 
0164 
0165 %----------------- Get frequency vector
0166 [f,r,m,L] = ltpda_compute_f(fsmax, lmin, Kdes, Kmin, Jdes, Nolap);
0167   
0168 %----------------- compute CPSD Estimates
0169 
0170 if nargout == 2
0171   disp('** Computing Coherence matrix with variance estimates')
0172   tic
0173   [Txy, ENBW, Sig] = mlcohere(iS, f, r, m, L, fsmax, Win, Order);
0174   toc
0175 else
0176   disp('** Computing Coherence matrix')
0177   tic
0178   [Txy, ENBW]  = mlcohere(iS, f, r, m, L, fsmax, Win, Order);
0179   toc
0180 end
0181   
0182 %----------------- Build output Matrix of AOs
0183 disp('** Building output AO matrix')
0184 bs  = cell(na);
0185 sig = cell(na);
0186 for j=1:na  % Loop over input signal
0187   bi = as(j);
0188   for k=1:na  % Loop over output signal
0189     bo = as(k);   
0190     % create new output fsdata
0191     fs = fsdata(f, squeeze(Txy(j,k,:)), fsmax);
0192     fs = set(fs, 'name', sprintf('LCOHERE(%s->%s)', bi.data.name, bo.data.name));
0193     fs = set(fs, 'yunits', [bo.data.yunits '/' bi.data.yunits]);
0194     fs = set(fs, 'enbw', squeeze(ENBW(j,k,:)));
0195     % create new output history
0196     h = history(ALGONAME, VERSION, pl, [bi.hist bo.hist]);
0197     h = set(h, 'invars', invars([j k]));
0198     % make output analysis object
0199     b = ao(fs, h);
0200     % set name
0201     if isempty(invars{k})
0202       nk = bo.name;
0203     else
0204       nk = invars{k};
0205     end
0206     if isempty(invars{j})
0207       nj = bi.name;
0208     else
0209       nj = invars{j};
0210     end
0211     b = set(b, 'name', sprintf('LCOHERE(%s->%s)', nj, nk));    
0212     % add to outputs
0213     bs(j,k) = {b};
0214     
0215     % Make variance outputs
0216     if nargout == 2
0217       % create new output fsdata
0218       fs = fsdata(f, squeeze(Sig(j,k,:)), fsmax);
0219       fs = set(fs, 'name', ['\sigma^2' sprintf('(%s->%s)', bi.data.name, bo.data.name)]);
0220       fs = set(fs, 'yunits', [bo.data.yunits '/' bi.data.yunits]);
0221       fs = set(fs, 'enbw', squeeze(ENBW(j,k,:)));
0222       % create new output history
0223       h = history(ALGONAME, VERSION, pl, [bi.hist bo.hist]);
0224       h = set(h, 'invars', invars{[j k]});
0225       % make output analysis object
0226       b = ao(fs, h);
0227       % set name
0228       if isempty(invars{k})
0229         nk = bo.name;
0230       else
0231         nk = invars{k};
0232       end
0233       if isempty(invars{j})
0234         nj = bi.name;
0235       else
0236         nj = invars{j};
0237       end
0238       b = set(b, 'name', ['\sigma^2' sprintf('(%s->%s)', nj, nk)]);
0239       % add to outputs
0240       sig(j,k) = {b};
0241     
0242     end
0243   end % End loop over output signal
0244 end % End loop over input signal
0245 
0246 % now we have a cell matrix of AOs but
0247 % we want a normal matrix
0248 b  = [bs{:}];
0249 bs = reshape(b, na, na);
0250 varargout{1} = bs;
0251 if nargout == 2
0252   b  = [sig{:}];
0253   s  = reshape(b, na, na);
0254   varargout{2} = s;
0255 end
0256 %--------------------------------------------------------------------------
0257 % Returns the length of the shortest vector in samples
0258 function lmin = findShortestVector(as)
0259 
0260 lmin = 1e20;
0261 for j=1:length(as)
0262   if len(as(j)) < lmin
0263     lmin = len(as(j));
0264   end
0265 end
0266 
0267 
0268 %--------------------------------------------------------------------------
0269 % Returns the max Fs of a set of AOs
0270 function fs = findFsMax(as)
0271 
0272 fs = 0;
0273 for j=1:length(as)
0274   a = as(j);
0275   if a.data.fs > fs
0276     fs = a.data.fs;
0277   end
0278 end
0279 
0280 %--------------------------------------------------------------------------
0281 % Get default params
0282 function plo = getDefaultPL()
0283 
0284 disp('* creating default plist...');
0285 plo = plist();
0286 plo = append(plo, param('Kdes', 100));
0287 plo = append(plo, param('Jdes', 1000));
0288 plo = append(plo, param('Kmin', 2));
0289 plo = append(plo, param('Win',  specwin('Kaiser', 1, 150)));
0290 plo = append(plo, param('Nolap', -1));
0291 plo = append(plo, param('Order', 0));
0292 disp('* done.');
0293 
0294 
0295 %--------------------------------------------------------------------------
0296 % Compute all TF combinations
0297 function varargout = mlcohere(varargin)
0298 
0299 % MLCOHERE
0300 %
0301 % Txy = mlcohere(X,f,r,m,L,fs,win)
0302 %
0303 % M Hewitson 19-05-07
0304 %
0305 %
0306 
0307 % Get inputs
0308 X     = varargin{1};
0309 f     = varargin{2};
0310 r     = varargin{3};
0311 m     = varargin{4};
0312 L     = varargin{5};
0313 fs    = varargin{6};
0314 win   = varargin{7};
0315 order = varargin{8};
0316 
0317 % do we need to compute variance?
0318 if nargout == 2
0319   computeVariance = 1;
0320 else
0321   computeVariance = 0;
0322 end
0323 
0324 twopi    = 2.0*pi;
0325 twopiofs = 2.0*pi/(1.0*fs);
0326 
0327 si   = size(X);
0328 nx   = si(2);
0329 nc   = si(1);
0330 nf   = length(f);
0331 Txy  = zeros(nc,nc,nf);
0332 Sig  = zeros(nc,nc,nf);
0333 ENBW = zeros(nc,nc,nf);
0334 
0335 disp_each = round(nf/100)*10;
0336 
0337 for fi=1:nf
0338   
0339   % compute DFT exponent and window
0340   l = L(fi);
0341   switch win.name
0342     case 'Kaiser'       
0343       win = specwin(win.name, l, win.psll);
0344     otherwise
0345       win = specwin(win.name, l);
0346   end    
0347       
0348   p     = 1i * twopi * m(fi)/l.*[0:l-1];
0349   C     = win.win .* exp(p);
0350   Xolap = (1-win.rov/100);
0351   
0352   if mod(fi,disp_each) == 0 || fi == 1 || fi == nf
0353     disp(sprintf('++ computing frequency %d of %d: %f Hz', fi, nf, f(fi)));
0354   end
0355   
0356   % Loop over input channels
0357   for j=1:nc
0358     % loop over output channels
0359     for k=1:nc      
0360       
0361       % Count N segments
0362       start = 1;
0363       nsegs = 0;
0364       while start + l - 1 <= nx
0365         % make next step
0366         start = floor(start + l*Xolap);
0367         nsegs = nsegs + 1;        
0368       end
0369       
0370       % do segments
0371       Axy = 0.0;
0372       Px  = 0.0;
0373       Py  = 0.0;
0374       T   = zeros(nsegs,1);
0375       i2  = zeros(nsegs,1);
0376       start = 1;
0377       nsegs = 0;
0378       p = 1;
0379       while start + l - 1 <= nx
0380 
0381         % get segments
0382         idx = start:start+l-1;
0383         xi    = X(j,idx).';
0384         xo    = X(k,idx).';
0385         
0386         % detrend segments
0387         switch order
0388           case -1
0389             % do nothing
0390           case 0
0391             xi = xi - mean(xi);
0392             xo = xo - mean(xo);
0393           case 1
0394             xi = detrend(xi);
0395             xo = detrend(xo);
0396           otherwise
0397             xi = polydetrend(xi, order);
0398             xo = polydetrend(xo, order);
0399         end
0400 
0401         % compute DFT
0402         ti    = C*xi;
0403         to    = C*xo;
0404 %         i2(p) = ti*conj(ti);
0405         T(p)  = to .* conj(ti);
0406         Axy   = Axy + T(p);
0407         Px    = Px + ti.*conj(ti);
0408         Py    = Py + to.*conj(to);
0409         
0410         % make next step
0411         start = floor(start + l*Xolap);
0412         nsegs = nsegs + 1;        
0413         p     = p + 1;
0414       end
0415 
0416 %       A2ns    = 2.0*A/nsegs;
0417       S1      = win.ws;
0418       S12     = S1*S1;
0419       S2      = win.ws2;
0420       ENBW(j,k,fi) = fs*S2/S12;
0421       
0422       sw = sum(i2);
0423       Txy(j,k,fi)  = abs(Axy).^2/Px/Py;
0424       
0425       % compute variance
0426       if computeVariance
0427         A     = 0;
0428         start = 1;
0429         p     = 1;
0430         scale = fs*S2;
0431         while start + l - 1 <= nx
0432 %         for p=1:ns
0433           S  = (T(p)/scale-Txy(j,k,fi))^2;
0434           A  = A + S * i2(p);
0435           % make next step
0436           start = floor(start + l*Xolap);
0437           p     = p + 1;
0438         end
0439         Sig(j,k,fi) = A/sw;        
0440       end
0441       
0442     end % End loop over output channels
0443   end % End loop over input channels
0444 end % End loop over frequencies
0445 
0446 varargout{1} = Txy;
0447 if nargout == 2
0448   varargout{2} = ENBW;
0449 elseif nargout == 3
0450   varargout{3} = Sig;
0451   varargout{2} = ENBW;
0452 end
0453 
0454 
0455 % END

Generated on Mon 03-Sep-2007 12:12:34 by m2html © 2003