function r = cross_corr ( n, x, y, m )
%*****************************************************************************80
%
%% cross_corr() computes the sample cross correlation between two signal samples.
%
% Discussion:
%
% The sample cross correlation is defined, for 0 <= i < N, as
%
% R(i) = 1/N * sum ( 0 <= j <= N - 1 - i ) X(i+j) * Y(j)
%
% The sample cross correlation is an estimate of the cross
% correlation function.
%
% It is usually the case that the signals X and Y are assumed to
% have zero mean. Here, we compute the means and adjust the
% calculation accordingly:
%
% R(i) = 1/N * sum ( 0 <= j <= N - 1 - i )
% ( X(i+j) - Xbar ) * ( Y(j) - Ybar )
%
% Experience suggests that only the first 5 or 10 percent of
% the lags are statistically reliable, so that one might choose
% M = N / 20 or M = N / 10, for instance.
%
% Licensing:
%
% This code is distributed under the MIT license.
%
% Modified:
%
% 20 June 2010
%
% Author:
%
% John Burkardt
%
% Reference:
%
% Sophocles Orfanidis,
% Introduction to Signal Processing,
% Prentice-Hall, 1995,
% ISBN: 0-13-209172-0,
% LC: TK5102.5.O246.
%
% Input:
%
% integer N, the number of equally spaced signal samples.
%
% real X(N), Y(N), the signal samples.
%
% integer M, the maximum lag to consider.
% 0 <= M < N.
%
% Output:
%
% real R(1:M+1), the sample correlations.
%
r(1:m+1) = 0.0;
xbar = sum ( x(1:n) ) / n;
ybar = sum ( y(1:n) ) / n;
for i = 0 : m
for j = 0 : n - i - 1
r(i+1) = r(i+1) + ( x(i+j+1) - xbar ) * ( y(j+1) - ybar );
end
end
r(1:m+1) = r(1:m+1) / n;
return
end