%function pb=log_map_simu(BPH,number_of_states,D,rho_in_dB)

BPH=1;
number_of_states=8;
rho_in_dB=7;
rho=10^(rho_in_dB/10);
N=10;
fanout=2^BPH;
L=floor(log(number_of_states)/log(fanout));
source=[randint(1,BPH*N),zeros(1,L)];
depth_of_trellis=length(source);

% derive the state transfer matrix and the former state matrix
nextstate=zeros(number_of_states,fanout); 
formerstate=number_of_states.*ones(number_of_states,fanout,fanout);
for i=0:number_of_states-1
    for j=0:fanout-1
        next_state=G_func(i,j,L,fanout);
        nextstate(i+1,j+1)=next_state;
       if(formerstate(next_state+1,1,j+1)==number_of_states)
           formerstate(next_state+1,1,j+1)=i; %input is zero
       else
           formerstate(next_state+1,2,j+1)=i; %input is one
       end
    end
end 

%G-function generates frequency sequence
f=zeros(1,depth_of_trellis);  
P=0;
for i=1:depth_of_trellis
    f(i)=nextstate(P+1,source(i)+1);
    P=f(i);
end

%simulate the FFT output
E=1;
sgma=sqrt(E/(BPH*2*rho));
demod_input=zeros(number_of_states,depth_of_trellis+1);
demod_input(:,1)=[1;zeros(number_of_states-1,1)];
for i=1:depth_of_trellis
    for j=0:number_of_states-1
        if(j~=f(i))
           rc=sgma*randn;
           rs=sgma*randn;
       else
           rc=sqrt(E)+sgma*randn;
           rs=sgma*randn;
       end
       demod_input(j+1,i+1)=sqrt(rc^2+rs^2);
   end
end
demod_input=demod_input/sgma^2;

% start max-log-map demodulation
alpha=zeros(number_of_states,depth_of_trellis);
alpha(:,1)=[0;-1e10*ones(number_of_states-1,1)];
gamma=zeros(2*fanout,depth_of_trellis);
max=-1e10*ones(1,depth_of_trellis-1);
beta=zeros(number_of_states,depth_of_trellis);
beta(:,depth_of_trellis)=[0;-1e10*ones(number_of_states-1,1)];
lu=zeros(1,depth_of_trellis); % decision variable
decis=zeros(1,depth_of_trellis);

%trace forward to compute alpha
for i=1:depth_of_trellis-1
    for j=1:number_of_states
       if(formerstate(j,1,1)~=number_of_states)
           gamma(1,i)=demod_input(formerstate(j,1,1)+1,i)+demod_input(j,i+1);
           gamma(2,i)=demod_input(formerstate(j,2,1)+1,i)+demod_input(j,i+1);
           alpha0=alpha(formerstate(j,1,1)+1,i)+gamma(1,i);
           alpha1=alpha(formerstate(j,2,1)+1,i)+gamma(2,i);
       else
           gamma(1,i)=demod_input(formerstate(j,1,2)+1,i)+demod_input(j,i+1);
           gamma(2,i)=demod_input(formerstate(j,2,2)+1,i)+demod_input(j,i+1);
           alpha0=alpha(formerstate(j,1,2)+1,i)+gamma(1,i);
           alpha1=alpha(formerstate(j,2,2)+1,i)+gamma(2,i);
       end
       if(alpha0<=-80)
           alpha0=0;
       else
           alpha0=exp(alpha0);
       end
       if(alpha1<=-80)
           alpha1=0;
       else
           alpha1=exp(alpha1);
       end
       if((alpha0+alpha1)>1e-30)
           alpha(j,i+1)=log(alpha0+alpha1);
       else
           alpha(j,i+1)=-1e10;
       end
       if(max(i)<alpha(j,i+1))
           max(i)=alpha(j,i+1);
       end
   end
   alpha(:,i+1)=alpha(:,i+1)-max(i);
end
for j=1:number_of_states
    if(formerstate(j,1,1)~=number_of_states)
        gamma(1,depth_of_trellis)=demod_input(formerstate(j,1,1)+1,depth_of_trellis)...
            +demod_input(j,depth_of_trellis+1);
        gamma(2,depth_of_trellis)=demod_input(formerstate(j,2,1)+1,depth_of_trellis)...
            +demod_input(j,depth_of_trellis+1);
    else
        gamma(1,depth_of_trellis)=demod_input(formerstate(j,1,2)+1,depth_of_trellis)...
            +demod_input(j,depth_of_trellis+1);
        gamma(2,depth_of_trellis)=demod_input(formerstate(j,2,2)+1,depth_of_trellis)...
            +demod_input(j,depth_of_trellis+1);
    end
end

%trace backward to compute beta and do error counting
num_of_err=0;
for i=depth_of_trellis:-1:2
    temp1=0;
    temp0=0;
    for j=1:number_of_states
        gamma(3,i)=demod_input(nextstate(j,1)+1,i+1)+demod_input(j,i);
        gamma(4,i)=demod_input(nextstate(j,2)+1,i+1)+demod_input(j,i);
        beta0=beta(nextstate(j,1)+1,i)+gamma(3,i);
        beta1=beta(nextstate(j,2)+1,i)+gamma(4,i);
        if(beta0<=-80)
            beta0=0;
        else
            beta0=exp(beta0);
        end
        if(beta1<=-80)
            beta1=0;
        else
            beta1=exp(beta1);
        end
        if((beta0+beta1)>1e-30)
            beta(j,i-1)=log(beta0+beta1)-max(i-1);
        else
            beta(j,i-1)=-1e10;
        end
        temp0=temp0+exp(alpha(j,i)+gamma(3,i)+beta(nextstate(j,1)+1,i));
        temp1=temp1+exp(alpha(j,i)+gamma(4,i)+beta(nextstate(j,2)+1,i));
    end
    
    %decision
    if(temp0==0)
        lu(i)=100;
    elseif(temp1==0)
        lu(i)=-100;
    else
        lu(i)=log(temp1/temp0);
    end
    if(lu(i)>80)
        lu(i)=80;
    elseif(lu(i)<-80)
        lu(i)=-80;
    end
    if(lu(i)>=0)
        decis(i)=1;
    else
        decis(i)=0;
    end
    if(source(i)~=decis(i))
        num_of_err=num_of_err+1;
    end
end
temp0=0;
temp1=0;
for j=1:number_of_states
    gamma(3,1)=demod_input(nextstate(j,1)+1,2)+demod_input(j,1);
    gamma(4,1)=demod_input(nextstate(j,2)+1,2)+demod_input(j,1);
    temp0=temp0+exp(alpha(j,1)+gamma(3,1)+beta(nextstate(j,1)+1,1));
    temp1=temp1+exp(alpha(j,1)+gamma(4,1)+beta(nextstate(j,2)+1,1));
end
if(temp0==0)
    lu(1)=100;
elseif(temp1==0)
    lu(1)=-100;
else
    lu(1)=log(temp1/temp0);
end
if(lu(1)>80)
    lu(1)=80;
elseif(lu(1)<-80)
    lu(1)=-80;
end
if(lu(1)>=0)
    decis(1)=1;
else
    decis(1)=0;
end
if(source(1)~=decis(1))
    num_of_err=num_of_err+1;
end
pb=num_of_err/depth_of_trellis;
sprintf('pb=%f',pb)