%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% .m script of web demo "Insert title here" 
% Written for GNU Octave
% 
% INSTITUTE OF TELECOMMUNICATIONS  
% University of Stuttgart 
% www.inue-uni-stuttgart.de 
% author: Xin Meng and Xiaojie Wang    
% date: 10.03.2017     
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% This example script just prints the values of the given parameters.
% Make sure that the function prototype matches the prototype defined
% in the XML config file. (Parameter names are arbitrary but should be the same
%for better readability)

function [] = ErgodicCapacityvsNrequalNt(
LengthN,
beta,
SNRSet,    
withWF,                          % Channel capacity with Waterfilling
withCorr,                        % Channel capacity with Correlation
Corr,
Corrtx,
Corrrx,
otherRatio,
beta2,
filename)	                 

% define formats
landscape = "-S930,350";
portrait  = "-S640,480";

% select format
output_format = portrait;

% Create an invisible figure.
fh = figure(1);
set(fh, "visible", "off");

% define style
set(0, "defaultlinelinewidth", 2);
set(0, "defaultaxesfontsize", 8);
set(0, "defaultaxesfontname", "Helvetica");
set(0, "defaulttextfontsize", 8);
set(0, "defaulttextfontname", "Helvetica");

% This to test and demonstrate the debug engine:
fprintf("This is a debug message. Checking format: %s", output_format);

% check input parameters 
if(~(1e-2<beta && beta<11 && isscalar(beta)))
   beta=1;
end

if(~(1e-2<beta2 && beta2<11 && isscalar(beta2)))
   beta2=1;
end

if(~(isscalar(SNRSet)))
   SNRSet = 15;
end


if(~(0<Corrtx && Corrtx<1 && isscalar(Corrtx) ))
   Corrtx=0.8;
end

if(~(0<Corrrx && Corrrx<1 && isscalar(Corrrx) ))
   Corrtx=0.8;
end

% END, set up graphical output
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
%Start main function


nsample = 10000;


% determine step size 
if(LengthN<21)
  stepAntennas = 2;
elseif(LengthN<41)
  stepAntennas = 4;
  elseif(LengthN<61)
  stepAntennas = 6;
  elseif(LengthN<81)
  stepAntennas = 8;
  else
  stepAntennas = 10;
end

idx = 1; 
for k = 1:stepAntennas:LengthN
    if(beta>=1) 
      Nt = k;
      Nr = ceil(beta*k);
    else
       Nt = ceil(k/beta);
       Nr = k;
    end
     if(withCorr)
       [ErgodicValue_corr(idx),ErgodicWaterfillingValue_corr(idx),~,~] = ErgodicCapacityCorr(Nr,Nt,SNRSet,Corrtx,Corrrx,Corr,withWF,nsample);
     end
    %Waterfillinggain_corr(k)=ErgodicWaterfillingValue_corr(k)-ErgodicValue_corr(k);
    [ErgodicValue_iid(idx),ErgodicWaterfillingValue_iid(idx),~,~] = ErgodicCapacityiid(Nr,Nt,SNRSet,withWF,nsample);
    %Waterfillinggain_iid(k)=ErgodicWaterfillingValue_iid(k)-ErgodicValue_iid(k);
    if(otherRatio)
    if(beta2>=1) 
      Nt_2 = k;
      Nr_2 = ceil(beta2*k);
    else
       Nt_2 = ceil(k/beta2);
       Nr_2 = k;
    end
       [ErgodicValue_iid2(idx),ErgodicWaterfillingValue_iid2(idx),~,~] = ErgodicCapacityiid(Nr_2,Nt_2,SNRSet,withWF,nsample);
       if(withCorr)
        [ErgodicValue_corr2(idx),ErgodicWaterfillingValue_corr2(idx),~,~] = ErgodicCapacityCorr(Nr_2,Nt_2,SNRSet,Corrtx,Corrrx,Corr,withWF,nsample);
       end
       x2(idx)=Nt_2;
    end
    x(idx)=Nt;
    idx = idx +1 ;
end

figure(1);
    xlabel('Number of Tx antennas N_{t}')
    ylabel('Capacity in bit/channel use')
    %title('MIMO Capacity vs Nt=Nr')
    init_legend = ['1';'2';'3';'4';'5';'6';'7';'8'];
    legend_str = cellstr(init_legend);
    if(otherRatio)
      xlim([max(x(1),x2(1)) min(x(end), x2(end))])
    else
    xlim([x(1) x(end)])
    end
    grid on  
    hold on;

    label_legend = cellstr(init_legend)
    
