function pe=map_simu(BPH,number_of_states,n,rho_in_dB)
% this function simulates the MAP demodulation
% progress of the DFH system
% n:correlation interval
%rho_in_dB=7;
rho=10^(rho_in_dB/10);
N=10000;
%BPH=1;
%number_of_states=8;
fanout=2^BPH;
L=floor(log(number_of_states)/log(fanout));
%n=5;
dsource=[randint(1,N,[0,fanout-1]),zeros(1,BPH*(n-1))];
depth_of_trellis=length(dsource);
f=zeros(1,depth_of_trellis);
D=0;
for i=1:depth_of_trellis
    f(i)=G_func1(D,dsource(i),L);
    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=sgma*randn;
           rs=sgma*randn;
       else
           rc=sqrt(Eb)+sgma*randn;
           rs=sgma*randn;
       end
       demod_input(j+1,i)=rc^2+rs^2;
   end
end
for i=0:number_of_states-1
    for j=0:fanout-1
        next_state=G_func1(i,j,L);
        nextstate(i+1,j+1)=next_state;
        %input(i+1,next_state+1)=j;
    end
end 
state_metric=zeros(fanout^(n),1);
decis=zeros(1,N);
symbol_matrix=zeros(fanout^(n),n);
for i=1:fanout^(n)
    symbol_matrix(i,:)=deci2bin(i-1,n);
end
latest_state=0;
mid_state=zeros(fanout^(n),n);
for i=1:depth_of_trellis-n+1
    for j=1:fanout^(n)   %求出第一次中間狀態
        middlestate=latest_state;
        for k=1:n
            mid_state(j,k)=nextstate(middlestate+1,symbol_matrix(j,k)+1);
            middlestate=mid_state(j,k);
        end
    end
    for j=1:fanout^(n)
        for k=i:i+n-1
            state_metric(j)=state_metric(j)+demod_input(mid_state(j,k-i+1)+1,k);
        end
    end
    [C,I]=max((sum(reshape(state_metric,fanout^(n-1),fanout)))');
    latest_state=nextstate(latest_state+1,I);
    decis(i)=I-1;
    state_metric=zeros(fanout^(n),1);
end
num_of_err=0;
for i=1:N
    if(dsource(i)~=decis(i))
       num_of_err=num_of_err+1;
   end
end
pe=num_of_err/N;
%sprintf('p=%f',p)