outer_decod_c.asv

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

function outer_out_c=outer_decod_c(outer_input,trellis)
% this function does the MAP decode process of the
% outer convolutional code

number_of_states=trellis.numStates;
nextstate=trellis.nextStates;
fanout=trellis.numInputSymbols;
output=trellis.outputs;
number_of_out=trellis.numOutputSymbols;
r=log(number_of_out)/log(2);
depth_of_trellis=length(outer_input)/r;
formerstate=number_of_states.*ones(number_of_states,fanout,fanout);
outseq=zeros(number_of_states,r*fanout);
for i=0:number_of_states-1
    for j=0:fanout-1
        y=deci2change(output(i+1,j+1),r,2);
        outseq(i+1,r*j+1:r*(j+1))=y;
        k=1;
        while(formerstate(nextstate(i+1,j+1)+1,k,j+1)~=number_of_states)
            k=k+1;
        end
        formerstate(nextstate(i+1,j+1)+1,k,j+1)=i; 
    end
end

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)];
lu=zeros(size(outer_input)); % decision variable

%trace forward to compute alpha
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
                    t=outseq(formerstate(j,h,k)+1,r*(k-1)+1:r*k);
                    gamma(h,i)=outer_input((i-1)*r+1:i*r)*(2*t-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
                t=outseq(formerstate(j,h,k)+1,r*(k-1)+1:r*k);
                gamma(h,depth_of_trellis)=outer_input((depth_of_trellis-1)*r+1...
                :depth_of_trellis*r)*(2*t-1)';
            end
        end
    end
end

for i=depth_of_trellis:-1:2
    temp=zeros(2,r);
    for j=1:number_of_states
        for k=1:fanout
            t=outseq(nextstate(j,k)+1,r*(k-1)+1:r*k);
            gamma(k+fanout,i)=outer_input((i-1)*r+1:i*r)*(2*t-1)';
            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 p=1:r
            for h=1:fanout
                if(outseq(j,r*(h-1)+p)==0)
                    temp(1,p)=temp(1,p)+exp(alpha(j,i)+gamma(h+fanout,i)+beta(nextstate(j,h)+1,i));
                else
                    temp(2,p)=temp(2,p)+exp(alpha(j,i)+gamma(h+fanout,i)+beta(nextstate(j,h)+1,i));
                end
            end
        end
    end
    for p=1:r
        if(temp(1,p)==0)
            lu(r*(i-1)+p)=100;
        elseif(temp(2,p)==0)
            lu(r*(i-1)+p)=-100;
        else
            lu(r*(i-1)+p)=log(temp(2,p)/temp(1,p));
        end
        if(lu(r*(i-1)+p)>80)
            lu(r*(i-1)+p)=80;
        elseif(lu(r*(i-1)+p)<-80)
            lu(r*(i-1)+p)=-80;
        end
    end
end
temp=zeros(2,r);
for j=1:number_of_states
    for k=1:fanout
        t=outseq(nextstate(j,k)+1,r*(k-1)+1:r*k);
        gamma(k+fanout,1)=outer_input((i-1)*r+1:i*r)*(2*t-1)';
    end
    for p=1:r
        for h=1:fanout
            if(outseq(j,r*(h-1)+p)==0)
                temp(1,p)=temp(1,p)+exp(alpha(j,1)+gamma(h+fanout,1)+beta(nextstate(j,h)+1,1));
            else
                temp(2,p)=temp(2,p)+exp(alpha(j,1)+gamma(h+fanout,1)+beta(nextstate(j,h)+1,1));
            end
        end
    end
end
for p=1:r
    if(temp(1,p)==0)
        lu(p)=100;
    elseif(temp(2,p)==0)
        lu(p)=-100;
    else
        lu(p)=log(temp(2,p)/temp(1,p));
    end
    if(lu(p)>80)
        lu(p)=80;
    elseif(lu(p)<-80)
        lu(p)=-80;
    end
end
outer_out_c=lu-outer_input;