function [mybit,decoder_output_table,BER,right_channel_output,channel_output_table,cumulated_metric_table,survivor_state]=viterbi_2_1_3
changdu=input('Please input the length of the code sequence\n\n');   
mybit=fix(rand(1,changdu)+0.5);

%  (2,1,3)卷積編碼
G=[1 1 1;1 0 1];                                                                          %  (2,1,3)生成矩陣G=[1 1 1;1 0 1]; 
k=1;
right_channel_output=cnv_encd(mybit);
lgth=length(right_channel_output);
BER=zeros(1,lgth);
cumulated_metric_table=zeros(1,lgth);
decoder_output_table=zeros(lgth,lgth/2-2);
channel_output_table=zeros(lgth,lgth);

%  模擬信道誤碼
for z=1:lgth                                                                             % 產生隨機誤碼                                                                           

%    channel=zeros(1,lgth);                                                              % 通過改變lgth的值 ,調整誤碼個數                                                 
%    for i=1:z:lgth
%        channel(i)=1;
%    end

if z==1
   channel=zeros(1,lgth);  
else
    channel=fix(rand(1,lgth)+0.2);                                                       % 0.5可改
end
channel_output=mod(right_channel_output+channel,2);
channel_output_table(z,:)=channel_output;                                                % 數據轉移
n=size(G,1);                                                                             % n=G的行數     

%  檢查大小
if rem(size(G,2),k)~=0
    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;                                                                           % 約束長度＝寄存器個數+1
number_of_states=2^((L-1)*k);                                                            % 狀態數

%  產生狀態轉移矩陣、輸出矩陣和輸入矩陣
for j=0:number_of_states-1                                                               % 對(2,1,3) j=0,1,2,3
    for l=0:2^k-1                                                                        % l=0,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);                                         % 狀態圖上根據寄存器狀態和input作出的編碼
        nextstate(j+1,l+1)=next_state;                                                   % 生成矩陣，含有十進制表示的下一狀態
        output(j+1,l+1)=bin2deci(branch_output);                                         % 生成矩陣，含有十進制表示的卷積編碼
    end
end
state_metric=zeros(number_of_states,2);                                                  % 4行，2列 (初始化)
depth_of_trellis=length(channel_output)/n;                                               % 需要譯碼的個數   設為6
channel_output_matrix=reshape(channel_output,n,depth_of_trellis);                        % 將經過信道的編碼序列轉換為以需要譯碼的個數為行數，以n為列數的矩陣
survivor_state=zeros(number_of_states,depth_of_trellis+1);                               % 4行，7列 (初始化)

%  開始無尾信道輸出解碼
for i=1:depth_of_trellis-L+1                                                             % i=1,2,3,4
    flag=zeros(1,number_of_states);                                                      % 1行，4列 (初始化)
    if i<=L
        step=2^((L-i)*k);                                                                % i=1,2,3  step= 4,2,1
    else
        step=1;                                                                          % i=4  step=1
    end
    for j=0:step:number_of_states-1                                                      % i=1 j=0;i=2 j=0, 2;i=3,4 j=0,1,2,3 
        for l=0:2^k-1                                                                    % l=0,1
            branch_metric=0;                                                             % 初始化
            binary_output=deci2bin(output(j+1,l+1),n);                                   % 將此時刻十進制的編碼結果用二進制表示
            for ll=1:n                                                                   % ll=1,2
                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(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

%  開始有尾信道輸出解碼
for i=depth_of_trellis-L+2:depth_of_trellis
    flag=zeros(1,number_of_states);
    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

%  從最佳路徑中產生解碼
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
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);
    decoder_output_matrix(:,i)=dec_output_bin(k:-1:1)';
end
decoder_output=reshape(decoder_output_matrix,1,k*(depth_of_trellis-L+1));
decoder_output_table(z,:)=decoder_output;                                                % 數據轉移
cumulated_metric=state_metric(1,1);
cumulated_metric_table(z)=cumulated_metric;                                              % 數據轉移
BER(z)=Bit_Error_Rate(mybit,decoder_output);
%  BER1(z)=Bit_Error_Rate(right_channel_output,channel_output)
end

%----------------------------------------subfunction Bit_Error_Rate-------------------------------------------------

function y=Bit_Error_Rate(a,b)
y=nnz(mod(a+b,2))/length(a);

%-------------------------------------------subfunction cnv_encd----------------------------------------------------

function the_output1=cnv_encd(mybit)       %  the_output1為完全正確譯碼，編碼輸出結果多出2*移位寄存器個數；  the_output2為不完全譯碼，譯碼結果與編碼之前的原始信息差移位寄存器的個數
g=[1 1 1 0 1 1];
n=length(mybit);
juzhen=zeros(n,(n-1)*2+6);                 %  多余個數＝2*約束長度-2=2*移位寄存器個數
for i=1:n
    juzhen(i,(i-1)*2+1:(i-1)*2+6)=g;
end
the_output=mybit*juzhen;
the_output=mod(the_output,2);
the_output1=the_output;                    %  需要時可以把分號去掉，看完全編碼和不完全編碼的結果
the_output2=the_output(1,1:n*2);




% survivor_state的每一行對應每一個狀態，每一列對應每一個時刻，它的數值為從前一時刻到這一時刻、進入到所在的行狀態所保留的幸存路徑的來源狀態   
% cumulated_metric為譯碼結果所對應的編碼序列與信道輸入含誤碼序列的漢明距離