function [beta,betap,sAt,Ht] = htilde2(At,Bt,L);
% [beta,betap,sAt,Ht] = htilde2(Atilde,Btilde,L);
% Purpose Developement of DMC, etc.
% For given Atilde and Btilde (see Di Ruscio 1994) 
% this function extract the (L*m x L*r) toepiliz 
% matrix H^d=beta of impulse responses.

[Lm,Lm]=size(At); % Lm=L*m 
[Lm,Lr]=size(Bt); % Lr = L*r
r=Lr/L;
% At a (Lm x Lm) matrix and Bt a (Lm x Lr) matrix

% Hankel matrix from impulse response matrices for the ESSM
w=zeros(Lm,Lr);
for i=1:L
  w = w + At^(i-1)*Bt;
  Ht((i-1)*Lm+1:i*Lm,:)=w;
end

% Construction of beta_L from upper part of Hankel matrix
beta=zeros(Lm,Lr);
sAt=zeros(Lm,Lm);
for j=1:L
   w=zeros(Lm,r);
   for i=1:L-j+1
     k=i+j-1;
     w=w+Ht((i-1)*Lm+1:i*Lm,(k-1)*r+1:k*r);
   end
   beta(:,(j-1)*r+1:j*r)=w;
   sAt=sAt+At^j;
end

% Construction of beta_p (past) from lower part of Hankel matrix
betap=zeros(Lm,(L-1)*r); bp=betap;
for j=2:L
   w=zeros(Lm,r);
   for i=1:L-j+1
     k=i+j-1;
     w=w+Ht((k-1)*Lm+1:k*Lm,(i-1)*r+1:i*r);
   end
   betap(:,(j-2)*r+1:(j-1)*r)=w;
   bp(:,(L-1)*r-(j-1)*r+1:(L-1)*r-(j-2)*r)=w;
end
betap=bp;
