


LTPDA_LXSPEC performs cross-spectral analysis of various forms.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
DESCRIPTION: LTPDA_LXSPEC performs cross-spectral analysis of various forms.
The function is a helper function for various higher level
functions. It is meant to be called from other functions
(e.g., ltpda_ltfe).
CALL: Y = ltpda_lxspec(X, f, r, m, L, fs, win, order, method);
[Y,Ys] = ltpda_lxspec(X, f, r, m, L, fs, win, order, method);
INPUTS: X - a matrix of input data samples
f,r,m,L - outputs of ltpda_compute_f
fs - sample rate of data vectors
win - specwin object used to window the data
order - order of detrending:
-1 - no detrending
0 - subtract mean
1 - subtract linear fit
N - subtract fit of polynomial, order N
method - A collection of 'Xpsd', 'TFE', 'Cohere'
OUTPUTS: Y - A cell array containing one NxN matrix of vectors for each
of the requested estimate types
Ys - A cell array of NxN matrices of variance estiamtes for Y
VERSION: $Id: ltpda_xspec.m,v 1.7 2007/10/08 11:28:29 hewitson Exp $
HISTORY: 28-10-2007 M Hewitson
Creation
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

0001 function varargout = ltpda_lxspec(varargin) 0002 % LTPDA_LXSPEC performs cross-spectral analysis of various forms. 0003 % 0004 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0005 % 0006 % DESCRIPTION: LTPDA_LXSPEC performs cross-spectral analysis of various forms. 0007 % The function is a helper function for various higher level 0008 % functions. It is meant to be called from other functions 0009 % (e.g., ltpda_ltfe). 0010 % 0011 % CALL: Y = ltpda_lxspec(X, f, r, m, L, fs, win, order, method); 0012 % [Y,Ys] = ltpda_lxspec(X, f, r, m, L, fs, win, order, method); 0013 % 0014 % INPUTS: X - a matrix of input data samples 0015 % f,r,m,L - outputs of ltpda_compute_f 0016 % fs - sample rate of data vectors 0017 % win - specwin object used to window the data 0018 % order - order of detrending: 0019 % -1 - no detrending 0020 % 0 - subtract mean 0021 % 1 - subtract linear fit 0022 % N - subtract fit of polynomial, order N 0023 % method - A collection of 'Xpsd', 'TFE', 'Cohere' 0024 % 0025 % OUTPUTS: Y - A cell array containing one NxN matrix of vectors for each 0026 % of the requested estimate types 0027 % Ys - A cell array of NxN matrices of variance estiamtes for Y 0028 % 0029 % VERSION: $Id: ltpda_xspec.m,v 1.7 2007/10/08 11:28:29 hewitson Exp $ 0030 % 0031 % HISTORY: 28-10-2007 M Hewitson 0032 % Creation 0033 % 0034 % 0035 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0036 0037 % Get inputs 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 % do we need to compute variance? 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 % compute DFT exponent and window 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 %-------------- Count segments 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 % now allocate indices for segments 0103 % idxs = repmat(struct('idx', zeros(1,l)), 1, nsegs); 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 %------------------ Process the channels 0118 0119 % Loop over input channels 0120 for j=1:nc 0121 % loop over output channels 0122 for k=1:nc 0123 % compute lower triangle of matrix 0124 if j<k 0125 0126 % do segments 0127 Sxx = 0; 0128 Syy = 0; 0129 Sxy(:) = 0; 0130 i2(:) = 0; 0131 0132 % Loop over segments 0133 for s=1:nsegs 0134 % get segments 0135 xi(:) = X(idxs(s,:),j); 0136 xo(:) = X(idxs(s,:),k); 0137 0138 % detrend segments 0139 eval(detrender); 0140 0141 % compute DFT 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 % End segment loop 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 %----------- compute variance 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 % End compute variance 0168 0169 end % End if j<k 0170 end % End loop over output channels 0171 end % End loop over input channels 0172 end % End loop over frequencies 0173 0174 % Fill the holes 0175 aoOnes = ones(1,nf); 0176 for j=1:nc 0177 for k=1:nc 0178 if j == k % the answer is 1 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