%%=====================================================
% Note : This file shows the different trials to 
%        validate the idea of trellis agglomeration.
%%=====================================================
% Date: 15/07/2012
% Author: Emmanuel Boutillon
%  Fast trellis computation method of a CC when there is no redundancy
%  In this program, we show in 4 examples that if there are only information
%  bits (no redundancy), then it is always possible to compact the N trellis
%  steps in N_state - 1 trellis stages (being 7 sections for the LTE).
%  To do this, it is necessary to pre-compute the N_state - 1 information bit pb
%  equivalent to the N information bits by doing the "soft logic or" of all the 
%  bits of the N sequence having the same residue modulo N_state -1.
% ===========================================
% This matlab program consists of 4 study cases
% 1) 4-state trellis in the probability domain
% 2) 4-state trellis in the log domain with max-log-map
% 3) 8-state trellis in the probability domain
% 4) 8-state trellis in the log domain with max-log-map
%==========================================
% Note: Since agglomeration bits is simple, it is not too expensive to implement
% it. It can be implemented with the log-map rather than the max-log-map in 
% order to improve performance.
%%==========================================
% Change history
%%==========================================
% Date: 15/07/2012
% Author: Emmanuel Boutillon
% Addition of the fifth study case, with the case of 8-state trellis in the log 
% domain with min-log, only 3 aggregated bits of minimum reliability are involved
% in the output.
% In other words, it is possible to saturate the 4 remaining agglomerated bits of 
% greater reliability: It is possible to replace, for each of the 4 corresponding
% trellis stages, the 8 ACS operations by just 8 selections (low consumption).
% This also allows not to keep more than 3 minimums (thus to only do 3 comparisons)
% to pre-compute on the fly these 3 trellis compression bits between two half iterations
%
% Thus, these are the 5 study cases of the program:
% 1) 4-state trellis in the probability domain
% 2) 4-state trellis in the log domain whit max-log-map
% 3) 8-state trellis in the probability domain
% 4) 8-state trellis in the log domain whit max-log-map
% 5) 8-state trellis in the log domain whit max-log-map only taking into account 
% the log2(8)=3 least reliable bits.
%%==========================================
% Date: 28/08/2013
% Author: Emmanuel Boutillon
% Addition of the sixth study case: code to generate the tables in paper 
% "Simplified Compression of Redundancy Free Trellis Sections in Turbo-Decoder"
% Accepted/Published in IEEE Communications Letters 2014.
%%=======================================================
% Study case 1: CC 4 recursive states Pol = 3. Study case in the probability domain
%=========================================================
N = 21;  %Should be a multiple of 3.

% Random selection of 18 probabilities (between 0 and 1). 
p = (cos(pi*(1+cos( pi*rand(1,N)))/2)+1)/2; 


% Initial state of the forward node metric (we can put in any number)
a(:,1) = [1 0 0 0 ];

% Recursive computation of the N trellis stages
for i=1:N
    G= [p(i)  (1-p(i)) 0 0; 0 0 (1 -p(i)) p(i) ;(1-p(i)) p(i) 0 0 ; 0 0 p(i) (1-p(i))];
    a(:,i+1) = G*a(:,i);
end;

% Agglomeration of the bits having the same modulo 3 remainder (soft XOR)
for i=1:3
    pb(i) = 1;
    for j=0:(N/3-1)
        pb(i) = pb(i)*p(i+3*j) + (1 - pb(i))*(1 - p(i+3*j));
    end
end;


% Implementation of the 3 agglomerated trellis stages
b(:,1) = a(:,1);
for i=1:3
    G= [pb(i)  (1-pb(i)) 0 0; 0 0 (1 -pb(i)) pb(i) ;(1-pb(i)) pb(i) 0 0 ; 0 0 pb(i) (1-pb(i))];
    b(:,i+1) = G*b(:,i);
end;

% We check that the final result is the same
b
a
% Conclusion: it is the same

%% ============================================
% Study case 2: CC 4 recursive states Pol = 3. Study case with Max-Log-MAP
%=========================================================
LLR = log( p./(1-p));

% bit's probabilities
a(:,1) = [9 1 1  1];
a(:,1) = a(:,1)/sum(a(:,1));

A(:,1) = log(a(:,1)/a(1,1));

G0 = [1 0 0 0 ; 0 0 0 1; 0 1 0 0; 0 0 1 0];
G1 = [0 1 0 0;  0 0 1 0; 1 0 0 0; 0 0 0 1];

