logmap.m

传统的卷积码和Turbo码直到今天仍有广泛的应用

function [L_all,start_state]=logmap(rec,state_map,L_a_1row,iter,start_state,L1)
%rec中包含A,B,parity1三部分,等长,因此网格图中经历的状态转移次数L_total=rec长度的1/3
L_total = length(rec)/3;
nstates = 8;          % number of states in the trellis
L_all=zeros(1,2*L_total);
% Set up the trellis
state_map = trellis;
%无穷大
Infty = 1e10;
%前项概率alpha,三维分别对应4种输入,网格图中L_total+1个时刻,状态总数
Alpha=zeros(4,L_total+1,nstates);
% Initialization of Alpha
%前L1次等概初始化
%L1+1次迭代开始,用L1次迭代估计的编码器初始状态start_state来初始化
if iter<=L1
    Alpha(:,1,:) = -log(nstates); 
else
    Alpha(:,1,:)=-Infty;
    Alpha(:,1,start_state+1)=0;
end
%后项概率beta,三维分别对应4种输入,网格图中L_total+1个时刻,状态总数
Beta=zeros(4,L_total,nstates);
% Initialization of Beta
if iter<=L1
    Beta(:,L_total,:) = -log(nstates); 
else
    Beta(:,L_total,:)=-Infty;
    Beta(:,L_total,start_state+1)=0;
end
% L_a is seperated into two rows:各代表A,B的边信息
L_a=zeros(2,L_total);
L_a(1,:)=L_a_1row(1:2:end);
L_a(2,:)=L_a_1row(2:2:end);
tempmax=zeros(1,L_total);
% Trace forward, compute Alpha
%计算状态转移概率gamma,对4种输入分别计算,将平方计算简化为乘积运算
%rec(3*k-5),rec(3*k-4),rec(3*k-3)对应3比特组,前两位是信息,末位是校验.
%gamma(1,i)=rec(3*k-5)*(-1)+rec(3*k-4)*(-1)+rec(3*k-3)*output(00,i)+log(1/(1+exp(L_a(1,k-1))))+log(1/(1+exp(L_a(2,k-1))))
%output(00,i)是第i个状态下输入为00时的输出,通过查表state_map得到
for k = 2:L_total+1
    gamma = -Infty*ones(4,nstates);
    for state2 = 1:nstates
        gamma(1,state2) = (-rec(3*k-5)-rec(3*k-4)+rec(3*k-3)*state_map(state2,11))...
            -log(1+exp(L_a(1,k-1)))-log(1+exp(L_a(2,k-1)));
        gamma(2,state2) = (-rec(3*k-5)+rec(3*k-4)+rec(3*k-3)*state_map(state2,14))...
            -log(1+exp(L_a(1,k-1)))+L_a(2,k-1)-log(1+exp(L_a(2,k-1)));
        gamma(3,state2) = (rec(3*k-5)-rec(3*k-4)+rec(3*k-3)*state_map(state2,17))...
            +L_a(1,k-1)-log(1+exp(L_a(1,k-1)))-log(1+exp(L_a(2,k-1)));
        gamma(4,state2) = (rec(3*k-5)+rec(3*k-4)+rec(3*k-3)*state_map(state2,20))...
            +L_a(1,k-1)-log(1+exp(L_a(1,k-1)))+L_a(2,k-1)-log(1+exp(L_a(2,k-1)));
    end
%计算Alpha:Alpha(i,k,j)=sum(exp(Alpha(m,k-1,nextstate(i,j))+gamma(i,j)))--对m求和,m=0,1,2,3
%nextstate(i,j)是j状态下输入为i时转移到的状态
    for i=1:4
        for j=1:nstates
            Alpha(i,k,j)=exp(Alpha(1,k-1,state_map(j,i))+gamma(i,j))+exp(Alpha(2,k-1,state_map(j,i))+gamma(i,j))...
                +exp(Alpha(3,k-1,state_map(j,i))+gamma(i,j))+exp(Alpha(4,k-1,state_map(j,i))+gamma(i,j));
            if(Alpha(i,k,j)<1e-300)
                Alpha(i,k,j)=-Infty;
            else
                Alpha(i,k,j) = log(Alpha(i,k,j));  
            end  
        end
    end
%归一化
    tempmax(k)=Alpha(1,k,1);
    temp(1) = max(Alpha(1,k,:));
    temp(2) = max(Alpha(2,k,:));
    temp(3) = max(Alpha(3,k,:));
    temp(4) = max(Alpha(4,k,:));
    temp=sort(temp);
    tempmax(k) = temp(4);
    Alpha(1,k,:) = Alpha(1,k,:) - tempmax(k);
    Alpha(2,k,:) = Alpha(2,k,:) - tempmax(k);
    Alpha(3,k,:) = Alpha(3,k,:) - tempmax(k);
    Alpha(4,k,:) = Alpha(4,k,:) - tempmax(k);
