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