% Generation of the max-log-map trellis
for i=1:N
    A(:,i+1) = max ( (G0*A(:,i) + LLR(i)*ones(4,1)) , G1*A(:,i) - LLR(i)*ones(4,1)); 
end;

A(:,N+1) = A(:,N+1) - A(1,N+1);

% Agglomeration of the bits having the same modulo 3 remainder (min-sum)
for i=1:3
    LLRb(i) = 10000;
    for j=0:(N/3-1)
        LLRb(i) = sign(LLRb(i))*sign(LLR(i+3*j))*min(abs(LLR(i+3*j)), abs(LLRb(i)));
    end
end;

% Implementation of the 3 agglomerated trellis stages
B(:,1) = A(:,1);
for i=1:3
    B(:,i+1) = max ( (G0*B(:,i) + LLRb(i)*ones(4,1)) , G1*B(:,i) - LLRb(i)*ones(4,1)); 
end;

B(:,4) = B(:,4) - B(1,4);

% We check that the final result is the same
A(:,end)
B(:,4)
% Conclusion: It is the same

%% ============================================
% Etude 3: CC 8 recursive states Pol = 15. Study case un the probability domain
%=========================================================
clear all;
N = 14;    % Multiple of 7, but it is not required.


% probability's generation (between 0 and 1). 
p = (cos(pi*(1+cos( pi*rand(1,N)))/2)+1)/2; 

% Probabilities of the initial node metrics 
a(:,1) = [1 0 0 0 0 0 0 0 ];

% N trellis stages initialisation 
for i=1:N
P = p(i); Q = 1 - p(i);
G= [P Q 0 0 0 0 0 0 
      0 0 Q P 0 0 0 0 
      0 0 0 0 P Q 0 0 
      0 0 0 0 0 0 Q P 
      Q P 0 0 0 0 0 0 
      0 0 P Q 0 0 0 0 
      0 0 0 0 Q P 0 0 
      0 0 0 0 0 0 P Q];
    a(:,i+1) = G*a(:,i);
end;


% Agglomeration of the bits having the same modulo 7 remainder (soft XOR). 
for i=1:7
    pb(i) = 1;
    for j=0:(N/7-1)
        pb(i) = pb(i)*p(i+7*j) + (1 - pb(i))*(1 - p(i+7*j));
    end
end;


% Implementation of the 7 agglomerated trellis stages
b(:,1) = a(:,1);
for i=1:7
P = pb(i); Q = 1 - pb(i);
G= [P Q 0 0 0 0 0 0 
      0 0 Q P 0 0 0 0 
      0 0 0 0 P Q 0 0 
      0 0 0 0 0 0 Q P 
      Q P 0 0 0 0 0 0 
      0 0 P Q 0 0 0 0 
      0 0 0 0 Q P 0 0 
      0 0 0 0 0 0 P Q];
    b(:,i+1) = G*b(:,i);
end;

% We check that the final result is the same
b(:,end)
a(:,end)
% Conclusion: It is the same.

%% ===================================
% Study case 4: CC 8 recursive states Pol = 7. Study case in the LLR domain
%              with Max-Log-Map
%==============================================
N = 14;    % Multiple of 7 but it is not necessary.

% Probability generation (between de 0 and 1). 
p = (cos(pi*(1+cos( pi*rand(1,N)))/2)+1)/2; 

LLR = log( p./(1-p));

% bit's probabilities
a(:,1) = [9 1 1  1 1 1 1 1 ]';
a(:,1) = a(:,1)/sum(a(:,1));

A(:,1) = log(a(:,1)/a(1,1));

G0 = [1 0 0 0 0 0 0 0 ; 0 0 0 1 0 0 0 0 ; 0 0 0 0 1 0 0 0 ; 0 0 0 0 0 0 0 1 ;
        0 1 0 0 0 0 0 0 ; 0 0 1 0 0 0 0 0 ; 0 0 0 0 0 1 0 0 ; 0 0 0 0 0 0 1 0];
G1 = [0 1 0 0 0 0 0 0 ; 0 0 1 0  0 0 0 0 ; 0 0 0 0 0 1 0 0; 0 0 0 0 0 0 1 0; 
        1 0 0 0 0 0 0 0 ;  0 0 0 1 0 0 0 0 ; 0 0 0 0 1 0 0 0; 0 0 0 0 0 0 0 1];

% Max-log-map trellis generation 
for i=1:N
    A(:,i+1) = max ( (G0*A(:,i) + LLR(i)*ones(8,1)) , G1*A(:,i) - LLR(i)*ones(8,1)); 
end;