end     

% Trace backward, compute Beta
%计算状态转移概率gamma,对4种输入分别计算,将平方计算简化为乘积运算
%rec(3*k+1),rec(3*k+2),rec(3*k+3)对应3比特组,前两位是信息,末位是校验.
%gamma(1,i)=rec(3*k+1)*(-1)+rec(3*k+2)*(-1)+rec(3*k+3)*output(00,i)+log(1/(1+exp(L_a(1,k-1))))+log(1/(1+exp(L_a(2,k-1))))
%output(00,i)是第i个状态下输入为00时的输出,通过查表state_map得到
for k = L_total-1:-1:1
    gamma = -Infty*ones(4,nstates);
    for state1 = 1:nstates
        gamma(1,state1) = (-rec(3*k+1)-rec(3*k+2)+rec(3*k+3)*state_map(state1,23))....
            -log(1+exp(L_a(1,k+1)))-log(1+exp(L_a(2,k+1)));
        gamma(2,state1) = (-rec(3*k+1)+rec(3*k+2)+rec(3*k+3)*state_map(state1,26))....
            -log(1+exp(L_a(1,k+1)))+L_a(2,k+1)-log(1+exp(L_a(2,k+1)));
        gamma(3,state1) = (rec(3*k+1)-rec(3*k+2)+rec(3*k+3)*state_map(state1,29))....
            +L_a(1,k+1)-log(1+exp(L_a(1,k+1)))-log(1+exp(L_a(2,k+1)));
        gamma(4,state1) = (rec(3*k+1)+rec(3*k+2)+rec(3*k+3)*state_map(state1,32))....
            +L_a(1,k+1)-log(1+exp(L_a(1,k+1)))+L_a(2,k+1)-log(1+exp(L_a(2,k+1)));
    end
%计算Beta:Beta(i,k,j)=sum(exp(Beta(m,k+1,laststate(i,j))+gamma(i,j)))--对m求和,m=0,1,2,3
%laststate(i,j)是输入为i时转移到j状态的那个状态
    for i=1:4
        for j=1:nstates
            Beta(i,k,j)=exp(Beta(1,k+1,state_map(j,i+4))+gamma(i,j))+exp(Beta(2,k+1,state_map(j,i+4))+gamma(i,j))...
                +exp(Beta(3,k+1,state_map(j,i+4))+gamma(i,j))+exp(Beta(4,k+1,state_map(j,i+4))+gamma(i,j));
            if(Beta(i,k,j)<1e-300)
                Beta(i,k,j)=-Infty;
            else
                Beta(i,k,j) = log(Beta(i,k,j));  
            end  
        end
    end
%归一化
    Beta(1,k,:) = Beta(1,k,:) - tempmax(k+1);
    Beta(2,k,:) = Beta(2,k,:) - tempmax(k+1);
    Beta(3,k,:) = Beta(3,k,:) - tempmax(k+1);
    Beta(4,k,:) = Beta(4,k,:) - tempmax(k+1);
end

% Compute the soft output, log-likelihood ratio of symbols in the frame
%//delta(1,i)~delta(4,i):分别代表第i个时刻下信息为00,01,10,11的概率
for k = 1:L_total
    for state2 = 1:nstates
        gamma0 = (-rec(3*k-2)-rec(3*k-1)+rec(3*k)*state_map(state2,11))....
            -log(1+exp(L_a(1,k)))-log(1+exp(L_a(2,k)));
        gamma1 = (-rec(3*k-2)+rec(3*k-1)+rec(3*k)*state_map(state2,14))....
            -log(1+exp(L_a(1,k)))+L_a(2,k)-log(1+exp(L_a(2,k)));
        gamma2 = (rec(3*k-2)-rec(3*k-1)+rec(3*k)*state_map(state2,17))....
            +L_a(1,k)-log(1+exp(L_a(1,k)))-log(1+exp(L_a(2,k)));
        gamma3 = (rec(3*k-2)+rec(3*k-1)+rec(3*k)*state_map(state2,20))....
            +L_a(1,k)-log(1+exp(L_a(1,k)))+L_a(2,k)-log(1+exp(L_a(2,k)));  
