function [Hy,Hya,I,Imin90,Imax90,dvec,F] = MI_MIMO_final(M,K,T,SNR)

% [Hy,Hya,I,dvec,F] = MI_MIMO_final(M,K,T,SNR)        Fredrik Rusek Nov. 2009
%
% This routine computes the MI for non-coherent MxK-MIMO setups with
% Gaussian inputs; so far only for K\geqM. Each channel
% coefficient is assumed i.i.d CN(0,1). The parameter 'T'
% is the coherence time of the channel. The input to the channel is a KxT
% block of i.i.d. CN(0,1). The noise is  i.i.d. CN(0,1).
% The SNR definition is as given in Jindal-Lozano-Marzetta ISIT'09:
% 
% Y=\sqrt(SNR/K)HA+N
%
% Outputs: 'Hy' and 'Hya' are (approximations of) the differential entropies
% h(Y) and h(Y|A), note that they are not normalized by T. 'I' is the normalized mutual information
% [h(Y)-h(Y|A)]/T, i.e. information per complex input vector.
% The rows of 'F' are the functions log2(f_k(d)), alive on the vector
% 'dvec'.
% 'Imin90' and 'Imax90' delineate the 90% confidence interval around I,
% calculated under the assmption of a Gaussian distribution with variance
% given by the mean-sample variance of the generated realizations.
%                                                                                                                                                          %   
% Variables specified inside the routine: 
% ---------------------------------------
% 'Blockl' is the number of channel realizations used in the MC-simulation.
% f_k(d) is computed on 0:'Dd':'dmax'.   'dmax' can be initialized with a
% small value since the routine increases it if needed
%                                                                                                                 


if T<K,
  error('T must at least equal K');
end
if M>K,
  error('K must be larger or equal to M');
end

s2=K/(2*SNR); % AWGN variance per real dimension

Blockl=500000; 
dmax=800;           % A small value can be used at first, the routine increases it if necessary
Dd=.25;

d=0:Dd:dmax;          % The scale must be linear due to interpolation later.

% --------- Compute all the functions log2(f_k(d)) -----------

%x=0:.025:100;                                                                                    % f_k(d) is found by integrating over x
x=0:.01:100;  

F=zeros(M,length(d));
for dpos=1:length(d),
  xm=-2*s2-(T+1-M)/2+sqrt((T+1-M)*(T+1-M)/4+d(dpos));                                            % finds x-value that maximizes argument inside exp( )
  if xm>0,
    vmx=d(dpos)*xm/2/s2/(xm+2*s2)-xm-(T-M+1)*log(xm+2*s2);                                       % computes corresponding value
  else
    vmx=-(T-M+1)*log(2*s2); 
  end
  for k=1:M,
    F(k,dpos)=(vmx)/log(2)+log2(x(2)*sum(x.^(K-M+k-1).*exp(d(dpos)*x./(2*s2*(x+2*s2))-x-(vmx)-(T-M+1)*log(x+2*s2))));  % computes log2 f_k(d)
  end
end
dvec=d;
disp('Pre-calculations done.')
pause(.1);

% ------------------ computation of log(f_k(d)) done -------------


% Compute Constant
Clog=-(M*T)*log2(pi);
for ii=1:M,
  Clog=Clog-log2(factorial(K-ii));
end

% ---------------- Compute h(Y)  and h(Y|A) ------------------

hy=0;Hya=0;
hy_stats=zeros(Blockl,1);
Hya_stats=zeros(Blockl,1);
for bl=1:Blockl
  H=randn(M,K)*sqrt(.5)+i*randn(M,K)*sqrt(.5);
  A=randn(K,T)*sqrt(.5)+i*randn(K,T)*sqrt(.5);
  W=randn(M,T)*sqrt(s2)+i*randn(M,T)*sqrt(s2);
  Y=H*A+W;
  
  d=real(flipud(eig(Y*Y')));
  Z=zeros(M,M);
 
  
  % ---------------- h(Y) ------------------
  
  for k=1:M,
    
    if d(k)<dmax,
      p1=floor(d(k)/Dd+1);                                         % d(k) lies between dvec(p1) and dvec(p2)
      p2=p1+1;
    else                                                           % Extend the functions f_k(d) to include d(k)
      disp('Needs to extend functions');
      pause(.1)
      newd=dmax+Dd:Dd:d(k)+Dd;
      Fnew=zeros(M,length(newd));
      for dpos=1:length(newd),
	for kk=1:M,
	  xm=-2*s2-(T+1-M)/2+sqrt((T+1-M)*(T+1-M)/4+newd(dpos));    
	  if xm>0,
	    vmx=newd(dpos)*xm/2/s2/(xm+2*s2)-xm-(T-M+1)*log(xm+2*s2);         
	  else
	    vmx=-(T-M+1)*log(2*s2);
	  end
	  Fnew(kk,dpos)=(vmx)/log(2)+log2(x(2)*sum(x.^(K-M+kk-1).*exp(newd(dpos)*x./(2*s2*(x+2*s2))-x-(vmx)-(T-M+1)*log(x+2*s2))));
	end
      end
      F=[F Fnew];
      dvec=[dvec newd];
      dmax=newd(end);
      disp('Extensions done')
      pause(.1);
      p1=floor(d(k)/Dd+1);        
      p2=p1+1;
    end
    
    
    for l=1:M    
      Z(l,k)=((F(l,p2)-F(l,p1))/(Dd))*(d(k)-(dvec(p1)))+F(l,p1);   % Interpolation
    end
  end
  
 
  maxZ=max(Z);
  py=sum(maxZ)+log2(det(2.^(Z-ones(M,1)*maxZ))/det(flipud(vander(d)')))+Clog-(sum(d))/2/s2/log(2);
  hy=hy-py;
  hy_stats(bl)=-py;
  
  % ------------------ h(Y|A) (by a similar MC-approach) --------------------
  % This can be replaced by a determinant based approach 
  
  AA=A*A';
  [u,Lambda,u]=svd(AA);
  la=diag(Lambda,0)';
  X=Y*A'*u;
  xt=sum(abs(X.^2),1);
  addon=M*K*log2(pi)-M*(K-T)*log2(2*pi*s2)+(sum(d))/2/s2/log(2)+sum(M*log2(la+2*s2)-xt./(2*s2*log(2).*(la+2*s2)));
  Hya=Hya+addon;
  Hya_stats(bl)=addon;
  
end
Hy=hy/Blockl;
Hya=Hya/Blockl;

% ------------- Shift differential entropies so they correspond with normalizations in the paper -----------

Hy=Hy+M*T*log2(SNR/K);
Hya=Hya+M*T*log2(SNR/K);

% ------------- Generate outputs -----------
I=Hy-Hya;
I=I/T;

I_stats=(hy_stats-Hya_stats)/T;
Delta=1.645*std(I_stats)/sqrt(Blockl);
Imin90=I-Delta;
Imax90=I+Delta;

  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  