% WEBDEMO .m-file (Octave), Institute of Telecommunications, University of Stuttgart, Germany
% AUTHOR:   Nan Zhao

function [] = mse(
    QAMorder,                   % bpsk or qpsk
    nChannel,                   % number of channel
    nSymbol,			% number of symbol
    nTx,                        % number of transmit antenna
    nRx,                        % number of receive antenna
    linear_theory_en,           % enable theory curv of ZF/MMSE
    ZF_en,                      % enable ZF
    MMSE_en,                    % enable MMSE
    ZF_SIC_en,                  % enable ZF with unordered SIC
    MMSE_SIC_en,                % enable MMSE with  SIC
    snrStart,                   % in dB, start value of SNR
    snrStep,                    % in dB, step value of SNR
    snrEnd,                     % in dB, end value of SNR
    svg_en,                     % if 1, then also output svg-file
    webdemo_en,                 % if 0, then test-mode in Octave/Matlab commandline, no need for webserver running
    filename                    % output file for webdemo
)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% BEGIN, set up graphical output

% Create an invisible figure.
if (webdemo_en)
fig = figure(1); set(fig, "visible", "off");
else
fig = figure(1); set(fig, "visible", "on");
end

landscape = "-S930,350";
portrait  = "-S640,480";

output_format = portrait;     % default

set(0, "defaultlinelinewidth", 2);

  if (webdemo_en)
    set(0, "defaultaxesfontsize", 8);
    set(0, "defaultaxesfontname", "Helvetica");
    set(0, "defaulttextfontsize", 8);
    set(0, "defaulttextfontname", "Helvetica");
    const_lw = 1;
    const_msize = 4;
  elseif (svg_en)
    set(0, "defaultaxesfontsize", 6);
    set(0, "defaultaxesfontname", "Helvetica");
    set(0, "defaulttextfontsize", 6);
    set(0, "defaulttextfontname", "Helvetica");
    const_lw = 1;
    const_msize = 3;
  else
    set(0, "defaultaxesfontsize", 10);
    set(0, "defaultaxesfontname", "Helvetica");
    set(0, "defaulttextfontsize", 8);
    set(0, "defaulttextfontname", "Helvetica");
    const_lw = 2;
    const_msize = 2;
  end

% Set default color order for multiline plots
set(0, 'DefaultAxesColorOrder', unicmap);

% END, set up graphical output
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% BEGIN check some input parameters for consistency

userWarning = 0;
warnStr=['Warning: ',"\n"];
if(nChannel>1000000)
    nChannel = 100000;
    %error('too much channel.');
    warnStr= [warnStr "reset the number of channel to 100000\n"];
    userWarning=1;
end

if(nRx<nTx)
    nRx = nTx;
    %error('receive antennas cannot be smaller than transmit antennas');
    warnStr= [warnStr "reset the number of receive antennas equal to the number of transmit antennas\n"];
    userWarning=1;
end

