% ###############################################################################
% ##  Loesung: Serielle Verkettung                                             ##
% ##---------------------------------------------------------------------------##
% ##  Aufgabe 9.2.3) Serielle Verkettung zweier Faltungs-Codes                 ##
% ##---------------------------------------------------------------------------##
% ##  benoetigte Matlab-Dateien: serialconcat.m, iowef_block_conv.m,           ##
% ##                             pb_unionbound_awgn.m                          ##
% ###############################################################################

clear all;
close all;

disp('Aufgrund der Simulation kann die Berechnung sehr lange dauern...');
ant=input('Soll ein vorhandener Arbeitsbereich für Aufgabe 9.2.3 geladen werden (j/n)? ','s');
if (ant=='j')
   disp('Arbeitsbereich a9_2_3.mat geladen');
   file=1;
   load results/a9_2_3;
elseif (ant=='n')
   disp('Arbeitsbereich nicht geladen - simuliere neu!');
   file=0;
end;


%allgemeine Parameter
% Faltungscodes (oktale Darstellung)
Lc   = 3;
n    = 2;
k    = 1;

% Nichtrekursiver innerer Code
G1   = [5;7];
R1   = [0;0];

% Rekursiver innerer Code mit rekursivem Polynom 7
G2 = [1;5];
R2 = [0;7];

% Rekursiver innerer Code mit rekursivem Polynom 5
G3 = [1;7];
R3 = [0;5];

P    = [3;3];

EbN0 = 0:0.5:10;

% ##############################################################
% #####Aufgabenteil a) IOWEF der beiden Teilcodes L_pi=20  #####
% ##############################################################

disp('Aufgabenteil a) IOWEF der beiden Teilcodes fuer L_pi=20');

if(~file)
   N    = 20;              % Interleaver-Laenge
   dmax = N*n;
   lmax = N*k/n;
   wmax = lmax-(Lc-1);
   
   A1   = iowef_block_conv(Lc,n,k,G1,R1,N,lmax,wmax,P);
   A2   = iowef_block_conv(Lc,n,k,G1,R1,dmax,N,N,P);
   As   = serialconcat(A1,A2,N);
   
   [w,d] = size(As);
   cd20  = (1:w-1)*As(2:w,:)/N;
   
   % ####################################################################
   % #####Aufgabenteil b) Bitfehlerwahrscheinlichkeit fuer L_pi=20  #####
   % ####################################################################
   
   disp(sprintf('\nAufgabenteil b) Bitfehlerwahrscheinlichkeit fuer L_pi=20'));
   
   pb20 = pb_unionbound_awgn(EbN0,As,dmax,wmax);
   
   % ################################################################################
   % #####Aufgabenteil c) systematischer, rekursiver Faltungscode fuer L_pi=20  #####
   % ################################################################################
   
   disp(sprintf('\nAufgabenteil c) systematischer, rekursiver Faltungscode fuer L_pi=20'));
   A2     = iowef_block_conv(Lc,n,k,G2,R2,dmax,N,N,P);
   As     = serialconcat(A1,A2,N);
   
   [w,d]  = size(As);
   cd20r7= (1:w-1)*As(2:w,:)/N;
   
   pb20r7 = pb_unionbound_awgn(EbN0,As,dmax,wmax);
   
   % ##################################################################
   % #####Aufgabenteil d) veraenderte Rueckkopplung fuer L_pi=20  #####
   % ##################################################################
   
   disp(sprintf('\nAufgabenteil d) veraenderte Rueckkopplung fuer L_pi=20'));
   A2     = iowef_block_conv(Lc,n,k,G3,R3,dmax,N,N,P);
   As     = serialconcat(A1,A2,N);
   
   [w,d]  = size(As);
   cd20r5 = (1:w-1)*As(2:w,:)/N;
   
   pb20r5 = pb_unionbound_awgn(EbN0,As,dmax,wmax);
   
   
% ##############################################################
% #####Aufgabenteil a) IOWEF der beiden Teilcodes L_pi=50  #####
% ##############################################################

