rho_in_dB=7;
rho=10^(rho_in_dB/10);
N=1000;
BPH=1;
number_of_states=8;
fanout=2^BPH;
L=floor(log(number_of_states)/log(fanout));
%source=[rem(unidrnd(BPH,1,N),BPH),zeros(1,BPH*L)];
%source1=reshape(source,BPH,length(source)/BPH);
%dsource=zeros(1,size(source1,2));
%dsource=source1(1,:)*2+source1(2,:);
dsource=[unidrnd(BPH,1,N),zeros(1,BPH*L)];
depth_of_trellis=length(dsource);
f=zeros(1,depth_of_trellis);
D=0;
for i=1:depth_of_trellis
    f(i)=G_func(D,dsource(i),L,fanout);
    D=f(i);
end
Eb=1;
sgma=sqrt(Eb/(2*rho));
demod_input=zeros(number_of_states,depth_of_trellis);
for i=1:depth_of_trellis
    for j=0:number_of_states-1
        if(j~=f(i))
           rc(j+1)=sgma*randn;
           rs(j+1)=sgma*randn;
       else
           rc(j+1)=sqrt(Eb)+sgma*randn;
           rs(j+1)=sgma*randn;
       end
       demod_input(j+1,i)=rc(j+1)^2+rs(j+1)^2;
   end
end
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;
        %output(i+1,j+1)=next_state;
        input(i+1,next_state+1)=j;
    end
end 
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(G_func(j,m,L,fanout)+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(G_func(j,0,L,fanout)+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_matrix=zeros(BPH,depth_of_trellis-L);
decoder_output=zeros(BPH,depth_of_trellis-L);
for i=1:depth_of_trellis-L
    %dec_output_quar=input(state_sequence(1,i)+1,state_sequence(1,i+1)+1);
    dec_output_bin=input(state_sequence(1,i)+1,state_sequence(1,i+1)+1);
    %dec_output_bin=quar2bin(dec_output_quar,BPH);
    %decoder_output_matrix(:,i)=dec_output_bin';
    decoder_output(i)=dec_output_bin;
end
%decoder_output=reshape(decoder_output_matrix,1,N);
num_of_err=0;
for i=1:N
    if(dsource(i)~=decoder_output(i))
       num_of_err=num_of_err+1;
   end
end
p=num_of_err/N;
sprintf('p=%f',p)