%function pb=osa_simu(BPH,number_of_states,D,rho_in_dB)
% this function simulates the MAP demodulation
% progress of the DFH system using optimum
% soft-output algorithm
% BPH:the number of bits per hop
% number_of_states:the number of frequency slots
% n:correlation interval
% rho_in_dB:SNR in dB
BPH=1;
number_of_states=8;
D=3;
rho_in_dB=10;
rho=10^(rho_in_dB/10);
N=1000;
fanout=2^BPH;
L=floor(log(number_of_states)/log(fanout));
source=randint(1,BPH*N);
dsource1=zeros(1,N);
if(BPH~=1)
    source1=reshape(source,BPH,N);
    for i=1:N
         dsource1(i)=change2deci(source1(:,i)',2);
    end
else
    dsource1=source;
end
dsource=zeros(1,N+D);
dsource=[dsource1,randint(1,D,fanout)]; %generate info source
depth_of_trellis=length(dsource);

% derive the state transfer matrix and the former state matrix
nextstate=zeros(number_of_states,fanout); 
formerstate=zeros(number_of_states,fanout);
input=zeros(number_of_states,number_of_states);
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;
        former_state=inv_G_func(i,j,L,fanout);
        formerstate(i+1,j+1)=former_state;
    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,dsource(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)=rc^2+rs^2;
   end
end

% start OSA demodulation
prob_xz=[1,zeros(1,number_of_states-1)];
alpha=zeros(number_of_states,fanout);
beta=zeros(number_of_states,fanout);
if(D>L)
   S_x=zeros(D-L,fanout,number_of_states); % state survivor
   S_x(:,:,1)=[(1/fanout)*ones(1,fanout);zeros(D-L-1,fanout)];
   S1_x=zeros(D-L,fanout,number_of_states); % matrix for computing the decis_u
end
decis_u=zeros(1,fanout); % soft outputs of each input symbol
decis=zeros(1,depth_of_trellis*BPH);  % decision results
for i=1:depth_of_trellis
    %if D=L
    if(D==L)
        % 1)calculate alpha for all branch
        for j=1:number_of_states
            for k=1:fanout
                %mm=demod_input(nextstate(j,k)+1,i+1)+demod_input(j,i);
                mm=demod_input(nextstate(j,k)+1,i+1);
                alpha(j,k)=prob_xz(j)*mm;
            end
        end
        alpha=alpha./sum(sum(alpha,1),2);
        % 2)for each state Xk+1,do
        for j=1:number_of_states
            % step a),calculate the innovation of prob_xz
            temp1=0;
            for k=1:fanout
                temp1=temp1+alpha(formerstate(j,k)+1,input_data(j,L,fanout)+1);
            end
            prob_xz(j)=temp1;
            % step b),calculate beta for fanout states of Xk
            if(prob_xz(j)==0)
            beta(j,:)=zeros(1,fanout);
            else
                for k=1:fanout
                    beta(j,k)=alpha(formerstate(j,k)+1,input_data(j,L,fanout)+1)/prob_xz(j);
                end
            end
        end
        % 3)calculate the information packet
        for k=1:fanout
            temp3=0;
            for j=1:number_of_states
                temp3=temp3+beta(j,k)*prob_xz(j);
            end
            decis_u(k)=temp3;
        end
        [C,I]=max(decis_u);
        if(BPH~=1)
           decis((i-1)*BPH+1:i*BPH)=deci2change(I-1,BPH,2);
        else
           decis(i)=I-1;
       end
    % if D>L
    else
        % 1)calculate alpha for all branch
        for j=1:number_of_states
            for k=1:fanout
                 mm=demod_input(nextstate(j,k)+1,i+1);
                alpha(j,k)=prob_xz(j)*mm;
            end
        end
        % 2)for each state Xk+1,do 
        for j=1:number_of_states
            % step a),calculate the innovation of prob_xz
            temp1=0;
            for k=1:fanout
                temp1=temp1+alpha(formerstate(j,k)+1,input_data(j,L,fanout)+1);
            end
            prob_xz(j)=temp1;
            % step b),calculate beta for fanout states of Xk
            if(prob_xz(j)==0)
               beta(j,:)=zeros(1,fanout);
            else
                for k=1:fanout
                    beta(j,k)=alpha(formerstate(j,k)+1,input_data(j,L,fanout)+1)/prob_xz(j);
                end
            end
            % step c),calculate S(xk+1) from fanout state survivors S(xk)            temp2=zeros(D-L,fanout);
            for k=1:fanout
                temp2=temp2+(beta(j,k).*S_x(:,:,formerstate(j,k)+1));
            end
            S1_x(:,:,j)=temp2;
        end
        % 3)calculate the information packet from the last row of all S1_x
        for k=1:fanout
            temp3=0;
            for j=1:number_of_states
                temp3=temp3+(S1_x(D-L,k,j)*prob_xz(j));
            end
            decis_u(k)=temp3;
        end
        [C,I]=max(decis_u,[],2);
        if(BPH~=1)
            decis((i-1)*BPH+1:i*BPH)=deci2change(I-1,BPH,2);
        else
            decis(i)=I-1;
        end
        % 4)shift the contents of all S1_x by one row
        if(D-L==1)
            for j=1:number_of_states
                S_x(1,:,j)=beta(j,:);
            end
        else
            S_x(2:D-L,:,:)=S1_x(1:D-L-1,:,:);
            for j=1:number_of_states
                S_x(1,:,j)=beta(j,:);
            end
        end
    end
end

% bit error rate calculation
num_of_err=0;
for i=1:length(source)
    if(source(i)~=decis(i+D*BPH))
       num_of_err=num_of_err+1;
   end
end
pb=num_of_err/length(source); % bit error probability
sprintf('pb=%f',pb)