disp('Aufgabenteil a) IOWEF der beiden Teilcodes fuer L_pi=50');
   N    = 50;
   dmax = N*n;
   lmax = N*k/n;
   wmax = lmax-(Lc-1);
   
   A1   = iowef_block_conv(Lc,n,k,G1,R1,N,lmax,wmax,P);
   A2   = iowef_block_conv(Lc,n,k,G1,R1,dmax,N,N,P);
   As   = serialconcat(A1,A2,N);
   
   [w,d] = size(As);
   cd50  = (1:w-1)*As(2:w,:)/N;
   
   % ####################################################################
   % #####Aufgabenteil b) Bitfehlerwahrscheinlichkeit fuer L_pi=50  #####
   % ####################################################################
   
   disp(sprintf('\nAufgabenteil b) Bitfehlerwahrscheinlichkeit fuer L_pi=50'));
   pb50 = pb_unionbound_awgn(EbN0,As,dmax,wmax);
   
   % ################################################################################
   % #####Aufgabenteil c) systematischer, rekursiver Faltungscode fuer L_pi=50  #####
   % ################################################################################
   
   disp(sprintf('\nAufgabenteil c) systematischer, rekursiver Faltungscode fuer L_pi=50'));
   A2    = iowef_block_conv(Lc,n,k,G2,R2,dmax,N,N,P);
   As    = serialconcat(A1,A2,N);
   
   [w,d]   = size(As);
   cd_50r7 = (1:w-1)*As(2:w,:)/N;
   
   pb50r7  = pb_unionbound_awgn(EbN0,As,dmax,wmax);
   
   % ##################################################################
   % #####Aufgabenteil d) veraenderte Rueckkopplung fuer L_pi=50  #####
   % ##################################################################
   
   disp(sprintf('\nAufgabenteil d) veraenderte Rueckkopplung fuer L_pi=50'));
   A2     = iowef_block_conv(Lc,n,k,G3,R3,dmax,N,N,P);
   As     = serialconcat(A1,A2,N);
   
   [w,d]  = size(As);
   cd50r5 = (1:w-1)*As(2:w,:)/N;
   
   pb50r5 = pb_unionbound_awgn(EbN0,As,dmax,wmax);
   
   
   
% ###############################################################
% #####Aufgabenteil a) IOWEF der beiden Teilcodes L_pi=100  #####
% ###############################################################

disp('Aufgabenteil a) IOWEF der beiden Teilcodes fuer L_pi=100');
   N    = 100;
   dmax = N*n;
   lmax = N*k/n;
   wmax = lmax-(Lc-1);
   
   A1   = iowef_block_conv(Lc,n,k,G1,R1,N,lmax,wmax,P);
   A2   = iowef_block_conv(Lc,n,k,G1,R1,dmax,N,N,P);
   As   = serialconcat(A1,A2,N);
   
   [w,d]  = size(As);
   cd100  = (1:w-1)*As(2:w,:)/N;
   
   % #####################################################################
   % #####Aufgabenteil b) Bitfehlerwahrscheinlichkeit fuer L_pi=100  #####
   % #####################################################################
   
   disp(sprintf('\nAufgabenteil b) Bitfehlerwahrscheinlichkeit fuer L_pi=100'));
   pb100 = pb_unionbound_awgn(EbN0,As,dmax,wmax);
   
   % #################################################################################
   % #####Aufgabenteil c) systematischer, rekursiver Faltungscode fuer L_pi=100  #####
   % #################################################################################
   
   disp(sprintf('\nAufgabenteil c) systematischer, rekursiver Faltungscode fuer L_pi=100'));
   A2    = iowef_block_conv(Lc,n,k,G2,R2,dmax,N,N,P);
   As    = serialconcat(A1,A2,N);
   
   [w,d]    = size(As);
   cd_100r7 = (1:w-1)*As(2:w,:)/N;
   
   pb100r7 = pb_unionbound_awgn(EbN0,As,dmax,wmax);
   
   % ##################################################################
   % #####Aufgabenteil d) veraenderte Rueckkopplung fuer L_pi=100  #####
   % ##################################################################
   
   disp(sprintf('\nAufgabenteil d) veraenderte Rueckkopplung fuer L_pi=100'));
   A2     = iowef_block_conv(Lc,n,k,G3,R3,dmax,N,N,P);
   As     = serialconcat(A1,A2,N);
   
   [w,d]   = size(As);
   cd100r5 = (1:w-1)*As(2:w,:)/N;
   
   pb100r5 = pb_unionbound_awgn(EbN0,As,dmax,wmax);
   
   
% ###############################################################
% #####Aufgabenteil a) IOWEF der beiden Teilcodes L_pi=200  #####
% ###############################################################

