Merge branch 'error-correction'

.gitignore

 autom4te.cache/
 build-aux/
 m4/
-src/*.o
-src/.deps/
-src/audiowmark
 aclocal.m4
 configure
 config.h

src/.gitignore

+*.o
+.deps/
+.libs/
+test/
+audiowmark
+testconvcode

src/Makefile.am

 bin_PROGRAMS = audiowmark
 
-audiowmark_SOURCES = audiowmark.cc wavdata.cc wavdata.hh fft.cc fft.hh
+COMMON_SRC = utils.hh utils.cc convcode.hh convcode.cc
+
+audiowmark_SOURCES = audiowmark.cc wavdata.cc wavdata.hh fft.cc fft.hh $(COMMON_SRC)
 audiowmark_LDFLAGS = $(SNDFILE_LIBS) $(FFTW_LIBS)
+
+noinst_PROGRAMS = testconvcode
+
+testconvcode_SOURCES = testconvcode.cc $(COMMON_SRC)

src/TODO

+possible improvements:
+- dynamic bit strength
+- mp3 support with libmad

src/audiowmark.cc

@@ -6,6 +6,8 @@
 
 #include "fft.hh"
 #include "wavdata.hh"
+#include "utils.hh"
+#include "convcode.hh"
 
 #include <assert.h>
 
@@ -19,15 +21,17 @@ using std::min;
 namespace Params
 {
   static size_t frame_size      = 1024;
-  static int    frames_per_bit  = 4;
+  static int    frames_per_bit  = 2;
   static size_t bands_per_frame = 30;
   static int max_band          = 100;
   static int min_band          = 20;
   static double water_delta    = 0.015; // strength of the watermark
   static double pre_scale      = 0.95;  // rescale the signal to avoid clipping after watermark is added
   static bool mix              = true;
+  static bool hard             = false; // hard decode bits? (soft decoding is better)
   static int block_size        = 32;    // block size for mix step (non-linear bit storage)
   static unsigned int seed     = 0;
+  static size_t payload_size   = 128;  // number of payload bits for the watermark
 }
 
 void
@@ -149,6 +153,10 @@ parse_options (int   *argc_p,
 	{
           Params::mix = false;
 	}
+      else if (check_arg (argc, argv, &i, "--hard"))
+	{
+          Params::hard = true;
+	}
       else if (check_arg (argc, argv, &i, "--seed", &opt_arg))
 	{
           Params::seed = atoi (opt_arg);
@@ -279,56 +287,25 @@ gen_shuffle (vector<T>& result, int seed)
     }
 }
 
-static unsigned char
-from_hex_nibble (char c)
+template<class T> vector<T>
+randomize_bit_order (const vector<T>& bit_vec, bool encode)
 {
-  int uc = (unsigned char)c;
+  vector<unsigned int> order;
 
-  if (uc >= '0' && uc <= '9') return uc - (unsigned char)'0';
-  if (uc >= 'a' && uc <= 'f') return uc + 10 - (unsigned char)'a';
-  if (uc >= 'A' && uc <= 'F') return uc + 10 - (unsigned char)'A';
+  for (size_t i = 0; i < bit_vec.size(); i++)
+    order.push_back (i);
 
-  return 16;	// error
-}
+  gen_shuffle (order, /* seed */ 0);
 
-vector<int>
-bit_str_to_vec (const string& bits)
-{
-  vector<int> bitvec;
-  for (auto nibble : bits)
+  vector<T> out_bits (bit_vec.size());
+  for (size_t i = 0; i < bit_vec.size(); i++)
     {
-      unsigned char c = from_hex_nibble (nibble);
-      if (c >= 16)
-        return vector<int>(); // error
-
-      bitvec.push_back ((c & 8) > 0);
-      bitvec.push_back ((c & 4) > 0);
-      bitvec.push_back ((c & 2) > 0);
-      bitvec.push_back ((c & 1) > 0);
+      if (encode)
+        out_bits[i] = bit_vec[order[i]];
+      else
+        out_bits[order[i]] = bit_vec[i];
     }
-  return bitvec;
-}
-
-string
-bit_vec_to_str (const vector<int>& bit_vec)
-{
-  string bit_str;
-
-  for (size_t pos = 0; pos + 3 < bit_vec.size(); pos += 4) // convert only groups of 4 bits
-    {
-      int nibble = 0;
-      for (int j = 0; j < 4; j++)
-        {
-          if (bit_vec[pos + j])
-            {
-              // j == 0 has the highest value, then 1, 2, 3 (lowest)
-              nibble |= 1 << (3 - j);
-            }
-        }
-      const char *to_hex = "0123456789abcdef";
-      bit_str += to_hex[nibble];
-    }
-  return bit_str;
+  return out_bits;
 }
 
 struct MixEntry
@@ -409,6 +386,21 @@ add_watermark (const string& infile, const string& outfile, const string& bits)
       fprintf (stderr, "audiowmark: cannot parse bits %s\n", bits.c_str());
       return 1;
     }
