log_map_simu4.asv

短波信道抗多音干扰的性能分析及其仿真

%function pb=log_map_simu4(rho_in_dB,BPH,number_of_states)

BPH=1;
number_of_states=8;
rho_in_dB=5;
N=1000;
fanout=2^BPH;
rho=10^(rho_in_dB/10);
L=floor(log(number_of_states)/log(fanout));
source=[randint(1,N,fanout),zeros(1,L)];
dsource=zeros(1,N*BPH);
if(BPH~=1)
    for i=1:N
        dsource((i-1)*BPH+1:i*BPH)=deci2change(source(i),BPH,2);
    end
else
    dsource=source(1:N);
end
depth_of_trellis=length(source);

nextstate=zeros(number_of_states,fanout); 
formerstate=number_of_states.*ones(number_of_states,fanout,fanout);
for i=0:number_of_states-1
    for j=0:fanout-1
        k=1;
        next_state=G_func(i,j,L,fanout);
        nextstate(i+1,j+1)=next_state;
        while(formerstate(next_state+1,k,j+1)~=number_of_states)
            k=k+1;
        end
        formerstate(next_state+1,k,j+1)=i; 
    end
end

f=zeros(1,depth_of_trellis);  
P=0;
for i=1:depth_of_trellis
    f(i)=nextstate(P+1,source(i)+1);
    P=f(i);
end

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)=sqrt(rc^2+rs^2);
   end
end
demod_input=demod_input/sgma^2;

alpha=zeros(number_of_states,depth_of_trellis);
alpha(:,1)=[0;-1e10*ones(number_of_states-1,1)];
alpha1=zeros(1,fanout);
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)];
beta1=zeros(1,fanout);
lu=zeros(depth_of_trellis,fanout); % decision variable
decis=zeros(1,depth_of_trellis);
decision=zeros(1,BPH*N);

for i=1:depth_of_trellis-1
    for j=1:number_of_states
        for k=1:fanout
            if(formerstate(j,1,k)~=number_of_states)
                for h=1:fanout
                    gamma(h,i)=demod_input(formerstate(j,h,k)+1,i)+demod_input(j,i+1);
                    alpha1(h)=alpha(formerstate(j,h,k)+1,i)+gamma(h,i);
                    if(alpha1(h)<=-80)
                        alpha1(h)=0;
                    else
                        alpha1(h)=exp(alpha1(h));
                    end
                end
                if(sum(alpha1)>1e-30)
                     alpha(j,i+1)=log(sum(alpha1));
                else
                     alpha(j,i+1)=-1e10;
                end
                if(max(i)<alpha(j,i+1))
                     max(i)=alpha(j,i+1);
                end
            end
         end
     end
     alpha(:,i+1)=alpha(:,i+1)-max(i);
 end
for j=1:number_of_states
    for k=1:fanout
        if(formerstate(j,1,k)~=number_of_states)
            for h=1:fanout
                gamma(h,depth_of_trellis)=demod_input(formerstate(j,h,k)+1,depth_of_trellis)...
                    +demod_input(j,depth_of_trellis+1);
            end
        end
    end
end

for i=depth_of_trellis:-1:2
    temp=zeros(1,fanout);
    for j=1:number_of_states
        for k=1:fanout
            gamma(k+fanout,i)=demod_input(nextstate(j,k)+1,i+1)+demod_input(j,i);
            beta1(k)=beta(nextstate(j,k)+1,i)+gamma(k+fanout,i);
            if(beta1(k)<=-80)
                beta1(k)=0;
            else
                beta1(k)=exp(beta1(k));
            end
            
        end
        if(sum(beta1)>1e-30)
            beta(j,i-1)=log(sum(beta1))-max(i-1);
        else
            beta(j,i-1)=-1e10;
        end
        for h=1:fanout
            
            if(temp(h)==0)
                temp(h)=exp(-100);
            end
        end
    end
    lu(i,:)=log(temp);
    [C,I]=max(lu(i,:));
    decis(i)=I-1;
    if(lu(i,I)>80)
        lu(i,I)=80;
    elseif(lu(i,I)<-80)
        lu(i,I)=-80;
    end
end
temp=zeros(1,fanout);
for j=1:number_of_states
    for k=1:fanout
        gamma(k+fanout,1)=demod_input(nextstate(j,k)+1,2)+demod_input(j,1);
        temp(k)=temp(k)+exp(alpha(j,1)+gamma(k+fanout,1)+beta(nextstate(j,k)+1,1));
    end
end
lu(1,:)=log(temp);
[C,I]=max(lu(1,:));
decis(1)=I-1;
if(lu(1,I)>80)
    lu(1,I)=80;
elseif(lu(1,I)<-80)
    lu(1,I)=-80;
end

if(BPH~=1)
    for i=1:N
        decision((i-1)*BPH+1:i*BPH)=deci2change(decis(i),BPH,2);
    end
else
    decision=decis(1:N);
end
num_of_err=0;
for i=1:BPH*N
    if(dsource(i)~=decision(i))
        num_of_err=num_of_err+1;
    end
end
pb=num_of_err/(BPH*N);
sprintf('pb=%f',pb)