0001
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 function varargout = xspec(varargin)
0031
0032 import utils.const.*
0033
0034
0035 as = varargin{1};
0036 pl = varargin{2};
0037 method = varargin{3};
0038 mi = varargin{4};
0039 invars = varargin{5};
0040
0041 VERSION = '$Id: xspec.m,v 1.15 2008/09/03 17:05:51 ingo Exp $';
0042
0043
0044 mi.setMversion([VERSION '-->' mi.mversion]);
0045
0046 na = numel(as);
0047 if na < 2
0048 error('### XSPEC needs at least two time-series AOs.');
0049 end
0050
0051
0052 copies = zeros(size(as));
0053
0054 fsmax = ao.findFsMax(as);
0055 fspl = plist(param('fsout', fsmax));
0056 for i=1:na
0057
0058 if ~isa(as(i).data, 'tsdata')
0059 error('### ao.xspec requires tsdata (time-series) inputs.');
0060 end
0061
0062 if as(i).data.fs ~= fsmax
0063 utils.helper.msg(msg.PROC1, 'resampling AO %s to %f Hz', as(i).name, fsmax);
0064
0065
0066 as(i) = copy(as(i), 1);
0067 copies(i) = 1;
0068 resample(as(i), fspl);
0069 end
0070 end
0071
0072
0073
0074
0075 lmin = ao.findShortestVector(as);
0076 nsecs = lmin / fsmax;
0077 for i=1:na
0078 if len(as(i)) ~= lmin
0079 utils.helper.msg(msg.PROC2, 'truncating AO %s to %d secs', as(i).name, nsecs);
0080
0081 if ~copies(i)
0082
0083
0084 as(i) = copy(as(i), 1);
0085 copies(i) = 1;
0086 end
0087 as(i).select(1:lmin);
0088 end
0089 end
0090
0091
0092
0093
0094 Nfft = find(pl, 'Nfft');
0095 if ~isempty(Nfft)
0096 if Nfft < 0
0097 Nfft = abs(Nfft) * lmin;
0098 utils.helper.msg(msg.PROC2, 'using default Nfft of %g', Nfft);
0099 pl.pset('Nfft', Nfft);
0100 end
0101 end
0102
0103
0104 Win = find(pl, 'Win');
0105 if ~isempty(Win)
0106 if length(Win.win)~=Nfft
0107 switch Win.type
0108 case 'Kaiser'
0109 Win = specwin(Win.type, Nfft, Win.psll);
0110 otherwise
0111 Win = specwin(Win.type, Nfft);
0112 end
0113 utils.helper.msg(msg.PROC2, 'reset window to %s(%d)', strrep(Win.type, '_', '\_'), length(Win.win));
0114 pl.pset('Win', Win);
0115 end
0116 end
0117
0118
0119 Olap = find(pl, 'Olap');
0120 if ~isempty(Olap)
0121 if Olap < 0
0122
0123 Olap = Win.rov;
0124 utils.helper.msg(msg.PROC2, 'using recommended overlap of %2.1f%%', Olap);
0125 pl.pset('Olap', Olap);
0126 end
0127 end
0128
0129
0130 order = find(pl, 'Order');
0131 if isempty(order)
0132 order = 0;
0133 utils.helper.msg(msg.PROC1, 'using default detrending order 0 (mean)');
0134 pl.append('Order', order);
0135 end
0136
0137
0138 bs(na,na) = ao;
0139 for j=1:na
0140 for k=1:na
0141 if k>=j
0142 utils.helper.msg(msg.PROC1, 'computing (%d,%d) %s -> %s', j, k, invars{j}, invars{k});
0143
0144
0145
0146
0147 [txy, f, info] = ao.welch(as(j), as(k), method, pl);
0148
0149
0150 if size(as(1).data.y,1) == 1
0151 txy = txy.';
0152 f = f.';
0153 end
0154
0155
0156 fs = fsdata(f, txy, fsmax);
0157 fs.setEnbw(info.enbw);
0158 fs.setXunits('Hz');
0159 fs.setYunits(info.units);
0160 fs.setNavs(info.navs);
0161
0162
0163 bs(j,k) = ao(fs);
0164
0165
0166
0167
0168
0169
0170
0171 if j<k
0172 bs(j,k).addHistory(mi, pl, [invars(j) invars(k)], [as(j).hist, as(k).hist]);
0173 else
0174 bs(j,k).addHistory(mi, pl, [invars(j) invars(k)], [as(k).hist, as(j).hist]);
0175 end
0176
0177
0178 bs(j,k).setName(sprintf('%s(%s->%s)', method, invars{j}, invars{k}), 'internal');
0179
0180 end
0181 end
0182 end
0183
0184
0185 for j=1:na
0186 for k=1:na
0187 if k<j
0188
0189
0190 bs(j,k) = bs(k,j);
0191 end
0192 end
0193 end
0194
0195
0196 varargout{1} = bs;
0197 end
0198