rho_in_dB=5;
rho=10^(rho_in_dB/10);
N=5000;
BPH=2;
number_of_states=64;
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,:);
depth_of_trellis=length(dsource);
ts=0.2/1024;
t=[0:ts:ts*(1024-1)];
f0=300;
df=10;
fset=[f0:df:f0+df*(number_of_states-1)];
f=zeros(1,depth_of_trellis);
D=0;
for i=1:depth_of_trellis
    interval=G_func(D,dsource(i),L);
    f(i)=fset(interval+1);
    D=interval;
end
pulse=cos(2*pi*t'*f);
Eb=norm(pulse(:,1))^2*ts/BPH;
sgma=sqrt(Eb/(2*rho));

demod_pro=ts*fft(data,1024);
channel_output=demod_pro+sgma
Pdemod_pro=demod_pro.*conj(demod_pro)/1024;
demod_input=Pdemod_pro(1:2:1023,:);
for i=0:number_of_states-1
    for j=0:fanout-1
        next_state=G_func(i,j,L);
        nextstate(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(fset(G_func(j,m,L)+1)/df+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(fset(G_func(j,0,L)+1)/df+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);
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=quar2bin(dec_output_quar,BPH);
    decoder_output_matrix(:,i)=dec_output_bin';
end
decoder_output=reshape(decoder_output_matrix,1,N);
num_of_err=0;
for i=1:N
    if(source(i)~=decoder_output(i))
       num_of_err=num_of_err+1;
   end
end
p=num_of_err/N;
sprintf('p=%f',p)