function [Hy,Hya,I,Imin90,Imax90,dvec,F] = MI_MIMO(Nr,Nt,T,SNR) % Author: Fredrik Rusek % 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)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;