function H = blkHankel(y,r,kk)
% H = blkHankel(y,r,k)
%  Return a block Hankel matrix filled with the data y containing k block-rows 
%     
%  INPUT       DESCRIPTION                                     DIMENSION
%  =====       ===========                                     ==========
%   y          a sequence of vectors, e.g., time sequences     m x N
%   r          number of columns in each block of y 
%   k          number of block rows of H
%
%  OUTPUT      DESCRIPTION                                     DIMENSION
%  ======      ===========                                     ==========
%   H          block Hankel matrix of the output               m*k x r*(N+1-k)

  [m,N] = size(y);
  if (N < m), y = y'; [m,N] = size(y); end % time sequences are the _rows_ of y
  N = N/r;

% Check dimensions
  if (kk < 0),       error('blkHankel: k must be positive');  end
  jj = N+1-kk;       % number of block columns in H
  if (m*kk > r*jj),  error('blkHankel: # rows > # columns,  reduce k');  end

% fill the block-Hankel matrix
  H = zeros(m*kk,r*jj);
  for ii=1:kk
    H(m*(ii-1)+1:m*ii,:) = y(:,r*(ii-1)+1:r*(ii+jj-1));
  end

% ----------------------------------------------------------------- blkHankel
