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.m,v 1.5 2007/07/16 12:52:21 ingo 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 = set(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 = set(pl, 'Win', Win);
0123 end
0124 end
0125
0126
0127 Nolap = find(pl, 'Nolap');
0128 if ~isempty(Nolap)
0129 if Nolap < 0
0130
0131 Nolap = Win.rov/100;
0132 disp(sprintf('! Using recommended overlap of %d', Nolap))
0133 pl = set(pl, 'Nolap', Nolap);
0134 end
0135 end
0136
0137
0138 bs = cell(na);
0139 for j=1:na
0140 aj = as(j);
0141 for k=1:na
0142 ak = as(k);
0143
0144
0145 switch method
0146 case 'TF'
0147
0148 [txy,f] = tfestimate(aj.data.x,ak.data.x,Win.win,round(Nolap*Nfft),Nfft,fsmax);
0149 case 'COHERE'
0150
0151 [txy,f] = mscohere(aj.data.x,ak.data.x,Win.win,round(Nolap*Nfft),Nfft,fsmax);
0152 case 'CPSD'
0153
0154 [txy,f] = cpsd(aj.data.x,ak.data.x,Win.win,round(Nolap*Nfft),Nfft,fsmax);
0155 otherwise
0156 error('### Unknown method.')
0157 end
0158
0159
0160 fs = fsdata(f, txy, fsmax);
0161 fs = set(fs, 'name', sprintf('%s(%s->%s)', method, aj.data.name, ak.data.name));
0162 fs = set(fs, 'enbw', Win.nenbw*fsmax/Nfft);
0163
0164
0165
0166
0167
0168
0169
0170 if j<k
0171 h = history(ALGONAME, VERSION, pl, [aj.hist ak.hist]);
0172 else
0173 h = history(ALGONAME, VERSION, pl, [ak.hist aj.hist]);
0174 end
0175 h = set(h, 'invars', [invars(j) invars(k)]);
0176
0177
0178 b = ao(fs, h);
0179
0180
0181 if isempty(invars{k})
0182 nk = ak.name;
0183 else
0184 nk = invars{k};
0185 end
0186 if isempty(invars{j})
0187 nj = aj.name;
0188 else
0189 nj = invars{j};
0190 end
0191 b = set(b, 'name', sprintf('%s(%s->%s)', method, nj, nk));
0192
0193
0194 bs(j,k) = {b};
0195 end
0196 end
0197
0198
0199
0200 b = [bs{:}];
0201 varargout{1} = reshape(b, na, na);
0202
0203
0204
0205
0206 function lmin = findShortestVector(as)
0207
0208 lmin = 1e20;
0209 for j=1:length(as)
0210 if len(as(j)) < lmin
0211 lmin = len(as(j));
0212 end
0213 end
0214
0215
0216
0217
0218 function fs = findFsMax(as)
0219
0220 fs = 0;
0221 for j=1:length(as)
0222 a = as(j);
0223 if a.data.fs > fs
0224 fs = a.data.fs;
0225 end
0226 end
0227
0228
0229