0001 function varargout = ltpda_lcohere(varargin)
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
0022
0023
0024
0025
0026
0027
0028
0029
0030
0031
0032
0033
0034
0035
0036
0037
0038
0039
0040
0041
0042
0043
0044
0045
0046
0047
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
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
0074 bs = [];
0075
0076
0077 if isempty(ps)
0078 pl = getDefaultPL();
0079 else
0080 pl = combine(ps, getDefaultPL);
0081 end
0082
0083
0084 na = length(as);
0085
0086
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
0095 if ~isa(a.data, 'tsdata')
0096 error('### ltpda_lcohere requires tsdata (time-series) inputs.');
0097 end
0098
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
0110
0111
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
0131 N = len(as(1));
0132 iS = zeros(N,na);
0133 for j=1:na
0134 a = as(j);
0135 iS(:,j) = a.data.x;
0136 end
0137
0138
0139
0140
0141 Kdes = find(pl, 'Kdes');
0142 Kmin = find(pl, 'Kmin');
0143 Jdes = find(pl, 'Jdes');
0144 Win = find(pl, 'Win');
0145 Nolap = find(pl, 'Nolap');
0146 if Nolap == -1
0147
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');
0153
0154
0155 debug = find(pl, 'DEBUG');
0156
0157
0158 [f,r,m,L] = ltpda_compute_f(fsmax, lmin, Kdes, Kmin, Jdes, Nolap, debug);
0159
0160
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
0175 disp('** Building output AO matrix')
0176 bs = cell(na);
0177 sig = cell(na);
0178 for j=1:na
0179 bi = as(j);
0180 for k=1:na
0181 bo = as(k);
0182
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
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
0198 b = ao(fs, h);
0199
0200 b = set(b, 'name', sprintf('LCOHERE(%s->%s)', nj, nk));
0201
0202 bs(j,k) = {b};
0203
0204
0205 if nargout == 2
0206
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
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
0217 b = ao(fs, h);
0218
0219 b = set(b, 'name', ['\sigma^2' sprintf('(%s->%s)', nj, nk)]);
0220
0221 sig(j,k) = {b};
0222
0223 end
0224 end
0225 end
0226
0227
0228
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
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
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
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
0278 function varargout = mlcohere(varargin)
0279
0280
0281
0282
0283
0284
0285
0286
0287
0288
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
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
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
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
0359
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
0374
0375
0376 for j=1:nc
0377
0378 for k=1:nc
0379
0380 if j<k
0381
0382
0383 Sxx = 0;
0384 Syy = 0;
0385 Sxy = 0;
0386 i2(:) = 0;
0387
0388
0389 for s=1:nsegs
0390
0391 xi(:) = X(idxs(s,:),j);
0392 xo(:) = X(idxs(s,:),k);
0393
0394
0395 eval(detrender);
0396
0397
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
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
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
0428
0429 end
0430 end
0431 end
0432 end
0433
0434
0435 aoOnes = ones(1,nf);
0436 for j=1:nc
0437 for k=1:nc
0438 if j == k
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