+  if (bitvec.size() > Params::payload_size)
+    {
+      fprintf (stderr, "audiowmark: number of bits in message '%s' larger than payload size\n", bits.c_str());
+      return 1;
+    }
+  if (bitvec.size() < Params::payload_size)
+    {
+      /* expand message automatically; good for testing, maybe not so good for the final product */
+      vector<int> expanded_bitvec;
+      for (size_t i = 0; i < Params::payload_size; i++)
+        expanded_bitvec.push_back (bitvec[i % bitvec.size()]);
+      bitvec = expanded_bitvec;
+    }
+  /* add forward error correction, bitvec will now be a lot larger */
+  bitvec = randomize_bit_order (conv_encode (bitvec), /* encode */ true);
 
   printf ("loading %s\n", infile.c_str());
 
@@ -520,34 +512,52 @@ add_watermark (const string& infile, const string& outfile, const string& bits)
     }
 
   /* generate synthesis window */
-  vector<float> synth_window (Params::frame_size);
-  for (size_t i = 0; i < Params::frame_size; i++)
+  // we want overlapping synthesis windows, so the window affects the last, the current and the next frame
+  vector<float> synth_window (Params::frame_size * 3);
+  for (size_t i = 0; i < synth_window.size(); i++)
     {
-      const double threshold = 0.2;
+      const double overlap = 0.1;
 
       // triangular basic window
-      const double tri = min (1.0 - fabs (double (2 * i)/Params::frame_size - 1.0), threshold) / threshold;
+      double tri;
+      double norm_pos = (double (i) - Params::frame_size) / Params::frame_size;
 
+      if (norm_pos > 0.5) /* symmetric window */
+        norm_pos = 1 - norm_pos;
+      if (norm_pos < -overlap)
+        {
+          tri = 0;
+        }
+      else if (norm_pos < overlap)
+        {
+          tri = 0.5 + norm_pos / (2 * overlap);
+        }
+      else
+        {
+          tri = 1;
+        }
       // cosine
       synth_window[i] = (cos (tri*M_PI+M_PI)+1) * 0.5;
     }
 
+  for (size_t pos = 0; pos < in_signal.size(); pos++)
+    out_signal[pos] = in_signal[pos] * Params::pre_scale;
+
   for (int f = 0; f < frame_count (wav_data); f++)
     {
       for (int ch = 0; ch < wav_data.n_channels(); ch++)
         {
-          /* add watermark to output frame */
-          vector<float> frame = get_frame (wav_data, f, ch);
-
           /* mix watermark signal to output frame */
           vector<float> fft_delta_out = ifft (fft_delta_spect[f * wav_data.n_channels() + ch]);
 
-          for (size_t i = 0; i < frame.size(); i++)
-            frame[i] += fft_delta_out[i] * synth_window[i];
+          int last_frame_start = (f - 1) * Params::frame_size;
+          for (int i = 0; i < int (synth_window.size()); i++)
+            {
+              int pos = (last_frame_start + i) * wav_data.n_channels() + ch;
 
-          /* modify out signal */
-          for (size_t i = 0; i < frame.size(); i++)
-            out_signal[(f * Params::frame_size + i) * wav_data.n_channels() + ch] = frame[i] * Params::pre_scale;
+              if (pos >= 0 && pos < int (out_signal.size()))
+                out_signal[pos] += fft_delta_out[i % Params::frame_size] * synth_window[i] * Params::pre_scale;
+            }
         }
     }
 
@@ -582,10 +592,37 @@ truncate_to_block_size (WavData& wav_data)
   wav_data.set_samples (in_signal);
 }
 
-vector<int>
+vector<float>
+normalize_soft_bits (const vector<float>& soft_bits)
+{
+  vector<float> norm_soft_bits;
+
+  /* soft decoding produces better error correction than hard decoding */
+  if (Params::hard)
+    {
+      for (auto value : soft_bits)
+        norm_soft_bits.push_back (value > 0 ? 1.0 : 0.0);
+    }
+  else
+    {
+      /* figure out average level of each bit */
+      double mean = 0;
+      for (auto value : soft_bits)
+        mean += fabs (value);
+      mean /= soft_bits.size();
+
+      /* rescale from [-mean,+mean] to [0.0,1.0] */
+      for (auto value : soft_bits)
+        norm_soft_bits.push_back (0.5 * (value / mean + 1));
+    }
+
+  return norm_soft_bits;
+}
+
+vector<float>
 mix_decode (const WavData& wav_data, vector<vector<complex<float>>>& fft_out, vector<vector<complex<float>>>& fft_orig_out)
 {
-  vector<int> bit_vec;
+  vector<float> soft_bit_vec;
 
   for (int block = 0; block < block_count (wav_data); block++)
     {
@@ -617,19 +654,19 @@ mix_decode (const WavData& wav_data, vector<vector<complex<float>>>& fft_out, ve
             }
           if ((f % Params::frames_per_bit) == (Params::frames_per_bit - 1))
             {
-              bit_vec.push_back ((umag > dmag) ? 1 : 0);
+              soft_bit_vec.push_back (umag - dmag);
               umag = 0;
               dmag = 0;
             }
         }
     }
-  return bit_vec;
+  return normalize_soft_bits (soft_bit_vec);
 }
 
