%Ponderea de eroare P(m,N) pentru modelele de canal tip tandem Bernoulli, tandem Gilbert-Elliot, tandem Bernoulli-Gilbert

pbr=0.01;%modelul tandem Bernoulli 
ptbr=2*pbr*(1-pbr);
m=12;
N=15;
m1=1:13;
a=1:m+1;
a1=1:m;
Nfact=prod(1:N);
mfact=gamma(m1+1);
Nmfact=gamma(N-m1+1);
CNm=Nfact./(mfact.*Nmfact);
T1=ptbr.^m1;
T2=(1-ptbr).^(N-m1);
T12=T1.*T2;
PtbrmN=CNm.*T12;

p1=0.240053*10^(-2);%modelul tandem Gilbert-Elliot
q1=0.100485;
p2=0.278662*10^(-2);
q2=0.119516;
pg1=0.220591*10^(-3);
pg2=0.563331*10^(-3);
pb1=0.419360;
pb2=0.417201;
h=(p1+q1)*(p2+q2);
z1=q1*q2/h;
z2=q1*p2/h;
z3=p1*q2/h;
z4=p1*p2/h;
P1=pg1+pg2-2*pg1*pg2;
P2=pg1+pb2-2*pg1*pb2;
P3=pb1+pg2-2*pb1*pg2;
P4=pb1+pb2-2*pb1*pb2;
p11=(1-p1)*(1-p2); p12=(1-p1)*p2; p13=p1*(1-p2); p14=p1*p2;
p21=(1-p1)*q2;     p22=(1-p1)*p2; p23=p1*(1-p2); p24=p1*p2;
p31=q1*(1-p2);     p32=q1*p2;     p33=p1*(1-p2); p34=(1-q1)*p2;
p41=q1*q2;         p42=q1*p2;     p43=p1*(1-p2); p44=(1-q1)*(1-q2);
c1=zeros(m+1,N+1);
c2=zeros(m+1,N+1);
c3=zeros(m+1,N+1);
c4=zeros(m+1,N+1);
c1(1,1)=0;    c2(1,1)=0;    c3(1,1)=0;    c4(1,1)=0;
c1(1,2)=1-P1; c2(1,2)=1-P2; c3(1,2)=1-P3; c4(1,2)=1-P4;
c1(2,2)=P1;   c2(2,2)=P2;   c3(2,2)=P3;   c4(2,2)=P4;

for k=3:N+1
  c1(1,k)=c1(1,k-1)*p11*(1-P1)+c2(1,k-1)*p12*(1-P1)+c3(1,k-1)*p13*(1-P1)+c4(1,k-1)*p14*(1-P1);
  c2(1,k)=c1(1,k-1)*p21*(1-P2)+c2(1,k-1)*p22*(1-P2)+c3(1,k-1)*p23*(1-P2)+c4(1,k-1)*p24*(1-P2);
  c3(1,k)=c1(1,k-1)*p31*(1-P3)+c2(1,k-1)*p32*(1-P3)+c3(1,k-1)*p33*(1-P3)+c4(1,k-1)*p34*(1-P3);
  c4(1,k)=c1(1,k-1)*p41*(1-P4)+c2(1,k-1)*p42*(1-P4)+c3(1,k-1)*p43*(1-P4)+c4(1,k-1)*p44*(1-P4);
end

for i=2:m+1,
   for j=2:N+1,
     c1(i,j)=c1(i,j-1)*p11*(1-P1)+c2(i,j-1)*p12*(1-P1)+c3(i,j-1)*p13*(1-P1)+c4(i,j-1)*p14*(1-P1)+...
        c1(i-1,j-1)*p11*P1+c2(i-1,j-1)*p12*P1+c3(i-1,j-1)*p13*P1+c4(i-1,j-1)*p14*P1;
     c2(i,j)=c1(i,j-1)*p21*(1-P2)+c2(i,j-1)*p22*(1-P2)+c3(i,j-1)*p23*(1-P2)+c4(i,j-1)*p24*(1-P2)+...
        c1(i-1,j-1)*p21*P2+c2(i-1,j-1)*p22*P2+c3(i-1,j-1)*p23*P2+c4(i-1,j-1)*p24*P2;
     c3(i,j)=c1(i,j-1)*p31*(1-P3)+c2(i,j-1)*p32*(1-P3)+c3(i,j-1)*p33*(1-P3)+c4(i,j-1)*p34*(1-P3)+...
        c1(i-1,j-1)*p31*P3+c2(i-1,j-1)*p32*P3+c3(i-1,j-1)*p33*P3+c4(i-1,j-1)*p34*P3;
     c4(i,j)=c1(i,j-1)*p41*(1-P4)+c2(i,j-1)*p42*(1-P4)+c3(i,j-1)*p43*(1-P4)+c4(i,j-1)*p44*(1-P4)+...
        c1(i-1,j-1)*p41*P4+c2(i-1,j-1)*p42*P4+c3(i-1,j-1)*p43*P4+c4(i-1,j-1)*p44*P4;
   end
end

PtgemN=z1.*c1(2:m+1,N+1)'+z2.*c2(2:m+1,N+1)'+z3.*c3(2:m+1,N+1)'+z4.*c4(2:m+1,N+1)';

p=0.278662*10^(-2);%modelul tandem Bernoulli-Gilbert
q=0.119516;
pg=0.105521*10^(-1);
pb=0.418857;
%crearea matricelor initiale
g=zeros(m+1,N);
b=zeros(m+1,N);
g(1,:)=(1-pg)*[1 0 0 0 0 0 0 0 0 0 0 0 0 0 0];
g(2,:)=pg*[1 0 0 0 0 0 0 0 0 0 0 0 0 0 0];
b(1,:)=(1-pb)*[1 0 0 0 0 0 0 0 0 0 0 0 0 0 0];
b(2,:)=pb*[1 0 0 0 0 0 0 0 0 0 0 0 0 0 0];

%calculul elementelor prin recurenta
for k=2:N
   g(1,k)=g(1,k-1)*(1-p)*(1-pg)+b(1,k-1)*p*(1-pg);
   b(1,k)=b(1,k-1)*(1-q)*(1-pb)+g(1,k-1)*q*(1-pb);
end

for i=2:m,
   for j=2:N,
     g(i,j)=g(i,j-1)*(1-p)*(1-pg)+g(i-1,j-1)*(1-p)*pg+b(i,j-1)*p*(1-pg)+b(i-1,j-1)*p*pg;
     b(i,j)=b(i,j-1)*(1-q)*(1-pb)+b(i-1,j-1)*(1-q)*pb+g(i,j-1)*q*(1-pb)+g(i-1,j-1)*q*pb;
   end
end
%calculul ponderii de eroare
PtbgmN=(q/(p+q))*g(2:13,N).'+(p/(p+q))*b(2:13,N).';

loglog(a,PtbrmN,a1,PtgemN,a1,PtbgmN)
axis([1 13 10^(-5) 1])
title('Ponderea de eroare P(m,N) pentru modelele de canal tip tandem Bernoulli, tandem Gilbert-Elliot si tandem Bernoulli-Gilbert');
grid
xlabel('m');
ylabel('Ponderea de eroare P(m,N)');
%text(2.4,0.00015,'TBR');
%text(6.5,0.005,'TGE');
%text(6.5,0.00015,'TBG'); 