home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l455 / 8.ddi / SIGNAL.DI$ / XCORR.M < prev    next >
Encoding:
Text File  |  1993-03-11  |  3.9 KB  |  146 lines
  1. function c = xcorr(a, b, option)
  2. %XCORR    Cross-correlation function estimates.
  3. %    XCORR(A,B), where A and B are length M vectors, returns the
  4. %    length 2*M-1 cross-correlation sequence in a column vector.
  5. %    XCORR(A), when A is a vector, is the auto-correlation sequence.
  6. %    XCORR(A), when A is an M-by-N matrix, is a large matrix with
  7. %    2*M-1 rows whose N^2 columns contain the cross-correlation
  8. %    sequences for all combinations of the columns of A.
  9. %    The zeroth lag of the output correlation is in the middle of the 
  10. %    sequence, at element or row M.
  11. %    By default, XCORR computes a raw correlation with no normalization.
  12. %    XCORR(A,'biased') or XCORR(A,B,'biased') returns the "biased"
  13. %    estimate of the cross-correlation function.  The biased estimate
  14. %    scales the raw cross-correlation by 1/M.
  15. %    XCORR(...,'unbiased') returns the "unbiased" estimate of the
  16. %    cross-correlation function.  The unbiased estimate scales the raw
  17. %    correlation by 1/(M-abs(k)), where k is the index into the result.
  18. %    XCORR(...,'coeff') normalizes the sequence so that the
  19. %    correlations at zero lag are identically 1.0.
  20. %    See also XCOV, CORRCOEF, CONV and XCORR2.
  21.  
  22. %    L. Shure 1-9-88, reviseed 4-13-92 by LS.
  23. %    Copyright (c) 1988-92 by the MathWorks, Inc.
  24.  
  25. %    References:
  26. %      [1] J.S. Bendat and A.G. Piersol, "Random Data:
  27. %          Analysis and Measurement Procedures", John Wiley
  28. %          and Sons, 1971, p.332.
  29. %      [2] A.V. Oppenheim and R.W. Schafer, Digital Signal 
  30. %          Processing, Prentice-Hall, 1975, pg 539.
  31.  
  32. onearray = 1;
  33. if nargin == 1
  34.     option = 'none';
  35.     if min(size(a)) == 1    % a is a vector
  36.         a = [a(:) a(:)];
  37.     else
  38.         onearray = 0;
  39.     end
  40. elseif nargin == 2
  41.     if isstr(b)
  42.         option = b; clear b
  43.         na = max(size(a));
  44.         if min(size(a)) == 1    % a is a vector
  45.             a = [a(:) a(:)];
  46.         else    % a is a matrix
  47.             onearray = 0;
  48.             [m,n] = size(a);
  49.         end
  50.     else    % b is truly a second arg
  51.         if min(size(a)) ~= 1 & min(size(b)) ~= 1
  52.             error('You may only specify 2 vector arrays.')
  53.         end
  54.         option = 'none';
  55.         onearray = 2;
  56.     end
  57. else
  58.     if max(size(a)) ~= max(size(b)) & ~strcmp(option,'none')
  59.         error('OPTION must be ''none'' for different length vectors A and B')
  60.     end
  61.     onearray = 2;
  62. end
  63. % check validity of option
  64. nopt = nan;
  65. if strcmp(option, 'none')
  66.     nopt = 0;
  67. elseif strcmp(option, 'coeff')
  68.     nopt = 1;
  69. elseif strcmp(option, 'biased')
  70.     nopt = 2;
  71. elseif strcmp(option, 'unbiased')
  72.     nopt = 3;
  73. end
  74. if isnan(nopt)
  75.     error('Unknown OPTION')
  76. end
  77. if onearray == 2
  78.     [ar,ac] = size(a);
  79.     na = max([ar ac]);
  80.     nb = max(size(b));
  81.     if na > nb
  82.         b(na) = 0;
  83.     elseif na < nb
  84.         a(nb) = 0;
  85.     end
  86.     a = [a(:) b(:)];
  87. end
  88. [nr, nc] = size(a);
  89. nsq  = nc^2;
  90. mr = 2 * nr - 1;
  91. c = zeros(mr,nsq);
  92. ci = zeros(1,nsq);
  93. cj = ci;
  94. nfft = 2^nextpow2(2*nr);
  95. for i = 1:nc
  96.     atmpi = a(:,i);
  97.     if ~any(any(atmpi))
  98.         real1 = 1;
  99.     else
  100.         real1 = 0;
  101.     end
  102.     atmpi = fft([atmpi(:); zeros(nfft-nr,1)]);
  103.     for j = 1:i
  104.         col = (i-1)*nc+j;
  105.         colaux = (j-1)*nc+i;
  106.         tmp = fft([a(:,j); zeros(nfft-nr,1)]); % pad with zeros for fft
  107.         tmp = fftshift(ifft(atmpi.*conj(tmp)));
  108.         c(:,colaux) = tmp((1:mr)+nfft/2-nr+1);
  109.         ci(col) = i;
  110.         cj(col) = j;
  111.         ci(colaux) = j;
  112.         cj(colaux) = i;
  113.         if ~any(any(imag(a(:,j)))) & real1
  114.             c(:,colaux) = real(c(:,colaux));
  115.         end
  116.         if i~= j
  117.             c(:,col) = conj(c(mr:-1:1,colaux));
  118.         end
  119.     end
  120. end
  121. if nopt == 1    % return normalized by sqrt of each autocorrelation at 0 lag
  122. % do column arithmetic to get correct autocorrelations
  123.     cdiv = ones(mr,1)*sqrt(c(nr,1+(ci-1)*(nc+1)).*c(nr,1+(cj-1)*(nc+1)));
  124.     c = c ./ cdiv;
  125. elseif nopt == 2    % biased result, i.e. divide by nr for each element
  126.     c = c / nr;
  127. elseif nopt == 3    % unbiased result, i.e. divide by nr-abs(lag)
  128.     c = c ./ ([1:nr (nr-1):-1:1]' * ones(1,nsq));
  129. end
  130. if onearray == 1
  131.     c = c(:,1);    % just want the autocorrelation
  132.     [am, an] = size(a);
  133.     if am == 1
  134.         c = c';
  135.     end
  136. elseif onearray == 2    % produce only cross-correlation
  137.     c = c(:,2);
  138.     if ar == 1
  139.         c = c';
  140.     end
  141. end
  142. if ~any(any(imag(a)))
  143.     c = real(c);
  144. end
  145.  
  146.