-vector<int>
+vector<float>
 linear_decode (const WavData& wav_data, vector<vector<complex<float>>>& fft_out, vector<vector<complex<float>>>& fft_orig_out)
 {
-  vector<int> bit_vec;
+  vector<float> soft_bit_vec;
 
   double umag = 0, dmag = 0;
   for (int f = 0; f < frame_count (wav_data); f++)
@@ -659,39 +696,38 @@ linear_decode (const WavData& wav_data, vector<vector<complex<float>>>& fft_out,
         }
       if ((f % Params::frames_per_bit) == (Params::frames_per_bit - 1))
         {
-          bit_vec.push_back ((umag > dmag) ? 1 : 0);
+          soft_bit_vec.push_back (umag - dmag);
           umag = 0;
           dmag = 0;
         }
     }
-  return bit_vec;
+  return normalize_soft_bits (soft_bit_vec);
 }
 
 int
-get_watermark (const string& infile, const string& orig_pattern)
+decode_and_report (const WavData& wav_data, const string& orig_pattern, vector<vector<complex<float>>>& fft_out, vector<vector<complex<float>>>& fft_orig_out)
 {
-  WavData wav_data;
-  if (!wav_data.load (infile))
-    {
-      fprintf (stderr, "audiowmark: error loading %s: %s\n", infile.c_str(), wav_data.error_blurb());
-      return 1;
-    }
-  vector<int> bit_vec;
-
-  // to keep the watermark detection code simpler, we truncate samples to avoid partial filled blocks
-  truncate_to_block_size (wav_data);
-
-  vector<vector<complex<float>>> fft_out = compute_frame_ffts (wav_data);
-  vector<vector<complex<float>>> fft_orig_out; /* no original data -> blind decode */
-
+  vector<float> soft_bit_vec;
   if (Params::mix)
     {
-      bit_vec = mix_decode (wav_data, fft_out, fft_orig_out);
+      soft_bit_vec = mix_decode (wav_data, fft_out, fft_orig_out);
     }
   else
     {
-      bit_vec = linear_decode (wav_data, fft_out, fft_orig_out);
+      soft_bit_vec = linear_decode (wav_data, fft_out, fft_orig_out);
     }
+  if (soft_bit_vec.size() < conv_code_size (Params::payload_size))
+    {
+      fprintf (stderr, "audiowmark: input file too short to retrieve watermark\n");
+      fprintf (stderr, " - number of recovered raw bits %zd\n", soft_bit_vec.size());
+      fprintf (stderr, " - need at least %zd raw bits to get watermark\n", conv_code_size (Params::payload_size));
+      return 1;
+    }
+  /* truncate to the required length */
+  soft_bit_vec.resize (conv_code_size (Params::payload_size));
+
+  vector<int> bit_vec = conv_decode_soft (randomize_bit_order (soft_bit_vec, /* encode */ false));
+
   printf ("pattern %s\n", bit_vec_to_str (bit_vec).c_str());
   if (!orig_pattern.empty())
     {
@@ -704,11 +740,31 @@ get_watermark (const string& infile, const string& orig_pattern)
           if (bit_vec[i] != orig_vec[i % orig_vec.size()])
             bit_errors++;
         }
+      printf ("bit_error_raw %d %d\n", bit_errors, bits);
       printf ("bit_error_rate %.5f %%\n", double (100.0 * bit_errors) / bits);
     }
   return 0;
 }
 
