% This script calculates numerically the constrained capacity of different
% modulation schemes considering equiprobable signals, incoherent detection
% and soft output.
% Envelope detectors are considered, in which the output follows a Ricean
% distribution.
% The following modulation schemes are supported:
% OOK, 4-ASK
% The results are saved in file 'mutualInf-MOD.mat'

%   Author: André Noll Barreto, Universidade de Brasília
%   E-mail: andrebarreto@ene.unb.br
%   History:
%       v1.0 Jan 2014 - created
%
%   Copyright (C) 2014, André Noll Barreto
%
%    This program is free software: you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation, either version 3 of the License, or
%    (at your option) any later version.
%
%    This program is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.
%
%    You should have received a copy of the GNU General Public License
%    along with this program.  If not, see <http://www.gnu.org/licenses/>.


modTypes.ook = 0;
modTypes.ask4 = 4;

% parameters
mod = modTypes.ask4; % modulation scheme
EsN0dB = -10:20; % SNR



integPoints = 4000;
integConf = .999; % integral limits are defined based on integConf ratio of
                  % noise samples

test = false;

                  
mutualInf = zeros(size(EsN0dB));

switch mod
    case modTypes.ook
        A = [0,1];
        modName = 'OOK';
        
    case modTypes.ask4
        A = [0,1,2,3]/ sqrt(3.5);
        modName = '4-ASK';
end


fileName = ['mutualInf-' modName '-inc.mat'];

%scatter(real(A), imag(A))

Es = mean(abs(A).^2);
EsN0vec = 10.^(EsN0dB/10);

M = length(A);

maxr = max(A);
minr = min(A);

fprintf('Calculating constrained capacity for %s:\n', modName)

for EsN0ind = 1:length(EsN0vec)
    
    fprintf('\tEs/N0 = %d dB\n', EsN0dB(EsN0ind));
    
    EsN0 = EsN0vec(EsN0ind);
    N0 = Es/EsN0;
    sigma = sqrt(N0/2);
    
    % supondo distribuição gaussiana = Rice
    integLim = icdf('normal',integConf, 0, sigma);
    
    integStep = (maxr + integLim) / (integPoints + 1);
    minrtmp = 0;
    maxrtmp = ceil((integLim+maxr)/integStep)*integStep;
    intervalr = minrtmp:integStep:maxrtmp;
    
    prA = 1/M;
              
    % pre-calculate p_y   
    fprintf('\t\tpre-calculating p_y\n');
    tic
    nelr = length(intervalr);

    py_given_Ak = zeros(M, nelr);
    for m = 1:M
        x = intervalr;
        mu = A(m);
        py_given_Ak(m,:) = (x / sigma^2) .* ...
                           exp(-(x.^2 + mu^2)/2/sigma^2) .* ...
                           besseli(0, x*mu/sigma^2);        
    end
    
    py = sum(py_given_Ak) / M;
        
    for m = 1:length(A)
              
        fprintf('\t\t\tcalculating integral\n');           
        y = py_given_Ak(m,:) .* log2(py_given_Ak(m,:) ./ py);
        y(isnan(y)) = 0;
        integ = trapz(intervalr, y);
               
        mutualInf(EsN0ind) = mutualInf(EsN0ind) + prA * integ ;
        
    end    
    
    if any(isnan(mutualInf)) || ...
           (EsN0ind > 1 && mutualInf(EsN0ind) < mutualInf(EsN0ind-1)) 
        mutualInf(EsN0ind:end) = log2(M);
        break
    end       
        
    fprintf('\t mutual information is %d\n', mutualInf(EsN0ind) )
    
end

if ~test
    save (fileName, 'EsN0dB', 'mutualInf')
end


