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
0031
0032
0033
0034
0035
0036
0037
0038
0039
0040
0041
0042
0043
0044
0045
0046
0047 function varargout = welch(varargin)
0048
0049 if nargin == 3
0050 a = varargin{1};
0051 esttype = varargin{2};
0052 pl = varargin{3};
0053 x = a.data.y;
0054 inunits = a.data.yunits;
0055 else
0056 a = varargin{1};
0057 b = varargin{2};
0058 esttype = varargin{3};
0059 pl = varargin{4};
0060 if a.data.fs ~= b.data.fs
0061 error('### Two time-series have different sample rates.');
0062 end
0063 inunits = b.data.yunits / a.data.yunits;
0064 x = {a.data.y, b.data.y};
0065 end
0066
0067
0068 win = find(pl, 'Win');
0069 nfft = find(pl, 'Nfft');
0070 olap = find(pl, 'Olap')/100;
0071 scale = find(pl, 'scale');
0072 Xolap = floor(olap*nfft);
0073 fs = a.data.fs;
0074 order = find(pl, 'order');
0075
0076 [x,M,isreal_x,y,Ly,win,winName,winParam,noverlap,k,L,options] = ...
0077 ao.welchparse(x,esttype,win.win, Xolap, nfft, fs);
0078
0079
0080 Sxx = zeros(options.nfft,1);
0081
0082
0083 freqVectorSpecified = false; nrow = 1;
0084 if length(options.nfft) > 1,
0085 freqVectorSpecified = true;
0086 [ncol,nrow] = size(options.nfft);
0087 end
0088
0089
0090
0091
0092
0093 if freqVectorSpecified,
0094 Sxx = zeros(length(options.nfft),1);
0095 else
0096 Sxx = zeros(options.nfft,1);
0097 end
0098 range = options.range;
0099
0100 LminusOverlap = L-noverlap;
0101 xStart = 1:LminusOverlap:k*LminusOverlap;
0102 xEnd = xStart+L-1;
0103 switch esttype,
0104 case {'ms','psd'},
0105 for i = 1:k,
0106 if order < 0
0107 seg = x(xStart(i):xEnd(i));
0108 else
0109 [seg,coeffs] = ltpda_polyreg(x(xStart(i):xEnd(i)), order);
0110 end
0111 [Sxxk,w] = ao.computeperiodogram(seg,win,...
0112 options.nfft,esttype,options.Fs);
0113 Sxx = Sxx + Sxxk;
0114 end
0115 case 'cpsd',
0116 for i = 1:k,
0117 if order < 0
0118 xseg = x(xStart(i):xEnd(i));
0119 else
0120 [xseg,coeffs] = ltpda_polyreg(x(xStart(i):xEnd(i)), order);
0121 end
0122 if order < 0
0123 yseg = y(xStart(i):xEnd(i));
0124 else
0125 [yseg,coeffs] = ltpda_polyreg(y(xStart(i):xEnd(i)), order);
0126 end
0127 [Sxxk,w] = ao.computeperiodogram({xseg,...
0128 yseg},win,options.nfft,esttype,options.Fs);
0129 Sxx = Sxx + Sxxk;
0130 end
0131 case 'tfe'
0132 Sxy = zeros(options.nfft,1);
0133 for i = 1:k,
0134 if order < 0
0135 xseg = x(xStart(i):xEnd(i));
0136 else
0137 [xseg,coeffs] = ltpda_polyreg(x(xStart(i):xEnd(i)), order);
0138 end
0139 if order < 0
0140 yseg = y(xStart(i):xEnd(i));
0141 else
0142 [yseg,coeffs] = ltpda_polyreg(y(xStart(i):xEnd(i)), order);
0143 end
0144 [Sxxk,w] = ao.computeperiodogram(xseg,...
0145 win,options.nfft,esttype,options.Fs);
0146 [Syxk,w] = ao.computeperiodogram({yseg,...
0147 x(xStart(i):xEnd(i))},win,options.nfft,esttype,options.Fs);
0148 Sxx = Sxx + Sxxk;
0149 Sxy = Sxy + Syxk;
0150 end
0151 case 'mscohere'
0152
0153
0154 Sxy = zeros(options.nfft,1);
0155 Syy = zeros(options.nfft,1);
0156 for i = 1:k,
0157 if order < 0
0158 xseg = x(xStart(i):xEnd(i));
0159 else
0160 [xseg,coeffs] = ltpda_polyreg(x(xStart(i):xEnd(i)), order);
0161 end
0162 if order < 0
0163 yseg = y(xStart(i):xEnd(i));
0164 else
0165 [yseg,coeffs] = ltpda_polyreg(y(xStart(i):xEnd(i)), order);
0166 end
0167 [Sxxk,w] = ao.computeperiodogram(xseg,...
0168 win,options.nfft,esttype,options.Fs);
0169 [Syyk,w] = ao.computeperiodogram(yseg,...
0170 win,options.nfft,esttype,options.Fs);
0171 [Sxyk,w] = ao.computeperiodogram({xseg,...
0172 yseg},win,options.nfft,esttype,options.Fs);
0173 Sxx = Sxx + Sxxk;
0174 Syy = Syy + Syyk;
0175 Sxy = Sxy + Sxyk;
0176 end
0177 end
0178
0179 Sxx = Sxx./k;
0180 if any(strcmpi(esttype,{'tfe','mscohere'})),
0181 Sxy = Sxy./k;
0182
0183 if strcmpi(esttype,'mscohere'),
0184 Syy = Syy./k;
0185 end
0186 end
0187
0188
0189
0190 if ~freqVectorSpecified,
0191 w = psdfreqvec('npts',options.nfft, 'Fs',options.Fs);
0192 else
0193 if strcmpi(options.range,'onesided')
0194 warning(generatemsgid('InconsistentRangeOption'),...
0195 'Ignoring ''onesided'' option. When a frequency vector is specified, a ''twosided'' PSD is computed');
0196 end
0197 options.range = 'twosided';
0198 end
0199
0200
0201
0202
0203 [Pxx,w,units] = computepsd(Sxx,w,options.range,options.nfft,options.Fs,esttype);
0204
0205 if any(strcmpi(esttype,{'tfe','mscohere'})),
0206
0207 [Pxy,f,xunits] = computepsd(Sxy,w,options.range,options.nfft,options.Fs,esttype);
0208
0209
0210 if strcmpi(esttype,'tfe'),
0211 Pxx = Pxy ./ Pxx;
0212 end
0213
0214
0215 if strcmpi(esttype,'mscohere'),
0216
0217
0218 [Pyy,f,xunits] = computepsd(Syy,w,options.range,options.nfft,options.Fs,esttype);
0219 Pxx = (abs(Pxy).^2)./(Pxx.*Pyy);
0220 end
0221 end
0222
0223
0224 [Pxx, info] = ao.welchscale(Pxx, win, fs, scale, inunits);
0225 info.navs = k;
0226
0227 varargout = {Pxx, w, info};
0228 end
0229