function mainBER(numOfMemElem1,g1_1,g2_1,numOfMemElem2,g1_2,g2_2,code1,code2,Uncoded,Punctured,UnPunctured,Improved_Results,filename) 
fh = figure(1); 

if code1 == 0
numOfMemElem1 = 1;
g1_1 = 3;
g2_1 = 2;

elseif code1 == 1
numOfMemElem1 = 2;

g1_1 = 7;
g2_1 = 5;

elseif code1 == 2
numOfMemElem1 = 3;

g1_1 = 13;
g2_1 = 15;

elseif code1 == 3
numOfMemElem1 = 4;

g1_1 = 23;
g2_1 = 37;

elseif code1 == 4
numOfMemElem1 = 5;

g1_1 = 67;
g2_1 = 45;

elseif code1 == 5
numOfMemElem1 = 6;

g1_1 = 147;
g2_1 = 117;

elseif code1 == 6
numOfMemElem1 = 8;
g1_1 = 435;
g2_1 = 777;

elseif code1 == 7
numOfMemElem1 = 4;
g1_1 = 37;
g2_1 = 21;
end

if code2 == 0
numOfMemElem2 = 1;

g1_2 = 3;
g2_2 = 2;

elseif code2 == 1
numOfMemElem2 = 2;

g1_2 = 7;
g2_2 = 5;

elseif code2 == 2
numOfMemElem2 = 3;

g1_2 = 13;
g2_2 = 15;

elseif code2 == 3
numOfMemElem2 = 4;

g1_2 = 23;
g2_2 = 37;

elseif code2 == 4
numOfMemElem2 = 5;

g1_2 = 67;
g2_2 = 45;

elseif code2 == 5
numOfMemElem2 = 6;

g1_2 = 147;
g2_2 = 117;

elseif code2 == 6
numOfMemElem2 = 8;
g1_2 = 435;
g2_2 = 777;

elseif code2 == 7
numOfMemElem2 = 4;
g1_2 = 37;
g2_2 = 21;
end

global NumberOfbits    
global Ia 
global BitsUncoded 
global histLimit 
NumImproved               = 30000       ; 
if Improved_Results 
        Numbits           = NumImproved  ;                                  
else 
        Numbits           = 10000        ;      
end 

if ~(Uncoded || Punctured || UnPunctured) 
Punctured = 1; 
end 
if numOfMemElem1 == numOfMemElem2 && g1_1 == g1_2 && g2_1 == g2_2 
    Identical_decoders    = 1            ;    
else 
    Identical_decoders    = 0            ; 
end 
if  Identical_decoders   == 0 
    Numbits               = 10000        ; 
end 
NumberOfbits              = Numbits      ; 
Ia                        = 0:0.1:1      ; 
histLimit                 = 32           ;                       % Maximum truncating value for histogram proper usage 
BitsUncoded               = randi(2,1,NumberOfbits)-1 ;           % generate a random sequence of '1's and '0's 
 

if Identical_decoders == 1 
    [trellis, template, ~, notPuncInd, ~, BitsBPSK]                             = Construct(numOfMemElem1, g1_1, g2_1,0,1); 
    [trellis_unPunc, template_unPunc, ~, notPuncInd_unPunc, ~, BitsBPSK_unPunc] = Construct(numOfMemElem1, g1_1, g2_1,0,0); 
else 
    [trellis1, template1, ~, notPuncInd1, ~, BitsBPSK1] = Construct(numOfMemElem1, g1_1, g2_1,0,1); 
    [trellis1_unPunc, template1_unPunc, ~, notPuncInd1_unPunc, ~, BitsBPSK1_unPunc] = Construct(numOfMemElem1, g1_1, g2_1,0,0); 
    [trellis2, template2, ~, notPuncInd2, ~, BitsBPSK2] = Construct(numOfMemElem2, g1_2, g2_2,0,1); 
    [trellis2_unPunc, template2_unPunc, ~, notPuncInd2_unPunc, ~, BitsBPSK2_unPunc] = Construct(numOfMemElem2, g1_2, g2_2,0,0); 
 
