viterbi_ber.m

实现（2

function [mybit,decoder_output,BER,channel_output,cumulated_metric,survivor_state,y]=viterbi_ber
%changdu=input('Please input the length of the code sequence\n\n');   
biaozhi=input('please choose the restrict length of the encoder\n3:restrict length=3\n7:restrict length=7\n');
save mybiaozhi biaozhi
%mybit=fix(rand(1,changdu)+0.5);                                                        %  产生随机信息
%pause
mybit=[fix(rand(1,84)) 0 1 1 1 1 1 1 1 1 1 1 1 fix(rand(1,104)) 0 1 0 0 fix(rand(1,12))];

%  卷积编码
if biaozhi==3
   G=[1 1 1;1 0 1];                                                                      %  (2,1,3)生成矩阵   
elseif biaozhi==7
   G=[1 1 1 1 0 0 1;1 0 1 1 0 1 1];                                                      %  (2,1,7)生成矩阵
end   
k=1;
right_channel_output=cnv_encd(mybit);
%pause
lgth=length(right_channel_output);
%xunhuan=input('please input the number of the cycle\n'); 
%BER=zeros(1,xunhuan);
%cumulated_metric_table=zeros(1,xunhuan);
%decoder_output_table=zeros(xunhuan,lgth/2-biaozhi+1);
%channel_output_table=zeros(xunhuan,lgth);
 
%  模拟信道误码
%for z=1:xunhuan                                                                                                                                       
   channel=zeros(1,lgth);
   %p=input('please input the number of the erroe bit\n'); 
   %for z=1:2
       %if z==1
           %channel_output=right_channel_output;
           %end
       %if z==2
          %for i=1:444
              %channel(i)=1;
              %end             
                                                                             %信道误码率＝1/z,即每循环一次，就有一个信道误码率p      
  
  %for yz=18:32
   %   channel(yz)=1;
   %end
   %end
   %channel_output=mod(right_channel_output+channel,2)
   %end
%channel_output=input('please input the erroe code sequence\n');
%channel_output_table(z,:)=channel_output;                                                % 数据转移
channel_output=right_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);
        myinput(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=myinput(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;
%decoder_output_table(z,:)=decoder_output;                                                % 数据转移
cumulated_metric=state_metric(1,1);
%cumulated_metric_table(z)=cumulated_metric;                                              % 数据转移
BER=Bit_Error_Rate(mybit,decoder_output);
%channel_BER(z)=Bit_Error_Rate(right_channel_output,channel_output);
y=mod((decoder_output+mybit),2);
%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为不完全译码，译码结果与编码之前的原始信息差移位寄存器的个数
load mybiaozhi
if biaozhi==3
   g=[1 1 1 0 1 1];                                                                         %  (2,1,3)生成矩阵G=[1 1 1;1 0 1];    
elseif biaozhi==7
   g=[1 1 1 0 1 1 1 1 0 0 0 1 1 1];
end
n=length(mybit);
juzhen=zeros(n,(n-1)*2+2*biaozhi);                 %  多余个数＝2*约束长度-2=2*移位寄存器个数
for i=1:n
    juzhen(i,(i-1)*2+1:(i-1)*2+2*biaozhi)=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为译码结果所对应的编码序列与信道输入含误码序列的汉明距离