A(:,N+1) = A(:,N+1) - A(1,N+1);

% Agglomeration of the bits having the same modulo 7 remainder (min-sum)
for i=1:7
    LLRb(i) = 10000;
    for j=0:(N/7-1)
        LLRb(i) = sign(LLRb(i))*sign(LLR(i+7*j))*min(abs(LLR(i+7*j)), abs(LLRb(i)));
    end
end;

% Implementation of the 7 agglomerated trellis stages
B(:,1) = A(:,1);
for i=1:7
    B(:,i+1) = max ( (G0*B(:,i) + LLRb(i)*ones(8,1)) , G1*B(:,i) - LLRb(i)*ones(8,1)); 
end;

% We check that the final result is the same
A(:,end)-A(1,end)
B(:,end)-B(1,end)
% Conclusion: It is the same.

%%
% ============================================
% Study case 5: CC 8 recursive states = 15. Study case in the LLR domain with 
%               max-log-map.
% Implementation of the 7 agglomerated trellis stages only taking the 3 minimums 
% (in absolute value) of the 7 agglomerated bits
%=========================================================
% Sort the agglomerated bits in ascending order of reliability
[dummy I] = sort(abs(LLRb));   %I = Index of the LLR in ascending order of reliability

% Generation of a LLRc vector by saturating the most 4 reliable LLRb.
LLRc = 10000*ones(1,7).*sign(LLRb);     
LLRc(I(1)) = LLRb(I(1));                       
LLRc(I(2)) = LLRb(I(2));
LLRc(I(3)) = LLRb(I(3));

% 7 trellis stages generation with LRRc
C(:,1) = A(:,1);
for i=1:7
    C(:,i+1) = max ( (G0*C(:,i) + LLRc(i)*ones(8,1)) , G1*C(:,i) - LLRc(i)*ones(8,1)); 
end;

% We check that the final result is the same
A(:,end)-A(1,end)
C(:,end)-C(1,end)
% Conclusion: It is the same

% ============================================
% Study case 6: Matlab code to generate array in paper
% "Simplified Compression of Redundancy Free Trellis Sections in
% Turbo-Decoder"
%=========================================================

%% version 28/08/2013 for article 3 min

clear all;

N = 16;    

LLR = floor(  15*(sign(rand(1,N)-0.5) + randn(1,N)));
LLR = min(LLR, 31);
LLR = max(LLR,-31);
LLR = [-14    31    24     12    31    20   6  -31     19    15   -19    -8    15    12     5   -11];

% Probability of the bits.
A(:,1) = floor(64*rand(8,1));
A(:,1) = A(:,1) - min(A(:,1));
A(:,1) = [35, 0, 20, 25, 17, 3, 16, 31]';
G0 = [1 0 0 0 0 0 0 0 ; 0 0 0 1 0 0 0 0 ; 0 0 0 0 1 0 0 0 ; 0 0 0 0 0 0 0 1 ;
        0 1 0 0 0 0 0 0 ; 0 0 1 0 0 0 0 0 ; 0 0 0 0 0 1 0 0 ; 0 0 0 0 0 0 1 0];
G1 = [0 1 0 0 0 0 0 0 ; 0 0 1 0  0 0 0 0 ; 0 0 0 0 0 1 0 0; 0 0 0 0 0 0 1 0; 
        1 0 0 0 0 0 0 0 ;  0 0 0 1 0 0 0 0 ; 0 0 0 0 1 0 0 0; 0 0 0 0 0 0 0 1];

% Generation trellis max-log-map
for i=1:N
    A(:,i+1) = max ( (G0*A(:,i) + LLR(i)*ones(8,1)) , G1*A(:,i) - LLR(i)*ones(8,1));
    A(:,i+1) = A(:,i+1) - min(A(:,i+1));
end;

% aglomeration of bits with same modulo 7 (min-sum).
for i=1:7
    LLRb(i) = 10000;
    for j=0:((N-i)/7)
        LLRb(i) = sign(LLRb(i))*sign(LLR(i+7*j))*min(abs(LLR(i+7*j)), abs(LLRb(i)));
    end
end;

% Permutation of (N mod 7).

q = mod(N,7);

% Realization of the q = N mod 7 shuffling steps.
B(:,1) = A(:,1);
for i=1:q
    B(:,i+1) = G0*B(:,i); 
end;

% Realization of the 7 trellis steps.
for i=q+1:q+7
    k = i- q - 1;
    ind = mod( k+q, 7) + 1;
    B(:,i + 1) = max ( (G0*B(:,i) + LLRb( ind )*ones(8,1)) , G1*B(:,i) - LLRb(ind)*ones(8,1)); 
    B(:,i+1) = B(:,i+1) - min(B(:,i+1));