plot(x,ErgodicValue_iid,"color",unicmap(1));
legend_str(1) = cellstr(['\beta =',num2str(beta),', w/o WF, w/o Corr']);
if (withWF)

    plot(x,ErgodicWaterfillingValue_iid,"color",unicmap(2),"linestyle","--");
    legend_str(2) = cellstr(['\beta =',num2str(beta),', w/ WF, w/o Corr']);
else
    legend_str(2) = cellstr('');
end

if (withCorr)
    plot(x,ErgodicValue_corr,"color",unicmap(3),"marker",'+');
    legend_str(3) = cellstr(['\beta =',num2str(beta),', w/o WF, w/ Corr']);
    hold on;
    if(withWF)
     plot(x,ErgodicWaterfillingValue_corr,"color",unicmap(4),"linestyle","--","marker",'+');
     legend_str(4) = cellstr(['\beta =',num2str(beta),', w/ WF, w/ Corr']);
    else
      legend_str(4) = cellstr('');
    end

else
    legend_str(3) = cellstr('');
    legend_str(4) = cellstr('');
end

if (otherRatio)

    plot(x2,ErgodicValue_iid2,"color",unicmap(5));
    legend_str(5) = cellstr(['\beta =',num2str(beta2),', w/o WF, w/o Corr']);
    if(withWF)
      plot(x2,ErgodicWaterfillingValue_iid2,"color",unicmap(6),"linestyle","--");
      legend_str(6) = cellstr(['\beta =',num2str(beta2),', w/ WF, w/o Corr']);
     else
      legend_str(6) = cellstr('');
    end
    if(withCorr)
        plot(x2,ErgodicValue_corr2,"color",unicmap(7),"marker",'+');
        legend_str(7) = cellstr(['\beta =',num2str(beta2),', w/o WF, w/ Corr']);
        hold on;
       if(withWF)
         plot(x2,ErgodicWaterfillingValue_corr2,"color",unicmap(8),"linestyle","--","marker",'+');
         legend_str(8) = cellstr(['\beta =',num2str(beta2),', w/ WF, w/ Corr']);
        else
          legend_str(8) = cellstr('');
       end
    else
       legend_str(7) = cellstr('');
    end
else
    legend_str(5) = cellstr('');
    legend_str(6) = cellstr('');
    legend_str(7) = cellstr('');
    legend_str(8) = cellstr('');
end

idx=1;
    for i = 1:8
    	if(~isempty(char(legend_str(i))))
    		label_legend(idx) = legend_str(i);
    		idx = idx+1;
    	end
end

    legend({char(label_legend(1)),char(label_legend(2)),char(label_legend(3)),char(label_legend(4)),char(label_legend(5)),char(label_legend(6)),char(label_legend(7)),char(label_legend(8))},"location", "northwest");
        
    hold off;

% Create an image file. This png image is displayed on the webpage. Its size 
% must be 640x480.
print(filename, "-dpng", output_format);



end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% End main function
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [ErgodicCapacity_iid,ErgodicCapacityWaterfilling_iid,ErgodicCapacityn_iid,ErgodicCapacityWaterfillingn_iid] = ErgodicCapacityiid(Nr,Nt,SNR_DB,withWF,nsample)
%nsample = 10000;        
SNR = 10^(SNR_DB/10);
ErgodicCapacityWaterfilling_iid = zeros(1,length(SNR_DB));

