viterbi.m

数字通信第四版原书的例程

function [msg, expen, codd] = viterbi(code, tran_func, leng, tran_prob, plot_flag);
%VITERBI Convolution decode using Viterbi algorithm.
%       MSG = VITERBI(CODE, TRAN_FUNC) uses Viterbi algorithm to recover the
%       message in CODE with transfer function matrix provided in TRAN_FUNC.
%       The transfer function matrix is provided in state system form as
%       described in the help file for SIM2TRUN. You can use SIM2TRAN to
%       generate the transfer function from a block diagram of convolution
%       code. CODE is received in a binary symmetric channel (BSC). All
%       elements in CODE should be binary. CODE is a matrix with the column
%       size 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 Viterbi algorithm
%       to recover the message in CODE with 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 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; second row specifies the probability
%       of zero signal transfer probability to the receiving range; 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 soft-decision decode.
%
%       [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 of PLOT_FLAG is a positive integer.
%       PLOT_FLAG is the figure number to be used to plot the diagram. In the
%       trellis diagram, yellow paths are the possible path; 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, SIM2TRAN, SIM2GEN, OCT2GEN.

%       Wes Wang 8/6/94, 10/11/95.
%       Copyright (c) 1995-96 by The MathWorks, Inc.
%       $Revision: 1.1 $  $Date: 1996/04/01 18:04:36 $

% routine check.
if nargin < 2
    error('Not enough input variable 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;
  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);
    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;