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