Home > m > sigproc > frequency_domain > ltpda_lxspec.m

ltpda_lxspec

PURPOSE ^

LTPDA_LXSPEC performs cross-spectral analysis of various forms.

SYNOPSIS ^

function varargout = ltpda_lxspec(varargin)

DESCRIPTION ^

 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


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

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

Generated on Fri 02-Nov-2007 19:39:27 by m2html © 2003