end;


% Starting 3-min algorithm.

[dummy I] = sort(abs(LLRb));   %I = Index of LLR in accent order


% Generation of the LLRc vector with saturation of the 4 most relialable LLRb.
LLRc = 10000*ones(1,7).*sign(LLRb);     
LLRc(I(1)) = LLRb(I(1));                       
LLRc(I(2)) = LLRb(I(2));
LLRc(I(3)) = LLRb(I(3));

% Realisation of the q = N mod 7 steps of shuffling
C(:,1) = A(:,1);
for i=1:q
    C(:,i+1) = G0*C(:,i); 
end;

% Realization of the 7 trellis steps. 
for i=q+1:q+7
    k = i- q - 1;
    ind = mod( k+q, 7) + 1;
    C(:,i + 1) = max ( (G0*C(:,i) + LLRc( ind )*ones(8,1)) , G1*C(:,i) - LLRc(ind)*ones(8,1)); 
    if max(C(:,i+1)) > 10000
        C(:,i+1) = C(:,i+1)-10000;
    end;
    C(:,i+1) = C(:,i+1) - min(C(:,i+1));
end;


N = 16;

fid = fopen('table_a.txt','w');

fprintf(fid,'\\begin{table} \n');
fprintf(fid, '\\centering \n');
fprintf(fid,'\\caption{Classical treillis for $R$=16} \n');
fprintf(fid,'\\begin{tabular}{|c|c|c|c|c|c|c|c|c|c|c|c|c|} \\hline \n'); 
fprintf(fid,'time & 0 ');
    for i=2:3 
        fprintf(fid,'& %d ',i-1);
    end
    fprintf(fid,'& ... ');
    for i=N-5:N+1 
        fprintf(fid,'& %d ',i-1);
    end
fprintf(fid,'\\\\  \\hline \n');
for j=1:8
        fprintf(fid,'$\\alpha_%d$ & ',j-1);
        for i = 1 : 3
            fprintf(fid, ' %d &', A(j,i));
        end;
        fprintf(fid, '... &');
        for i = N-5 : N
            fprintf(fid, ' %d &', A(j,i));
        end;
        fprintf(fid, ' %d \\\\ \\hline \n ', A(j,N+1));
end;
fprintf(fid, '\\end{tabular}\n');
fprintf(fid, '	\\label{tab:full_trellis}\n');
fprintf(fid, '\\end{table}\n');

N = 7 + mod(N,7);

fprintf(fid,'\\begin{table} \n');
fprintf(fid, '\\centering \n');
fprintf(fid,'\\caption{compressed trellis (7-minimums)}\n');
fprintf(fid,'\\begin{tabular}{|c|c|c|c|c|c|c|c|c|c|c|} \\hline \n'); 
fprintf(fid,' time ');
for i=1:N+1 
         fprintf(fid,'& %d ',i-1);
 end
fprintf(fid,'\\\\  \\hline \n');
for j=1:8
         fprintf(fid,'$\\alpha_%d$ & ',j-1);
         for i = 1 : N
             fprintf(fid, ' %d &', B(j,i));
         end;
         fprintf(fid, ' %d \\\\ \\hline \n ', B(j,N+1));
end;
fprintf(fid, '\\end{tabular}\n');
fprintf(fid, '	\\label{tab:7min_trellis}\n');
fprintf(fid, '\\end{table}\n');

fprintf(fid,'\\begin{table} \n');
fprintf(fid, '\\centering \n');
fprintf(fid,'\\caption{compressed trellis (3-minimums)}\n');
fprintf(fid,'\\begin{tabular}{|c|c|c|c|c|c|c|c|c|c|c|} \\hline \n'); 
fprintf(fid,' time ');
for i=1:N+1 
         fprintf(fid,'& %d ',i-1);
 end
fprintf(fid,'\\\\  \\hline \n');
for j=1:8
         fprintf(fid,'$\\alpha_%d$ & ',j-1);
         for i = 1 : N
             fprintf(fid, ' %d &', C(j,i));
         end;
         fprintf(fid, ' %d \\\\ \\hline \n ', C(j,N+1));
end;
fprintf(fid, '\\end{tabular}\n');
fprintf(fid, '	\\label{tab:3min_trellis}\n');
fprintf(fid, '\\end{table}\n');

fclose(fid);
%     
  
max(max(abs(B-C)))
