www.pudn.com > ctc-matlab.rar > logmap.m, change:2005-03-21,size:9272b

```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
```