clear;
clc;
BPH=2;
number_of_states=16;
fanout=2^BPH;
L=floor(log(number_of_states)/log(fanout));
nextstate=zeros(number_of_states,fanout); 
output=zeros(number_of_states,fanout); 
for i=0:number_of_states-1
    for j=0:fanout-1
        [next_state,out_put]=G_func(i,j,L,fanout);
        nextstate(i+1,j+1)=next_state;
        output(i+1,j+1)=out_put;
        input(i+1,next_state+1)=j;
    end
end
dfhtrel.numinput=fanout;
dfhtrel.numoutput=fanout*number_of_states;
dfhtrel.numstates=number_of_states;
dfhtrel.nextstates=nextstate;
dfhtrel.outputs=output;
dfhtrel.inputs=input;

K=3;
trel=poly2trellis(K,[7,5]);

tim=20;
rho_in_dB=11;
p=zeros(1,length(tim));
pb=0;
for i=1:tim
    p(i)=mlsd(rho_in_dB,dfhtrel,trel,K);
    pb=pb+p(i);
    i
end
pb=pb/tim;
sprintf('pb=%f',pb)