+int
+get_watermark (const string& infile, const string& orig_pattern)
+{
+  WavData wav_data;
+  if (!wav_data.load (infile))
+    {
+      fprintf (stderr, "audiowmark: error loading %s: %s\n", infile.c_str(), wav_data.error_blurb());
+      return 1;
+    }
+
+  // to keep the watermark detection code simpler, we truncate samples to avoid partial filled blocks
+  truncate_to_block_size (wav_data);
+
+  vector<vector<complex<float>>> fft_out = compute_frame_ffts (wav_data);
+  vector<vector<complex<float>>> fft_orig_out; /* no original data -> blind decode */
+
+  return decode_and_report (wav_data, orig_pattern, fft_out, fft_orig_out);
+}
+
 int
 get_watermark_delta (const string& origfile, const string& infile, const string& orig_pattern)
 {
@@ -733,30 +789,7 @@ get_watermark_delta (const string& origfile, const string& infile, const string&
   vector<vector<complex<float>>> fft_out = compute_frame_ffts (wav_data);
   vector<vector<complex<float>>> fft_orig_out = compute_frame_ffts (orig_wav_data);
 
-  vector<int> bit_vec;
-  if (Params::mix)
-    {
-      bit_vec = mix_decode (wav_data, fft_out, fft_orig_out);
-    }
-  else
-    {
-      bit_vec = linear_decode (wav_data, fft_out, fft_orig_out);
-    }
-  printf ("pattern %s\n", bit_vec_to_str (bit_vec).c_str());
-  if (!orig_pattern.empty())
-    {
-      int bits = 0, bit_errors = 0;
-
-      vector<int> orig_vec = bit_str_to_vec (orig_pattern);
-      for (size_t i = 0; i < bit_vec.size(); i++)
-        {
-          bits++;
-          if (bit_vec[i] != orig_vec[i % orig_vec.size()])
-            bit_errors++;
-        }
-      printf ("bit_error_rate %.5f %%\n", double (100.0 * bit_errors) / bits);
-    }
-  return 0;
+  return decode_and_report (wav_data, orig_pattern, fft_out, fft_orig_out);
 }
 
 int
@@ -773,9 +806,10 @@ gentest (const string& infile, const string& outfile)
   const vector<float>& in_signal = wav_data.samples();
   vector<float> out_signal;
 
-  /* 10 seconds of audio - starting at 30 seconds of the original track */
+  /* 42 seconds of audio - starting at 30 seconds of the original track */
+  /* this is approximately the minimal amount of audio data required for storing a 128-bit encoded message */
   const size_t offset = 30 * wav_data.n_channels() * int (wav_data.mix_freq());
-  const size_t n_samples = 10 * wav_data.n_channels() * int (wav_data.mix_freq());
+  const size_t n_samples = 42 * wav_data.n_channels() * int (wav_data.mix_freq());
   if (in_signal.size() < (offset + n_samples))
     {
       fprintf (stderr, "audiowmark: input file %s too short\n", infile.c_str());

src/ber-test.sh

-# pseudo random pattern
+#!/bin/bash
 
-PATTERN=4e1243bd22c66e76c2ba9eddc1f91394e57f9f83
 TRANSFORM=$1
 if [ "x$AWM_SET" == "x" ]; then
   AWM_SET=small
@@ -8,12 +7,20 @@ fi
 if [ "x$AWM_SEEDS" == "x" ]; then
   AWM_SEEDS=0
 fi
+if [ "x$AWM_REPORT" == "x" ]; then
+  AWM_REPORT=ber
+fi
+if [ "x$AWM_FILE" == "x" ]; then
+  AWM_FILE=t
+fi
 
 {
   if [ "x$AWM_SET" == "xsmall" ]; then
     ls test/T*
   elif [ "x$AWM_SET" == "xbig" ]; then
     cat test_list
+  elif [ "x$AWM_SET" == "xhuge" ]; then
+    ls huge/T*
   else
     echo "bad AWM_SET $AWM_SET" >&2
     exit 1
@@ -23,20 +30,34 @@ do
   for SEED in $AWM_SEEDS
   do
     echo $i
-    audiowmark add "$i" t.wav $PATTERN $AWM_PARAMS --seed $SEED >/dev/null
+
+    if [ "x$AWM_RAND_PATTERN" != "x" ]; then
+      # random pattern, 128 bit
+      PATTERN=$(
+        for i in $(seq 16)
+        do
+          printf "%02x" $((RANDOM % 256))
+        done
+      )
+    else
+      # pseudo random pattern, 128 bit
+      PATTERN=4e1243bd22c66e76c2ba9eddc1f91394
+    fi
+
+    audiowmark add "$i" ${AWM_FILE}.wav $PATTERN $AWM_PARAMS --seed $SEED >/dev/null
     if [ "x$TRANSFORM" == "xmp3" ]; then
       if [ "x$2" == "x" ]; then
         echo "need mp3 bitrate" >&2
         exit 1
       fi
-      lame -b $2 t.wav t.mp3 --quiet
-      rm t.wav
-      ffmpeg -i t.mp3 t.wav -v quiet -nostdin
+      lame -b $2 ${AWM_FILE}.wav ${AWM_FILE}.mp3 --quiet
+      rm ${AWM_FILE}.wav
+      ffmpeg -i ${AWM_FILE}.mp3 ${AWM_FILE}.wav -v quiet -nostdin
 
       # some (low) mpeg quality settings use a lower sample rate
-      if [ "x$(soxi -r t.wav)" != "x44100" ]; then
-        sox t.wav tr.wav rate 44100
-        mv tr.wav t.wav
+      if [ "x$(soxi -r ${AWM_FILE}.wav)" != "x44100" ]; then
+        sox ${AWM_FILE}.wav ${AWM_FILE}r.wav rate 44100
+        mv ${AWM_FILE}r.wav ${AWM_FILE}.wav
       fi
     elif [ "x$TRANSFORM" == "xdouble-mp3" ]; then
       if [ "x$2" == "x" ]; then
@@ -44,27 +65,27 @@ do
         exit 1
       fi
       # first mp3 step (fixed bitrate)
-      lame -b 128 t.wav t.mp3 --quiet
-      rm t.wav
-      ffmpeg -i t.mp3 t.wav -v quiet -nostdin
+      lame -b 128 ${AWM_FILE}.wav ${AWM_FILE}.mp3 --quiet
+      rm ${AWM_FILE}.wav
+      ffmpeg -i ${AWM_FILE}.mp3 ${AWM_FILE}.wav -v quiet -nostdin
 
       # second mp3 step
-      lame -b $2 t.wav t.mp3 --quiet
-      rm t.wav
-      ffmpeg -i t.mp3 t.wav -v quiet -nostdin
+      lame -b $2 ${AWM_FILE}.wav ${AWM_FILE}.mp3 --quiet
+      rm ${AWM_FILE}.wav
+      ffmpeg -i ${AWM_FILE}.mp3 ${AWM_FILE}.wav -v quiet -nostdin
 
       # some (low) mpeg quality settings use a lower sample rate
-      if [ "x$(soxi -r t.wav)" != "x44100" ]; then
-        sox t.wav tr.wav rate 44100
-        mv tr.wav t.wav
+      if [ "x$(soxi -r ${AWM_FILE}.wav)" != "x44100" ]; then
+        sox ${AWM_FILE}.wav ${AWM_FILE}r.wav rate 44100
+        mv ${AWM_FILE}r.wav ${AWM_FILE}.wav
       fi
     elif [ "x$TRANSFORM" == "xogg" ]; then
       if [ "x$2" == "x" ]; then
         echo "need ogg bitrate" >&2
         exit 1
       fi
-      oggenc -b $2 t.wav -o t.ogg --quiet
-      oggdec t.ogg -o t.wav --quiet
+      oggenc -b $2 ${AWM_FILE}.wav -o ${AWM_FILE}.ogg --quiet
+      oggdec ${AWM_FILE}.ogg -o ${AWM_FILE}.wav --quiet
     elif [ "x$TRANSFORM" == "x" ]; then
       :
     else
@@ -72,8 +93,17 @@ do
       exit 1
     fi
     # blind decoding
-    audiowmark cmp t.wav $PATTERN $AWM_PARAMS --seed $SEED
+    audiowmark cmp ${AWM_FILE}.wav $PATTERN $AWM_PARAMS --seed $SEED
     # decoding with original
     # audiowmark cmp-delta "$i" t.wav $PATTERN $AWM_PARAMS --seed $SEED
   done
-done | grep bit_error_rate | awk 'BEGIN { max_er = er = n = 0 } { er += $2; n++; if ($2 > max_er) max_er = $2;} END { print er / n, max_er; }'
+done | grep bit_error_rate | {
+  if [ "x$AWM_REPORT" == "xber" ]; then
+    awk 'BEGIN { max_er = er = n = 0 } { er += $2; n++; if ($2 > max_er) max_er = $2;} END { print er / n, max_er; }'
+  elif [ "x$AWM_REPORT" == "xfer" ]; then
+    awk 'BEGIN { bad = n = 0 } { if ($2 > 0) bad++; n++; } END { print bad, n, bad * 100.0 / n; }'
+  else
+    echo "unknown report $AWM_REPORT" >&2
+    exit 1
+  fi
+}

src/convcode.cc

+#include "utils.hh"
+
+#include <array>
+#include <algorithm>
+#include <assert.h>
+
+using std::vector;
+using std::string;
+using std::min;
+
+int
+parity (unsigned int v)
+{
+  int p = 0;
+
+  while (v)
+    {
+      p ^= (v & 1);
+      v >>= 1;
+    }
+  return p;
+}
+
+// rate 1/6 code generator poynomial from "In search of a 2dB Coding Gain", Yuen and Vo
+// minimum free distance 56
+constexpr  unsigned int rate        = 6;
+constexpr  unsigned int order       = 15;
+constexpr  auto         generators  = std::array<unsigned,6> { 046321, 051271, 070535, 063667, 073277, 076531 };
+
+/*
+constexpr  unsigned int order       = 9;
+constexpr  auto         generators  = std::array<unsigned,3> { 0557, 0663, 0711 };
+
+constexpr  unsigned int order       = 14;
+constexpr  auto         generators  = std::array<unsigned,3> { 021645, 035661, 037133 };
+
+constexpr  unsigned int order       = 18;
+constexpr  auto         generators  = std::array<unsigned,3> { 0552137, 0614671, 0772233 };
+*/
+
+constexpr  unsigned int state_count = (1 << order);
+constexpr  unsigned int state_mask  = (1 << order) - 1;
+
+size_t
+conv_code_size (size_t msg_size)
+{
+  return (msg_size + order) * rate;
+}
+
+vector<int>
+conv_encode (const vector<int>& in_bits)
+{
+  vector<int> out_vec;
+  vector<int> vec = in_bits;
+
+  /* termination: bring encoder into all-zero state */
+  for (unsigned int i = 0; i < order; i++)
+    vec.push_back (0);
+
+  unsigned int reg = 0;
+  for (auto b : vec)
+    {
+      reg = (reg << 1) | b;
+
+      for (auto poly : generators)
+        {
+          int out_bit = parity (reg & poly);
+
+          out_vec.push_back (out_bit);
+        }
+    }
+  return out_vec;
+}
+
+/* decode using viterbi algorithm */
+vector<int>
+conv_decode_soft (const vector<float>& coded_bits)
+{
+  vector<int> decoded_bits;
+
+  assert (coded_bits.size() % rate == 0);
+
+  struct StateEntry
+  {
+    int   last_state;
+    float delta;
+    int   bit;
+  };
+  vector<vector<StateEntry>> error_count;
+  for (size_t i = 0; i < coded_bits.size() + rate; i += rate) /* 1 extra element */
+    error_count.emplace_back (state_count, StateEntry {0, -1, 0});
+
+  error_count[0][0].delta = 0; /* start state */
+
+  /* precompute state -> output bits table */
+  vector<float> state2bits;
+  for (unsigned int state = 0; state < state_count; state++)
+    {
+      for (size_t p = 0; p < generators.size(); p++)
+        {
+          int out_bit = parity (state & generators[p]);
+          state2bits.push_back (out_bit);
+        }
+    }
+
+  for (size_t i = 0; i < coded_bits.size(); i += rate)
+    {
+      vector<StateEntry>& old_table = error_count[i / rate];
+      vector<StateEntry>& new_table = error_count[i / rate + 1];
+
+      for (unsigned int state = 0; state < state_count; state++)
+        {
+          /* this check enforces that we only consider states reachable from state=0 at time=0*/
+          if (old_table[state].delta >= 0)
+            {
+              for (int bit = 0; bit < 2; bit++)
+                {
+                  int   new_state = ((state << 1) | bit) & state_mask;
+
+                  float delta = old_table[state].delta;
+                  int   sbit_pos = new_state * rate;
+
+                  for (size_t p = 0; p < generators.size(); p++)
+                    {
+                      const float cbit = coded_bits[i + p];
+                      const float sbit = state2bits[sbit_pos + p];
+
+                      /* decoding error weight for this bit; if input is only 0.0 and 1.0, this is the hamming distance */
+                      delta += (cbit - sbit) * (cbit - sbit);
+                    }
+
+                  if (delta < new_table[new_state].delta || new_table[new_state].delta < 0) /* better match with this link? replace entry */
+                    {
+                      new_table[new_state].delta      = delta;
+                      new_table[new_state].last_state = state;
+                      new_table[new_state].bit        = bit;
+                    }
+                }
+            }
+        }
+    }
+
+  unsigned int state = 0;
+  for (size_t idx = error_count.size() - 1; idx > 0; idx--)
+    {
+      decoded_bits.push_back (error_count[idx][state].bit);
+
+      state = error_count[idx][state].last_state;
+    }
+  std::reverse (decoded_bits.begin(), decoded_bits.end());
+
+  /* remove termination */
+  assert (decoded_bits.size() >= order);
+  decoded_bits.resize (decoded_bits.size() - order);
+
+  return decoded_bits;
+}
+
+vector<int>
+conv_decode_hard (const vector<int>& coded_bits)
+{
+  /* for the final application, we always want soft decoding, so we don't
+   * special case hard decoding here, so this will be a little slower than
+   * necessary
+   */
+  vector<float> soft_bits;
+  for (auto b : coded_bits)
+    soft_bits.push_back (b ? 1.0f : 0.0f);
+
+  return conv_decode_soft (soft_bits);
+}

src/convcode.hh

+#ifndef AUDIOWMARK_CONV_CODE_HH
+#define AUDIOWMARK_CONV_CODE_HH
+
+#include <vector>
+#include <string>
+
+size_t           conv_code_size (size_t msg_size);
+std::vector<int> conv_encode (const std::vector<int>& in_bits);
+std::vector<int> conv_decode_hard (const std::vector<int>& coded_bits);
+std::vector<int> conv_decode_soft (const std::vector<float>& coded_bits);
+
+#endif /* AUDIOWMARK_CONV_CODE_HH */

src/fer-test.sh

+#!/bin/bash
+SEEDS="$1"
+MAX_SEED=$(($SEEDS - 1))
+P="$2"
+shift 2
+echo "n seeds       : $SEEDS"
+echo "ber-test args : $@"
+echo "params        : $P"
+for seed in $(seq 0 $MAX_SEED)
+do
+  echo $(AWM_SEEDS=$seed AWM_PARAMS="$P" AWM_REPORT="fer" ber-test.sh "$@")
+done | awk '{bad += $1; files += $2; print bad, files, bad * 100. / files }'

src/gen-fer-adoc.sh

+#!/bin/bash
+
+for TEST in mp3 double-mp3 ogg
+do
+  echo ".$TEST"
+  echo '[frame="topbot",options="header",cols="12*>"]'
+  echo '|=========================='
+  echo -n "| "
+  for D in $(seq 15 -1 5)
+  do
+    DELTA=$(printf "0.0%02d\n" $D)
+    echo -n "| $DELTA"
+  done
+  echo
+  for BITRATE in 512 256 196 128 96 64
+  do
+    echo -n "| $BITRATE" | sed 's/512/wav/g'
+    for D in $(seq 15 -1 5)
+    do
+      DELTA=$(printf "0.0%02d\n" $D)
+      cat $DELTA-$TEST-* | grep ^$BITRATE | awk '{bad += $2; n += $3} END {fer=100.0*bad/n; bold=fer>0?"*":" ";printf ("| %s%.2f%s", bold, fer, bold)}'
+    done
+    echo
+  done
+  echo '|=========================='
+  echo
+done

src/gen-fer-mk.sh

+#!/bin/bash
+
+DELTA_RANGE="0.005 0.006 0.007 0.008 0.009 0.010 0.011 0.012 0.013 0.014 0.015"
+SEEDS="$(seq 0 19)"
+
+echo -n "all: "
+for SEED in $SEEDS
+do
+  for DELTA in $DELTA_RANGE
+  do
+    echo -n "$DELTA-ogg-$SEED $DELTA-mp3-$SEED $DELTA-double-mp3-$SEED "
+  done
+done
+
+echo
+echo
+
+for SEED in $SEEDS
+do
+  for DELTA in $DELTA_RANGE
+  do
+    FILE="$DELTA-ogg-$SEED"
+    echo "$FILE:"
+    echo -e "\t( cd ..; AWM_SET=huge AWM_PARAMS='--water-delta $DELTA' AWM_REPORT=fer AWM_SEEDS=$SEED AWM_FILE='t-$FILE' ber-ogg.sh ) >x$FILE"
+    echo -e "\tmv x$FILE $FILE"
+    echo
+
+    FILE="$DELTA-mp3-$SEED"
+    echo "$FILE:"
+    echo -e "\t( cd ..; AWM_SET=huge AWM_PARAMS='--water-delta $DELTA' AWM_REPORT=fer AWM_SEEDS=$SEED AWM_FILE='t-$FILE' ber-mp3.sh ) >x$FILE"
+    echo -e "\tmv x$FILE $FILE"
+    echo
+
+    FILE="$DELTA-double-mp3-$SEED"
+    echo "$FILE:"
+    echo -e "\t( cd ..; AWM_SET=huge AWM_PARAMS='--water-delta $DELTA' AWM_REPORT=fer AWM_SEEDS=$SEED AWM_FILE='t-$FILE' ber-double-mp3.sh ) >x$FILE"
+    echo -e "\tmv x$FILE $FILE"
+    echo
+  done
+done

src/seed-test.sh

-MAX_SEED=$1
+SEEDS="$1"
+MAX_SEED=$(($SEEDS - 1))
 P1="$2"
 P2="$3"
 shift 3
-echo "max seed      : $MAX_SEED"
+echo "n seeds       : $SEEDS"
 echo "ber-test args : $@"
 echo "left          : $P1"
 echo "right         : $P2"

src/testconvcode.cc

+#include "utils.hh"
+#include "convcode.hh"
+
+#include <random>
+
+#include <assert.h>
+#include <sys/time.h>
+
+using std::vector;
+using std::string;
+
+static double
+gettime()
+{
+  timeval tv;
+  gettimeofday (&tv, 0);
+
+  return tv.tv_sec + tv.tv_usec / 1000000.0;
+}
+
+vector<int>
+generate_error_vector (size_t n, int errors)
+{
+  vector<int> ev (n);
+
+  while (errors)
+    {
+      size_t pos = rand() % ev.size();
+      if (ev[pos] != 1)
+        {
+          ev[pos] = 1;
+          errors--;
+        }
+    }
+  return ev;
+}
+int
+main (int argc, char **argv)
+{
+  if (argc == 1)
+    {
+      vector<int> in_bits = bit_str_to_vec ("80f12381");
+
+      printf ("input vector (k=%zd):  ", in_bits.size());
+      for (auto b : in_bits)
+        printf ("%d", b);
+      printf ("\n");
+
+      vector<int> coded_bits = conv_encode (in_bits);
+      printf ("coded vector (n=%zd): ", coded_bits.size());
+      for (auto b : coded_bits)
+        printf ("%d", b);
+      printf ("\n");
+      printf ("coded hex: %s\n", bit_vec_to_str (coded_bits).c_str());
+
+      assert (coded_bits.size() == conv_code_size (in_bits.size()));
+
+      vector<int> decoded_bits = conv_decode_hard (coded_bits);
+      printf ("output vector (k=%zd): ", decoded_bits.size());
+      for (auto b : decoded_bits)
+        printf ("%d", b);
+      printf ("\n");
+
+      assert (decoded_bits.size() == in_bits.size());
+      int errors = 0;
+      for (size_t i = 0; i < decoded_bits.size(); i++)
+        if (decoded_bits[i] != in_bits[i])
+          errors++;
+      printf ("decoding errors: %d\n", errors);
+    }
+  if (argc == 2 && string (argv[1]) == "error")
+    {
+      size_t max_bit_errors = conv_code_size (128) * 0.5;
+
+      for (size_t bit_errors = 0; bit_errors < max_bit_errors; bit_errors++)
+        {
+          size_t coded_bit_count = 0;
+          int bad_decode = 0;
+          constexpr int test_size = 20;
+
+          for (int i = 0; i < test_size; i++)
+            {
+              vector<int> in_bits;
+              while (in_bits.size() != 128)
+                in_bits.push_back (rand() & 1);
+
+              vector<int> coded_bits = conv_encode (in_bits);
+              coded_bit_count = coded_bits.size();
+
+              vector<int> error_bits = generate_error_vector (coded_bits.size(), bit_errors);
+              for (size_t pos = 0; pos < coded_bits.size(); pos++)
+                coded_bits[pos] ^= error_bits[pos];
+
+              vector<int> decoded_bits = conv_decode_hard (coded_bits);
+
+              assert (decoded_bits.size() == 128);
+
+              int errors = 0;
+              for (size_t i = 0; i < 128; i++)
+                if (decoded_bits[i] != in_bits[i])
+                  errors++;
+              if (errors > 0)
+                bad_decode++;
+            }
+          printf ("%f %f\n", (100.0 * bit_errors) / coded_bit_count, (100.0 * bad_decode) / test_size);
+        }
+    }
+  if (argc == 2 && string (argv[1]) == "soft-error")
+    {
+      for (double stddev = 0; stddev < 1.5; stddev += 0.01)
+        {
+          size_t coded_bit_count = 0;
+          int bad_decode1 = 0, bad_decode2 = 0;
+          constexpr int test_size = 20;
+
+          int local_be = 0;
+          for (int i = 0; i < test_size; i++)
+            {
+              vector<int> in_bits;
+              while (in_bits.size() != 128)
+                in_bits.push_back (rand() & 1);
+
+              vector<int> coded_bits = conv_encode (in_bits);
+              coded_bit_count = coded_bits.size();
+
+              std::default_random_engine generator;
+              std::normal_distribution<double> dist (0, stddev);
+
+              vector<float> recv_bits;
+              for (auto b : coded_bits)
+                recv_bits.push_back (b + dist (generator));
+
+              vector<int> decoded_bits1 = conv_decode_soft (recv_bits);
+
+              vector<int> recv_hard_bits;
+              for (auto b : recv_bits)
+                recv_hard_bits.push_back ((b > 0.5) ? 1 : 0);
+
+              for (size_t x = 0; x < recv_hard_bits.size(); x++)
+                local_be += coded_bits[x] ^ recv_hard_bits[x];
+
+              vector<int> decoded_bits2 = conv_decode_hard (recv_hard_bits);
+
+              assert (decoded_bits1.size() == 128);
+              assert (decoded_bits2.size() == 128);
+
+              int e1 = 0;
+              int e2 = 0;
+              for (size_t i = 0; i < 128; i++)
+                {
+                  if (decoded_bits1[i] != in_bits[i])
+                    e1++;
+                  if (decoded_bits2[i] != in_bits[i])
+                    e2++;
+                }
+              if (e1)
+                bad_decode1++;
+              if (e2)
+                bad_decode2++;
+            }
+          printf ("%f %f %f\n", double (100 * local_be) / test_size / coded_bit_count, (100.0 * bad_decode1) / test_size, (100.0 * bad_decode2) / test_size);
+        }
+    }
+  if (argc == 2 && string (argv[1]) == "perf")
+    {
+      vector<int> in_bits;
+      while (in_bits.size() != 128)
+        in_bits.push_back (rand() & 1);
+
+      const double start_t = gettime();
+      const size_t runs = 20;
+      for (size_t i = 0; i < runs; i++)
+        {
+          vector<int> out_bits = conv_decode_hard (conv_encode (in_bits));
+          assert (out_bits == in_bits);
+        }
+      printf ("%.1f ms/block\n", (gettime() - start_t) / runs * 1000.0);
+    }
+}

src/utils.cc

+#include "utils.hh"
+
+using std::vector;
+using std::string;
+
+static unsigned char
+from_hex_nibble (char c)
+{
+  int uc = (unsigned char)c;
+
+  if (uc >= '0' && uc <= '9') return uc - (unsigned char)'0';
+  if (uc >= 'a' && uc <= 'f') return uc + 10 - (unsigned char)'a';
+  if (uc >= 'A' && uc <= 'F') return uc + 10 - (unsigned char)'A';
+
+  return 16;	// error
+}
+
+vector<int>
+bit_str_to_vec (const string& bits)
+{
+  vector<int> bitvec;
+  for (auto nibble : bits)
+    {
+      unsigned char c = from_hex_nibble (nibble);
+      if (c >= 16)
+        return vector<int>(); // error
+
+      bitvec.push_back ((c & 8) > 0);
+      bitvec.push_back ((c & 4) > 0);
+      bitvec.push_back ((c & 2) > 0);
+      bitvec.push_back ((c & 1) > 0);
+    }
+  return bitvec;
+}
+
+string
+bit_vec_to_str (const vector<int>& bit_vec)
+{
+  string bit_str;
+
+  for (size_t pos = 0; pos + 3 < bit_vec.size(); pos += 4) // convert only groups of 4 bits
+    {
+      int nibble = 0;
+      for (int j = 0; j < 4; j++)
+        {
+          if (bit_vec[pos + j])
+            {
+              // j == 0 has the highest value, then 1, 2, 3 (lowest)
+              nibble |= 1 << (3 - j);
+            }
+        }
+      const char *to_hex = "0123456789abcdef";
+      bit_str += to_hex[nibble];
+    }
+  return bit_str;
+}
+

src/utils.hh

+#ifndef AUDIOWMARK_UTILS_HH
+#define AUDIOWMARK_UTILS_HH
+
+#include <vector>
+#include <string>
+
+std::vector<int> bit_str_to_vec (const std::string& bits);
+std::string      bit_vec_to_str (const std::vector<int>& bit_vec);
+
+#endif /* AUDIOWMARK_UTILS_HH */