stc/exp/ml_stc_linux_make_v1.0/matlab/STC matlab implementation/dual_viterbi.m

   1 function [ y min_cost paths ] = dual_viterbi( code, x, w, m )\r
   2 %DUAL_VITERBI pure Matlab implementation of the Viterbi algorithm over \r
   3 % syndrome trellis. For given code, binary cover vector x, weight \r
   4 % vector w, finds the closest binary 'stego' vector y that produces given \r
   5 % syndrome m. Use dual_viterbi() for faster C++ implementation.\r
   6 %\r
   7 % Example:\r
   8 %   code = create_code_from_submatrix([7 5 1],3);\r
   9 %   w = ones(code.n,1);\r
  10 %   x = double(rand(code.n,1)<0.5);\r
  11 %   m = double(rand(sum(code.shift),1)<0.5);\r
  12 %   [y min_cost] = dual_viterbi(code, x, w, m);\r
  13 %   x = x'; x\r
  14 %   y = y'; y\r
  15 %   m = [m' ; calc_syndrome(code,y)']; m\r
  16 %   min_cost\r
  17 %\r
  18 % Tomas Filler (tomas.filler@binghamton.edu)\r
  19 % http://dde.binghamton.edu/filler\r
  20 \r
  21 %% check input parameters\r
  22 if length(m)~=sum(code.shift)\r
  23     error('dualViterbi:wrongInputArg', 'Syndrome must be of the same length as the number of shifts in code.');\r
  24 end\r
  25 if length(x)~=length(w)\r
  26     error('dualViterbi:wrongInputArg', 'Cover vector must be of the same size as weight vector.');\r
  27 end\r
  28 if length(x)~=code.n\r
  29     error('dualViterbi:wrongInputArg', 'Cover vector must be of the same length as the code.');\r
  30 end\r
  31 %% initialize datastructures\r
  32 C = zeros(2^code.l,code.n);\r
  33 costs = inf*ones(2^code.l, 1);\r
  34 costs(1) = 0;\r
  35 paths = zeros(2^code.l, code.n);\r
  36 m_id = 1; % message bit id\r
  37 y = zeros(size(x));\r
  38 %% run forward algorithm\r
  39 for i = 1:code.n  % for every bit in cover\r
  40     costs_old = costs;\r
  41     hi = uint32(code.h(i)); % column represented as uint32\r
  42     ji = uint32(1);\r
  43     for j = 1:2^code.l\r
  44         c1 = costs_old(ji) + x(i)*w(i);\r
  45         c2 = costs_old(bitxor(ji-1,hi)+1) + (1-x(i))*w(i);\r
  46         if c1<c2\r
  47             costs(j) = c1;\r
  48             paths(j,i) = ji; % store index of the previous path\r
  49         else\r
  50             costs(j) = c2;\r
  51             paths(j,i) = bitxor(ji-1,hi)+1; % store index of the previous path\r
  52         end\r
  53         ji = ji + 1;\r
  54     end\r
  55     % shift the trellis by shift(i)\r
  56     for j = 1:code.shift(i)\r
  57         if isinf(min(costs( m(m_id)+1:2:end )))\r
  58             error('dualViterbi:noSolution', 'No possible solution exists.');\r
  59         end\r
  60         if m(m_id)==0\r
  61             costs = [costs(1:2:end) ; inf*ones(2^(code.l-1),1)];\r
  62         else\r
  63             costs = [costs(2:2:end) ; inf*ones(2^(code.l-1),1)];\r
  64         end\r
  65         m_id = m_id + 1;\r
  66     end\r
  67     C(:,i) = costs;\r
  68 end\r
  69 %% backward run\r
  70 [min_cost ind] = min(costs);\r
  71 m_id = m_id - 1;\r
  72 for i = code.n:-1:1\r
  73     for j = 1:code.shift(i)\r
  74         ind = 2*ind + m(m_id) - 1; % invert the shift in syndrome trellis\r
  75         m_id = m_id - 1;\r
  76     end\r
  77      y(i) = paths(ind,i)~=ind;\r
  78      ind = paths(ind,i);\r
  79 end\r
  80 end\r