viterbi_decoder.m

实现GSM通信的快速接入控制信道（FACCH）的信道编码的MATLAB仿真

%VITERBI _decorder   卷积码的维特比硬判决解码器 

%输入参数：             G：生成矩阵，G是一个n×LK矩阵，该矩阵的每一行确定了从移位计错器到第n个输出间的连接
%                       k：并行输入比特数
%          channel_output：需要解码的数据
%输出参数：decoder_output：解码数据输出（已经截掉尾部的零）
%survivor_state为记录最佳路径的矩阵,cumulated_metric为累积最小汉明距
%程序中调用的外部函数：nxt_stat.m：利用当前状态和输入数据求网格图的下一转移状态
%                      bin2deci.m：二进制转十进制
%                      deci2bin.m：十进制转二进制
%                      metric.m  ：计算汉明距离
%该译码思想是按照最小汉明距准则确定幸存路径，然后根据幸存路径由后到前得到各级的状态，最后根据各级状态和状态转移图来确定各级的输入
function [decoder_output,survivor_state,cumulated_metric]=viterbi_decoder(G,k,channel_output)

n=size(G,1); %算出n的大小
% 检查输入各矩阵格式是否合格
if rem(size(G,2),k)~=0           %寄存器个数不是k的整数倍，出错
   error('Size of G and k do not agree')
end
if rem(size(channel_output,2),n)~=0 %输入序列与编码器不匹配，出错
   error('Channel out not of the right size')
end
%确定L，用十进制数表示寄存器状态
L=size(G,2)/k;

