viterbi.m

实现了卷积码的viterbi译码

function [msg, expen, codd] = viterbi2(code, tran_func, leng, tran_prob, plot_flag);
%VITERBI
%
%WARNING: This is an obsolete function and may be removed in the future.
%         Please use VITDEC instead.

%VITERBI Convolution decode using Viterbi algorithm.
%       MSG = VITERBI(CODE, TRAN_FUNC) uses the Viterbi algorithm to recover the
%       message in CODE with the transfer function matrix provided in TRAN_FUNC.
%       The transfer function matrix is provided in state system form.
%       CODE is received in a binary symmetric channel (BSC). All
%       elements in CODE should be binary. CODE is a matrix with the column
%       size the same as the code length. See function OCT2GEN for the format
%       of an octal form of transfer function matrix. 
%
%       MSG = VITERBI(CODE, TRAN_FUNC, LEN) keeps the trace memory length
%       no larger than LEN.
%       
%       MSG = VITERBI(CODE, TRAN_FUNC, LEN, TRAN_PROB) uses the Viterbi algorithm
%       to recover the message in CODE with the transfer function matrix provided
%       in TRAN_FUNC. TRAN_FUNC is in octal form. CODE is assumed to be
%       received in a discrete memoryless channel (DMC). The element in CODE
%       may not be in a binary form. TRAN_PROB is a parameter specifying the
%       transfer probability. TRAN_PROB is a three row matrix specifying the
%       transfer probability in soft decision case. The first row of TRAN_PROB
%       specifies the receiving range; the second row specifies the probability
%       of "zero" signal transfer probability to the receiving range; the third
%       row specifies the probability of "one" signal transfer probability to
%       the receiving range. For example,
%       TRAN_PROB = [-Inf 1/10, 1/2, 9/10; 
%                     2/5, 1/3, 1/5, 1/15; 
%                    1/16, 1/8, 3/8, 7/16]
%       means that transferred "zero" signal has 2/5 probability receiving in
%       range (-Inf, 1/10], 1/3 in range (1/10, 1/2], 1/5 in range 
%       (1/2, 9/10], and 1/15 in range (9/10, Inf). The "one" signal has
%       1/16 probability in range (-Inf, 1/10], 1/8 in range (1/10, 1/2], 
%       3/8 in range (1/2, 9/10], 7/16 in range (9/10, Inf). When TRAN_PROB is
%       assigned to be zero, this function processes hard-decision decode.
%       In order to do soft-decision decode,  you must use a three row matrix.
%       For more information, see the Tutorial chapter of the Communications 
%       Toolbox User's Guide. 
%
%       [MSG, DIST] = VITERBI(CODE, TRAN_FUNC) returns the Hamming distance
%       in each step of recovering the message.
%
%       [MSG, REC_PROB, SURVIVOR] = VITERBI(...) outputs the survivors of the
%       code.
%
%       [...] = VITERBI(CODE, TRANS_FUNC, LEN, TRAN_PROB, PLOT_FLAG) plots the
%       trellis diagram for the decoding.  PLOT_FLAG, a positive integer,
%       is the figure number to be used to plot the diagram. In the
%       trellis diagram, yellow paths are the possible paths; red paths are
%       the survivors. The number in the state position is the metric or the
%       accumulated Hamming distance in the BSC case. When PLOT_FLAG is not
%       an integer or it is a negative number, the function does not plot
%       the diagram. Note that the function computation is very slow if the
%       plot is required.
%
%       See also CONVENCO, OCT2GEN.

%       Wes Wang 8/6/94, 10/11/95.
%       Copyright 1996-2002 The MathWorks, Inc.
%       $Revision: 1.19 $

% routine check.
if nargin < 2
   error('Not enough input variables for VITERBI.');
elseif nargin < 3
   leng = 0;
end;