Hw = (randn(Nr, Nt, nsample) + sqrt(-1) * randn(Nr, Nt, nsample)) / sqrt(2);    
for i = 1:nsample
    H_iid =  Hw(:,:,i);
    ErgodicCapacityn_iid(i) = (log2(abs(det(eye(Nr)+H_iid*H_iid'*SNR/Nt)))); 
    if(withWF)
    lambda_iid = eig(H_iid * H_iid');
    [mu_iid pl_iid] = Waterfilling(lambda_iid, SNR);
    ErgodicCapacityWaterfillingn_iid(i) = sum(log2(1 + pl_iid .* lambda_iid));
    end
end
ErgodicCapacity_iid = mean(ErgodicCapacityn_iid);
 if(withWF)
ErgodicCapacityWaterfilling_iid = mean(ErgodicCapacityWaterfillingn_iid);
end
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [ErgodicCapacity_corr,ErgodicCapacityWaterfilling_corr,ErgodicCapacityn_corr,ErgodicCapacityWaterfillingn_corr] = ErgodicCapacityCorr(Nr,Nt,SNR_DB,Corrtx,Corrrx,Corr,withWF,nsample)
%nsample = 10000;        
SNR = 10^(SNR_DB/10);
ErgodicCapacityn_corr = zeros(1,nsample);
ErgodicCapacityWaterfillingn_corr = zeros(1,nsample);
ErgodicCapacityWaterfilling_corr = zeros(1,length(SNR_DB));

R_t_corr = eye(Nt);

R_r_corr = eye(Nr);

rtx =Corrtx;
rrx =Corrrx;

if Corr == 1
for i = 1:Nt
    R_t_corr_firstrow(i) = rtx^((i-1)^2);
end
R_t_corr=toeplitz(R_t_corr_firstrow);
end


if Corr == 2
for i = 1:Nr
    R_r_corr_firstrow(i) = rrx^((i-1)^2);
end
R_r_corr=toeplitz(R_r_corr_firstrow);
end

if Corr == 3
for i = 1:Nt
    R_t_corr_firstrow(i) = rtx^((i-1)^2);
end
R_t_corr=toeplitz(R_t_corr_firstrow);

for i = 1:Nr
    R_r_corr_firstrow(i) = rrx^((i-1)^2);
end
R_r_corr=toeplitz(R_r_corr_firstrow);

end


R_r_sqrt_corr = sqrtm(R_r_corr);
R_t_sqrt_corr= sqrtm(R_t_corr);


Hw = (randn(Nr, Nt, nsample) + sqrt(-1) * randn(Nr, Nt, nsample)) / sqrt(2);    
for i = 1:nsample

    H_corr = R_r_sqrt_corr' * Hw(:,:,i) * R_t_sqrt_corr;
    ErgodicCapacityn_corr(i) = (log2(abs(det(eye(Nr)+H_corr*H_corr'*SNR/Nt)))); 
    if(withWF)
     lambda_corr = eig(H_corr * H_corr');
     [mu_corr pl_corr] = Waterfilling(lambda_corr, SNR);
     ErgodicCapacityWaterfillingn_corr(i) = sum(log2(1 + pl_corr .* lambda_corr));
    end
end
ErgodicCapacity_corr = mean(ErgodicCapacityn_corr);
if(withWF)
ErgodicCapacityWaterfilling_corr = mean(ErgodicCapacityWaterfillingn_corr);
end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [mu, pl] = Waterfilling(lambda, SNR)

[lambda idx] = sort(lambda, 'descend');
lambda = lambda(find(lambda > 0)); 
pl = -1;
while (min(pl) < 0)
       mu = (SNR + sum(1 ./ lambda)) / length(lambda);
       pl = mu - 1 ./ lambda;
       lambda = lambda(1:end-1);
end
pl = [pl; zeros(length(idx) - length(pl), 1)]; 
pl(idx) = pl;
end



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [ uni_colors ] = unicmap(n)
%UNICMAP Get a matrix of colors defined by the University of Stuttgart
%   unicmap returns a matrix of 9 colors defined by the University of 
%           Stuttgart.
%   unicmap(n) returns the University color with index n. n must be in
%              [1..9].

% uni_hellgrau = [224/255 224/255 224/255];
% uni_mittelgrau = [128/255 128/255 128/255];
% uni_dunkelgrau = [80/255 80/255 80/255];    

narginchk(0, 1);

uni_colors = [0 126/255 198/255]; % mittelblau
uni_colors = [uni_colors; 238/255 28/255 35/255]; % rot
uni_colors = [uni_colors; 140/255 198/255 62/255]; % apfelgruen
uni_colors = [uni_colors; 244/255 111/255 33/255]; % orange
uni_colors = [uni_colors; 237/255 0 140/255]; % pink
uni_colors = [uni_colors; 19/255 31/255 31/255]; % anthrazit
uni_colors = [uni_colors; 0/255 170/255 173/255]; % turquoise
uni_colors = [uni_colors; 128/255 40/255 145/255]; % lila
uni_colors = [uni_colors; 1 221/255 0]; % gelb

if nargin == 1
    if (~isscalar(n) || ~isnumeric(n))
        error('First argument must be a numeric scalar.')
    end
    if (n < 1 || n > size(uni_colors, 1))
        error(['First argument mus be postive and less than or equal to '...
            num2str(size(uni_colors, 1))]);
    end
    uni_colors = uni_colors(n,:);
end

end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



% END, function definitions
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



%%%%
%function [ErgodicCapacity_iid,ErgodicCapacityWaterfilling_iid,ErgodicCapacityn_iid,ErgodicCapacityWaterfillingn_iid] = ErgodicCapacityiid(Nr,Nt,SNRSet)
%nsample = 1000;        
%SNR = 10^(SNRSet/10);
%
%Hw = (randn(Nr, Nt, nsample) + sqrt(-1) * randn(Nr, Nt, nsample)) / sqrt(2);    
%for i = 1:nsample
%    H_iid =  Hw(:,:,i);
%    ErgodicCapacityn_iid(i) = (log2(abs(det(eye(Nr)+H_iid*H_iid'*SNR/Nt)))); 
%    
%    lambda_iid = eig(H_iid * H_iid');
%    [mu_iid pl_iid] = Waterfilling(lambda_iid, SNR);
%    ErgodicCapacityWaterfillingn_iid(i) = sum(log2(1 + pl_iid .* lambda_iid));
%end
%ErgodicCapacity_iid = mean(ErgodicCapacityn_iid);
%ErgodicCapacityWaterfilling_iid = mean(ErgodicCapacityWaterfillingn_iid);
%end
%%
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%function [ErgodicCapacity_corr,ErgodicCapacityWaterfilling_corr,ErgodicCapacityn_corr,ErgodicCapacityWaterfillingn_corr] = ErgodicCapacityCorr(Nr,Nt,SNRSet)
%nsample = 1000;        
%SNR = 10^(SNRSet/10);
%ErgodicCapacityn_corr = zeros(1,nsample);
%ErgodicCapacityWaterfillingn_corr = zeros(1,nsample);
%
%R_t_corr = eye(Nt);
%R_r_corr = eye(Nr);
%
%rtx =0.9;
%rrx =0.9;
%
%for i = 1:Nt
%    R_t_corr_firstrow(i) = rtx^(i-1);
%end
%R_t_corr=toeplitz(R_t_corr_firstrow);
%for i = 1:Nr
%    R_r_corr_firstrow(i) = rrx^(i-1);
%end
%R_r_corr=toeplitz(R_r_corr_firstrow);
%
%
%R_r_sqrt_corr = sqrtm(R_r_corr);
%R_t_sqrt_corr= sqrtm(R_t_corr);
%
%Hw = (randn(Nr, Nt, nsample) + sqrt(-1) * randn(Nr, Nt, nsample)) / sqrt(2);    
%for i = 1:nsample
%
%    H_corr = R_r_sqrt_corr' * Hw(:,:,i) * R_t_sqrt_corr;
%    ErgodicCapacityn_corr(i) = (log2(abs(det(eye(Nr)+H_corr*H_corr'*SNR/Nt)))); 
%    
%    lambda_corr = eig(H_corr * H_corr');
%    [mu_corr pl_corr] = Waterfilling(lambda_corr, SNR);
%    ErgodicCapacityWaterfillingn_corr(i) = sum(log2(1 + pl_corr .* lambda_corr));
%end
%ErgodicCapacity_corr = mean(ErgodicCapacityn_corr);
%ErgodicCapacityWaterfilling_corr = mean(ErgodicCapacityWaterfillingn_corr);
%end
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%function [mu, pl] = Waterfilling(lambda, SNR)
%
%[lambda idx] = sort(lambda, 'descend');
%lambda = lambda(find(lambda > 0)); 
%pl = -1;
%while (min(pl) < 0)
%       mu = (SNR + sum(1 ./ lambda)) / length(lambda);
%       pl = mu - 1 ./ lambda;
%       lambda = lambda(1:end-1);
%end
%pl = [pl; zeros(length(idx) - length(pl), 1)]; 
%pl(idx) = pl;
%end



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
