viterbi_cpm.m

无线通信中连续相位调制CPM的维特比解调

function [decoder_output,survivor_state,cumulated_metric]=viterbi_cpm(channel_output)
%   CPM信号的Viterbi接收，参数：进制数M、调制指数h（分数形式，互质）、相位状态数phase_num、记忆长度L、基带采样率Ns，	
%   默认的情况下：头尾各有两个零     
 
L=3;h=1/2;Ns=4;M=2;number_of_states=16;
for j=0:number_of_states-1
  for l=0:1
    [next_state,branch_output]=nxt_stat_cpm(j,l);     %   四个参数:状态序号和输入值，后面两个是状态参数  
    %   已验证 input 矩阵是对的。
    input(j+1,next_state+1)=l;
    %   已验证 nextstate 矩阵是对的。
    nextstate(j+1,l+1)=next_state;
    output(j+1,l+1,1:4)=branch_output;
  end
end
state_metric=zeros(number_of_states,2);     %   状态度量是16*2的矩阵
depth_of_trellis=length(channel_output)/4;  %   网格深度（数据长度）
channel_output_matrix=reshape(channel_output,4,depth_of_trellis);   %   信道输出矩阵，即译码输入值
survivor_state=zeros(number_of_states,depth_of_trellis+1);  %   幸存路径就是幸存状态。
%  start decoding of non-tail channel outputs   %   non-tail ??
for i=1:depth_of_trellis-5+1    %   i 是译码数据的序号         L=5;
  flag=zeros(1,number_of_states);   %   flag 是1*16的矩阵
  %     下面怎么改？
  %     定义初始状态转移
  if i < 4
      switch i
          case 1
              step_end=1;stepC=[0];
          case 2
              step_end=2;stepC=[12,13];
          case 3
              step_end=4;stepC=[8,9,10,11];
      end
  elseif i>3 && rem(i,2)==1
      % 每次只有8个可能状态
      step_end=8;stepC=[0,1,2,3,8,9,10,11];
  else
      step_end=8;stepC=[4,5,6,7,12,13,14,15];
  end
  for jjj=0:step_end-1   %   状态数
    j=stepC(jjj+1);
    for l=0:1   %   输入或者0或者1，考虑2进制
      branch_metric=0;
      for index=1:4     %   4 是采样率
          binary_output(index)=output(j+1,l+1,index);
      end
      %     计算度量值，基带复数相减，取模平方，再积分求和。这个metric和采样率有关系！
      branch_metric=branch_metric+sum(abs(transpose(channel_output_matrix(:,i))-binary_output));
      
      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);
end
%  start decoding of the tail channel-outputs
for i=depth_of_trellis-5+2:depth_of_trellis     %   L=5
  %     下面主要改尾巴部分。让它怎么收缩，只考虑输入0
  flag=zeros(1,number_of_states);
  switch i
      case 997
          step_end=8;stepC=[0,1,2,3,8,9,10,11];
      case 998
          step_end=4;stepC=[4,6,12,14];
      case 999
          step_end=1;stepC=[0];
      case 1000
          step_end=1;stepC=[12];
  end
  for jjj=0:step_end-1   %   状态数
      j=stepC(jjj+1);
      branch_metric=0;
      for index=1:4     %   4 是采样率
          binary_output(index)=output(j+1,1,index);
      end
      %     计算度量值，基带复数相减，取模平方，再积分求和。这个metric和采样率有关系！
      branch_metric=branch_metric+sum(abs(transpose(channel_output_matrix(:,i))-binary_output));
      
      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);
%   根据状态度量metric在状态的末尾end找到幸存状态
%   这里的最后状态还有待处理，好像不对。
index_survivor=survivor_end(stepC,state_metric(:,1));
state_sequence(1,depth_of_trellis+1)=index_survivor;
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
decodeder_output_matrix=zeros(1,depth_of_trellis-4);
for i=1:depth_of_trellis-4
  dec_output=input(state_sequence(1,i)+1,state_sequence(1,i+1)+1);
  decoder_output_matrix(1,i)=dec_output;
end
decoder_output=decoder_output_matrix;
cumulated_metric=state_metric(index_survivor+1,1);