viterbi.m

超宽带程序

%卷积码的译码
function decoder_output=viterbi(G,k,channel_output)
%VITERBI		卷积码的维特比解码器
%[decoder_ouput,survivor_state,cumulated_metric]=viterbi(G,k,channel_output)
%				G是一个n*Lk矩阵，该矩阵的每一行确定了从移位记错器到第n个输出间的连接，
%				是码速率。
%				survivor_state是表示通过网络的最佳路径的矩阵。
%				量度在另一个函数metric（x,y）中给出，而且可根据
%				硬判决和软判决来指定。
%				该算法最小化了量度而不是最大化似然

n=size(G,1);	%取出矩阵G的一维大小，即得出输出端口
	%	检查大小
if rem(size(G,2),k)~=0		%当G列数不是k的整数倍时
   error('Size of G and k do not agree')		%发出出错信息
end
if rem(size(channel_output,2),n)~=0			%当输出量元素个数不是输出端口的整数倍时
   error('Channel output not of the right size')		%发出出错信息
end

L=size(G,2)/k;			%得出移位数，即寄存器组的个数

	%		由于L-1个寄存器组的状态即可表示出输出状态，
	%		所以总的状态数number_of_states可由前L-1个
	%		寄存器的状态组合来确定
number_of_states=2^((L-1)*k);

	%	产生状态转移矩阵、输出矩阵和输入矩阵
for j=0:number_of_states-1		%j表示当前寄存器组的状态因为状态是从零
   									%开始的,所以循环从0到number_of_states-1
                              
      for l=0:2^k-1				%l为从k个输入端的信号组成的状态，总的状
                                 	%态数为2^k，所以循环从0到2^k-1
                                    
      %	nxt_stat完成从当前的状态和输入的矢量得出下寄存器组的一个状态
      [next_state,memory_contents]=nxt_stat(j,l,L,k);
      
      %	input数组值是用于记录当前状态到下一个状态所要的输入信号矢量
      %  input数组的维数：	一维坐标x=j+1指当前状态的值
      %					  	二维坐标y=next_state+1指下一个状态的值
      %	由于Matlab中数组的下标是从1开始的，而状态值
      %	是从0开始的，所以以上坐标值为:状态值+1
      input(j+1,next_state+1)=l;	
      %	branch_output用于记录在状态j下输入l时的输出
      branch_output=rem(memory_contents*G',2);
      %	nextstate数组记录了当前状态j下输入l时的下一个状态
      nextstate(j+1,l+1)=next_state;
      %	output数组记录了当前状态j下输入l时的输出（十进制）
      output(j+1,l+1)=bin2deci(branch_output);
   end
end

%		state_metric数组用于记录译码过程在每状态时的汉明距离
%		state_metric大小为number_of_states×2，（:,1）当前
%		状态位置的汉明距离，为确定值，而（:,2）为当前状态加输入
%		得到的下一个状态汉明距离，为临时值
state_metric=zeros(number_of_states,2);
%		depth_of_trellis用于记录网格图的深度
depth_of_trellis=length(channel_output)/n;
%		输出矩阵，每一列为一个输出状态
channel_output_matrix=reshape(channel_output,n,depth_of_trellis);
%		survivor_state描述译码过程中在网格图中的路径
survivor_state=zeros(number_of_states,depth_of_trellis+1);
%开始无尾信道输出的解码
for i=1:depth_of_trellis-L+1	%i指示网格图的深度
   %		flag矩阵用于记录网格图中的某一列是否被访问过
   flag=zeros(1,number_of_states);
   if i<=L
      step=2^((L-i)*k);		%在网格图的开始处，并不是所有的状态都取
   else							%用step来说明这个变化
      step=1;					%状态数从1→2→4→...→number_of_states
   end
   for j=0:step:number_of_states-1		%j表示寄存器的当前状态
      for l=0:2^k-1					%l为当前的输入
         branch_metric=0;			%用于记录码间距离
         %		将当前状态下输入状态l时的输出output转为n位二进制，以便
         %		计算码间距离（说明：数组坐标大小变化同上）。
         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
         %		选择码间距离较小的那条路径
         %		选择方法：
         %				当下一个状态没有被访问时就直接赋值，否则，用比它小的将其覆盖
         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数组的一维坐标为下一个状态值，二维坐标为此状态
            %	在网格图中的列位置，记录的数值为当前状态，这样就可以从网格图中
            %	某位置的某个状态得出其对应上一个列位置的状态，从而能很方便的完
            %	成译码过程。
            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);		%移动state_metric，将临时值移为确定值
end

%开始尾部信道输出解码
for i=depth_of_trellis-L+2:depth_of_trellis
   flag=zeros(1,number_of_states);
   %  状态数从number_of_states→number_of_states/2→...→2→1
   %	程序说明同上，只不过输入矢量只为0
   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;
         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
%	从最佳路径中产生解码
%	译码过程可从数组survivor_state的最后一个位置向前逐级译码
state_sequence=zeros(1,depth_of_trellis+1);
%	survivor_state数组的最后的输出状态肯定是“0”
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
decorder_output_matrix=zeros(k,depth_of_trellis-L+1);
for i=1:depth_of_trellis-L%+1
   %		根据数组input的定义来得出从当前状态到下一个状态的输入信号矢量
   dec_output_deci=input(state_sequence(1,i)+1,state_sequence(1,i+1)+1);
   %		转成二进制信号
   dec_output_bin=deci2bin(dec_output_deci,k);
   %		将一次译码存入译码输出矩阵decoder_output_matrix相应的位置
   decoder_output_matrix(:,i)=dec_output_bin(k:-1:1)';
end
%		按照一维序列形式重新组织输出
decoder_output=reshape(decoder_output_matrix,1,k*(depth_of_trellis-L));%+1
%		state_metric为网格图最后一个列位置中“0”状态位置的汉明距
%		离，这个值就是整个译码过程中的汉明距离。
cumulated_metric=state_metric(1,1);