%%============================================================
%% This file related experimentation to find simple method to 
%% Quantize the LLR using 8-PSK constellation. 
%%E. Boutillon, 13/08/2011
%%============================================================
function [llr0 llr1 llr2 llr0q llr1q llr2q] = Quantif_LRR(R, alpha)

% Number of bit to quantized rau : 5 between 0 and 3.1 with step of 0.1
% Number of bit to quantized theta : 7 (between 0 and 2pi).

% Values of sigma for the different code rate using 8-PSK
% Rate 3/5, sigma = sqrt(0.278705347
% Rate 2/3, sigma = sqrt(0.215662519)
% Rate 3/4, sigma = sqrt(0.160234798)
% Rate 5/6, sigma = sqrt(0.115148896)
% Rate 8/9, sigma = sqrt(0.084649243)
% Rate 9/10, sigma = sqrt(0.079184924)

%% Simulation with rate 2/3
close
sigma = sqrt(0.215662519);

%%===================================================================
% We search the scaling factor for LLR quantization so that 1/31 of the
% value are equal to +15 after quantization of LLR.
% The formula used to compute quantized LLR_Q from LLR is given by:
% LLR_Q = floor( LLR*14.5/R + 0.5);
%===================================================================
%R = loi_quantif_LLR(sigma)    

%===================================================================
% Generation of the 8-PSK modulation
%===================================================================
coord = (1:8);
constellation = exp(2*pi*1i*coord/8);
plot(constellation);

%===================================================================
% Quantization of the Polar_quantized_LLR value.
%===================================================================
for rau = 1:32
    for theta = 1:128
        y     = rau*0.1*exp( 1i*2*pi*(theta-1)/128);
        LLR   = LLR_compute_8PSK( constellation, y, sigma);
        llr_Q = floor(LLR*14.5/R + 0.5);
        llr0(rau, theta) = max(-15, min( 15, llr_Q(1)));
        llr1(rau, theta) = max(-15, min( 15, llr_Q(2)));
        llr2(rau, theta) = max(-15, min( 15, llr_Q(3)));
        llr0q(rau, theta)= QuantizationLLR0( rau, theta,alpha);
        llr1q(rau, theta)= QuantizationLLR1( rau, theta,alpha);
        llr2q(rau, theta)= QuantizationLLR2( rau, theta,alpha);
    end
end

figure(1)
set(gcf,'DefaultTextFontName','Time New Roman')
set(gcf,'DefaultTextFontSize',14)
set(gcf,'DefaultTextInterpreter','latex')
plot(0:128, [llr0(10,:) llr0(10,1)],'b')
hold on
plot(0:128, [llr1(10,:) llr1(10,1)],'k')
plot(0:128, [llr2(10,:) llr2(10,1)],'r')
xlabel(1,'string','$\bar{\theta}$ on $n_{\theta}$ = 7 bits','interpreter','latex', 'fontsize',12,'units','norm');
ylabel(1,'string','$\bar{L}^i(\rho,\theta)$ on $n_b = 5$ bits','interpreter','latex', 'fontsize',12);
title(1,'string','Variation of LLR pour $\rho$ = 1 as a function of $\bar{\theta}$','interpreter','latex', 'fontsize',12);
text('string','$\bar{L}^0$','interpreter','latex', 'fontsize',12, 'units','norm', 'pos', [0.8 0.88]);
text('string','$\bar{L}^1$','interpreter','latex', 'fontsize',12, 'units','norm', 'pos', [0.55 0.88]);
text('string','$\bar{L}^2$','interpreter','latex', 'fontsize',12, 'units','norm', 'pos', [0.30 0.65]);
axis([ 0 128 -15 15])
grid on


figure(2)
set(gcf,'DefaultTextFontName','Time New Roman')
set(gcf,'DefaultTextFontSize',14)
set(gcf,'DefaultTextInterpreter','latex')
plot(0:128, [llr0(5,:) llr0(5,1)],'b')
hold on
plot(0:128, [llr0(10,:) llr0(10,1)],'k')
plot(0:128, [llr0(20,:) llr0(20,1)],'r')
%plot(0:128, [llr0(30,:) llr0(30,1)],'m')
xlabel(1,'string','$\bar{\theta}$ on $n_\theta = 7$ bits','interpreter','latex', 'fontsize',12,'units','norm');
ylabel(1,'string','$\bar{L}^0(\rho,\theta)$','interpreter','latex', 'fontsize',12);
title(1,'string','$\bar{L}^0$ as function of $\theta$ for several values of $\rho$, ','interpreter','latex', 'fontsize',12);
text('string','$\rho=0.5$','interpreter','latex', 'fontsize',12, 'units','norm', 'pos', [0.25 0.4]);
text('string','$\rho=1.0$','interpreter','latex', 'fontsize',12, 'units','norm', 'pos', [0.25 0.25]);
text('string','$\rho=2.0$','interpreter','latex', 'fontsize',12, 'units','norm', 'pos', [0.25 0.05]);
%text('string','$\rho=3.0$','interpreter','latex', 'fontsize',12, 'units','norm', 'pos', [0.25 0.05]);
axis([ 0 128 -15 15])
grid on

figure(3)
set(gcf,'DefaultTextFontName','Time New Roman')
set(gcf,'DefaultTextFontSize',14)
set(gcf,'DefaultTextInterpreter','latex')
plot(0:128, [llr0(5,:) llr0(5,1)],'b')
hold on
plot(0:128, [llr0(10,:) llr0(10,1)],'k')
plot(0:128, [llr0(20,:) llr0(20,1)],'r')
plot(0:128, [llr0(30,:) llr0(30,1)],'m')

plot(0:128, [llr0q(5,:) llr0q(5,1)],'.b')
plot(0:128, [llr0q(10,:) llr0q(10,1)],'.k')
plot(0:128, [llr0q(20,:) llr0q(20,1)],'.r')
plot(0:128, [llr0q(30,:) llr0q(30,1)],'.m')
xlabel(1,'string','$\bar{\theta}$ coded on $n_\theta = 7$ bits','interpreter','latex', 'fontsize',12,'units','norm');
ylabel(1,'string','$\bar{L}^0(\rho,\theta)$','interpreter','latex', 'fontsize',12);
title(1,'string','$\bar{L}^0$ as a function of $\theta$ for several values of $\rho$','interpreter','latex', 'fontsize',12);
text('string','$\rho=0.5$','interpreter','latex', 'fontsize',12, 'units','norm', 'pos', [0.25 0.4]);
text('string','$\rho=1.0$','interpreter','latex', 'fontsize',12, 'units','norm', 'pos', [0.25 0.25]);
text('string','$\rho=2.0$','interpreter','latex', 'fontsize',12, 'units','norm', 'pos', [0.25 0.05]);
axis([ 0 128 -15 15])
grid on

k=1;
for i=1:32
   for j=1:128 
       xx(k) = i*cos((j-1)*2*pi/128);
       yy(k) = i*sin((j-1)*2*pi/128);
       kk(k) = llr0q(i,j) - llr0(i,j);
       k = k+1;
   end
end
figure(4)
plot3(xx,yy,kk);
grid on;
sum( abs(kk))
max( abs(kk))
figure(5)
plot(llr0(1,1:128) - llr0q(1,1:128));
hold on;
plot(llr0(1,1:128),'*');
plot(llr0(10,1:128) - llr0q(10,1:128),'r');
plot(llr0(20,1:128) - llr0q(20,1:128),'g');
R
%===================================================================
% Function "bit thrue" approximation quantization LLR0.
%    S is the scaling factors. 
%===================================================================
function llrQ = QuantizationLLR0( rau, theta, S);

    % rotation of theta of pi/8 to have a simetrical distribution
    theta = theta - 1;
    theta = theta - 8;
    if theta < 0;
        theta = theta + 128;
    end;
    
    % quantization of LLR0 : first step: use symetries
    sign = -1;
    if theta > 64 
        theta = theta - 64;
        sign = 1;
    else
        theta = theta;
    end;
    if theta > 32
        theta = 64 - theta;
    end;
    llrQ = sign*min( 15, floor(rau*min(theta, 30)*S + 0.5));
end;

%===================================================================
% Function "bit thrue" approximation quantization LLR0.
%    S is the scaling factors. 
%===================================================================
function llrQ = QuantizationLLR1( rau, theta, S);

    % rotation of theta of pi/8 to have a simetrical distribution
    % + rotation of theta of pi/2 to have the same computation.
    theta = theta - 1;
    theta = theta - 8;
    if theta < 0;
        theta = theta + 128;
    end;
    % + rotation of theta of pi/2 to have the same computation.
    theta = theta + 32;
    if theta > 128;
        theta = theta - 128;
    end;
    
     % quantization of LLR1 : first step: use symetries
    sign = -1;
    if theta > 64 
        theta = theta - 64;
        sign = 1;
    else
        theta = theta;
    end;
    if theta > 32
        theta = 64 - theta;
    end;
    llrQ = sign*min( 15, floor(rau*min(theta, 30)*S + 0.5));
end;

%===================================================================
% Function "bit thrue" approximation quantization LLR2.
%    S is the scaling factors. 
%===================================================================
function llrQ = QuantizationLLR2( rau, theta, S);

    % rotation of theta of pi/8 to have a simetrical distribution
    theta = theta - 1;
    theta = theta - 8;
    if theta < 0;
        theta = theta + 128;
    end;
 
    % quantization of LLR2 : first step: use symetries
      if theta > 64 
         theta = theta-64;
      else
         theta = theta;
      end;
      sign = 1;
      if theta > 32
          theta = theta - 32;
          sign = -1;
      end;
      llrQ = sign*min( 15, floor(rau*min(theta-16, 15)*S + 0.5));

end;    
    

function R = loi_quantif_LLR(sigma)
coord = 1:8;
constellation = exp(2*pi*1i*coord/8);

%% First step : quantization of input symbols.

N = 1000;

y = [];
for j = 1:8;
   y = [y (constellation(j) + sigma*(randn(1,N) + 1i*randn(1,N)))];
end;

llr_acc = [];
for i=1:length(y)
    llr = LLR_compute_8PSK( constellation, y(i), sigma);
    llr_acc = [llr_acc llr];
end;    
 
 R = 0; delta = 2;
hist( llr_acc, 30);

   nb_above_threahold = sum( (sign(abs(llr_acc) - R)+1)/2)
   while nb_above_threahold ~= floor(8*N/31);
       if nb_above_threahold > floor(8*N/31);
           R = R + delta;
       else
           delta = delta/2;
           R = R - delta;
       end;
       nb_above_threahold = sum( (sign(llr_acc - R)+1)/2);
   end
end;


    
    function LLR = LLR_compute_8PSK( const, y, sigma);
    
    dist = (const - y).*conj(const-y);
    
    dist = -dist/(2*sigma^2);
    
    P1 = exp(dist(5))+exp(dist(6))+exp(dist(7))+exp(dist(8));
    P0 = exp(dist(1))+exp(dist(2))+exp(dist(3))+exp(dist(4));
    LLR(1) = log( P1/P0);
    
    P1 = exp(dist(3))+exp(dist(4))+exp(dist(5))+exp(dist(6));
    P0 = exp(dist(1))+exp(dist(2))+exp(dist(7))+exp(dist(8));
    LLR(2) = log( P1/P0);
   
    P1 = exp(dist(2))+exp(dist(3))+exp(dist(6))+exp(dist(7));
    P0 = exp(dist(1))+exp(dist(4))+exp(dist(5))+exp(dist(8));
    LLR(3) = log( P1/P0);

    end;
    
    end%     
%     
% %     
% function ent = Quantif_input
% coord = 1:8;
% constellation = exp(2*pi*1i*coord/8);
% 
% %% First step : quantization of input symbols.
% 
% N = 100000;
% 
% sigma = sqrt(0.3);
% y = [];
%     for j = 1:8;
%             y = [y (constellation(j) + sigma*(randn(1,N) + 1i*randn(1,N)))];
%     end;
% 
% y = [abs(real(y)) abs(imag(y))];
%     
% for i = 1:20;
%         delta = i*0.02;
%         for k=1:15
%             seuil = k*delta - delta/2;
%             f(k)  = -sum( sign(y - seuil) -1)/2;
%             if k==1 
%                 p(k) = f(k);
%             else
%                 p(k) = f(k) - f(k-1);
%             end
%         end
%     p(16) = 16*N - f(15);
%     p = p/(16*N);
%     ent(i,:) = p;
%     entropy(i) = sum( -p.*log(p));
% end;
% end;    
%  
%     