0001 function varargout = ltpda_lxspec(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
0036
0037
0038 X = varargin{1};
0039 f = varargin{2};
0040 r = varargin{3};
0041 m = varargin{4};
0042 L = varargin{5};
0043 fs = varargin{6};
0044 win = varargin{7};
0045 order = varargin{8};
0046
0047
0048 if nargout == 2
0049 computeVariance = 1;
0050 else
0051 computeVariance = 0;
0052 end
0053
0054 twopi = 2.0*pi;
0055
0056 si = size(X);
0057 nx = si(1);
0058 nc = si(2);
0059 nf = length(f);
0060 Txy = zeros(nc,nc,nf);
0061 Sig = zeros(nc,nc,nf);
0062 ENBW = zeros(nc,nc,nf);
0063
0064 switch order
0065 case -1
0066 detrender = '';
0067 case 0
0068 detrender = 'xi = xi - mean(xi);xo = xo - mean(xo);';
0069 case 1
0070 detrender = 'xi = detrend(xi);xo = detrend(xo);';
0071 otherwise
0072 detrender = 'xi = polydetrend(xi, order); xo = polydetrend(xo, order);';
0073 end
0074
0075 disp_each = round(nf/100)*10;
0076 for fi=1:nf
0077
0078
0079 l = L(fi);
0080 switch win.name
0081 case 'Kaiser'
0082 win = specwin(win.name, l, win.psll);
0083 otherwise
0084 win = specwin(win.name, l);
0085 end
0086
0087 p = 1i * twopi * m(fi)/l.*[0:l-1];
0088 C = win.win .* exp(p);
0089 Xolap = (1-win.rov/100);
0090 if mod(fi,disp_each) == 0 || fi == 1 || fi == nf
0091 disp(sprintf('++ computing frequency %03d of %03d: %f Hz', fi, nf, f(fi)));
0092 end
0093
0094
0095 nsegs = 0;
0096 start = 1;
0097 lm1 = l-1;
0098 while start + l - 1 <= nx
0099 start = floor(start + l*Xolap);
0100 nsegs = nsegs + 1;
0101 end
0102
0103
0104 idxs = zeros(nsegs,l);
0105 nsegs = 0;
0106 start = 1;
0107 while start + l - 1 <= nx
0108 idxs(1+nsegs,:) = start:start+lm1;
0109 start = floor(start + l*Xolap);
0110 nsegs = nsegs + 1;
0111 end
0112 Sxy = zeros(nsegs,1);
0113 i2 = zeros(nsegs,1);
0114 xi = zeros(l,1);
0115 xo = zeros(l,1);
0116
0117
0118
0119
0120 for j=1:nc
0121
0122 for k=1:nc
0123
0124 if j<k
0125
0126
0127 Sxx = 0;
0128 Syy = 0;
0129 Sxy(:) = 0;
0130 i2(:) = 0;
0131
0132
0133 for s=1:nsegs
0134
0135 xi(:) = X(idxs(s,:),j);
0136 xo(:) = X(idxs(s,:),k);
0137
0138
0139 eval(detrender);
0140
0141
0142 ti = C*xi;
0143 to = C*xo;
0144 cto = conj(to);
0145 Sxy(s) = cto*ti;
0146 Sxx = Sxx + ti*conj(ti);
0147 Syy = Syy + to*cto;
0148
0149 end
0150 disp(['-processed ' num2str(nsegs) ' segments']);
0151
0152 S1 = sum(win.ws);
0153 S12 = S1*S1;
0154 S2 = win.win*win.win';
0155 ENBW(j,k,fi) = fs*S2/S12;
0156 Txy(j,k,fi) = (abs(sum(Sxy))^2)/(Sxx*Syy);
0157
0158
0159 if computeVariance
0160 A = 0;
0161 scale = fs*S2;
0162 for s=1:nsegs
0163 S = (Sxy(s)/scale-Txy(j,k,fi))^2;
0164 A = A + S * i2(s);
0165 end
0166 Sig(j,k,fi) = A/sum(i2);
0167 end
0168
0169 end
0170 end
0171 end
0172 end
0173
0174
0175 aoOnes = ones(1,nf);
0176 for j=1:nc
0177 for k=1:nc
0178 if j == k
0179 Txy(j,k,:) = aoOnes;
0180 end
0181 if j>k
0182 Txy(j,k,:) = Txy(k,j,:);
0183 if computeVariance
0184 Sig(j,k,:) = Sig(k,j,:);
0185 end
0186 end
0187 end
0188 end
0189
0190
0191 varargout{1} = Txy;
0192 if nargout == 2
0193 varargout{2} = ENBW;
0194 elseif nargout == 3
0195 varargout{3} = Sig;
0196 varargout{2} = ENBW;
0197 end