function mi=mi_2d_awgn(signal, SNRdB, N)
%MI_2D_AWGN Mutual information of a complex signal set over AWGN.
%
%   MI_2D_AWGN(SIGNAL, SNRdB) computes the mutual information
%   (information rate, in bits per channel use) achievable by the 
%   complex signal set SIGNAL at the SNR points specified in input 
%   vector SNRdB. A uniform distribution of the signal points and a 
%   complex additive white gaussian noise (AWGN) channel are assumed.
%
%   MI_2D_AWGN(SIGNAL, SNRdB, N) uses Gauss-Hermite quadrature with
%   N points to compute the mutual information. Default is N=10.

%   Raphael Le Bidan.
%   Created 12 Apr 2011.
%
%   References:
%   [1] P. E. McIllree, Channel capacity calculation for M-ary
%       N-dimensional signal sets, EE Ms Thesis, University of South
%       Australia, Feb. 1995
%   [2] Abramovitz & Stegun, Handbook of Mathematical Functions, 10th
%       printing, Dover, 1972. Formula 25.4.46
%       http://math.sfu.ca/~cbm/aands/page_890.htm
%   [3] Wikipedia, "Gauss-Hermite quadrature",
%       http://en.wikipedia.org/wiki/Gauss-Hermite_quadrature

if nargin==2
	N=10;
end

% Compute abscissas and weights for order-N Gauss-Hermite integration
Hn=HermitePoly(N);
Hnm1=HermitePoly(N-1);
xn=roots(Hn);
wn=2^(N-1)*factorial(N)*sqrt(pi)./(N^2*polyval(Hnm1,xn).^2);

% Normalize signal set to average unit energy
M=length(signal);
Es=sum(abs(signal).^2)/M;
signal=signal/sqrt(Es);

% Compute mutual information
mi=zeros(1,length(SNRdB));
for s=1:length(SNRdB)
	SNR=10.^(SNRdB(s)/10);
	z=0;
	for x1=signal
		Ew=0;
		for i=1:N
			inner_sum=0;
			for j=1:N
				sum_exp=0;
				for x2=signal
					w=xn(i)+1i*xn(j);
					sum_exp=sum_exp+exp(-abs(sqrt(SNR)*(x1-x2)+w)^2);
				end
				inner_sum=inner_sum+wn(j)*log2(sum_exp/pi);
			end
			Ew=Ew+wn(i)*inner_sum;
		end
		Ew=Ew/pi;
		z=z+Ew;
	end
	z=z/M;
	mi(s)=log2(M)-log2(pi*exp(1))-z;
end

%================================================================
% HermitePoly.m by David Terr, Raytheon, 5-10-04

% Given nonnegative integer n, compute the 
% Hermite polynomial H_n. Return the result as a vector whose mth
% element is the coefficient of x^(n+1-m).
% polyval(HermitePoly(n),x) evaluates H_n(x).


function hk = HermitePoly(n)

if n==0 
    hk = 1;
elseif n==1
    hk = [2 0];
else
    
    hkm2 = zeros(1,n+1);
    hkm2(n+1) = 1;
    hkm1 = zeros(1,n+1);
    hkm1(n) = 2;

    for k=2:n
        
        hk = zeros(1,n+1);

        for e=n-k+1:2:n
            hk(e) = 2*(hkm1(e+1) - (k-1)*hkm2(e));
        end
        
        hk(n+1) = -2*(k-1)*hkm2(n+1);
        
        if k<n
            hkm2 = hkm1;
            hkm1 = hk;
        end
        
    end
    
end