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 VERSION  = '$Id:$';
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

Generated on Mon 31-Mar-2008 13:54:54 by m2html © 2003