disp('Aufgabenteil a) IOWEF der beiden Teilcodes fuer L_pi=200');
   N    = 200;
   dmax = N*n;
   lmax = N*k/n;
   wmax = lmax-(Lc-1);
   
   A1   = iowef_block_conv(Lc,n,k,G1,R1,N,lmax,wmax,P);
   A2   = iowef_block_conv(Lc,n,k,G1,R1,dmax,N,N,P);
   As   = serialconcat(A1,A2,N);
   
   [w,d] = size(A1);
   cd1   = (1:w-1)*A1(2:w,:)/N;
   [w,d] = size(A2);
   cd2   = (1:w-1)*A2(2:w,:)/N;
   [w,d] = size(As);
   cd200 = (1:w-1)*As(2:w,:)/N;
   
   % #####################################################################
   % #####Aufgabenteil b) Bitfehlerwahrscheinlichkeit fuer L_pi=200  #####
   % #####################################################################
   
   disp(sprintf('\nAufgabenteil b) Bitfehlerwahrscheinlichkeit fuer L_pi=200'));
   pb200 = pb_unionbound_awgn(EbN0,As,dmax,wmax);
   
   % #################################################################################
   % #####Aufgabenteil c) systematischer, rekursiver Faltungscode fuer L_pi=200  #####
   % #################################################################################
   
   disp(sprintf('\nAufgabenteil c) systematischer, rekursiver Faltungscode fuer L_pi=200'));
   A2      = iowef_block_conv(Lc,n,k,G2,R2,dmax,N,N,P);
   As      = serialconcat(A1,A2,N);
   
   [w,d]   = size(A2);
   cd2r7   = (1:w-1)*A2(2:w,:)/N;
   [w,d]   = size(As);
   cd200r7 = (1:w-1)*As(2:w,:)/N;
   
   pb200r7 = pb_unionbound_awgn(EbN0,As,dmax,wmax);
   
   % ###################################################################
   % #####Aufgabenteil d) veraenderte Rueckkopplung fuer L_pi=200  #####
   % ###################################################################
   
   disp(sprintf('\nAufgabenteil d) veraenderte Rueckkopplung fuer L_pi=200'));
   A2     = iowef_block_conv(Lc,n,k,G3,R3,dmax,N,N,P);
   As     = serialconcat(A1,A2,N);
   
   [w,d]   = size(A2);
   cd2r5   = (1:w-1)*A2(2:w,:)/N;
   [w,d]   = size(As);
   cd200r5 = (1:w-1)*As(2:w,:)/N;
   
   pb200r5 = pb_unionbound_awgn(EbN0,As,dmax,wmax)
end



% Ausgabe der Ergebnisse 
figure('name', 'Koeffizienten und Union Bound', 'numbertitle', 'off'); 
%subplot(121);
gh=semilogy(0:length(cd1)-1,cd1,'k-',...
   0:length(cd20)-1,cd20,'k--',0:length(cd50)-1,cd50,'k-.',...
   0:length(cd100)-1,cd100,'k.-',0:length(cd200)-1,cd200,'k:');
%set(gh,'linewidth',1.5,'markersize',8);
grid on;
axis([0 40 10^-10 1e15]);
set(gca,'fontname','times','fontsize',15);
xlabel('x');
ylabel('y');
set(gcf,'PaperPosition',[2 5 4.7 3.4]);
title('t');
set(gca,'ytick',[10^-10 10^-5 10^-0 10^5 10^10 10^15 ],'xtick',[0 5 10 15 20 25 30 35 40]);
%xlabel('d \rightarrow'); 
%ylabel('c_d \rightarrow');
legend('Teilcode','Lpi=20','Lpi=50','Lpi=100','Lpi=200',2);
%title('a) Koeffizienten c_d fuer Verkettung nichtrekursiver Codes');
print -deps b9_2_6a

%subplot(122)
figure;
gh=semilogy(EbN0,pb20,'k-s',EbN0,pb50,'k-x',EbN0,pb100,'k-o',EbN0,pb200,'k-d');
%set(gh,'linewidth',1.5,'markersize',8);
grid on;
set(gca,'fontname','times','fontsize',15);
xlabel('x');
ylabel('y');
set(gcf,'PaperPosition',[2 5 4.7 3.4]);
title('t');
set(gca,'ytick',[10^-10 10^-8 10^-6 10^-4 10^-2 10^0 ],'xtick',[0 2 4 6 8 10]);
%xlabel('E_b/N_0 in dB \rightarrow');
%ylabel('P_b \rightarrow');
%title('b) Union Bound fuer serielle Verkettung, NSC innen');
legend('Lpi=20','Lpi=50','Lpi=100','Lpi=200');
axis([0 10 10^-10 1]);
print -deps b9_2_6b

