function [phiar,sigmaar]=m5346blockar(Xblocked,norder,errorcheck)
%
% function [phiar,sigmaar]=m5346blockar(Xblocked,norder,errorcheck)
% generates an ar model from a block of data runs.  This routine can be used
% for the standard situation, based on a single run, by using a vector for Xblocked.
%
% An important feature of this implementation is the removal of means.  The
% mean is removed separately from each block, the covariance estimates are made
% for each block, and then the covariance estimates are averaged.  A standard
% Yule-Walker equation is then applied.  Note that the subtracted means are NOT returned
% as a function output.
%
% The model is X(n)=phiar(1)X(n-1)+phiar(2)X(n-2)+...+phiar(norder)X(n-norder) + Z(n)
%
% where norder=order of the model, Z(n) are IID normal N(0,sigmaar^2).
%
% Inputs:
%
% Xblocked is (length signal) x (number of runs)
%
% norder is the desired model order
%
% errorcheck set errorcheck=1 for simple error checking
%
% Outputs:
%
% phiarma:  A vector of length (order of  model) containing coefficients
%
% sigmaarma: standard deviation of white noise term in model
%


[lengthx nblocks]=size(Xblocked);

gamma=zeros([(norder+1) 1]);

% Note that we remove the mean SEPARATELY from each block.

%  Use a gamma0=gamma(1) gamma1=gamma(2) ... shift in notation form the standard
%  algorithm.

for kk=1:nblocks  % loop over blocks
   xg=Xblocked(:,kk)-mean(Xblocked(:,kk)); % remove mean
for k=1:(norder+1)
   gamma(k)=gamma(k)+mean( xg(1:(length(xg)-k+1)) .* xg(k:(length(xg)))  );
end  % end of k loop for gammas
end % end of kk loop over blocks
gamma=gamma/nblocks;

capgamma=zeros([norder norder]);  % Define capgamma matrix
for k1=1:norder
   for k2=1:norder
      capgamma(k1,k2)=gamma(abs(k1-k2)+1);
   end % end of k2
end % end of k1

phiar=inv(capgamma)*gamma(2:length(gamma));

sigmaar=sqrt(gamma(1)-phiar'*gamma(2:length(gamma)));