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 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
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(na, N);
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 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');
0150 Nolap = find(pl, 'Nolap');
0151 if Nolap == -1
0152
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');
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
0166 [f,r,m,L] = ltpda_compute_f(fsmax, lmin, Kdes, Kmin, Jdes, Nolap);
0167
0168
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
0183 disp('** Building output AO matrix')
0184 bs = cell(na);
0185 sig = cell(na);
0186 for j=1:na
0187 bi = as(j);
0188 for k=1:na
0189 bo = as(k);
0190
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
0196 h = history(ALGONAME, VERSION, pl, [bi.hist bo.hist]);
0197 h = set(h, 'invars', invars([j k]));
0198
0199 b = ao(fs, h);
0200
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
0213 bs(j,k) = {b};
0214
0215
0216 if nargout == 2
0217
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
0223 h = history(ALGONAME, VERSION, pl, [bi.hist bo.hist]);
0224 h = set(h, 'invars', invars{[j k]});
0225
0226 b = ao(fs, h);
0227
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
0240 sig(j,k) = {b};
0241
0242 end
0243 end
0244 end
0245
0246
0247
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
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
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
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
0297 function varargout = mlcohere(varargin)
0298
0299
0300
0301
0302
0303
0304
0305
0306
0307
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
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
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
0357 for j=1:nc
0358
0359 for k=1:nc
0360
0361
0362 start = 1;
0363 nsegs = 0;
0364 while start + l - 1 <= nx
0365
0366 start = floor(start + l*Xolap);
0367 nsegs = nsegs + 1;
0368 end
0369
0370
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
0382 idx = start:start+l-1;
0383 xi = X(j,idx).';
0384 xo = X(k,idx).';
0385
0386
0387 switch order
0388 case -1
0389
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
0402 ti = C*xi;
0403 to = C*xo;
0404
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
0411 start = floor(start + l*Xolap);
0412 nsegs = nsegs + 1;
0413 p = p + 1;
0414 end
0415
0416
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
0426 if computeVariance
0427 A = 0;
0428 start = 1;
0429 p = 1;
0430 scale = fs*S2;
0431 while start + l - 1 <= nx
0432
0433 S = (T(p)/scale-Txy(j,k,fi))^2;
0434 A = A + S * i2(p);
0435
0436 start = floor(start + l*Xolap);
0437 p = p + 1;
0438 end
0439 Sig(j,k,fi) = A/sw;
0440 end
0441
0442 end
0443 end
0444 end
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