%卷積碼的譯碼
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);