0001 function varargout = ltpda_ltfe(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.9 2007/09/01 10:18:06 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_ltfe 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 TF matrix with variance estimates')
0172 tic
0173 [Txy, ENBW, Sig] = mltfe(iS, f, r, m, L, fsmax, Win, Order);
0174 toc
0175 else
0176 disp('** Computing TF matrix')
0177 tic
0178 [Txy, ENBW] = mltfe(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('LTFe(%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 try
0198 h = set(h, 'invars', invars([j k]));
0199 catch
0200 warning('!!! Unable to set input variable names for history. Perhaps you input a vector of AOs.');
0201 end
0202
0203 b = ao(fs, h);
0204
0205 try
0206 if isempty(invars{k})
0207 nk = bo.name;
0208 else
0209 nk = invars{k};
0210 end
0211 if isempty(invars{j})
0212 nj = bi.name;
0213 else
0214 nj = invars{j};
0215 end
0216 catch
0217 warning('!!! Unable to set input variable names for history. Perhaps you input a vector of AOs.');
0218 nk = bo.name;
0219 nj = bi.name;
0220 end
0221
0222 b = set(b, 'name', sprintf('LTFe(%s->%s)', nj, nk));
0223
0224 bs(j,k) = {b};
0225
0226
0227 if nargout == 2
0228
0229 fs = fsdata(f, squeeze(Sig(j,k,:)), fsmax);
0230 fs = set(fs, 'name', ['\sigma^2' sprintf('(%s->%s)', bi.data.name, bo.data.name)]);
0231 fs = set(fs, 'yunits', [bo.data.yunits '/' bi.data.yunits]);
0232 fs = set(fs, 'enbw', squeeze(ENBW(j,k,:)));
0233
0234 h = history(ALGONAME, VERSION, pl, [bi.hist bo.hist]);
0235 try
0236 h = set(h, 'invars', invars([j k]));
0237 catch
0238 warning('!!! Unable to set input variable names for history. Perhaps you input a vector of AOs.');
0239 end
0240
0241 b = ao(fs, h);
0242
0243 try
0244 if isempty(invars{k})
0245 nk = bo.name;
0246 else
0247 nk = invars{k};
0248 end
0249 if isempty(invars{j})
0250 nj = bi.name;
0251 else
0252 nj = invars{j};
0253 end
0254 catch
0255 nk = bo.name;
0256 nj = bi.name;
0257 end
0258 b = set(b, 'name', ['\sigma^2' sprintf('(%s->%s)', nj, nk)]);
0259
0260 sig(j,k) = {b};
0261
0262 end
0263 end
0264 end
0265
0266
0267
0268 b = [bs{:}];
0269 bs = reshape(b, na, na);
0270 varargout{1} = bs;
0271 if nargout == 2
0272 b = [sig{:}];
0273 s = reshape(b, na, na);
0274 varargout{2} = s;
0275 end
0276
0277
0278 function lmin = findShortestVector(as)
0279
0280 lmin = 1e20;
0281 for j=1:length(as)
0282 if len(as(j)) < lmin
0283 lmin = len(as(j));
0284 end
0285 end
0286
0287
0288
0289
0290 function fs = findFsMax(as)
0291
0292 fs = 0;
0293 for j=1:length(as)
0294 a = as(j);
0295 if a.data.fs > fs
0296 fs = a.data.fs;
0297 end
0298 end
0299
0300
0301
0302 function plo = getDefaultPL()
0303
0304 disp('* creating default plist...');
0305 plo = plist();
0306 plo = append(plo, param('Kdes', 100));
0307 plo = append(plo, param('Jdes', 1000));
0308 plo = append(plo, param('Kmin', 2));
0309 plo = append(plo, param('Win', specwin('Kaiser', 1, 150)));
0310 plo = append(plo, param('Nolap', -1));
0311 plo = append(plo, param('Order', 0));
0312 disp('* done.');
0313
0314
0315
0316
0317 function varargout = mltfe(varargin)
0318
0319
0320
0321
0322
0323
0324
0325
0326
0327
0328 X = varargin{1};
0329 f = varargin{2};
0330 r = varargin{3};
0331 m = varargin{4};
0332 L = varargin{5};
0333 fs = varargin{6};
0334 win = varargin{7};
0335 order = varargin{8};
0336
0337
0338 if nargout == 2
0339 computeVariance = 1;
0340 else
0341 computeVariance = 0;
0342 end
0343
0344 twopi = 2.0*pi;
0345 twopiofs = 2.0*pi/(1.0*fs);
0346
0347 si = size(X);
0348 nx = si(2);
0349 nc = si(1);
0350 nf = length(f);
0351 Txy = zeros(nc,nc,nf);
0352 Sig = zeros(nc,nc,nf);
0353 ENBW = zeros(nc,nc,nf);
0354
0355 disp_each = round(nf/100)*10;
0356
0357 for fi=1:nf
0358
0359
0360 l = L(fi);
0361 switch win.name
0362 case 'Kaiser'
0363 win = specwin(win.name, l, win.psll);
0364 otherwise
0365 win = specwin(win.name, l);
0366 end
0367
0368 p = 1i * twopi * m(fi)/l.*[0:l-1];
0369 C = win.win .* exp(p);
0370 Xolap = (1-win.rov/100);
0371
0372 if mod(fi,disp_each) == 0 || fi == 1 || fi == nf
0373 disp(sprintf('++ computing frequency %d of %d: %f Hz', fi, nf, f(fi)));
0374 end
0375
0376
0377 for j=1:nc
0378
0379 for k=1:nc
0380
0381
0382 start = 1;
0383 nsegs = 0;
0384 while start + l - 1 <= nx
0385
0386 start = floor(start + l*Xolap);
0387 nsegs = nsegs + 1;
0388 end
0389
0390
0391 A = 0.0;
0392 T = zeros(nsegs,1);
0393 i2 = zeros(nsegs,1);
0394
0395 start = 1;
0396 nsegs = 0;
0397 p = 1;
0398 while start + l - 1 <= nx
0399
0400
0401 idx = start:start+l-1;
0402 xi = X(j,idx).';
0403 xo = X(k,idx).';
0404
0405
0406 switch order
0407 case -1
0408
0409 case 0
0410 xi = xi - mean(xi);
0411 xo = xo - mean(xo);
0412 case 1
0413 xi = detrend(xi);
0414 xo = detrend(xo);
0415 otherwise
0416 xi = polydetrend(xi, order);
0417 xo = polydetrend(xo, order);
0418 end
0419
0420
0421 ti = C*xi;
0422 to = C*xo;
0423 i2(p) = (ti*conj(ti));
0424 T(p) = to/ti;
0425 A = A + T(p)*i2(p);
0426
0427
0428 start = floor(start + l*Xolap);
0429 nsegs = nsegs + 1;
0430 p = p + 1;
0431 end
0432
0433 sw = sum(i2);
0434 Txy(j,k,fi) = conj(A/sw);
0435
0436 S1 = win.ws;
0437 S12 = S1*S1;
0438 S2 = win.ws2;
0439 ENBW(j,k,fi) = fs*S2/S12;
0440
0441
0442 if computeVariance
0443 A = 0;
0444 start = 1;
0445 p = 1;
0446 while start + l - 1 <= nx
0447
0448 S = (T(p)-Txy(j,k,fi))^2;
0449 A = A + S * i2(p);
0450
0451 start = floor(start + l*Xolap);
0452 p = p + 1;
0453 end
0454 Sig(j,k,fi) = abs(A/sw);
0455 end
0456
0457 end
0458 end
0459 end
0460
0461 varargout{1} = Txy;
0462 if nargout >= 2
0463 varargout{2} = ENBW;
0464 end
0465 if nargout == 3
0466 varargout{3} = Sig;
0467 end
0468
0469