function [Hy,Hya,I,Imin90,Imax90,dvec,F] = MI_MIMO(Nr,Nt,T,SNR)

% Author: Fredrik Rusek <fredrik.rusek@eit.lth.se>
% Electrical and Information Technology, Lund University, Sweden
% This version: May 2010
%
% [Hy,Hya,I,Imin90,Imax90,dvec,F] = MI_MIMO(Nr,Nt,T,SNR)
%
% This routine computes the MI for non-coherent block faded NtxNr-MIMO setups for
% Gaussian inputs. 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 NtxT
% block of i.i.d. CN(0,1). The noise is a NrxT block of i.i.d. CN(0,1).
% The SNR definition is as given in Jindal-Lozano-Marzetta ISIT'09:
% 
% Y=\sqrt(SNR/Nt)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 (bits) 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.
%
% The routine may fail and return NaN for scenarios of simultaneously large T, large Nr and Nt,
% and high SNR, but it delivers reliable results for most scenarios of interest, and in particular
% for small-to-moderate T.
%                                                                                                                                                          %   
% 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
%                                                


gamma=SNR/Nt;

s2=Nt/(2*SNR);               % AWGN variance per real dimension, this is used for h(Y|X)

Blockl=30000; 
dmax=1000;           % A small value can be used at first, the routine increases it if necessary
Dd=1;

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

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

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

F=zeros(min(Nr,Nt),length(d));
for dpos=1:length(d),
 
  xm=(sqrt(d(dpos)*gamma)-1)/gamma;
  if xm>0,
    vmx=d(dpos)*xm*gamma/(xm*gamma+1)-xm;
  else
    vmx=0; 
  end
  for k=1:min(Nr,Nt),
    F(k,dpos)=(vmx)/log(2)+log2(x(2)*sum((x.^(max(0,Nt-Nr))).*(x/gamma).^(k-1).*exp(d(dpos)*gamma*x./(x*gamma+1)-x-(vmx)-(T-Nr+1)*log(x*gamma+1))));  % 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=-(Nr*T)*log2(pi);
for ii=1:Nt,
  Clog=Clog-log2(factorial(Nt-ii));
end
Clog=Clog-max(0,Nr-Nt)*log2(gamma);


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

hy=0;Hya=0;
hy_stats=zeros(Blockl,1);hy_stats2=zeros(Blockl,1);
Hya_stats=zeros(Blockl,1);
for bl=1:Blockl
  H=randn(Nr,Nt)*sqrt(.5)+i*randn(Nr,Nt)*sqrt(.5);
  A=randn(Nt,T)*sqrt(.5)+i*randn(Nt,T)*sqrt(.5);
  W=randn(Nr,T)*sqrt(.5)+i*randn(Nr,T)*sqrt(.5);
  Y=sqrt(SNR/Nt)*H*A+W;
  
  d=real(flipud(eig(Y*Y')));
  Z=zeros(Nr,Nr);
 
  
  % ---------------- h(Y) ------------------
  
  for k=1:Nr,
    
    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(min(Nr,Nt),length(newd));
      for dpos=1:length(newd),
	for kk=1:min(Nr,Nt),
	  xm=(sqrt(newd(dpos)*gamma)-1)/gamma;
	  if xm>0,
	    vmx=newd(dpos)*xm*gamma/(xm*gamma+1)-xm;
	  else
	    vmx=0;
	  end
	  Fnew(kk,dpos)=(vmx)/log(2)+log2(x(2)*sum(x.^(max(0,Nt-Nr)).*(x/gamma).^(k-1).*exp(newd(dpos)*gamma*x./(x*gamma+1)-x-(vmx)-(T-Nr+1)*log(x*gamma+1))));
	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:min(Nt,Nr)    
      Z(k,l)=((F(l,p2)-F(l,p1))/(Dd))*(d(k)-(dvec(p1)))+F(l,p1);   % Interpolation
    end
    for l=Nt+1:Nr
      Z(k,l)=log2(d(k))*(l-Nt-1);
    end
  end
  
  maxZv=zeros(1,Nr);
  maxZ=max(Z);
  Z=(Z-ones(Nr,1)*maxZ);
  
  % ---------- An attempt to fix ill-conditioned Z-matrix ----------
  if  isfinite(log2(abs(det(2.^Z))/det((vander(d)))))==0,
    if Nr>Nt,
      for kk=2:Nr,
	maxZv(kk)=min(Z(kk,:))/2;
	Z(kk,:)=Z(kk,:)+abs(min(Z(kk,:)))/2;
      end
    else
       for kk=2:Nr,
	 maxZv(kk)=min(Z(kk,:));
	 Z(kk,:)=Z(kk,:)+abs(min(Z(kk,:)));
      end
    end
  end
  
  % ----------------------------------------------------------------
  
  py=sum(maxZ)+sum(maxZv)+log2(abs(det(2.^Z))/det((vander(d))))+Clog-(sum(d))/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/sqrt(gamma);
  xt=sum(abs(X.^2),1);
  addon=Nr*Nt*log2(pi)-Nr*(Nt-T)*log2(2*pi*s2)+(sum(d)/gamma)/2/s2/log(2)+sum(Nr*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 conditional differential entropy so it correspond with normalizations in the paper -----------

Hya=Hya+Nr*T*log2(SNR/Nt);

% ------------- 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;
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  