0001 function varargout = ltpda_lcpsd(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_lcpsd 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 CPSD matrix with variance estimates')
0172 tic
0173 [Txy, ENBW, Sig] = mlcpsd(iS, f, r, m, L, fsmax, Win, Order);
0174 toc
0175 else
0176 disp('** Computing CPSD matrix')
0177 tic
0178 [Txy, ENBW] = mlcpsd(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('LCPSD(%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('LCPSD(%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 = mlcpsd(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 A = 0.0;
0372 T = zeros(nsegs,1);
0373 i2 = zeros(nsegs,1);
0374 start = 1;
0375 nsegs = 0;
0376 p = 1;
0377 while start + l - 1 <= nx
0378
0379
0380 idx = start:start+l-1;
0381 xi = X(j,idx).';
0382 xo = X(k,idx).';
0383
0384
0385 switch order
0386 case -1
0387
0388 case 0
0389 xi = xi - mean(xi);
0390 xo = xo - mean(xo);
0391 case 1
0392 xi = detrend(xi);
0393 xo = detrend(xo);
0394 otherwise
0395 xi = polydetrend(xi, order);
0396 xo = polydetrend(xo, order);
0397 end
0398
0399
0400 ti = C*xi;
0401 to = C*xo;
0402
0403 T(p) = (to .* conj(ti));
0404
0405 A = A + T(p);
0406
0407
0408 start = floor(start + l*Xolap);
0409 nsegs = nsegs + 1;
0410 p = p + 1;
0411 end
0412
0413 A2ns = 2.0*A/nsegs;
0414 S1 = win.ws;
0415 S12 = S1*S1;
0416 S2 = win.ws2;
0417 ENBW(j,k,fi) = fs*S2/S12;
0418
0419 sw = sum(i2);
0420
0421 Txy(j,k,fi) = 2.0*A/nsegs/fs/S2;
0422
0423
0424 if computeVariance
0425 A = 0;
0426 start = 1;
0427 p = 1;
0428 scale = fs*S2;
0429 while start + l - 1 <= nx
0430
0431 S = (T(p)/scale-Txy(j,k,fi))^2;
0432 A = A + S * i2(p);
0433
0434 start = floor(start + l*Xolap);
0435 p = p + 1;
0436 end
0437 Sig(j,k,fi) = A/sw;
0438 end
0439
0440 end
0441 end
0442 end
0443
0444 varargout{1} = Txy;
0445 if nargout == 2
0446 varargout{2} = ENBW;
0447 elseif nargout == 3
0448 varargout{3} = Sig;
0449 varargout{2} = ENBW;
0450 end
0451
0452
0453