end 

legBER={}; 

 
    if Identical_decoders == 1 
        xLeft  = zeros(1,3); 
        xRight = zeros(1,3); 
        yDown  = zeros(1,3); 
        yTop = zeros(1,3); 
 
        %figure; 
        if Uncoded 
        disp('Working on BER curve (Uncoded)') 
        [xlimit, ylimit] = drawBERcurve(template, trellis, 0, BitsBPSK, notPuncInd,3,Improved_Results,1,filename); %uncoded
        legBER = {'uncoded'}; 
        xLeft (1) = xlimit(1); 
        xRight(1) = xlimit(2); 
        yDown (1) = ylimit(1); 
        yTop(1) = ylimit(2); 
        hold on 

        end 
        if Punctured 
        disp('Working on BER curve (Coded 2/3 punctured)') 
        [xlimit, ylimit] = drawBERcurve(template, trellis, 0, BitsBPSK, notPuncInd,1,Improved_Results,1,filename); 
        legBER = [legBER, 'Coded 2/3 punctured']; 
        xLeft (2) = xlimit(1); 
        xRight(2) = xlimit(2); 
        yDown (2) = ylimit(1); 
        yTop(2) = ylimit(2); 
        hold on

        end 
        if UnPunctured 
        disp('Working on BER curve (Coded without puncturation)') 
        [xlimit, ylimit] = drawBERcurve(template_unPunc, trellis_unPunc, 0, BitsBPSK_unPunc, notPuncInd_unPunc,2,Improved_Results,1,filename); 
        legBER = [legBER, 'Coded 1/2 without puncturation']; 
        xLeft (3) = xlimit(1); 
        xRight(3) = xlimit(2); 
        yDown (3) = ylimit(1); 
        yTop(3) = ylimit(2); 
        hold off 
        end 
 
    else 
        %figure; 
        xLeft  = zeros(1,5); 
        xRight = zeros(1,5); 
        yDown  = zeros(1,5); 
        yTop = zeros(1,5); 
        
        if Uncoded 
        EncInd=3; %Uncoded 
        [xlimit, ylimit] = drawBERcurve(template1, trellis1, 0, BitsBPSK1, notPuncInd1,EncInd,Improved_Results,1,filename); %uncoded 
        legBER = [legBER, 'uncoded']; 
        xLeft (1) = xlimit(1); 
        xRight(1) = xlimit(2); 
        yDown (1) = ylimit(1); 
        yTop(1) = ylimit(2); 
        hold on 
        end 
        if Punctured 
        EncInd=1; %first encoder alone 
        [xlimit, ylimit] = drawBERcurve(template1, trellis1, 0, BitsBPSK1, notPuncInd1,EncInd,Improved_Results,1,filename); 
          legBER = [legBER, 'Encoder 1 punctured'];        
        xLeft (2) = xlimit(1); 
        xRight(2) = xlimit(2); 
        yDown (2) = ylimit(1); 
        yTop(2) = ylimit(2); 
        hold on 
        end 
        if UnPunctured 
        EncInd=11; %first encoder alone WITHOUT puncturing 
        [xlimit, ylimit] = drawBERcurve(template1_unPunc, trellis1_unPunc, 0, BitsBPSK1_unPunc, notPuncInd1_unPunc, EncInd,Improved_Results,1,filename); 
          legBER = [legBER, 'Encoder 1 without puncturation'];         
        xLeft(3) = xlimit(1); 
        xRight(3) = xlimit(2); 
        yDown (3) = ylimit(1); 
        yTop(3) = ylimit(2); 
        hold on 
        end 
        if Punctured 
        EncInd=2; %second encoder alone 
        [xlimit, ylimit] = drawBERcurve(template2, trellis2, 0, BitsBPSK2, notPuncInd2,EncInd,Improved_Results,1,filename); 
        legBER = [legBER, 'Encoder 2 punctured'];           
        xLeft (4) = xlimit(1); 
        xRight(4) = xlimit(2); 
        yDown (4) = ylimit(1); 
        yTop  (4) = ylimit(2); 
        hold on 
        end 
        if UnPunctured 
        EncInd=12; %second encoder alone WITHOUT puncturing 
        drawBERcurve(template2_unPunc, trellis2_unPunc, 0, BitsBPSK2_unPunc, notPuncInd2_unPunc, EncInd,Improved_Results,1,filename); 
          legBER = [legBER, 'Encoder 2 without puncturation'];        
        xLeft (5) = xlimit(1); 
        xRight(5) = xlimit(2); 
        yDown (5) = ylimit(1); 
        yTop  (5) = ylimit(2); 
        hold off 
        end 
 
        
        
    end 

    %xlimit=([min(xLeft) max(xRight) ]); 
    %ylimit=([min(yDown) max(yTop) ]); 
    xlabel('E_b/N_o in dB') 
    ylabel('BER') 
    title ('BER curve') 
    legend (legBER) 
    %print(filename, '-dgif','-S640,480', '-r0');
 print(fh, '-dpng','-r72', filename);  
   %xlim([min(xLeft) max(xRight)]); 
    %ylim([min(yDown) max(yTop) ]); 
     
 end 
 
 function [ output] = LOGMAP( template, trellis, BitsBPSK, Ia, noise, variance, L_a_input, ref, notPuncInd)                      
        global NumberOfbits    
        global BitsUncoded 
        global histLimit         
        if ref == 3 %Uncoded 
            received                                 = BitsBPSK + noise                                      ; 
            BitsReceived                             = zeros(1, 2*NumberOfbits)                              ;             
            BitsReceived( notPuncInd )               = received                                              ; 
            BitsReceivedUncoded                      = BitsReceived(1:2:end)                                 ; 
            BitsCorrected( BitsReceivedUncoded >= 0) = 1                                                     ; 
            BitsCorrected( BitsReceivedUncoded <=0 ) = 0                                                     ; 
            BER = sum(xor(BitsCorrected, BitsUncoded(1:length(BitsCorrected))))/length(BitsCorrected);
            %[Numb,BER]                               = biterr(BitsCorrected, BitsUncoded(1:length(BitsCorrected)))                 ;         
        else %Coded 
                received                             = BitsBPSK + noise                                      ; 
                BitsReceived                         = zeros(1, 2*NumberOfbits)                              ; 
                BitsReceived( notPuncInd )           = received                                              ; 
        end 
       				                    
        L_c                    = 2 / variance                                                                ;       
        L_a                    = L_a_input                                                                   ; 
        gammas                 = constructGammas(template, trellis, L_a,BitsReceived, variance)              ; 
        [Alphas, Betas]        = constructAlphasBetas(trellis, gammas)                                       ; 
        Betas_LLR              = Betas(reshape(trellis.nextStates,2*trellis.numStates,1)+1,:)                ; 
        Alphas_LLR             = Alphas((1:size(Alphas,1))' * ones(1,2), :)                                  ; 
        Gammas_LLR             = [ gammas(1:2:length(gammas(:,1)),:) ; gammas(2:2:length(gammas(:,1)),:) ]   ; 
        LLR_pre                = Betas_LLR + Alphas_LLR + Gammas_LLR                                         ; 
        x1                     = LLR_pre(trellis.numStates+1:2*trellis.numStates,:)                          ; 
        x0                     = LLR_pre(1:trellis.numStates,:)                                              ; 
        LLR_log1               = max(x1) + log(sum(exp(x1 - repmat(max(x1),trellis.numStates,1)),1))         ; 
        LLR_log0               = max(x0) + log(sum(exp(x0 - repmat(max(x0),trellis.numStates,1)),1))         ; 
        LLR_log                = LLR_log1 - LLR_log0                                                         ; 
        L_e                    = LLR_log - L_a  -L_c *(BitsReceived(1:2:end))                                ;   
        Corrected(LLR_log>=0)  = 1                                                                           ; 
        Corrected(LLR_log<0 )  = 0                                                                           ; 
        if ref ~= 3 
				BER = sum(xor(Corrected, BitsUncoded(1:length(Corrected))))/length(Corrected);           
           % [Numb, BER]            = biterr( Corrected, BitsUncoded(1:length(Corrected) ) )                  ; 
        end 
        L_e(L_e >=  histLimit) =  histLimit                                                                  ; 
        L_e(L_e <= -histLimit) = -histLimit                                                                  ; 
        Ie                     = measure_mutual_information_histogram(L_e, Ia)                               ; 
        output                 = BER                                                                         ; 
end 
 
 function [trellis, template, Encoded, notPuncInd, toBeSent, BitsBPSK] = Construct(numOfMemElem, g1, g2, int,punc) 
    global BitsUncoded; 
    trellis        = poly2trellis( numOfMemElem + 1, [g1 g2])             ;     %Encoding 
    trellis        = rcs(trellis)                
    template       = constructTemplate( trellis )                          ;     %Constructing template 
    Encoded        = convenc( BitsUncoded, trellis )                       ; 
    if punc == 1 
        [tBS,nP]   = puncture( Encoded )                                   ;     %Puncturing 
    else 
        tBS        = Encoded                                               ; 
        nP         = 1:length( Encoded )                                   ; 
    end 
    BitsBPSK       = 2*tBS -1                                              ;      %BPSK Modulation 
    notPuncInd     = nP                                                    ; 
    toBeSent       = tBS                                                   ; 
 end 
    
 function [trell]    = rcs(t)
    x                         = de2bi( t.outputs,log2(max(max(t.outputs))+1));
    y                         = reshape( x(:, 2), t.numStates, 2 )           ;
    for i = 1:length( t.nextStates(:, 1) ) %nextStates
        if y(i,1) == 1
            temp              = t.nextStates(i,1)                            ;
            t.nextStates(i,1) = t.nextStates(i,2)                            ;
            t.nextStates(i,2) = temp                                         ;
            temp              = t.outputs(i,1)                               ;
            t.outputs(i,1)    = t.outputs(i,2)                               ;
            t.outputs(i,2)    = temp                                         ;
        end
    end
            trell             = t                                            ;
    end
 
 function [Template] = constructTemplate(trell) 
        global NumberOfbits ;   
        % Constructing Output matrix 
        output_Matrix         = nan(trell.numStates, 2*trell.numStates)                                       ; 
        Input_Matrix          = nan(trell.numStates, trell.numStates)                                         ; 
        for i=1:trell.numStates 
            for j=1:2 
                output_Matrix( i, 2*(trell.nextStates(i,j))+1:2*(trell.nextStates(i,j))+2 )... 
                              = 2*de2bi(trell.outputs(i,j)', 2, 'left-msb') - 1                               ; 
                Input_Matrix ( i, trell.nextStates(i,j)+1) ... 
                              = find( trell.nextStates(i,:) == trell.nextStates(i,j) ) - 1                    ; 
            end 
        end 
        trellisBranchOutputs  = 2*de2bi( trell.outputs', log2( max(max(trell.outputs)) + 1 ), 'left-msb') - 1 ; 
        %Constructin previous states matrix 
        prevStates            = [ ]                                                                           ; 
        for i=1:trell.numStates 
            prevStates        = [ prevStates; ( find( sum( (trell.nextStates == i - 1)' )' == 1 ) - 1 )' ]    ; 
        end 
        previous              = prevStates                                                                    ; 
        for i = 2:trell.numStates-1 
            temp              = previous(i, 1)                                                                ; 
            previous(i,1)     = previous(i, 2)                                                                ; 
            previous(i,2)     = temp                                                                          ; 
        end 
        [a_prev, b_prev]      = sort( reshape( previous', 2*trell.numStates, 1 ) + 1, 'ascend' )              ; 
        for i = 3:2:(2*trell.numStates) - 2 
            temp              = b_prev(i  , 1)                                                                ; 
            b_prev(i, 1)      = b_prev(i+1, 1)                                                                ; 
            b_prev(i+1, 1)    = temp                                                                          ; 
        end 
        % Entrance stage for the trellis 
        activeStates          = 0                                                                             ; 
        count                 = 0                                                                             ; 
        output_Conc           = [ ]                                                                           ; 
        while length( activeStates ) < trell.numStates 
            output            = [ ]                                                                           ; 
            tN                = trell.nextStates( activeStates+1, : )                                         ; 
            for i = 1:trell.numStates 
                if ismember( i, activeStates + 1 ) 
                    output    = [ output; 2*de2bi( trell.outputs(i, :)', log2( max( max(trell.outputs))+1)... 
                                 , 'left-msb') - 1 ]                                                          ; 
                else 
                    output    = [ output; nan nan; nan nan ]                                                  ; 
                end 
            end 
            output_Conc       = [ output_Conc output]                                                         ; 
            [ newStates, not ]    = unique( reshape( tN, 1, numel(tN) ) )                                         ; 
            for i = 1:length(newStates) 
              if ~ismember( newStates(i), activeStates ) 
                 activeStates = [ activeStates newStates(i) ]                                                 ; 
              end 
            end 
            count             = count + 1                                                                     ; 
        end 
        %Exit stage for the trellis 
        activeStates          = 0                                                                             ; 
        countPrev             = 0                                                                             ; 
        output_Conc_Prev      = [ ]                                                                           ; 
        while length(activeStates) < trell.numStates 
            output            = [ ]                                                                           ; 
            output_temp       = nan( 2*trell.numStates, 2 )                                                   ; 
            for i = 1:trell.numStates 
                if ismember( i, activeStates + 1 ) 
                    xx        = previous(i, :)                                                                ; 
                    for k = 1:2 
                        offset                           = Input_Matrix (xx(k)+1, i        )                  ; 
                        output_temp(2*xx(k)+1+offset,:)  = output_Matrix(xx(k)+1, 2*i-1:2*i)                  ; 
                    end 
                end 
            end 
            output            = [ output_temp output ]                                                        ; 
            output_Conc_Prev  =[output output_Conc_Prev ]                                                     ; 
            tP                = prevStates(activeStates + 1 , :)                                              ; 
            [newPrevStates, not ] = unique( reshape(tP, 1, numel(tP) ) )                                          ; 
            for i = 1:length(newPrevStates) 
                if ~ismember( newPrevStates(i), activeStates) 
                    activeStates = [ activeStates newPrevStates(i) ]                                          ; 
                end 
            end 
            countPrev            = countPrev + 1                                                              ; 
        end 
        steady                   = repmat( trellisBranchOutputs, 1, NumberOfbits - count - countPrev )        ;    
        Template                 = [ output_Conc steady output_Conc_Prev ]                                    ; 
 end 
    
 function [gammas]   = constructGammas(template, trellis, L_a,BitsReceived,variance) 
     inf                   = 1000                                                                                 ; 
     L_a_Matrix            = repmat( L_a, 2*trellis.numStates, 1 )                                                ;         
     BitsReceived_Matrix   = repmat( BitsReceived, 2*trellis.numStates, 1 )                                       ; 
     gammas                = 0.5*(((2/variance)*BitsReceived_Matrix(:,1:2:end)+L_a_Matrix).*template(:,1:2:end)... 
                           +       (2/variance)*BitsReceived_Matrix(:,2:2:end)            .*template(:,2:2:end))  ;               
     gammas(isnan(gammas)) = -inf                                                                                 ; 
    end 
 
 function [Alphas_log,Betas_log] = constructAlphasBetas(trellis, gammas) 
    global NumberOfbits ;   
    numStates         = trellis.numStates; 
    inf               = 1000; 
    % gammas_ratio = exp(gammas); 
    [a,b]             = sort( reshape( trellis.nextStates', 2*trellis.numStates, 1 ) + 1, 'ascend') ; 
    bBeta             = reshape( trellis.nextStates', 2*trellis.numStates, 1 ) +1                   ; 
    gammas1           = gammas(b, :)                                                                ; 
    %Alpha elements 
    alphas            = zeros( 2*trellis.numStates,1)                                               ; 
    alphas_final      = zeros( trellis.numStates,NumberOfbits)                                      ; 
    alphas_final(:,1) = [ 0; -1000*ones( trellis.numStates-1, 1 ) ]                                 ; 
    %Beta elements 
    iBeta             = NumberOfbits - 1                                                            ; 
    Betas             = zeros( 2*trellis.numStates, 1 )                                             ; 
    Betas_final       = zeros( trellis.numStates, NumberOfbits)                                     ; 
    nn                = NumberOfbits                                                                ; 
    Betas_final(:,nn) = [0; -1000*ones( trellis.numStates-1, 1 ) ]                                  ; 
    
    for i=2:NumberOfbits 
        xx            = alphas_final(:, i-1)'                                                       ; 
        alphas_temp   = reshape( xx([1:size(xx,1)]' * ones(1,2), :),2*numel(alphas_final(:,i-1)),1) ; 
        alphas_temp   = alphas_temp (b,:)                                                           ; 
        alphas        = (alphas_temp + gammas1(:,i-1))                                              ; 
        x             = Betas_final(:,iBeta+1)                                                      ; 
        Betas_temp    = x(bBeta,:)                                                                  ; 
        Betas         = (Betas_temp + gammas(:,iBeta+1))                                            ; 
        ii            = 1                                                                           ; 
    for j = 1:trellis.numStates 
             alphas_final(j,i)    = max( alphas(ii), alphas(ii+1) ) ... 
                                  + log( 1 + exp( -abs( alphas(ii) - alphas(ii+1) ) ) )             ; 
             Betas_final(j,iBeta) = max( Betas (ii), Betas (ii+1) ) ... 
                                  + log( 1 + exp( -abs( Betas (ii) - Betas (ii+1) ) ) )             ; 
             ii                   = ii + 2                                                          ; 
        end 
        
        if ( i <= (log2(numStates)) || i >= NumberOfbits - (log2( numStates ) ) + 1 ) 
            alphas_final(isnan(alphas_final))   = -inf                                              ; 
            Betas_final(isnan(Betas_final))     = -inf                                              ; 
        end 
        iBeta                                   = iBeta - 1                                         ; 
    end     
    Betas_log   = Betas_final                                                                       ; 
    Alphas_log  = alphas_final                                                                      ; 
    end 
 
 function [xlim, ylim] = drawBERcurve(template, trellis, Ia_ip, BitsBPSK, notPuncInd, EncInd, Improved_Results, repeat,filename)    
     
        global NumberOfbits; 
        %cstr = 'rkbgm'; 
        cstr = [1 0 0; 0 0 0; 0 0 1; 0 1 0; 0 1 1];
        ref=2; 
        SNR_BER = 0:1:9; 
        if Improved_Results 
            SNR_BER = 0:0.5:9; 
        end 
        if EncInd == 3 %uncoded 
            ref=3; 
        end 

        L_a_input         = generate_La( Ia_ip )                                  ; 
                      
        BER               = zeros(1, length(SNR_BER) )                            ; 
        for iSNR_BER = 1:length(SNR_BER) 
            disp ( [ 'SNR = ', num2str( SNR_BER(iSNR_BER) ), ' dB' ] )            ; 
            SNR_BER_ratio = 10 ^ (SNR_BER (iSNR_BER) / 10 )                       ; 
            variance_BER  = 1 / SNR_BER_ratio                                     ; 
            noise_BER     = normrnd( 0, sqrt(variance_BER), 1, length(BitsBPSK) ) ;            
            BER(iSNR_BER) = LOGMAP(template, trellis, BitsBPSK, Ia_ip, noise_BER, variance_BER, L_a_input, ref, notPuncInd) ;              
        end   

            BER           = smooth( BER )                                         ; 
            BER ( BER == 0 )  = 1/NumberOfbits ;                                             ; 
        %BER ( BER < 1e-4 )= 1e-4                                                                    ; 
            if EncInd == 11 
                semilogy ( SNR_BER, BER, '--' , 'markersize', 2 ,'linewidth',2, 'color', cstr(1,:))                       ; 
            elseif EncInd == 12 
                semilogy ( SNR_BER, BER, '--' , 'markersize', 2 ,'linewidth',2, 'color', cstr(2,:))                       ; 
            else 
                semilogy ( SNR_BER, BER, 'o-' , 'markersize', 2 ,'linewidth',2, 'color', cstr(EncInd,:))                       ; 
            end 
            xlim          = [ SNR_BER(1) SNR_BER(end) + 1 ]                       ; 
            ylim          = [ min(BER) max(BER) ]                                 ; 
            grid on; 
               
 end 
    
 function [LLR] = generate_La(Ia) 
        global BitsUncoded; 
        global histLimit; 
        sigma_a = JInverseFunction( Ia )                         ; 
        mu_a    = ( sigma_a^2 ) / 2                              ; 
        u       = 2*BitsUncoded - 1                              ; 
        if     Ia==1  
        LLR     = u*histLimit                                    ; 
        elseif Ia==0                                             ; 
        LLR     =  zeros(1, length(u) )                          ; 
        else 
        LLR     =  normrnd( u*(mu_a), sigma_a, 1, length(u) )    ; 
        end 
    
    end 
 
  function [mutual_information] = measure_mutual_information_histogram(llrs,Ia_sent) 
        global BitsUncoded; 
        global histLimit; 
        indices1           = BitsUncoded == 1                                                                       ;            
        indices0           = BitsUncoded == 0                                                                       ; 
        Num_HistBins       = 400                                                                                    ;                       
        xAxis              = linspace( -histLimit,histLimit, Num_HistBins )                                         ;   
        hist1              = hist( llrs(indices1), xAxis )                                                          ; 
        hist0              = hist( llrs(indices0), xAxis )                                                          ; 
        count1             = sum( BitsUncoded )                                                                     ; 
        count0             = length(BitsUncoded) - count1                                                           ; 
        if count1~=0;        hist1=hist1/(count1); end                                                              ; 
        if count0~=0;        hist0=hist0/(count0); end                                                              ; 
        t0                 = sum( hist0(hist0~=0).*log2( (2*hist0(hist0~=0))./(hist1(hist0~=0)+hist0(hist0~=0))) )  ;   
        t1                 = sum( hist1(hist1~=0).*log2( (2*hist1(hist1~=0))./(hist1(hist1~=0)+hist0(hist1~=0))) )  ; 
        mutual_information = 0.5 * ( t0 + t1 )                                                                      ; 
        if Ia_sent == 1; mutual_information = 1; end                                                                ; 
  end 
    
  function [ JI ] = JInverseFunction( I ) 
        I_tp     =  0.3646                                                               ; 
        a_alpha1 =  1.09542                                                              ; 
        b_alpha1 =  0.214217                                                             ; 
        c_alpha1 =  2.33727                                                              ; 
        a_alpha2 =  0.706692                                                             ; 
        b_alpha2 =  0.386013                                                             ; 
        c_alpha2 = -1.75017                                                              ; 
        JI_temp  = zeros(1,length(I))                                                    ; 
        JI_temp(I >= 0   & I <= I_tp)  = (a_alpha1*(I(I >= 0 & I <= I_tp).^2)) ... 
                                       + (b_alpha1*(I(I >= 0 & I <= I_tp)))    ... 
                                       + (c_alpha1*(I(I >= 0 & I <= I_tp).^(1/2)))       ; 
        JI_temp(I > I_tp & I <= 1   )  = (-a_alpha2*log(-b_alpha2*((I(I > I_tp ... 
                                       & I <= 1))-1)))-(c_alpha2*(I(I > I_tp & ... 
                                       I <= 1 ) ) )                                      ; 
        JI      = JI_temp                                                                ; 
  end 
 
   function [punc, notPuncInd] = puncture(enc) 
    puncInd            = 2:4:length(enc)                   ; 
    notPuncInd         = setdiff( 1:length(enc), puncInd ) ; 
    punc_temp          = enc                               ; 
    punc_temp(puncInd) = [ ]                               ; 
    punc               = punc_temp                         ; 
    end 