%temp=delta*Alpha(k-1)*Beta(k)   
        temp0(state2) = sum(exp(gamma0 + Alpha(:,k,state_map(state2,1)) + Beta(1,k,state2)))...
            +sum(exp(gamma0 + Alpha(:,k,state_map(state2,1)) + Beta(2,k,state2)))...
            +sum(exp(gamma0 + Alpha(:,k,state_map(state2,1)) + Beta(3,k,state2)))...
            +sum(exp(gamma0 + Alpha(:,k,state_map(state2,1)) + Beta(4,k,state2)));
        temp1(state2) = sum(exp(gamma1 + Alpha(:,k,state_map(state2,2)) + Beta(2,k,state2)))...
            +sum(exp(gamma1 + Alpha(:,k,state_map(state2,2)) + Beta(1,k,state2)))...
            +sum(exp(gamma1 + Alpha(:,k,state_map(state2,2)) + Beta(3,k,state2)))...
            +sum(exp(gamma1 + Alpha(:,k,state_map(state2,2)) + Beta(4,k,state2)));
        temp2(state2) = sum(exp(gamma2 + Alpha(:,k,state_map(state2,3)) + Beta(3,k,state2)))...
            +sum(exp(gamma2 + Alpha(:,k,state_map(state2,3)) + Beta(1,k,state2)))...
            +sum(exp(gamma2 + Alpha(:,k,state_map(state2,3)) + Beta(2,k,state2)))...
            +sum(exp(gamma2 + Alpha(:,k,state_map(state2,3)) + Beta(4,k,state2)));
        temp3(state2) = sum(exp(gamma3 + Alpha(:,k,state_map(state2,4)) + Beta(4,k,state2)))...
            +sum(exp(gamma3 + Alpha(:,k,state_map(state2,4)) + Beta(1,k,state2)))...
            +sum(exp(gamma3 + Alpha(:,k,state_map(state2,4)) + Beta(2,k,state2)))...
            +sum(exp(gamma3 + Alpha(:,k,state_map(state2,4)) + Beta(3,k,state2)));
    end
    L_all(2*k-1)=log(sum(temp2)+sum(temp3))-log(sum(temp0)+sum(temp1));
    L_all(2*k)  =log(sum(temp1)+sum(temp3))-log(sum(temp0)+sum(temp2));
end
%状态估计,用Alpha的N时刻状态概率和Beta的0时刻状态概率联合估计
if iter==L1
        % determine the start and end state,and return as start_state
    p_endstate=zeros(1,nstates);
    p_startstate=zeros(1,nstates);
    % p_endstate
%用Alpha的最后一个时刻的状态估计编码器的终止状态,计算每个状态概率
    p2=zeros(1,nstates);
    for i=1:nstates
        p2(i)=sum(exp(Alpha(:,L_total+1,i)));
    end
    p1=sum(p2);
    p_endstate=log(p2)-log(p1);
    % p_startstate
%用Beta的第0个时刻的状态估计编码器的起始状态,计算每个状态概率
    p2=zeros(1,nstates);
    beta=zeros(4,nstates);
    delta = -Infty*ones(4,nstates);
    for state1 = 1:nstates
        delta(1,state1) = (-rec(1)-rec(2)+rec(3)*state_map(state1,23))-log(1+exp(L_a(1,1)))-log(1+exp(L_a(2,1)));
        delta(2,state1) = (-rec(1)+rec(2)+rec(3)*state_map(state1,26))-log(1+exp(L_a(1,1)))+L_a(2,1)-log(1+exp(L_a(2,1)));
        delta(3,state1) = (rec(1)-rec(2)+rec(3)*state_map(state1,29))+L_a(1,1)-log(1+exp(L_a(1,1)))-log(1+exp(L_a(2,1)));
        delta(4,state1) = (rec(1)+rec(2)+rec(3)*state_map(state1,32))+L_a(1,1)-log(1+exp(L_a(1,1)))+L_a(2,1)-log(1+exp(L_a(2,1)));
    end
    for i=1:4
        for j=1:nstates
            beta(i,j)=exp(Beta(1,1,state_map(j,i+4))+delta(i,j))+exp(Beta(2,1,state_map(j,i+4))+delta(i,j))...
                +exp(Beta(3,1,state_map(j,i+4))+delta(i,j))+exp(Beta(4,1,state_map(j,i+4))+delta(i,j));
            if(beta(i,j)<1e-300)
                beta(i,j)=-Infty;
            else
                beta(i,j) = log(beta(i,j));  
            end  
        end
    end
    p2=zeros(1,nstates);
    for i=1:nstates
        p2(i)=sum(exp(beta(:,i)));
    end
    p1=sum(p2);
    p_startstate=log(p2)-log(p1);
    % Sc
%用Alpha和Beta联合估计编码器的起始和终止状态,选出概率最大的一个
    p=zeros(1,nstates);
    start_state=0;
    for i=1:nstates
        p(i)=p_endstate(i)+p_startstate(i);
        if p(i)>p(start_state+1)
            start_state=i-1;
        end
    end
end