function pb=mlsd(rho_in_dB,dfhtrel,trel,K)
% this function simulates the viterbi hard dicision demodulation 
% of the DFH system with convolutional code

% information source
rho=10^(rho_in_dB/10);
seg=1;
len=1000;
source=[randint(seg,len-K),zeros(seg,K)];

% outer encoder 
code1=zeros(seg,2*len);
for i=1:seg
    code1(i,:)=convenc(source(i,:),trel,[]);
end

% interleave the output of the outer encoder
inte_in=zeros(seg,len);
inte_out=zeros(seg,len); 
for i=1:seg
    for j=1:len
        inte_in(i,j)=change2deci(code1(i,2*j-1:2*j),2);
    end
end
posit=zeros(size(inte_out));
for i=1:seg
    [inte_out(i,:),posit(i,:)]=interleave(inte_in(i,:));
end

% generate transmitted frequency sequence and inner encoder output
nextstate=dfhtrel.nextstates;
fanout=dfhtrel.numinput;
number_of_out=dfhtrel.numoutput;
number_of_states=dfhtrel.numstates;
input=dfhtrel.inputs;
output=dfhtrel.outputs;
E=1;
sgma=sqrt(E/(2*rho));
demod_input=zeros(number_of_out,len,seg);
for i=1:seg
    f=zeros(1,len);  
    P=0;
    for j=1:len
        f(j)=output(P+1,inte_out(i,j)+1);
        for k=0:number_of_out-1
            if(k~=f(j))
                rc=sgma*randn;
                rs=sgma*randn;
            else
                %rc=sqrt(E)+sgma*randn;
                rc=raylrnd(1/sqrt(2))+sgma*randn;
                rs=sgma*randn;
            end
            demod_input(k+1,j,i)=sqrt(rc^2+rs^2);
        end
        P=nextstate(P+1,inte_out(i,j)+1);
    end
end
%demod_input=demod_input./sgma^2;
sgma21=2*sgma^2*(1+rho);
sgma22=2*sgma^2;

% demodulation and decode
num_of_err=0;
ncode=zeros(size(code1));
for i=1:seg
    inner_output=vitdemod(demod_input(:,:,i),nextstate,fanout,number_of_states,input,output,sgma21,sgma22);
    outer_input=deinterleave(inner_output,posit(i,:));
    for j=1:len
        ncode(i,j*2-1:j*2)=deci2change(outer_input(j),2,2);
    end
    outer_out=vitdec(ncode,trel,50,'trunc','hard');
    for j=1:len
        if(source(i,j)~=outer_out(j))
            num_of_err=num_of_err+1;
        end
    end
end
pb=num_of_err/(seg*len);