log_map_simu1.asv

来自「短波信道抗多音干扰的性能分析及其仿真」· ASV 代码 · 共 139 行

function inner_out=log_map_simu1(source,rho_in_dB,BPH,number_of_states)
% this function simulates the MAP demodulation of the
% DFH system

fanout=2^BPH;
rho=10^(rho_in_dB/10);
depth_of_trellis=length(source);
dsource=zeros(1,depth_of_trellis*BPH);
if(BPH~=1)
    for i=1:depth_of_trellis
        dsource((i-1)*BPH+1:i*BPH)=deci2change(source(i),BPH,2);
    end
else
    dsource=source(1:depth_of_trellis);
end

alpha=zeros(number_of_states,depth_of_trellis);
alpha(:,1)=[0;-1e10*ones(number_of_states-1,1)];
gamma=zeros(2*fanout,depth_of_trellis);
max=-1e10*ones(1,depth_of_trellis-1);
beta=zeros(number_of_states,depth_of_trellis);
beta(:,depth_of_trellis)=[0;-1e10*ones(number_of_states-1,1)];
lu=zeros(1,depth_of_trellis); % decision variable
decis=zeros(1,depth_of_trellis);

%trace forward to compute alpha
for i=1:depth_of_trellis-1
    for j=1:number_of_states
       if(formerstate(j,1,1)~=number_of_states)
           gamma(1,i)=demod_input(formerstate(j,1,1)+1,i)+demod_input(j,i+1);
           gamma(2,i)=demod_input(formerstate(j,2,1)+1,i)+demod_input(j,i+1);
           alpha0=alpha(formerstate(j,1,1)+1,i)+gamma(1,i);
           alpha1=alpha(formerstate(j,2,1)+1,i)+gamma(2,i);
       else
           gamma(1,i)=demod_input(formerstate(j,1,2)+1,i)+demod_input(j,i+1);
           gamma(2,i)=demod_input(formerstate(j,2,2)+1,i)+demod_input(j,i+1);
           alpha0=alpha(formerstate(j,1,2)+1,i)+gamma(1,i);
           alpha1=alpha(formerstate(j,2,2)+1,i)+gamma(2,i);
       end
       if(alpha0<=-80)
           alpha0=0;
       else
           alpha0=exp(alpha0);
       end
       if(alpha1<=-80)
           alpha1=0;
       else
           alpha1=exp(alpha1);
       end
       if((alpha0+alpha1)>1e-30)
           alpha(j,i+1)=log(alpha0+alpha1);
       else
           alpha(j,i+1)=-1e10;
       end
       if(max(i)<alpha(j,i+1))
           max(i)=alpha(j,i+1);
       end
   end
   alpha(:,i+1)=alpha(:,i+1)-max(i);
end
for j=1:number_of_states
    if(formerstate(j,1,1)~=number_of_states)
        gamma(1,depth_of_trellis)=demod_input(formerstate(j,1,1)+1,depth_of_trellis)...
            +demod_input(j,depth_of_trellis+1);
        gamma(2,depth_of_trellis)=demod_input(formerstate(j,2,1)+1,depth_of_trellis)...
            +demod_input(j,depth_of_trellis+1);
    else
        gamma(1,depth_of_trellis)=demod_input(formerstate(j,1,2)+1,depth_of_trellis)...
            +demod_input(j,depth_of_trellis+1);
        gamma(2,depth_of_trellis)=demod_input(formerstate(j,2,2)+1,depth_of_trellis)...
            +demod_input(j,depth_of_trellis+1);
    end
end

%trace backward to compute beta and do error counting
num_of_err=0;
for i=depth_of_trellis:-1:2
    temp1=0;
    temp0=0;
    for j=1:number_of_states
        gamma(3,i)=demod_input(nextstate(j,1)+1,i+1)+demod_input(j,i);
        gamma(4,i)=demod_input(nextstate(j,2)+1,i+1)+demod_input(j,i);
        beta0=beta(nextstate(j,1)+1,i)+gamma(3,i);
        beta1=beta(nextstate(j,2)+1,i)+gamma(4,i);
        if(beta0<=-80)
            beta0=0;
        else
            beta0=exp(beta0);
        end
        if(beta1<=-80)
            beta1=0;
        else
            beta1=exp(beta1);
        end
        if((beta0+beta1)>1e-30)
            beta(j,i-1)=log(beta0+beta1)-max(i-1);
        else
            beta(j,i-1)=-1e10;
        end
        temp0=temp0+exp(alpha(j,i)+gamma(3,i)+beta(nextstate(j,1)+1,i));
        temp1=temp1+exp(alpha(j,i)+gamma(4,i)+beta(nextstate(j,2)+1,i));
    end
    
    %decision
    if(temp0==0)
        lu(i)=100;
    elseif(temp1==0)
        lu(i)=-100;
    else
        lu(i)=log(temp1/temp0);
    end
    if(lu(i)>80)
        lu(i)=80;
    elseif(lu(i)<-80)
        lu(i)=-80;
    end
end
temp0=0;
temp1=0;
for j=1:number_of_states
    gamma(3,1)=demod_input(nextstate(j,1)+1,2)+demod_input(j,1);
    gamma(4,1)=demod_input(nextstate(j,2)+1,2)+demod_input(j,1);
    temp0=temp0+exp(alpha(j,1)+gamma(3,1)+beta(nextstate(j,1)+1,1));
    temp1=temp1+exp(alpha(j,1)+gamma(4,1)+beta(nextstate(j,2)+1,1));
end
if(temp0==0)
    lu(1)=100;
elseif(temp1==0)
    lu(1)=-100;
else
    lu(1)=log(temp1/temp0);
end
if(lu(1)>80)
    lu(1)=80;
elseif(lu(1)<-80)
    lu(1)=-80;
end