


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_lxspec.html,v 1.9 2008/03/26 18:02:09 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_lxspec.html,v 1.9 2008/03/26 18:02:09 hewitson Exp $ 0030 % 0031 % HISTORY: 28-10-2007 M Hewitson 0032 % Creation 0033 % 0034 % 0035 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0036 0037 VERSION = '$Id: ltpda_lxspec.html,v 1.9 2008/03/26 18:02:09 hewitson Exp $'; 0038 CATEGORY = 'Internal'; 0039 0040 0041 % Get inputs 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 % do we need to compute variance? 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 % compute DFT exponent and window 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 %-------------- Count segments 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 % now allocate indices for segments 0107 % idxs = repmat(struct('idx', zeros(1,l)), 1, nsegs); 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 %------------------ Process the channels 0122 0123 % Loop over input channels 0124 for j=1:nc 0125 % loop over output channels 0126 for k=1:nc 0127 % compute lower triangle of matrix 0128 if j<k 0129 0130 % do segments 0131 Sxx = 0; 0132 Syy = 0; 0133 Sxy(:) = 0; 0134 i2(:) = 0; 0135 0136 % Loop over segments 0137 for s=1:nsegs 0138 % get segments 0139 xi(:) = X(idxs(s,:),j); 0140 xo(:) = X(idxs(s,:),k); 0141 0142 % detrend segments 0143 eval(detrender); 0144 0145 % compute DFT 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 % End segment loop 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 %----------- compute variance 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 % End compute variance 0172 0173 end % End if j<k 0174 end % End loop over output channels 0175 end % End loop over input channels 0176 end % End loop over frequencies 0177 0178 % Fill the holes 0179 aoOnes = ones(1,nf); 0180 for j=1:nc 0181 for k=1:nc 0182 if j == k % the answer is 1 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