figure('name', 'unterschiedliche Rueckkopplungspolynome', 'numbertitle', 'off'); 
%subplot(121)
gh=semilogy(EbN0,pb20r7,'k-s',EbN0,pb50r7,'k-x',EbN0,pb100r7,'k-o',EbN0,pb200r7,'k-d');
%set(gh,'linewidth',1.5,'markersize',8);
grid on;
set(gca,'fontname','times','fontsize',15);
xlabel('x');
ylabel('y');
set(gcf,'PaperPosition',[2 5 4.7 3.4]);
title('t');
set(gca,'ytick',[10^-10 10^-8 10^-6 10^-4 10^-2 10^0 ],'xtick',[0 2 4 6 8 10]);
%xlabel('E_b/N_0 in dB \rightarrow');
%ylabel('P_b \rightarrow');
%title('a) Serielle Verkettung, innen RSC (7)');
legend('Lpi=20','Lpi=50','Lpi=100','Lpi=200');
axis([0 10 10^-10 1]);
print -deps b9_2_7a

%subplot(122)
figure
gh=semilogy(EbN0,pb20r5,'k-s',EbN0,pb50r5,'k-x',EbN0,pb100r5,'k-o',EbN0,pb200r5,'k-d');
%set(gh,'linewidth',1.5,'markersize',8);
grid on;
set(gca,'fontname','times','fontsize',15);
xlabel('x');
ylabel('y');
set(gcf,'PaperPosition',[2 5 4.7 3.4]);
title('t');
set(gca,'ytick',[10^-10 10^-8 10^-6 10^-4 10^-2 10^0 ],'xtick',[0 2 4 6 8 10]);
%xlabel('E_b/N_0 in dB \rightarrow');
%ylabel('P_b \rightarrow');
%title('b) Serielle Verkettung, innen RSC (5)');
legend('Lpi=20','Lpi=50','Lpi=100','Lpi=200');
axis([0 10 10^-10 1]);
print -deps b9_2_7b

disp('Weiter mit beliebiger Taste');
%pause;

% #########################################################################
% #####Aufgabenteil e) Punktierung des verketteten Codes auf R_c=1/3  #####
% #########################################################################

disp(sprintf('\nAufgabenteil e) Punktierung des verketteten Codes auf R_c=1/3'));

N    = 100;
dmax = N*n;
lmax = N*k/n;
wmax = lmax-(Lc-1);

% Punktierung des inneren Codes  
A1   = iowef_block_conv(Lc,n,k,G1,R1,N,lmax,wmax,P);
P    = [3;1];
A2   = iowef_block_conv(Lc,n,k,G3,R3,dmax,N,N,P);
As1  = serialconcat(A1,A2,N);

% Punktierung des aeusseren Codes
A1   = iowef_block_conv(Lc,n,k,G1,R1,N,lmax,wmax,P);
N    = 75;
dmax = N*n;
P    = [3;3];
A2   = iowef_block_conv(Lc,n,k,G3,R3,dmax,N,N,P);
As2  = serialconcat(A1,A2,N);

% Abschaetzung der Fehlerrate mit Union Bound
pb1 = pb_unionbound_awgn(EbN0,As1,150,48);
pb2 = pb_unionbound_awgn(EbN0,As2,150,48);

figure('name', 'Union Bound-Abschaetzung fuer verschiedene Punktierungen', 'numbertitle', 'off'); 
gh=semilogy(EbN0,pb1,'k-s',EbN0,pb2,'k-x',EbN0,pb100r5,'k-o');
%set(gh,'linewidth',1.5,'markersize',8);
grid on;
set(gca,'fontname','times','fontsize',13);
xlabel('x');
ylabel('y');
set(gcf,'PaperPosition',[2 5 4.7 3.4]);
title('t');
set(gca,'ytick',[10^-10 10^-8 10^-6 10^-4 10^-2 10^0 ],'xtick',[0 2 4 6 8 10]);
%xlabel('E_b/N_0 in dB \rightarrow');
%ylabel('P_b \rightarrow');
%title('Punktierung der Teilcodes');
legend('innen','aussen','Rc=1/3');
axis([0 10 10^-10 1]);
print -deps b9_2_8