function [lag r]= ser_corr_fcn(x, n, cor, bi);
% r = scf(x, n, cor, bi);
% if cor = 1, calculates the serial correlation function
% i.e., the sample autocorrelation function, of x, for
% up to n lags. n should be < N/4, where N is the length
% of x.  If cor = 0, the serial covariance function is
% calculated: bi = 1 means the biased estimate is 
% calculated bi = 0 (the default) means the unbiased
% estimate is calculated.  Ref: articles by O.D.Anderson
% in the Encyclopedia of Statistics.  
% Written by Gerry Middleton, November 1996
N = length(x);
xbar = mean(x);
xd = x - xbar;
ss = xd'*xd;
for i = 1:n+1
   xdlag = x(1:N+1-i) - xbar;
   if cor == 1   
      r(i) = xd(i:N)'*xdlag / ss;
   elseif bi == 1
      r(i) = xd(i:N)'*xdlag / N;
   else
      r(i) = xd(i:N)'*xdlag / (N-n);
   end
end
lag = [0:n];
    