0001 function varargout = ltpda_xspec(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 as = varargin{1};
0036 pl = varargin{2};
0037 method = varargin{3};
0038 iALGO = varargin{4};
0039 iVER = varargin{5};
0040 invars = varargin{6};
0041
0042 ALGONAME = iALGO;
0043 VERSION = [iVER '/' '$Id: ltpda_xspec.html,v 1.14 2008/03/31 10:27:39 hewitson Exp $'];
0044
0045
0046
0047 bs = [];
0048
0049 na = length(as);
0050 if na < 2
0051 error('### LTPDA_XSPEC needs at least two AOs to make a transfer function.');
0052 end
0053
0054
0055
0056
0057 disp('*** Resampling AOs...');
0058 fsmax = findFsMax(as);
0059 fspl = plist(param('fsout', fsmax));
0060 bs = [];
0061 for i=1:na
0062 a = as(i);
0063
0064 if ~isa(a.data, 'tsdata')
0065 error('### ltpda_xspec requires tsdata (time-series) inputs.');
0066 end
0067
0068 if a.data.fs ~= fsmax
0069 warning(sprintf('!!! Resampling AO %s/%s to %f Hz', a.name, a.data.name, fsmax));
0070 a = resample(a, fspl);
0071 end
0072 bs = [bs a];
0073 end
0074 as = bs;
0075 clear bs;
0076 disp('*** Done.');
0077
0078
0079
0080
0081
0082 disp('*** Truncating all vectors...');
0083 lmin = findShortestVector(as);
0084 nsecs = lmin / fsmax;
0085 bs = [];
0086 for i=1:na
0087 a = as(i);
0088 if len(a) ~= lmin
0089 warning(sprintf('!!! Truncating AO %s/%s to %d secs', a.name, a.data.name, nsecs));
0090 bs = [bs select(a, 1:lmin)];
0091 else
0092 bs = [bs a];
0093 end
0094 end
0095 as = bs;
0096 clear bs;
0097 disp('*** Done.');
0098
0099
0100
0101
0102 Nfft = find(pl, 'Nfft');
0103 if ~isempty(Nfft)
0104 if Nfft < 0
0105 Nfft = abs(Nfft) * lmin;
0106 disp(sprintf('! Using default Nfft of %g', Nfft))
0107 pl = pset(pl, 'Nfft', Nfft);
0108 end
0109 end
0110
0111
0112 Win = find(pl, 'Win');
0113 if ~isempty(Win)
0114 if length(Win.win)~=Nfft
0115 switch Win.name
0116 case 'Kaiser'
0117 Win = specwin(Win.name, Nfft, Win.psll);
0118 otherwise
0119 Win = specwin(Win.name, Nfft);
0120 end
0121 disp(sprintf('! Reset window to %s(%d)', strrep(Win.name, '_', '\_'), length(Win.win)))
0122 pl = pset(pl, 'Win', Win);
0123 end
0124 end
0125
0126
0127 Olap = find(pl, 'Olap');
0128 if ~isempty(Olap)
0129 if Olap < 0
0130
0131 Olap = Win.rov;
0132 disp(sprintf('! Using recommended overlap of %2.1f%%', Olap))
0133 pl = pset(pl, 'Olap', Olap);
0134 end
0135 end
0136
0137
0138 order = find(pl, 'Order');
0139 if isempty(order)
0140 order = 0;
0141 warning('! using default detrending order 0 (mean)');
0142 pl = append(pl, 'Order', order);
0143 end
0144
0145
0146 bs = cell(na);
0147 for j=1:na
0148 aj = as(j);
0149 for k=1:na
0150 ak = as(k);
0151
0152
0153
0154 x = aj.data.y;
0155 y = ak.data.y;
0156 switch order
0157 case -1
0158
0159 case 0
0160 x = x - mean(x);
0161 y = y - mean(y);
0162 case 1
0163 x = detrend(x);
0164 y = detrend(y);
0165 otherwise
0166 y = polydetrend(y, order);
0167 y = polydetrend(y, order);
0168 end
0169
0170 switch method
0171 case 'TF'
0172
0173 [txy,f] = tfestimate(x,y,Win.win,round(Olap/100*Nfft),Nfft,fsmax);
0174 case 'COHERE'
0175
0176 [txy,f] = mscohere(x,y,Win.win,round(Olap/100*Nfft),Nfft,fsmax);
0177 case 'CPSD'
0178
0179 [txy,f] = cpsd(x,y,Win.win,round(Olap/100*Nfft),Nfft,fsmax);
0180 otherwise
0181 error('### Unknown method.')
0182 end
0183
0184
0185 fs = fsdata(f, txy, fsmax);
0186 fs = set(fs, 'name', sprintf('%s(%s->%s)', method, aj.data.name, ak.data.name));
0187 fs = set(fs, 'enbw', Win.nenbw*fsmax/Nfft);
0188
0189
0190
0191
0192
0193
0194
0195 if j<k
0196 h = history(ALGONAME, VERSION, pl, [aj.hist ak.hist]);
0197 else
0198 h = history(ALGONAME, VERSION, pl, [ak.hist aj.hist]);
0199 end
0200 if na == length(invars)
0201 h = set(h, 'invars', [invars(j) invars(k)]);
0202 nk = invars{k};
0203 nj = invars{j};
0204 else
0205 nk = ak.name;
0206 nj = aj.name;
0207 end
0208
0209
0210 b = ao(fs, h);
0211
0212
0213 b = setnh(b, 'name', sprintf('%s(%s->%s)', method, nj, nk));
0214
0215
0216 bs(j,k) = {b};
0217 end
0218 end
0219
0220
0221
0222 b = [bs{:}];
0223 varargout{1} = reshape(b, na, na);
0224
0225
0226
0227
0228 function lmin = findShortestVector(as)
0229
0230 lmin = 1e20;
0231 for j=1:length(as)
0232 if len(as(j)) < lmin
0233 lmin = len(as(j));
0234 end
0235 end
0236
0237
0238
0239
0240 function fs = findFsMax(as)
0241
0242 fs = 0;
0243 for j=1:length(as)
0244 a = as(j);
0245 if a.data.fs > fs
0246 fs = a.data.fs;
0247 end
0248 end
0249
0250
0251