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