% Fulvio Babich, University of Trieste, babich@phd.units.it
% 19/12/2016
% ============================================================================================
% Script containing some examples about parity decoding:
clear;
% SEZIONE DI INPUT
%i dati sono un ipercubo di dimensione dim e lato R
R=5;
dim=3; %<=4
iterations=3; %numero massimo iterazioni
% intervallo di test
ebn0dB=[0.0:0.5:5.5]; 
%numero di simulazioni
max_test=1000000;
%condizone di arresto
max_err=100;
%intervallo visualizzazione
interval=20000;
% FINE SEZIONE DI INPUT

res=sprintf('R=%d, dim=%d, iterations=%d', R,dim,iterations); disp(res);

%INIZIALIZZAZIONE
%dimensione blocoo
N=(R+1)^dim;
%tasso del codice
Rc=(R/(R+1))^dim;
%distanza minima
dmin=2^dim;
%numero di parole a distanza minima
Ndmin=R^dim*(R+1)^dim/dmin;
switch dim
    case 2
        word=ones(R+1,R+1);
    case 3
        word=ones(R+1,R+1,R+1);
    case 4
        word=ones(R+1,R+1,R+1,R+1);
    otherwise
        disp('errore dim')
        return;
end
snrv=10.^(ebn0dB/10)*Rc;
% approssimazione probabilità d'errore di parola teorica
pwt=1/2*Ndmin*erfc((snrv*dmin).^.5);
% approssimazione probabilità d'errore di bit teorica
pbt=pwt/N;
%probabilità d'errore di bit senza codice
pbnc=1/2*erfc((10.^(ebn0dB/10)).^.5);
%inizializzazione contatori
snr_l=length(snrv);
errors_s=zeros(1,snr_l);
errors_h=zeros(1,snr_l);
errors_s2=zeros(1,snr_l);
errors_b=zeros(1,snr_l);
test_done=zeros(1,snr_l);

%SIMULAZIONE
for il=1:snr_l
    snr=snrv(il);
    while test_done(il)<max_test & errors_s(il)<max_err
        test_done(il)=test_done(il)+1;
        noise=randn(size(word))/(2*snr)^.5;
        rec=word+noise;
        llr=4*rec*snr^.5;
        if dim==2
            errors=sum(sum(llr<0));
        %hard decoding
            if errors > 1
                errors_h(il)=errors_h(il)+1;
            end
        end
if rem(test_done(il),interval)==0, res=sprintf('ebn0dB %.1f, test %6d, errors %2d', ebn0dB(il),test_done(il),errors_s(il)); disp(res); end
%soft decoding
        llrs=llr;
        errors=1; ii=0;
        while errors & ii < iterations
            ii=ii+1;
            switch dim
                case 2
                    for i=1:R+1, temp=tanh(llrs/2); temp(i,:)=1; px(i,:)=2*atanh(prod(temp,1)); end, llrs=llrs+px;
                    for i=1:R+1, temp=tanh(llrs/2); temp(:,i)=1; px(:,i)=2*atanh(prod(temp,2)); end, llrs=llrs+px;
                    errors=sum(sum(llrs<0));
                case 3
                    for i=1:R+1, temp=tanh(llrs/2); temp(i,:,:)=1; px(i,:,:)=2*atanh(prod(temp,1)); end, llrs=llrs+px;
                    for i=1:R+1, temp=tanh(llrs/2); temp(:,i,:)=1; px(:,i,:)=2*atanh(prod(temp,2)); end, llrs=llrs+px;
                    for i=1:R+1, temp=tanh(llrs/2); temp(:,:,i)=1; px(:,:,i)=2*atanh(prod(temp,3)); end, llrs=llrs+px;
                    errors=sum(sum(sum(llrs<0)));
                case 4
                    for i=1:R+1, temp=tanh(llrs/2); temp(i,:,:,:)=1; px(i,:,:,:)=2*atanh(prod(temp,1)); end, llrs=llrs+px;
                    for i=1:R+1, temp=tanh(llrs/2); temp(:,i,:,:)=1; px(:,i,:,:)=2*atanh(prod(temp,2)); end, llrs=llrs+px;
                    for i=1:R+1, temp=tanh(llrs/2); temp(:,:,i,:)=1; px(:,:,i,:)=2*atanh(prod(temp,3)); end, llrs=llrs+px;
                    for i=1:R+1, temp=tanh(llrs/2); temp(:,:,:,i)=1; px(:,:,:,i)=2*atanh(prod(temp,4)); end, llrs=llrs+px;
                    errors=sum(sum(sum(sum(llrs<0))));
            end
        end
        if errors > 0
            errors_s(il)=errors_s(il)+1;
            errors_b(il)=errors_b(il)+errors;
            if errors > floor((dmin-1)/2)
               errors_s2(il)=errors_s2(il)+1;
            end
        end
    end
    res=sprintf('ebn0dB %.1f, test %6d, errors %2d', ebn0dB(il),test_done(il),errors_s(il)); disp(res);
end
pws=errors_s./test_done;
pbs=errors_b./test_done/N;
phs=errors_h./test_done;
figure(1)
semilogy(ebn0dB,pws,'r:o');
set(gca,'FontName','TimesNewRoman');
set(gca,'FontSize',12);
set(gca,'YLim',[1e-4,1]);
xlabel('E_b/N_0 [dB]');
ylabel('BlER');
hold on
semilogy(ebn0dB,pwt,'r--');
if dim==2
    semilogy(ebn0dB,phs,'b--*');
end
figure(2)
semilogy(ebn0dB,pbs,'r:o');
set(gca,'FontName','TimesNewRoman');
set(gca,'FontSize',12);
set(gca,'YLim',[1e-4,1]);
xlabel('E_b/N_0 [dB]');
ylabel('BER');
hold on
semilogy(ebn0dB,pbt,'r--');
%semilogy(ebn0dB,pbnc,'k');