if (exist('vitcore') == 3) 
  if isstr(tran_func)
    tran_func = sim2tran(tran_func);
  end;
  
  [len_m, len_n] = size(tran_func);
  if tran_func(len_m, len_n)>=0 & max(max(tran_func)) ~= Inf
    [tran_func, code_param] = oct2gen(tran_func);
    N = code_param(1);
    K = code_param(2);
    M = code_param(3);
    G = [];
    for i = 1 : M+1
      G = [G ; tran_func(:, i:M+1:N*(M+1))];
    end;
    
    a = eye(K * (M-1));
    a = [zeros(K, K), zeros(K, length(a)); a, zeros(length(a), K)]; 
    
    % B matrix:
    b = [eye(K); zeros(length(a)-K, K)];
    
    % C matrix
    c = [];
    for i = 1 : M
      c = [ c, G(i*K + 1 : (i + 1) * K, :)'];
    end;
      
    % D matrix
    d = G(1:K, :)';
    M = length(a);
    
    tran_func = [a b; c d];
    [len_m, len_n] = size(tran_func);
    if len_m < 4
      tran_func = [tran_func; zeros(4-len_m, len_n)];
      [len_m, len_n] = size(tran_func);
    end;
    tran_func = [tran_func zeros(len_m, 1)];
    len_n = len_n + 1;
    tran_func(1, len_n) = N;
    tran_func(2, len_n) = K;
    tran_func(3, len_n) = M;
    tran_func(len_m, len_n) = -Inf;
  elseif tran_func(len_m, len_n) == -Inf
    N = tran_func(1, len_n);
    K = tran_func(1, len_n);
    M = tran_func(1, len_n);      
  end;
   
  % This part is only for detecting the size of CODE, see if its size matches the vitcore's reqirement.
  % This is only because there isn't input detection inside VITCORE.C and not encourage users to directly use it.

  [n_code, m_code] = size(code);
  if max (n_code, m_code) <= 1
    msg = floor(code);
    if nargout > 1
      expen = [];
      if nargout == 3
	codd = [];
      end;
    end;
    return;
  end;
  
  if min (n_code, m_code) == 1
    code = code(:);
    code = vec2mat(code, N);
  elseif m_code ~= N
    error('CODE length does not match the TRAN_FUNC dimension in VITERBI.')
  end;
  
  [n_code, m_code] = size(code);
  if n_code == 1
    code = [code; code*0];
  end;
  
  if nargin == 5
    eval('[msg, expen, codd]=vitcore(code, tran_func, leng, tran_prob, plot_flag);');
  elseif nargin == 4
    eval('[msg, expen, codd]=vitcore(code, tran_func, leng, tran_prob);');
  else
    eval('[msg, expen, codd]=vitcore(code, tran_func, leng);');
  end;
else
  % the case of unlimited delay and memory.
  [n_code, m_code] = size(code);
  if (leng <= 0) | (n_code <= leng)
    leng = n_code;
  end;
  if leng <= 1
    leng = 2;
  end;
  
  if nargin < 4
    expen_flag = 0; 			% expense flag == 0; BSC code. == 1, with TRAN_PROB.
  else
    [N, K] = size(tran_prob);
    if N < 3
      expen_flag = 0;
    else
      expen_flag = 1;
      expen_work = zeros(1, N); 	% work space.
      tran_prob(2:3, :) = log10(tran_prob(2:3, :));
    end;
  end;
  
  if nargin < 5
    plot_flag(1) = 0;
  else
    if (plot_flag(1) <= 0)
      plot_flag(1) = 0;
    elseif floor(plot_flag(1)) ~= plot_flag(1)
      plot_flag(1) = 0;
    end;
  end;
  
  if isstr(tran_func)
    tran_func = sim2tran(tran_func);
  end;
  % initial parameters.
  [A, B, C, D, N, K, M] = gen2abcd(tran_func);
  
  if m_code ~= N
    if n_code == 1
      code = code';
      [n_code, m_code] = size(code);
    else
      error('CODE length does not match the TRAN_FUNC dimension in VITERBI.')
    end;
  end;
  
  % state size
  if isempty(C)
    states = zeros(tran_func(2,1), 1);
  else
    states = zeros(length(A), 1);
  end;
  num_state = length(states);
  
  % The STD state number reachiable is
  n_std_sta = 2^num_state;
  
  PowPowM = n_std_sta^2;
  
  % at the very begining, the current state number is 1, only the complete zeros.
  cur_sta_num = 1;
  
  % the variable expense is divided into n_std_sta section of size n_std_sta
  % length vector for each row. expense(k, (i-1)*n_std_sta+j) is the expense of
  % state j transfering to be state i
  expense = ones(leng, PowPowM) * NaN;
  expense(leng, [1:n_std_sta]) = zeros(1, n_std_sta);
  solu = zeros(leng, PowPowM);
  solution = expense(1,:);
  
  K2 = 2^K;
  
  if plot_flag(1)
    plot_flag(1) = figure(plot_flag(1));
    drawnow;
    delete(get(plot_flag(1),'child'))
    set(plot_flag(1), 'NextPlot','add', 'Clipping','off');
    if length(plot_flag) >= 5
      set(plot_flag(1),'Position',plot_flag(2:5));
    end;
    axis([0, n_code, 0, n_std_sta-1]);
    handle_axis = get(plot_flag(1),'Child');
    xx=text(0, 0, '0');
    set(xx,'Color',[0 1 1]);
    hold on
    set(handle_axis, 'Visible', 'off');
    plot_w = plot(0, 0, 'w--');
    plot_y = plot(0, 0, 'y-');
    plot_r = plot(0, 0, 'r-');
    for i = 0 : n_std_sta-1
      text(-n_code/8, i, ['State ', num2str(bi2de(fliplr(de2bi(i,num_state))))]);
    end;
    text(n_code/2.1, -n_std_sta/15, 'Time');
    set(handle_axis, 'Visible', 'off');
    flipped = zeros(1, n_std_sta);