function H = blkhankel(C,R)
%_________________________________________________________________________________________
%
%   function H = blkhankel(C,R)
%
%   Construct a block Hankel matrix, H, from the data matrices, C and R.   
%   The first block column is C and the last block row is R. 
%
%   Note that this means that the last block of C and the first block of 
%   R specify the same block in H.  If these are not equal the block from
%   C takes precedence.  This is also the case in the matlab hankel command.
%   Unlike the matlab hankel command this function also issues a warning.
%   
%   Both arguments are required in order to determine the row dimension of
%   the blocks.
%
%   Roy Smith,   8 November 2012.   Modified to use a calling syntax similar
%                to the matlab hankel formulation.
%
%_________________________________________________________________________________________

if nargin == 0,
  disp('Syntax: H = blkhankel(C,R)');
  return
elseif nargin ~= 2,
  error('Syntax: H = blkhankel(C,R)');
  return
end

[nrC,nc] = size(C);     % determine block column dimension (nc)
[nr,ncR] = size(R);     % determine block row dimension (nr)

p = nrC/nr;             % number of block rows
q = ncR/nc;             % number of block columns

%  Check for matrix dimension consistency.

if mod(p,1) ~= 0,
  error('Number of rows in column data is not an integer multiple of block row dimension')
end
if mod(q,1) ~= 0,
  error('Number of columns in row data is not an integer multiple of block column dimension')
end

if norm(C((p-1)*nr+1:p*nr,:) - R(:,1:nc),'fro') > sqrt(eps),
  warning('Data conflict: bottom left block - block row data is overwritten');
end

if p == 1,         % we have only a row
  H = R;
  H(:,1:nc) = C;   % but C must overwrite the first element by definition,
elseif q == 1,      % we have only a column
  H = C;
else
  H = zeros(nr*p,nc*q);
  H(:,1:nc) = C;
  for j = 1:q-1,
    H(1:(p-1)*nr,j*nc+1:(j+1)*nc) = H(nr+1:p*nr,(j-1)*nc+1:j*nc);
    H((p-1)*nr+1:p*nr,j*nc+1:(j+1)*nc) = R(:,j*nc+1:(j+1)*nc);
  end
end

%_________________________________________________________________________________________




