function pb=symbysym(rho_in_dB,BPH,number_of_states)

N=1000;
fanout=2^BPH;
rho=10^(rho_in_dB/10);
source=randint(1,N,fanout);
L=floor(log(number_of_states)/log(fanout));
depth_of_trellis=length(source);

nextstate=zeros(number_of_states,fanout); 
formerstate=number_of_states.*ones(number_of_states,fanout,fanout);
output=zeros(number_of_states,fanout); 
number_of_out=number_of_states*fanout;
for i=0:number_of_states-1
    for j=0:fanout-1
        k=1;
        [next_state,out_put]=G_func1(i,j,L,fanout);
        nextstate(i+1,j+1)=next_state;
        output(i+1,j+1)=out_put;
        while(formerstate(next_state+1,k,j+1)~=number_of_states)
            k=k+1;
        end
        formerstate(next_state+1,k,j+1)=i; 
    end
end

E=1;
sgma=sqrt(E/(BPH*2*rho));
demod_input=zeros(number_of_out,depth_of_trellis);
f=zeros(1,depth_of_trellis);  
P=0;
num=0;
for i=1:depth_of_trellis
    f(i)=output(P+1,source(i)+1);
    for j=0:number_of_out-1
        if(j~=f(i))
           rc=sgma*randn;
           rs=sgma*randn;
       else
           rc=raylrnd(1/sqrt(2))+sgma*randn;
           rs=sgma*randn;
       end
       demod_input(j+1,i)=sqrt(rc^2+rs^2);
   end
   P=nextstate(P+1,source(i)+1);
   [C,I]=max(demod_input(:,i));
   deci(i)=I-1;
   if(f(i)~=deci(i))
       num=num+1;
   end
end
ps=num/N;
t=fanout/(2*(fanout-1));
pb=t*ps;