% END
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% main function

    %  init parameters 
    nSNR = [snrStart:snrStep:snrEnd];

    % pre allocate space for parameters
    % scalar values
    mse_ZF_theory = zeros(size(nChannel));
    mse_MMSE_theory = zeros(size(nChannel));
    mse_ZF = zeros(size(nChannel));
    mse_MMSE = zeros(size(nChannel));
    mse_ZF_UOSIC = zeros(size(nChannel));
    mse_ZF_OSIC = zeros(size(nChannel));
    mse_MMSE_UOSIC = zeros(size(nChannel));
    mse_MMSE_OSIC = zeros(size(nChannel));

    % dB values
    MSE_ZF_theory_dB = zeros(size(nSNR));
    MSE_MMSE_theory_dB = zeros(size(nSNR));
    MSE_ZF_dB = zeros(size(nSNR));
    MSE_MMSE_dB = zeros(size(nSNR));
    MSE_ZF_UOSIC_dB = zeros(size(nSNR));
    MSE_ZF_OSIC_dB = zeros(size(nSNR));
    MSE_MMSE_UOSIC_dB = zeros(size(nSNR));
    MSE_MMSE_OSIC_dB = zeros(size(nSNR));

    
    % loop caculation
    for ii = 1:length(nSNR)

        parfor jj = 1:nChannel
            % y = W*r = W*h*s+W*n
            % r is the received undetected symbol
            % y is the detected, but not hard decided symbol, which is soft value
            % s_tx is the transmitted symbol
            [s_tx,sigma_n,h,r] = generate_MIMO_signal(nSNR(ii),nTx,nRx,nSymbol,QAMorder);
            [s_est_soft_ZF,w_ZF] = ZF(h,r);
            [s_est_soft_MMSE,w_MMSE] = MMSE(h,r,nTx,sigma_n);
            if(linear_theory_en)
                mse_ZF_theory(jj) = sigma_n*(real(trace(w_ZF*w_ZF')));

                d_MMSE = diag(diag(w_MMSE*h));
                z_MMSE = w_MMSE*h-d_MMSE;
                mse_MMSE_theory(jj) = real(trace(z_MMSE'*z_MMSE))+sigma_n*(real(trace(w_MMSE*w_MMSE')));
            end

            if(ZF_en)
                mse_ZF(jj) = (norm(s_est_soft_ZF-s_tx))^2/nSymbol;
            end

            if(MMSE_en)
                mse_MMSE(jj) = (norm(s_est_soft_MMSE-s_tx))^2/nSymbol;
            end

            if(ZF_SIC_en)
                s_est_soft_ZF_UOSIC = SIC(h,r,sigma_n,nTx,0,0,QAMorder);
                s_est_soft_ZF_OSIC = SIC(h,r,sigma_n,nTx,0,1,QAMorder);
                mse_ZF_UOSIC(jj) = (norm(s_est_soft_ZF_UOSIC-s_tx))^2/nSymbol;
                mse_ZF_OSIC(jj) = (norm(s_est_soft_ZF_OSIC-s_tx))^2/nSymbol;
            end    
            
            if(MMSE_SIC_en)
                s_est_soft_MMSE_UOSIC = SIC(h,r,sigma_n,nTx,1,0,QAMorder);
                s_est_soft_MMSE_OSIC = SIC(h,r,sigma_n,nTx,1,1,QAMorder);
                mse_MMSE_UOSIC(jj) = (norm(s_est_soft_MMSE_UOSIC-s_tx))^2/nSymbol;
                mse_MMSE_OSIC(jj) = (norm(s_est_soft_MMSE_OSIC-s_tx))^2/nSymbol;
            end

        end

        % conver the scalar value to dB value
        if(linear_theory_en)
            mse_ZF_theory_dB(ii) = 10*log10(mean(mse_ZF_theory/nTx));
            mse_MMSE_theory_dB(ii) = 10*log10(mean(mse_MMSE_theory/nTx));
        end
        
        if(ZF_en)
            mse_ZF_dB(ii) = 10*log10(mean(mse_ZF/nTx));
        end
        
        if(MMSE_en)
            mse_MMSE_dB(ii) = 10*log10(mean(mse_MMSE/nTx));
        end
	    
        if(ZF_SIC_en)
            mse_ZF_UOSIC_dB(ii) = 10*log10(mean(mse_ZF_UOSIC/nTx));
            mse_ZF_OSIC_dB(ii) = 10*log10(mean(mse_ZF_OSIC/nTx));
        end

        if(MMSE_SIC_en)
            mse_MMSE_UOSIC_dB(ii) = 10*log10(mean(mse_MMSE_UOSIC/nTx));
            mse_MMSE_OSIC_dB(ii) = 10*log10(mean(mse_MMSE_OSIC/nTx));
        end

    end

    % set the properity of figure
    clf
    figure(1);
    hold on;
    grid on;
    xlabel('SNR in dB');
    ylabel('MSE in dB');
    title(['MSE with ' num2str(nTx) ' x ' num2str(nRx) ' MIMO (Rayleigh channel)']);

    % init legend
    init_legend = ['1';'2';'3';'4';'5';'6';'7';'8'];
    legend_str = cellstr(init_legend);
    label_legend = cellstr(init_legend);

    % plot the checked curve
    if(linear_theory_en)
        plot(nSNR,mse_ZF_theory_dB,'--*','lineWidth',const_lw,'markersize',const_msize,'color',unicmap(9));
        plot(nSNR,mse_MMSE_theory_dB,'--*','lineWidth',const_lw,'markersize',const_msize,'color',unicmap(4));
        legend_str(1) = cellstr('ZF theory');
        legend_str(2) = cellstr('MMSE theory');
    else
        legend_str(1) = cellstr('');
        legend_str(2) = cellstr('');
    end

    if(ZF_en)
        plot(nSNR,mse_ZF_dB,'lineWidth',const_lw,'markersize',const_msize,'color',unicmap(1));    
        legend_str(3) = cellstr('ZF');
    else
        legend_str(3) = cellstr('');
    end

    if(MMSE_en)
        plot(nSNR,mse_MMSE_dB,'lineWidth',const_lw,'markersize',const_msize,'color',unicmap(2));    
        legend_str(4) = cellstr('MMSE');
    else
        legend_str(4) = cellstr('');
    end

    if(ZF_SIC_en)
        plot(nSNR,mse_ZF_UOSIC_dB,'-o','lineWidth',const_lw,'markersize',const_msize,'color',unicmap(5));
        plot(nSNR,mse_ZF_OSIC_dB,'--x','lineWidth',const_lw,'markersize',const_msize,'color',unicmap(6));
        legend_str(5) = cellstr('ZF with unordered SIC');
        legend_str(6) = cellstr('ZF with ordered SIC');
    else
        legend_str(5) = cellstr('');
        legend_str(6) = cellstr('');
    end

    if(MMSE_SIC_en)
        plot(nSNR,mse_MMSE_UOSIC_dB,'-o','lineWidth',const_lw,'markersize',const_msize,'color',unicmap(7));
        plot(nSNR,mse_MMSE_OSIC_dB,'--x','lineWidth',const_lw,'markersize',const_msize,'color',unicmap(8));
        legend_str(7) = cellstr('MMSE with unordered SIC');
        legend_str(8) = cellstr('MMSE with ordered SIC');
    else
        legend_str(7) = cellstr('');
        legend_str(8) = cellstr('');
    end

    % plot the legend if it's checked
    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))});

    hold off;

    % Warning message of invalid parameter settings
    if (userWarning)
     text(0.05, 0.5, warnStr, 'Color', 'red', 'FontSize', 7,...
            'Units', 'Normalized', 'HorizontalAlignment', 'left');
    end 
    
    % Create an image file. This image is displayed on the webpage. Its size
    % must be 640x480.
    % print(filename, "-dgif","-S640,480", "-r0");
    if (webdemo_en)
        print(filename, "-dpng",output_format);
    end

    %print(filename, "-djpg","-S1280,960", "-r0");
    if (svg_en)
    print(["fig_",filename], "-dsvg",output_format);
    end

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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% channel simulator
    function [ s_tx,sigma_n,h,r ] = generate_MIMO_signal(snr,nTx,nRx,nSymbol,QAMorder)

        sigma_n = 10^(-snr/10);

        if(QAMorder == 1)
            s_constellation = [1 -1];
            s_tx = randsrc(nTx,nSymbol,s_constellation);
        end

        if(QAMorder == 2)
    	    s_constellation = [(1+1i*1)/sqrt(2) (-1+1i*1)/sqrt(2) (1-1i*1)/sqrt(2) (-1-1i*1)/sqrt(2)];
    	    s_tx = randsrc(nTx,nSymbol,s_constellation);
        end

        if(QAMorder == 3)
            s_constellation = [(1+1i*1)/sqrt(10) (1-1i*1)/sqrt(10) (1+1i*3)/sqrt(10) (1-1i*3)/sqrt(10)...
                   (-1+1i*1)/sqrt(10) (-1-1i*1)/sqrt(10) (-1+1i*3)/sqrt(10) (-1-1i*3)/sqrt(10)...
                   (3+1i*1)/sqrt(10) (3-1i*1)/sqrt(10) (3+1i*3)/sqrt(10) (3-1i*3)/sqrt(10)...
                   (-3+1i*1)/sqrt(10) (-3-1i*1)/sqrt(10) (-3+1i*3)/sqrt(10) (-3-1i*3)/sqrt(10)];
            s_tx = randsrc(nTx,nSymbol,s_constellation);
        end

        h = (normrnd(0,1/sqrt(2),nRx,nTx) + 1i.*normrnd(0,1/sqrt(2),nRx,nTx))/sqrt(nTx);
        n = normrnd(0,sqrt(sigma_n./2),nRx,nSymbol) + 1i.*normrnd(0,sqrt(sigma_n./2),nRx,nSymbol);
        r = h*s_tx+n;
        
    end

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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% ZF, get the soft value after ZF detection
    function [ s_est_soft_ZF, w_ZF ] = ZF(h,r)
            w_ZF = (h'*h)\h';
            s_est_soft_ZF = w_ZF*r;
    end

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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% MMSE, get the soft value after MMSE detection
    function [ s_est_soft_MMSE, w_MMSE ] = MMSE(h,r,nTx,sigma_n)
            w_MMSE = (h'*h+sigma_n*eye(nTx))\h';
            s_est_soft_MMSE = w_MMSE*r;
    end

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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SIC
    function [ s_est_soft_SIC ] = SIC(h,r,sigma_n,nTx,useMMSE,activeOrdering,QAMorder)

      % h is the channel matrix
      % r is the received undemapped signal
      
      mmse_enable = useMMSE;
      idxes_tobeDetected = 1:nTx;
      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      % first iteration
      % calculate the product of receive matrix W's first row and the signal
      % vector r, obtaining the soft-detected value of first layer 
      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

      w_1 = (h'*h+sigma_n*eye(nTx)*mmse_enable)\h';
      
      % select strongst layer to detect
      if(activeOrdering)
	  streamIdx_1 = get_first_SIC_stream_todetect(w_1,h,sigma_n);
      else
	  streamIdx_1 = 1;
      end
      w1_1 = w_1(streamIdx_1,:);
      y_1 = w1_1*r;

      detected_stream_index = streamIdx_1;
      idxes_tobeDetected(streamIdx_1) = [];

      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      % second iteration
      % h2 = [h2,h3,...] h2,h3 are the ranks of H
      % substracte the detected signal and its corresponding parts from the
      % channel matrix h and the signal vector r.
      % if the <activeOrdering> is actived, the channel matrix h should be sorted
      % according to the SINR of different layers.
      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

      h1 = h(:,streamIdx_1);
      h_2 = h;
      h_2(:,streamIdx_1) = [];
      w_2 = (h_2'*h_2+sigma_n*eye(nTx-1)*mmse_enable)\h_2';
      
      % select strongst layer to detect
      if(activeOrdering)
	  streamIdx_2 = get_first_SIC_stream_todetect(w_2,h_2,sigma_n);
      else
	  streamIdx_2 = 1;
      end
      
      s_est_hard_1 = HardDecisionUnit(y_1,QAMorder);
      w1_2 = w_2(streamIdx_2,:);
      r_2 = r-h1*s_est_hard_1;
      y_2 = w1_2*r_2;

      s_est_soft_SIC = [y_1;y_2]
      
      detected_stream_index = [detected_stream_index;idxes_tobeDetected(streamIdx_2)];
      idxes_tobeDetected(streamIdx_2) = [];

      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      % third iteration,repeat the process of 2nd iteration
      % H3 = [h3,h4...]
      % w_3 = inv(H3+sigma_n*I(2))
      % r_3 = r_2-h2*s^_2  h2 is the 2nd rank of channel matrix
      % y_3 = w1_3*r_3     w1_3 is the first row of w_3
      % s^_3 is the result of hard decision from y_4
      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

      % nTx >= 3
      if size(h,2) == 3 || size(h,2) == 4
	  h2 = h_2(:,streamIdx_2);
	  h_3 = h_2;
	  h_3(:,streamIdx_2) = [];
	  w_3 = (h_3'*h_3+sigma_n*eye(nTx-2)*mmse_enable)\h_3';
	  
	  % select strongst layer to detect
	  if(activeOrdering)
	      streamIdx_3 = get_first_SIC_stream_todetect(w_3,h_3,sigma_n);
	  else
	      streamIdx_3 = 1;
	  end
	  
	  s_est_hard_2 = HardDecisionUnit(y_2,QAMorder);
	  w1_3 = w_3(streamIdx_3,:);
	  r_3 = r_2-h2*s_est_hard_2;
	  y_3 = w1_3*r_3;
	  

	  s_est_soft_SIC = [s_est_soft_SIC;y_3];
	  detected_stream_index = [detected_stream_index;idxes_tobeDetected(streamIdx_3)];
	  idxes_tobeDetected(streamIdx_3) = [];     
      end

      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      % fourth iteration, simply repeat previous step
      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

      if size(h,2) == 4
	  h3 = h_3(:,streamIdx_3);
	  h_4 = h_3;
	  h_4(:,streamIdx_3) = [];
	  w_4 = (h_4'*h_4+sigma_n*eye(nTx-3)*mmse_enable)\h_4';
	  
	  % select strongst layer to detect
	  if(activeOrdering)
	      streamIdx_4 = get_first_SIC_stream_todetect(w_4,h_4,sigma_n);
	  else
	      streamIdx_4 = 1;
	  end
	  
	  s_est_hard_3 = HardDecisionUnit(y_3,QAMorder);
	  w1_4 = w_4(streamIdx_4,:);
	  r_4 = r_3-h3*s_est_hard_3;
	  y_4 = w1_4*r_4;

	  s_est_soft_SIC = [s_est_soft_SIC;y_4];
	  detected_stream_index = [detected_stream_index;idxes_tobeDetected(streamIdx_4)];
	  idxes_tobeDetected(streamIdx_4) = [];
      end

      if(activeOrdering)
	  s_est_soft_SIC(detected_stream_index,:) = s_est_soft_SIC;
      end
    end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Subfunktion of SIC algorithm
function streamIdx = get_first_SIC_stream_todetect(w_rx,channelH,noise_var)
    G = w_rx*channelH;
    P_Sig = abs(diag(G)).^2;
    P_ICI = (abs(diag((G-diag(diag(G)))*((G-diag(diag(G)))')))).^2;
    P_noise = noise_var*diag((w_rx*(w_rx')));
    SINR = P_Sig./(P_ICI+P_noise);
    [~,streamIdxSorted] = sort(SINR, 'descend');
    streamIdx = streamIdxSorted(1);
end

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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Subfunktion of SIC algorithm
function [s_est_hard] = HardDecisionUnit(s_est_soft,QAMorder)
    if(QAMorder == 1)
      s_est_hard = sign(real(s_est_soft));
    end
    
    if(QAMorder == 2)
      s_est_hard = 1/sqrt(2)*(sign(real(s_est_soft))+1i*sign(imag(s_est_soft)));
    end
end

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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% open source code from octave
      function b = randsrc (n, m, alphabet, seed)

        switch (nargin)
          case 0,
        m = 1;
        n = 1;
        alphabet = [-1,1];
        seed = Inf;
          case 1,
        m = n;
        alphabet = [-1,1];
        seed = Inf;
          case 2,
        alphabet = [-1,1];
        seed = Inf;
          case 3,
        seed = Inf;      
          case 4,
          otherwise
        usage ("b = randsrc (n, [m, [alphabet, [seed]]])");
        endswitch

        ## Check alphabet
        [ar,ac] = size (alphabet);
        if (ac == 1)
          b = alphabet (1, 1) * ones (n, m);
          return;
        endif

        if (ar == 1)
          prob = [1:ac] / ac;
        elseif (ar == 2)
          prob = alphabet(2,:);
          alphabet = alphabet(1,:);
          if (abs(1-sum(prob)) > eps)
        error ("randsrc: probabilities must added up to one");
          endif
          prob = cumsum(prob);
        else
          error ("randsrc: alphabet must have 1 or 2 rows");
        endif
        
        ## Check seed;
        if (!isinf (seed))
          old_seed = rand ("seed");
          rand ("seed", seed);
        endif
        
        ## Create indexes with the right probabilities
        tmp = rand (n, m);
        b = ones (n, m);
        for i = 1:ac-1 
          b = b + (tmp > prob(i));
        end

        ## Map the indexes to the symbols
        b = alphabet(b);
        
        ## BUG: the above gives a row vector for some reason. Force what we want
        b = reshape(b, n, m);
          
        ## Get back to the old
        if (!isinf (seed))
          rand ("seed", old_seed);
        endif

      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; 1 221/255 0]; % gelb
        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

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