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