function pb=va_demod(rho_in_dB,BPH,number_of_states)
% this function simulates the viterbi demodulation of
% the DFH system

% rho_in_dB=5;
% BPH=2;
% number_of_states=16;

N=1000;
fanout=2^BPH;
rho=10^(rho_in_dB/10);
L=floor(log(number_of_states)/log(fanout));
source=[randint(1,N,fanout),zeros(1,L)];
dsource=zeros(1,N*BPH);
if(BPH~=1)
    for i=1:N
        dsource((i-1)*BPH+1:i*BPH)=deci2change(source(i),BPH,2);
    end
else
    dsource=source(1:N);
end

nextstate=zeros(number_of_states,fanout); 
output=zeros(number_of_states,fanout);
input=zeros(number_of_states,number_of_states);
number_of_out=number_of_states*fanout;
for i=0:number_of_states-1
    for j=0:fanout-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;
        input(i+1,next_state+1)=j;
    end
end 

depth_of_trellis=length(source);
E=1;
sgma=sqrt(E/(BPH*2*rho));
demod_input=zeros(number_of_out,depth_of_trellis);
f=zeros(1,depth_of_trellis);
D=0;
for i=1:depth_of_trellis
    f(i)=output(D+1,source(i)+1);
    for j=0:number_of_out-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)=sqrt(rc^2+rs^2);
   end
    D=nextstate(D+1,source(i)+1);
end
% demod_input=demod_input/sgma^2;

state_metric=zeros(number_of_states,2);
survivor_state=zeros(number_of_states,depth_of_trellis+1);
for i=1:depth_of_trellis-L
    flag=zeros(1,number_of_states);
    if i<=L+1
        step=2^((L+1-i)*BPH);
    else
        step=1;
    end
    for j=0:step:number_of_states-1
        for m=0:fanout-1
            branch_metric=demod_input(output(j+1,m+1)+1,i);
            if((state_metric(nextstate(j+1,m+1)+1,2)<state_metric(j+1,1)...
                    +branch_metric)|flag(nextstate(j+1,m+1)+1)==0)
                state_metric(nextstate(j+1,m+1)+1,2)=state_metric(j+1,1)+branch_metric;
                survivor_state(nextstate(j+1,m+1)+1,i+1)=j;
                flag(nextstate(j+1,m+1)+1)=1;
            end
        end
    end
    state_metric=state_metric(:,2:-1:1);
end
for i=depth_of_trellis-L+1:depth_of_trellis
    flag=zeros(1,number_of_states);
    last_stop=number_of_states/(2^((i-depth_of_trellis+L-1)*BPH));
    for j=0:last_stop-1
        branch_metric=demod_input(output(j+1,m+1)+1,i);
        if((state_metric(nextstate(j+1,1)+1,2)<state_metric(j+1,1)...
            +branch_metric)|flag(nextstate(j+1,1)+1)==0)
            state_metric(nextstate(j+1,1)+1,2)=state_metric(j+1,1)+branch_metric;
            survivor_state(nextstate(j+1,1)+1,i+1)=j;
            flag(nextstate(j+1,1)+1)=1;
        end
    end
    state_metric=state_metric(:,2:-1:1);
end
state_sequence=zeros(1,depth_of_trellis+1);
for i=1:depth_of_trellis
    state_sequence(1,depth_of_trellis-i+1)=survivor_state((state_sequence(1,depth_of_trellis+2-i)...
        +1),depth_of_trellis-i+2);
end
decoder_output=zeros(1,BPH*(depth_of_trellis-L));
for i=1:depth_of_trellis-L
    dec_output_deci=input(state_sequence(1,i)+1,state_sequence(1,i+1)+1);
    if(BPH~=1)
        dec_output_bin=deci2change(dec_output_deci,BPH,2);
        decoder_output((i-1)*BPH+1:i*BPH)=dec_output_bin;
    else
        decoder_output(1,i)=dec_output_deci;
    end
end

num_of_err=0;
for i=1:N*BPH
    if(dsource(i)~=decoder_output(i))
       num_of_err=num_of_err+1;
   end
end
pb=num_of_err/(N*BPH);