function [decoder_output,survivor_state,cumulated_metric]=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); 
%卷積碼的維特比譯碼函數 
 
%nxt_stat.m 記錄狀態函數 
%		next_state用于記錄下一個狀態的值 
%		memory_contents用于記錄 
function [next_state,memory_contents]=nxt_stat(current_state,input,L,k) 
	%	將當前狀態值（十進制）轉成位數為k*(L-1)的二進制 
binary_state=deci2bin(current_state,k*(L-1)); 
	%	將輸入狀態值（十進制）轉成位數為k的二進制序列 
binary_input=deci2bin(input,k); 
	%	寄存器組的下一個狀態值（二進制） 
next_state_binary=[binary_input,binary_state(1:(L-2)*k)]; 
	%	將寄存器組的下一個狀態值（二進制）轉成十進制 
next_state=bin2deci(next_state_binary); 
	%	用memory_contents來記錄各個寄存器在下一個狀態下的信息（二進制） 
	%	以便與生成矩陣相乘得出輸出 
memory_contents=[binary_input,binary_state]; 
%nxt_stat.m 記錄狀態函數