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 if nargin == 1
0044 in = char(varargin{1});
0045 if strcmp(in, 'Params')
0046 varargout{1} = getDefaultPL();
0047 return
0048 end
0049 end
0050
0051
0052 invars = {};
0053 as = [];
0054 ps = [];
0055 for j=1:nargin
0056 invars = [invars cellstr(inputname(j))];
0057 if isa(varargin{j}, 'ao')
0058 as = [as varargin{j}];
0059 end
0060 if isa(varargin{j}, 'plist')
0061 ps = [ps varargin{j}];
0062 end
0063 end
0064
0065 ALGONAME = mfilename;
0066 VERSION = '$Id: ltpda_ltfe.html,v 1.2 2007/07/10 05:37:14 hewitson Exp $';
0067
0068
0069 bs = [];
0070
0071
0072 if isempty(ps)
0073 pl = getDefaultPL();
0074 else
0075 pl = combine(ps, getDefaultPL);
0076 end
0077
0078
0079 na = length(as);
0080
0081
0082
0083 disp('*** Resampling AOs...');
0084 fsmax = findFsMax(as);
0085 fspl = plist(param('fsout', fsmax));
0086 bs = [];
0087 for i=1:na
0088 a = as(i);
0089
0090 if ~isa(a.data, 'tsdata')
0091 error('### ltpda_ltfe requires tsdata (time-series) inputs.');
0092 end
0093
0094 if a.data.fs ~= fsmax
0095 warning(sprintf('!!! Resampling AO %s/%s to %f Hz', a.name, a.data.name, fsmax));
0096 a = resample(a, fspl);
0097 end
0098 bs = [bs a];
0099 end
0100 as = bs;
0101 clear bs;
0102 disp('*** Done.');
0103
0104
0105
0106
0107 disp('*** Truncating all vectors...');
0108 lmin = findShortestVector(as);
0109 nsecs = lmin / fsmax;
0110 bs = [];
0111 for i=1:na
0112 a = as(i);
0113 if len(a) ~= lmin
0114 warning(sprintf('!!! Truncating AO %s/%s to %d secs', a.name, a.data.name, nsecs));
0115 bs = [bs select(a, 1:lmin)];
0116 else
0117 bs = [bs a];
0118 end
0119 end
0120 as = bs;
0121 clear bs;
0122 disp('*** Done.');
0123
0124
0125
0126 N = len(as(1));
0127 iS = zeros(na, N);
0128 for j=1:na
0129 a = as(j);
0130 iS(j,:) = a.data.x;
0131 end
0132
0133
0134
0135
0136 Kdes = find(pl, 'Kdes');
0137 Kmin = find(pl, 'Kmin');
0138 Jdes = find(pl, 'Jdes');
0139 if Jdes == -1
0140 Jdes = round(fsmax/4);
0141 disp(sprintf('! Using default Jdes of %d', Jdes))
0142 pl = set(pl, 'Jdes', Jdes);
0143 end
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
0153
0154 [f,r,m,L] = ltpda_compute_f(fsmax, lmin, Kdes, Kmin, Jdes, Nolap);
0155
0156
0157
0158 if nargout == 2
0159 disp('** Computing TF matrix with variance estimates')
0160 tic
0161 [Txy, Sig] = mltfe(iS, f, r, m, L, fsmax, Win);
0162 toc
0163 else
0164 disp('** Computing TF matrix')
0165 tic
0166 Txy = mltfe(iS, f, r, m, L, fsmax, Win);
0167 toc
0168 end
0169
0170
0171 disp('** Building output AO matrix')
0172 bs = cell(na);
0173 sig = cell(na);
0174 for j=1:na
0175 bi = as(j);
0176 for k=1:na
0177 bo = as(k);
0178
0179 fs = fsdata(f, squeeze(Txy(j,k,:)), fsmax);
0180 fs = set(fs, 'name', sprintf('LTFe(%s->%s)', bi.data.name, bo.data.name));
0181 fs = set(fs, 'yunits', [bo.data.yunits '/' bi.data.yunits]);
0182 fs = set(fs, 'enbw', 0);
0183
0184 h = history(ALGONAME, VERSION, pl, [bi.hist bo.hist]);
0185 h = set(h, 'invars', invars([j k]));
0186
0187 b = ao(fs, h);
0188
0189 if isempty(invars{k})
0190 nk = bo.name;
0191 else
0192 nk = invars{k};
0193 end
0194 if isempty(invars{j})
0195 nj = bi.name;
0196 else
0197 nj = invars{j};
0198 end
0199 b = set(b, 'name', sprintf('LTFe(%s->%s)', nj, nk));
0200
0201 bs(j,k) = {b};
0202
0203
0204 if nargout == 2
0205
0206 fs = fsdata(f, squeeze(Sig(j,k,:)), fsmax);
0207 fs = set(fs, 'name', ['\sigma^2' sprintf('(%s->%s)', bi.data.name, bo.data.name)]);
0208 fs = set(fs, 'yunits', [bo.data.yunits '/' bi.data.yunits]);
0209 fs = set(fs, 'enbw', 0);
0210
0211 h = history(ALGONAME, VERSION, pl, [bi.hist bo.hist]);
0212 h = set(h, 'invars', invars{[j k]});
0213
0214 b = ao(fs, h);
0215
0216 if isempty(invars{k})
0217 nk = bo.name;
0218 else
0219 nk = invars{k};
0220 end
0221 if isempty(invars{j})
0222 nj = bi.name;
0223 else
0224 nj = invars{j};
0225 end
0226 b = set(b, 'name', ['\sigma^2' sprintf('(%s->%s)', nj, nk)]);
0227
0228 sig(j,k) = {b};
0229
0230 end
0231 end
0232 end
0233
0234
0235
0236 b = [bs{:}];
0237 bs = reshape(b, na, na);
0238 varargout{1} = bs;
0239 if nargout == 2
0240 b = [sig{:}];
0241 s = reshape(b, na, na);
0242 varargout{2} = s;
0243 end
0244
0245
0246 function lmin = findShortestVector(as)
0247
0248 lmin = 1e20;
0249 for j=1:length(as)
0250 if len(as(j)) < lmin
0251 lmin = len(as(j));
0252 end
0253 end
0254
0255
0256
0257
0258 function fs = findFsMax(as)
0259
0260 fs = 0;
0261 for j=1:length(as)
0262 a = as(j);
0263 if a.data.fs > fs
0264 fs = a.data.fs;
0265 end
0266 end
0267
0268
0269
0270 function plo = getDefaultPL()
0271
0272 disp('* creating default plist...');
0273 plo = plist();
0274 plo = append(plo, param('Kdes', 100));
0275 plo = append(plo, param('Jdes', 1000));
0276 plo = append(plo, param('Kmin', 2));
0277 plo = append(plo, param('Win', specwin('Kaiser', 1, 150)));
0278 plo = append(plo, param('Nolap', -1));
0279 disp('* done.');
0280
0281
0282
0283
0284 function varargout = mltfe(varargin)
0285
0286
0287
0288
0289
0290
0291
0292
0293
0294
0295 X = varargin{1};
0296 f = varargin{2};
0297 r = varargin{3};
0298 m = varargin{4};
0299 L = varargin{5};
0300 fs = varargin{6};
0301 win = varargin{7};
0302
0303
0304 if nargout == 2
0305 computeVariance = 1;
0306 else
0307 computeVariance = 0;
0308 end
0309
0310 twopi = 2.0*pi;
0311 twopiofs = 2.0*pi/(1.0*fs);
0312
0313 si = size(X);
0314 nx = si(2);
0315 nc = si(1);
0316 nf = length(f);
0317 Txy = zeros(nc,nc,nf);
0318 Sig = zeros(nc,nc,nf);
0319
0320 disp_each = round(nf/100)*10;
0321
0322 for fi=1:nf
0323
0324
0325 l = L(fi);
0326 switch win.name
0327 case 'Kaiser'
0328 win = specwin(win.name, l, win.psll);
0329 otherwise
0330 win = specwin(win.name, l);
0331 end
0332
0333 p = 1i * twopi * m(fi)/l.*[0:l-1];
0334 C = win.win .* exp(p);
0335 Xolap = (1-win.rov/100);
0336
0337 if mod(fi,disp_each) == 0 || fi == 1 || fi == nf
0338 disp(sprintf('++ computing frequency %d of %d: %f Hz', fi, nf, f(fi)));
0339 end
0340
0341
0342 for j=1:nc
0343
0344 for k=1:nc
0345
0346
0347 start = 1;
0348 nsegs = 0;
0349 while start + l - 1 <= nx
0350
0351 start = floor(start + l*Xolap);
0352 nsegs = nsegs + 1;
0353 end
0354
0355
0356 A = 0.0;
0357 T = zeros(nsegs,1);
0358 i2 = zeros(nsegs,1);
0359 start = 1;
0360 nsegs = 0;
0361 p = 1;
0362 while start + l - 1 <= nx
0363
0364
0365
0366
0367
0368 idx = start:start+l-1;
0369 xi = X(j,idx).';
0370 xo = X(k,idx).';
0371
0372
0373 ti = C*xi;
0374 to = C*xo;
0375 i2(p) = ti*conj(ti);
0376 T(p) = (to/ti);
0377 A = A + T(p)*i2(p);
0378
0379
0380 start = floor(start + l*Xolap);
0381 nsegs = nsegs + 1;
0382 p = p + 1;
0383 end
0384
0385 sw = sum(i2);
0386 Txy(j,k,fi) = conj(A/sw);
0387
0388
0389 if computeVariance
0390 A = 0;
0391 start = 1;
0392 p = 1;
0393 while start + l - 1 <= nx
0394
0395 S = (T(p)-Txy(j,k,fi))^2;
0396 A = A + S * i2(p);
0397
0398 start = floor(start + l*Xolap);
0399 p = p + 1;
0400 end
0401 Sig(j,k,fi) = abs(A/sw);
0402 end
0403
0404 end
0405 end
0406 end
0407
0408 varargout{1} = Txy;
0409 if nargout == 2
0410 varargout{2} = Sig;
0411 end
0412
0413
0414