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
0044
0045
0046
0047
0048
0049
0050
0051
0052
0053 ALGONAME = mfilename;
0054 VERSION = '$Id: ltpda_ltfe.html,v 1.14 2008/03/31 10:27:39 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_ltfe 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 na = length(as);
0113 clear bs;
0114 disp('*** Done.');
0115
0116
0117
0118
0119 disp('*** Truncating all vectors...');
0120 lmin = findShortestVector(as);
0121 nsecs = lmin / fsmax;
0122 bs = [];
0123 for i=1:na
0124 a = as(i);
0125 if len(a) ~= lmin
0126 warning(sprintf('!!! Truncating AO %s/%s to %d secs', a.name, a.data.name, nsecs));
0127 bs = [bs select(a, 1:lmin)];
0128 else
0129 bs = [bs a];
0130 end
0131 end
0132 as = bs;
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 if find(pl, 'OLD SCHEDULER')
0177 [f,r,m,L,rr,rrr] = ltpda_compute_f(fsmax, lmin, Kdes, Kmin, Jdes, Nolap, 0);
0178 else
0179 [f, r, m, L, K] = ltpda_ltf_plan(lmin, fsmax, Nolap, 1, Lmin, Jdes, Kdes);
0180 end
0181
0182
0183
0184 Txy = mltfe(iS, f, r, m, L, fsmax, Win, Order, Nolap*100, Lmin);
0185
0186
0187 disp('** Building output AO matrix')
0188 bs = cell(na);
0189 for j=1:na
0190 bi = as(j);
0191 for k=1:na
0192 bo = as(k);
0193
0194 fs = fsdata(f, squeeze(Txy(j,k,:)), fsmax);
0195 fs = set(fs, 'name', sprintf('LTFE(%s->%s)', bi.data.name, bo.data.name));
0196 fs = set(fs, 'yunits', [bo.data.yunits '/' bi.data.yunits]);
0197
0198 h = history(ALGONAME, VERSION, pl, [bi.hist bo.hist]);
0199 try
0200 h = set(h, 'invars', invars([j k]));
0201 catch
0202 warning('!!! Unable to set input variable names for history. Perhaps you input a vector of AOs.');
0203 end
0204
0205 b = ao(fs, h);
0206
0207 try
0208 if isempty(invars{k})
0209 nk = bo.name;
0210 else
0211 nk = invars{k};
0212 end
0213 if isempty(invars{j})
0214 nj = bi.name;
0215 else
0216 nj = invars{j};
0217 end
0218 catch
0219 warning('!!! Unable to set input variable names for history. Perhaps you input a vector of AOs.');
0220 nk = bo.name;
0221 nj = bi.name;
0222 end
0223
0224 b = setnh(b, 'name', sprintf('LTFE(%s->%s)', nj, nk));
0225
0226 bs(j,k) = {b};
0227
0228 end
0229 end
0230
0231
0232
0233 b = [bs{:}];
0234 bs = reshape(b, na, na);
0235 varargout{1} = bs;
0236
0237
0238
0239 function lmin = findShortestVector(as)
0240
0241 lmin = 1e20;
0242 for j=1:length(as)
0243 if len(as(j)) < lmin
0244 lmin = len(as(j));
0245 end
0246 end
0247
0248
0249
0250
0251 function fs = findFsMax(as)
0252
0253 fs = 0;
0254 for j=1:length(as)
0255 a = as(j);
0256 if a.data.fs > fs
0257 fs = a.data.fs;
0258 end
0259 end
0260
0261
0262
0263 function plo = getDefaultPL()
0264
0265 plo = plist();
0266 plo = append(plo, param('Kdes', 100));
0267 plo = append(plo, param('Jdes', 1000));
0268 plo = append(plo, param('Kmin', 2));
0269 plo = append(plo, param('Lmin', 0));
0270 plo = append(plo, param('Win', specwin('Kaiser', 1, 150)));
0271 plo = append(plo, param('Nolap', -1));
0272 plo = append(plo, param('Order', 0));
0273
0274
0275
0276
0277 function varargout = mltfe(varargin)
0278
0279
0280
0281
0282
0283
0284
0285
0286
0287
0288 X = varargin{1};
0289 f = varargin{2};
0290 r = varargin{3};
0291 m = varargin{4};
0292 L = varargin{5};
0293 fs = varargin{6};
0294 win = varargin{7};
0295 order = varargin{8};
0296 olap = varargin{9};
0297 Lmin = varargin{10};
0298
0299
0300
0301 twopi = 2.0*pi;
0302 si = size(X);
0303 nc = si(1);
0304 nf = length(f);
0305 Txy = zeros(nc,nc,nf);
0306 disp_each = round(nf/100)*10;
0307 minReached = 0;
0308
0309
0310 for fi=1:nf
0311
0312
0313 l = L(fi);
0314
0315
0316
0317 if ~minReached
0318 switch win.name
0319 case 'Kaiser'
0320 win = specwin(win.name, l, win.psll);
0321 otherwise
0322 win = specwin(win.name, l);
0323 end
0324 if l==Lmin
0325 minReached = 1;
0326 end
0327 end
0328
0329
0330 p = 1i * twopi * m(fi)/l.*[0:l-1];
0331 C = win.win .* exp(p);
0332
0333 if mod(fi,disp_each) == 0 || fi == 1 || fi == nf
0334 disp(sprintf('++ computing frequency %d of %d: %f Hz', fi, nf, f(fi)));
0335 end
0336
0337
0338 for j=1:nc
0339
0340 for k=1:nc
0341 if k~=j
0342
0343 xi = X(j,:);
0344 xo = X(k,:);
0345
0346 [XY, XX, YY, nsegs] = ltpda_dft(xi, xo, l, C, olap, order);
0347 Txy(j,k,fi) = XY/XX;
0348 else
0349
0350 Txy(j,k,fi) = 1;
0351 end
0352 end
0353 end
0354 end
0355
0356
0357 varargout{1} = Txy;
0358
0359