viterbi.m

实现了卷积码的viterbi译码

    for i = 1 : n_std_sta 
      flipped(i) =  bi2de(fliplr(de2bi(i-1, M)));
    end;
  end;
  
  pre_state = 1;
  inp_pre = de2bi([0:K2-1]', K);
  cur_sta_pre = de2bi([0:n_std_sta-1], M);
  starter = 0;
  msg = [];
  expen = [];
  codd = [];
  
  for i = 1 : n_code
    % make room for one more storage.
    trace_pre = rem(i-2+leng, leng) + 1; % previous line of the trace.
    trace_num = rem(i-1, leng) + 1; 	% current line of the trace.
    expense(trace_num,:) = solution;
    
    if expen_flag
      for j = 1 : length(pre_state)
	jj = pre_state(j) - 1; 		% index number - 1 is the state.
	cur_sta = cur_sta_pre(pre_state(j),:)';
	indx_j = (pre_state(j) - 1) * n_std_sta;
	for num_N = 1 : N
	  expen_work(num_N) = max(find(tran_prob(1,:) <= code(i, num_N)));
	end;
	for num_k = 1 : K2
               inp = inp_pre(num_k, :)';
               if isempty(C)
                  tran_indx = pre_state(j) + n_std_sta * (num_k - 1);
                  nex_sta = A(tran_indx, :)';
                  out = B(tran_indx, :)';
               else
                  out = rem(C * cur_sta + D * inp,2);
                  nex_sta = rem(A * cur_sta + B * inp, 2);
               end;
               nex_sta_de = bi2de(nex_sta') + 1;
               % find the expense by the transfer probability
               expen_0 = find(out' <= 0.5);
               expen_1 = find(out' > 0.5);
               loca_exp = sum([tran_prob(2,expen_work(expen_0)) 0])...
                  +sum([tran_prob(3,expen_work(expen_1)) 0]);
               tmp = (nex_sta_de-1)*n_std_sta + pre_state(j);
               if isnan(expense(trace_num, tmp))
                  expense(trace_num, tmp) = loca_exp;
                  solu(trace_num, nex_sta_de + indx_j) = num_k;
               elseif expense(trace_num, tmp) < loca_exp
                  expense(trace_num, tmp) = loca_exp;
                  solu(trace_num, nex_sta_de + indx_j) = num_k;
               end;
            end;
         end;
      else
         for j = 1 : length(pre_state)
            jj = pre_state(j) - 1;           % index number - 1 is the state.
            cur_sta = cur_sta_pre(pre_state(j),:)';
            indx_j = (pre_state(j) - 1) * n_std_sta;
            for num_k = 1 : K2
               inp = inp_pre(num_k, :)';
               if isempty(C)
                  tran_indx = (num_k - 1) * n_std_sta + pre_state(j);
                  nex_sta = A(tran_indx, :)';
                  out = B(tran_indx, :)';
               else
                  out = rem(C * cur_sta + D * inp,2);
                  nex_sta = rem(A * cur_sta + B * inp, 2);
               end;
               nex_sta_de = bi2de(nex_sta') + 1;
               loca_exp = sum(rem(code(i, :) + out', 2));
               tmp = (nex_sta_de-1)*n_std_sta + pre_state(j);
               if isnan(expense(trace_num, tmp))
                  expense(trace_num, tmp) = loca_exp;
                  solu(trace_num, nex_sta_de + indx_j) = num_k;
               elseif expense(trace_num, tmp) > loca_exp
                  expense(trace_num, tmp) = loca_exp;
                  solu(trace_num, nex_sta_de + indx_j) = num_k;
               end;
            end;
         end;
      end;
      
      aft_state = [];
      for j2 = 1 : n_std_sta
         if max(~isnan(expense(trace_num, [1-n_std_sta : 0] + j2*n_std_sta)))
            aft_state = [aft_state, j2];
         end
      end;
      
      % go back one step to re-arrange the lines.
      if (plot_flag(1))
         plot_curv_x = [];
         plot_curv_y = [];
         for j = 1 : n_std_sta
            tmp = expense(trace_num, (j-1)*n_std_sta+1: j*n_std_sta);
            if expen_flag
               expen_rec(j) = max(tmp(find(~isnan(tmp))));
            else
               expen_rec(j) = min(tmp(find(~isnan(tmp))));
            end;
            for j2 = 1 : length(tmp)
               if ~isnan(tmp(j2))
                  plot_curv_x = [plot_curv_x, i-1+.1, i, NaN];
                  plot_curv_y = [plot_curv_y, ...
                        flipped(j2), ...
                        flipped(j), NaN];
               end;
            end;
         end;
      end;
      
      %find out the plot lines:
      
      if plot_flag(1)
         % plot only; no computation function here.
         plot_curv_x = [get(plot_y, 'XData'), NaN, plot_curv_x];
         plot_curv_y = [get(plot_y, 'YData'), NaN, plot_curv_y];
         for j = 1 : length(aft_state)
            xx=[xx; text(i, flipped(aft_state(j)), num2str(expen_rec(aft_state(j))))];
            set(xx,'Color',[0 1 1]);
         end;
         %        plot(plot_curv_x, plot_curv_y, 'y-');
         set(plot_y, 'XData', plot_curv_x, 'YData', plot_curv_y);
         if i < leng
            drawnow
         end;
      end;
      
      sol = expense(trace_num,:);
      % decision making.
      if i >= leng
         trace_eli = rem(trace_num, leng) + 1;
         % strike out the unnecessary.
         for j_k = 1 : leng - 2
            j_pre = rem(trace_num - j_k - 1 + leng, leng) + 1;
            sol = vitshort(expense(j_pre, :), sol, n_std_sta, expen_flag);
         end;
         
         tmp = (ones(n_std_sta,1) * expense(trace_eli, [starter+1:n_std_sta:PowPowM]))';
         sol = sol + tmp(:)';
         
         if expen_flag
            loc_exp =  max(sol(find(~isnan(sol))));   
         else
            loc_exp =  min(sol(find(~isnan(sol))));
         end
         dec = find(sol == loc_exp);
         dec = dec(1);
         dec = rem((dec - 1), n_std_sta);
         
         inp = de2bi(solu(trace_eli, starter*n_std_sta+dec+1)-1, K);
         if isempty(C)
            tran_indx = pre_state(j) + n_std_sta * (num_k - 1);
            out = B(tran_indx, :)';
         else
            cur_sta = cur_sta_pre(starter+1, :)';
            out = rem(C * cur_sta + D * inp(:),2);
         end;
         
         msg = [msg; inp];
         codd = [codd; out'];
         [n_msg, m_msg] = size(msg);
         if nargout > 1
            if expen_flag
               % find the expense by the transfer probability
               expen_0 = find(out' <= 0.5);
               expen_1 = find(out' > 0.5);
               loc_exp = sum([tran_prob(2,expen_work(expen_0)) 0])...
                  +sum([tran_prob(3,expen_work(expen_1)) 0]);
            else
               loc_exp = sum(rem(code(n_msg, :) + out', 2));
            end
            expen = [expen; loc_exp];
         end;
         
         if plot_flag(1)
            plot_curv_x = [get(plot_r, 'XData'), n_msg-1+.1, n_msg, NaN];
            plot_curv_y = [get(plot_r, 'YData'), flipped(starter+1), flipped(dec+1), NaN];
            set(plot_r, 'XData', plot_curv_x, 'YData', plot_curv_y,'LineWidth',2);
            drawnow
         end;
         
         starter = dec;
      end; %(if i >= leng)
      pre_state = aft_state;
      aft_state = [];
   end;
   
   for i =  1 : leng-1
      sol = solution;
      trace_eli = rem(trace_num + i, leng) + 1;
      if i < leng-1
         sol(1:n_std_sta) = expense(trace_num, 1:n_std_sta);
         for j_k = 1 : leng - 2 - i
            j_pre = rem(trace_num - j_k - 1 + leng, leng) + 1;
            sol = vitshort(expense(j_pre, :), sol, n_std_sta, expen_flag);
         end;
         tmp = (ones(n_std_sta,1) * expense(trace_eli, [starter+1:n_std_sta:PowPowM]))';
         sol = sol + tmp(:)';
         if expen_flag
            loc_exp =  max(sol(find(~isnan(sol))));   
         else
            loc_exp =  min(sol(find(~isnan(sol))));
         end
         dec = find(sol == loc_exp);
         dec = dec(1);
         dec = rem((dec - 1), n_std_sta);
      else
         dec = 0;
      end;
      
      inp = de2bi(solu(trace_eli, starter*n_std_sta+dec+1)-1, K);
      cur_sta = de2bi(starter, num_state);
      out = rem(C*cur_sta' + D * inp', 2);
      msg = [msg; inp];
      codd = [codd; out'];
      [n_msg, m_msg] = size(msg);
      
      if nargout > 1
         if expen_flag
            % find the expense by the transfer probability
            expen_0 = find(out' <= 0.5);
            expen_1 = find(out' > 0.5);
            loc_exp = sum([tran_prob(2,expen_work(expen_0)) 0])...
               +sum([tran_prob(3,expen_work(expen_1)) 0]);
         else
            loc_exp = sum(rem(code(n_msg, :) + out', 2));
         end
         expen = [expen; loc_exp];
      end;
      
      if plot_flag(1)
         plot_curv_x = [get(plot_r, 'XData'), n_msg-1+.1, n_msg, NaN];
         plot_curv_y = [get(plot_r, 'YData'), flipped(starter+1), flipped(dec+1), NaN];
         set(plot_r, 'XData', plot_curv_x, 'YData', plot_curv_y,'LineWidth',2);
         drawnow
      end;
      starter = dec;
   end;
   
   
   % cut the extra message length
   [n_msg, m_msg] = size(msg);
   msg(n_msg-M+1:n_msg, :) = [];
   if plot_flag(1)
      set(xx,'Color',[0 1 1]);
      hold off
   end;
end;
% end of VITERBI.M