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 ALGONAME = mfilename;
0048 VERSION = '$Id: ltpda_lcohere.html,v 1.12 2008/03/31 10:27:39 hewitson Exp $';
0049 CATEGORY = 'Signal Processing';
0050
0051
0052 if nargin == 1
0053 in = char(varargin{1});
0054 if strcmp(in, 'Params')
0055 varargout{1} = getDefaultPL();
0056 return
0057 elseif strcmp(in, 'Version')
0058 varargout{1} = VERSION;
0059 return
0060 elseif strcmp(in, 'Category')
0061 varargout{1} = CATEGORY;
0062 return
0063 end
0064 end
0065
0066
0067
0068 in_names = {};
0069 for ii = 1:nargin
0070 in_names{end+1} = inputname(ii);
0071 end
0072 [as, ps, invars] = collect_inputs(varargin, in_names);
0073
0074
0075 bs = [];
0076
0077
0078 if isempty(ps)
0079 pl = getDefaultPL();
0080 else
0081 pl = combine(ps, getDefaultPL);
0082 end
0083
0084
0085 na = length(as);
0086
0087
0088
0089 disp('*** Resampling AOs...');
0090 fsmax = findFsMax(as);
0091 fspl = plist(param('fsout', fsmax));
0092 bs = [];
0093 for i=1:na
0094 a = as(i);
0095
0096 if ~isa(a.data, 'tsdata')
0097 error('### ltpda_lcohere requires tsdata (time-series) inputs.');
0098 end
0099
0100 if a.data.fs ~= fsmax
0101 warning(sprintf('!!! Resampling AO %s/%s to %f Hz', a.name, a.data.name, fsmax));
0102 a = resample(a, fspl);
0103 end
0104 bs = [bs a];
0105 end
0106 as = bs;
0107 na = length(as);
0108 clear bs;
0109 disp('*** Done.');
0110
0111
0112
0113
0114 disp('*** Truncating all vectors...');
0115 lmin = findShortestVector(as);
0116 nsecs = lmin / fsmax;
0117 bs = [];
0118 for i=1:na
0119 a = as(i);
0120 if len(a) ~= lmin
0121 warning(sprintf('!!! Truncating AO %s/%s to %d secs', a.name, a.data.name, nsecs));
0122 bs = [bs select(a, 1:lmin)];
0123 else
0124 bs = [bs a];
0125 end
0126 end
0127 as = bs;
0128 clear bs;
0129 disp('*** Done.');
0130
0131
0132
0133 N = len(as(1));
0134 iS = zeros(na, N);
0135 for j=1:na
0136 a = as(j);
0137 iS(j,:) = a.data.y;
0138 end
0139
0140
0141
0142
0143 Kdes = find(pl, 'Kdes');
0144 Kmin = find(pl, 'Kmin');
0145 Jdes = find(pl, 'Jdes');
0146 Win = find(pl, 'Win');
0147 Nolap = find(pl, 'Nolap');
0148 if Nolap == -1
0149
0150 Nolap = Win.rov/100;
0151 disp(sprintf('! Using recommended overlap of %d', Nolap))
0152 pl = pset(pl, 'Nolap', Nolap);
0153 end
0154 Order = find(pl, 'Order');
0155 Lmin = find(pl, 'Lmin');
0156 if isempty(Lmin)
0157 Lmin = 0;
0158 end
0159
0160
0161 debug = find(pl, 'DEBUG');
0162
0163
0164 [f, r, m, L, K] = ltpda_ltf_plan(lmin, fsmax, Nolap, 1, Lmin, Jdes, Kdes);
0165
0166
0167
0168
0169 Txy = mlcohere(iS, f, r, m, L, fsmax, Win, Order, Nolap*100, Lmin);
0170
0171
0172
0173
0174
0175
0176
0177
0178
0179
0180
0181
0182
0183
0184
0185 disp('** Building output AO matrix')
0186 bs = cell(na);
0187 for j=1:na
0188 bi = as(j);
0189 for k=1:na
0190 bo = as(k);
0191
0192 fs = fsdata(f, squeeze(Txy(j,k,:)), fsmax);
0193 fs = set(fs, 'name', sprintf('LCOHERE(%s->%s)', bi.data.name, bo.data.name));
0194 fs = set(fs, 'yunits', [bo.data.yunits '/' bi.data.yunits]);
0195
0196 h = history(ALGONAME, VERSION, pl, [bi.hist bo.hist]);
0197 if na == length(invars)
0198 h = set(h, 'invars', invars([j k]));
0199 nk = invars{k};
0200 nj = invars{j};
0201 else
0202 nk = bo.name;
0203 nj = bi.name;
0204 end
0205
0206 b = ao(fs, h);
0207
0208 b = setnh(b, 'name', sprintf('LCOHERE(%s->%s)', nj, nk));
0209
0210 bs(j,k) = {b};
0211 end
0212 end
0213
0214
0215
0216 b = [bs{:}];
0217 bs = reshape(b, na, na);
0218 varargout{1} = bs;
0219
0220
0221
0222 function lmin = findShortestVector(as)
0223
0224 lmin = 1e20;
0225 for j=1:length(as)
0226 if len(as(j)) < lmin
0227 lmin = len(as(j));
0228 end
0229 end
0230
0231
0232
0233
0234 function fs = findFsMax(as)
0235
0236 fs = 0;
0237 for j=1:length(as)
0238 a = as(j);
0239 if a.data.fs > fs
0240 fs = a.data.fs;
0241 end
0242 end
0243
0244
0245
0246 function plo = getDefaultPL()
0247
0248 disp('* creating default plist...');
0249 plo = plist();
0250 plo = append(plo, param('Kdes', 100));
0251 plo = append(plo, param('Jdes', 1000));
0252 plo = append(plo, param('Kmin', 2));
0253 plo = append(plo, param('Lmin', 0));
0254 plo = append(plo, param('Win', specwin('Kaiser', 1, 200)));
0255 plo = append(plo, param('Nolap', -1));
0256 plo = append(plo, param('Order', 0));
0257 disp('* done.');
0258
0259
0260
0261
0262 function varargout = mlcohere(varargin)
0263
0264
0265
0266
0267
0268
0269
0270
0271
0272
0273 X = varargin{1};
0274 f = varargin{2};
0275 r = varargin{3};
0276 m = varargin{4};
0277 L = varargin{5};
0278 fs = varargin{6};
0279 win = varargin{7};
0280 order = varargin{8};
0281 olap = varargin{9};
0282 Lmin = varargin{10};
0283
0284
0285
0286 twopi = 2.0*pi;
0287 si = size(X);
0288 nc = si(1);
0289 nf = length(f);
0290 Txy = zeros(nc,nc,nf);
0291 disp_each = round(nf/100)*10;
0292 minReached = 0;
0293
0294
0295 for fi=1:nf
0296
0297
0298 l = L(fi);
0299
0300
0301
0302 if ~minReached
0303 switch win.name
0304 case 'Kaiser'
0305 win = specwin(win.name, l, win.psll);
0306 otherwise
0307 win = specwin(win.name, l);
0308 end
0309 if l==Lmin
0310 minReached = 1;
0311 end
0312 end
0313
0314
0315 p = 1i * twopi * m(fi)/l.*[0:l-1];
0316 C = win.win .* exp(p);
0317
0318 if mod(fi,disp_each) == 0 || fi == 1 || fi == nf
0319 disp(sprintf('++ computing frequency %04d of %04d: %f Hz', fi, nf, f(fi)));
0320 end
0321
0322
0323 for j=1:nc
0324
0325 for k=1:nc
0326 if k>j
0327
0328 xi = X(j,:);
0329 xo = X(k,:);
0330
0331 [XY, XX, YY, nsegs] = ltpda_dft(xi, xo, l, C, olap, order);
0332 Txy(j,k,fi) = (abs(XY).^2)/(XX.*YY);
0333 else
0334
0335 Txy(j,k,fi) = 1;
0336 end
0337 end
0338 end
0339 end
0340
0341
0342 for j=1:nc
0343 for k=1:nc
0344 if k<j
0345 Txy(j,k,:) = Txy(k,j,:);
0346 end
0347 end
0348 end
0349
0350
0351 varargout{1} = Txy;
0352
0353
0354