siso.m

采用 DVB反向链路（DVB-RCS）标准中的双二元Turbo码matlab仿真程序和仿真结果。支持1/2和1/3码率。支持MAP

function [OutLa ,OutLe]=SISO(InLi,Ar,Br,Wr,Yr,InVar,N,Sc)
%---[OutLa ,OutLe]=SISO(InLi,Ar,Br,Wr,Yr,in_var,N,Sc)----
%Turbo码解码器的核心模块：软入软出解码器
%-----------输出--------------
%InLi：输入的先验对数似然比信息,3*N维
%Ar，Br：接收到的加入噪声的信息比特部分,1*N维
%Wr，Yr：接收到的加入噪声的校验比特部分,1*N维
%InVar：整个系统噪声方差
%N：编码帧长
%Sc：循环状态，由发端传过来，因此该程序现在没有预译码部分
%----------输入-------------
%OutLa：解码得到的全部似然对数比,3*N维
%OutLe：解码得到的似然对数比外信息部分,3*N维
%
trellis=struct('NextStates',[0 7 4 3;4 3 0 7;1 6 5 2;5 2 1 6;...
    6 1 2 5;2 5 6 1;7 0 3 4;3 4 7 0],...
    'Outputs',[ [-1 -1] [1 1] [1 1] [-1 -1];...
    [-1 -1] [1 1] [1 1] [-1 -1];...
    [-1 1] [1 -1] [1 -1] [-1 1];...
    [-1 1] [1 -1] [1 -1] [-1 1];...
    [1 1] [-1 -1] [-1 -1] [1 1];...
    [1 1] [-1 -1] [-1 -1] [1 1];...
    [1 -1] [-1 1] [-1 1] [1 -1];...
    [1 -1] [-1 1] [-1 1] [1 -1] ]);
Lc=2/InVar;
%初始化，初始状态是Sc
alpha(1,1:8)=zeros(1,8);
alpha(1,Sc+1)=1;
beta(N,1:8)=zeros(1,8);
beta(N,Sc+1)=1;
% beta(N+1,1:8)=zeros(1,8);
% beta(N+1,Sc+1)=1;

%
AB_0123=[0 1 2 3];   %分别对应于ab='00''01''10''11'四种情况
AB=[[-1 -1] [-1 +1] [+1 -1] [+1 +1]];
for k=1:N   %帧长
    %输入AB=0的先验概率
    P0=1/(1+exp(InLi(1,k))+exp(InLi(2,k))+exp(InLi(3,k)));
    %输入AB=[0 1 2 3]的先验概率
    P0123=[P0 P0*exp(InLi(1,k)) P0*exp(InLi(2,k)) P0*exp(InLi(3,k))];
    %求gama
    for i=1:8   %状态数，对应于m'(S')
        out=trellis.Outputs(i,:);   %实际要的就是每个状态S'对应的输出和下状态
        NextS=trellis.NextStates(i,:);   %也就是上面矩阵中的1~8行
        gama(k,i,NextS+1)=exp(1/2*Lc*(Ar(k)*AB(1:2:7)+Br(k)*AB(2:2:8)+Wr(k)*out(1:2:7)+Yr(k)*out(2:2:8))).*P0123;
    end
    %normalize gama
    if max(gama(k,:))~=0
        gama(k,:,:)=gama(k,:,:)/max(gama(k,:));
    end
    %alpha(k)
    for i=1:8
        alpha(k+1,i)=sum(gama(k,:,i).*alpha(k,:));
    end
end
for k=N:-1:2
    % beta(k)
    for i=1:8
        tmp=0;
        for j=1:8
            tmp=gama(k,i,j)*beta(k,j)+tmp;
        end
        beta(k-1,i)=tmp;
    end
end
tmp0=zeros(N,1);
tmp1=zeros(N,1);
tmp2=zeros(N,1);
tmp3=zeros(N,1);
for i=1:8
    %AB为00时
    NextS=trellis.NextStates(i,1);
    tmp0=tmp0+alpha(1:N,i).*gama(1:N,i,NextS+1).*beta(1:N,NextS+1);
    %AB为01时
    NextS=trellis.NextStates(i,2);
    tmp1=tmp1+alpha(1:N,i).*gama(1:N,i,NextS+1).*beta(1:N,NextS+1);
    %AB为10时
    NextS=trellis.NextStates(i,3);
    tmp2=tmp2+alpha(1:N,i).*gama(1:N,i,NextS+1).*beta(1:N,NextS+1);
    %AB为11时
    NextS=trellis.NextStates(i,4);
    tmp3=tmp3+alpha(1:N,i).*gama(1:N,i,NextS+1).*beta(1:N,NextS+1);
end
%对双二元Turbo码，总共有3个后验对数似然比
OutLa(1,:)=log(tmp1./tmp0);
OutLa(2,:)=log(tmp2./tmp0);
OutLa(3,:)=log(tmp3./tmp0);
OutLe(1,:)=OutLa(1,:)-Lc*Br-InLi(1,:);
OutLe(2,:)=OutLa(2,:)-Lc*Ar-InLi(2,:);
OutLe(3,:)=OutLa(3,:)-Lc*Ar-Lc*Br-InLi(3,:);