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

pbr=0.01;%modelul Bernoulli 
N=15;
m1=1:13;
Nfact=prod(1:N);
mfact=gamma(m1+1);
Nmfact=gamma(N-m1+1);
CNm=Nfact./(mfact.*Nmfact);
T1=pbr.^m1;
T2=(1-pbr).^(N-m1);
T12=T1.*T2;
PbrmN=CNm.*T12;

p=0.187966*10^(-2);%modelul Gilbert-Elliot
q=0.080617;
pg=0.563331*10^(-3);
pb=0.417201;
m=12;
a1=1:12;
a=1:13;

%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
PgemN=(q/(p+q))*g(2:13,N).'+(p/(p+q))*b(2:13,N).';

p11=0.998511; %modelul Fritchman cu trei stari
p31=0.142542;
p13=0.001489;
p32=0.396651;
p23=0.135132;
p22=0.864868;
p33=1-p31-p32;
z1=p31/p13/(1+p31/p13+p32/p23);
z2=p32/p23/(1+p31/p13+p32/p23);
z3=1/(1+p31/p13+p32/p23);
c1=zeros(m+1,N+1);
c2=zeros(m+1,N+1);
c3=zeros(m+1,N+1);
c1(1,1)=1; c2(1,1)=1; c3(1,1)=1;
c1(1,2)=1; c1(1,2)=1; c3(1,2)=0;
c1(2,2)=0; c2(2,2)=0; c3(2,2)=1;

for k=3:N+1
  c3(1,k)=0;
  c1(1,k)=c1(1,k-1)*p11+c3(1,k-1)*p13;
  c2(1,k)=c2(1,k-1)*p22+c3(1,k-1)*p23;
end

for i=2:m+1
   for j=2:N+1
      c1(i,j)=c1(i,j-1)*p11+c3(i,j-1)*p13;
      c2(i,j)=c2(i,j-1)*p22+c3(i,j-1)*p23;
      c3(i,j)=c1(i-1,j-1)*p31+c2(i-1,j-1)*p32+c3(i-1,j-1)*p33;
   end
end

PfrmN=z1.*c1(2:13,N+1)'+z2.*c2(2:13,N+1)'+z3.*c3(2:13,N+1)';

loglog(a,PbrmN,a1,PgemN,a1,PfrmN)

axis([1 13 10^(-5) 0.15])  
title('Ponderea de eroare P(m,N) pentru modelele de canal tip Bernoulli, Gilbert-Elliot si Fritchman');
grid
xlabel('m');
ylabel('Ponderea de eroare P(m,N)');
text(2.4,0.00015,'Bernoulli');
text(6.5,0.005,'Gilbert-Elliot');
text(6.5,0.00015,'Fritchman'); 