number_of_states=2^((L-1)*k);
%产生网格图
%产生状态转移矩阵，输出矩阵和输入矩阵
%在nextstate(状态数x2)(即状态转移矩阵)中， 两列分别为i状态下输入0，1所得到的下一状态（为十进制数）
%在input(状态数x状态数)矩阵中保存状态的转移关系，input(i,j)是状态i到状态j的输入
%在output(状态数x2)(即输出矩阵)中，两列分别为i状态下输入0，1理想输出（为十进制数）
for j=0:number_of_states-1
   for l=0:2^k-1
      [next_state,memory_contents]=nxt_stat(j,l,L,k);
      input(j+1,next_state+1)=l;   %输入矩阵，由当前状态和下一状态确定输入的十进制表示
      branch_output=rem(memory_contents*G',2);
      nextstate(j+1,l+1)=next_state;  %转移矩阵，下一状态的十进制表示
      output(j+1,l+1)=bin2deci(branch_output);  %输出矩阵，输出的十进制表示
   end
end
%state_metric状态距离矩阵（状态数x2）,第二列表明到达下一状态的最小距离，第一列表明当前状态的最小距离.
state_metric=zeros(number_of_states,2);  %记录当前级各节点的最小汉明距和下一级各节点的汉明距
depth_of_trellis=length(channel_output)/n; %级数
%channel_output_matric(输出位数x网格深度)：将输入按输出位数分为矩阵，每列为一个时刻的输出。
channel_output_matrix=reshape(channel_output,n,depth_of_trellis); %矩阵的每一行代表编码器的一个输出
%survivor_state(状态数x结点深度)幸存状态矩阵：（i,j）表示第j个时刻时第i个状态是由谁过来的。
survivor_state=zeros(number_of_states,depth_of_trellis+1);  %通过网格的最佳路径的矩阵

%开始无尾信道输出的解码（即去掉结尾(L-1)个零）
%body encoder
for i=1:depth_of_trellis-L+1
   flag=zeros(1,number_of_states); %标志各状态有否被访问过
                                   %flag(1x状态数)：第j位表示在一个时刻中第j个状态是否到达过，1表示到达过。
%确定状态转移的步长，起始阶段的步长与中间态的步长不同
   if i<=L
      step=2^((L-i)*k); %表示第i级各状态的最小差值
   else
      step=1;
   end
   %进行加－比－选
   for j=0:step:number_of_states-1
      for l=0:2^k-1
         branch_metric=0;
         binary_output=deci2bin(output(j+1,l+1),n);  %将理想输出转化为二进制
         for ll=1:n  %计算汉明距
            branch_metric=branch_metric+metric(channel_output_matrix(ll,i),binary_output(ll));
         end
       %在AWGN信道下，最大似然估计转化为求最小汉明距
       %如果下一状态度量距离大于当前距离加汉明距，或是下一状态未被遍历过则设为当前状态下一状态的幸存状态，当前距离加汉明距设为下一状态的距离  
         if((state_metric(nextstate(j+1,l+1)+1,2)>state_metric(j+1,1)+branch_metric)|flag(nextstate(j+1,l+1)+1)==0)
            state_metric(nextstate(j+1,l+1)+1,2)=state_metric(j+1,1)+branch_metric; %更改汉明距
            survivor_state(nextstate(j+1,l+1)+1,i+1)=j; %更改幸存路径
            flag(nextstate(j+1,l+1)+1)=1;
         end
      end
   end
   state_metric=state_metric(:,2:-1:1); %将两列互换，目的是上一次计算的下一状态距离变为下一次计算的上一状态距离
end

%尾部信道输出的解码，最后的那些级没有包含所有的状态
%开始尾部信道输出解码
%归零计算
%tailer encoder
for i=depth_of_trellis-L+2:depth_of_trellis
   flag=zeros(1,number_of_states);        %flag(1x状态数)：第j位表示在一个时刻中第j个状态是否到达过，1表示到达过。
   last_stop=number_of_states/(2^((i-depth_of_trellis+L-2)*k));%所剩余的状态数
  %认为尾部的数据为零，每次状态转移，只有输入为零一条路径    
   for j=0:last_stop-1
      branch_metric=0;
      binary_output=deci2bin(output(j+1,1),n);
      for ll=1:n
         branch_metric=branch_metric+metric(channel_output_matrix(ll,i),binary_output(ll));
      end
      
      if((state_metric(nextstate(j+1,1)+1,2)>state_metric(j+1,1)+branch_metric)|flag(nextstate(j+1,1)+1)==0)
         state_metric(nextstate(j+1,1)+1,2)=state_metric(j+1,1)+branch_metric;%在第2列保存下一状态的距离
         survivor_state(nextstate(j+1,1)+1,i+1)=j; %更改幸存路径
         flag(nextstate(j+1,1)+1)=1;
      end
   end
   state_metric=state_metric(:,2:-1:1);%将两列互换，目的是上一次计算的下一状态距离变为下一次计算的上一状态距离
end

%开始回溯最佳路径，从最佳路径中找出解码
%state_sequence(1x结点深度)矩阵，1～dep分别记载各个阶段的路径（即前一状态数）。1～dep分别记载各个阶段的路径（即前一状态数）。
%从最佳路径中产生解码
%由后到前得到各级的状态
state_sequence=zeros(1,depth_of_trellis+1);
state_sequence(1,depth_of_trellis)=survivor_state(1,depth_of_trellis+1);%开始回溯最佳路径
for i=1:depth_of_trellis
   state_sequence(1,depth_of_trellis-i+1)=survivor_state((state_sequence(1,depth_of_trellis+2-i)+1),depth_of_trellis-i+2);
end

%解码输出矩阵，输出码元个数为结点深度减L
decoder_output_matrix=zeros(k,depth_of_trellis-L+1); %由各级的状态和输入矩阵得到各级的输入
for i=1:depth_of_trellis-L+1
   dec_output_deci=input(state_sequence(1,i)+1,state_sequence(1,i+1)+1); %解出码元
   dec_output_bin=deci2bin(dec_output_deci,k);%按k将输出码元转换
   decoder_output_matrix(:,i)=dec_output_bin(k:-1:1)';
end
%由各级的输入的得到输入序列
decoder_output=reshape(decoder_output_matrix,1,k*(depth_of_trellis-L+1));
cumulated_metric=state_metric(1,1);