core/conv: Add some generic code for convolutional coding/decoding

Far from perfect but suits our need thus far.

The viterbi with softbit input is quite cpu-intensive. Since
most received bursts are often mostly error free, you could
use a less cpu intensive algorithm (Fano ?) and with hard bit
input. Then only switch to viterbi soft bit input if the channel
is bad enough to justify it.

Soft output is not implemented as its usefulness for the block
coding is limited.

Signed-off-by: Sylvain Munaut <tnt@246tNt.com>

4 files changed, 596 insertions, 2 deletions

diff --git a/include/osmocom/core/Makefile.am b/include/osmocom/core/Makefile.am
index 6109f478..b21e0476 100644
--- a/include/osmocom/core/Makefile.am
+++ b/include/osmocom/core/Makefile.am
@@ -3,7 +3,7 @@ osmocore_HEADERS = signal.h linuxlist.h timer.h select.h msgb.h bits.h \
 		   gsmtap.h write_queue.h \
 		   logging.h rate_ctr.h gsmtap_util.h \
 		   plugin.h crc16.h panic.h process.h msgfile.h \
-		   backtrace.h
+		   backtrace.h conv.h
 
 if ENABLE_TALLOC
 osmocore_HEADERS += talloc.h
diff --git a/include/osmocom/core/conv.h b/include/osmocom/core/conv.h
new file mode 100644
index 00000000..af676eed
--- /dev/null
+++ b/include/osmocom/core/conv.h
@@ -0,0 +1,101 @@
+/*
+ * conv.h
+ *
+ * Copyright (C) 2011  Sylvain Munaut <tnt@246tNt.com>
+ *
+ * All Rights Reserved
+ *
+ * This program is free software; you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation; either version 2 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License along
+ * with this program; if not, write to the Free Software Foundation, Inc.,
+ * 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+ */
+
+#ifndef __OSMO_CONV_H__
+#define __OSMO_CONV_H__
+
+#include <stdint.h>
+
+#include <osmocom/core/bits.h>
+
+struct osmo_conv_code {
+	int N;
+	int K;
+	int len;
+
+	const uint8_t (*next_output)[2];
+	const uint8_t (*next_state)[2];
+
+	const uint8_t *next_term_output;
+	const uint8_t *next_term_state;
+
+	const int *puncture;
+};
+
+
+/* Encoding */
+
+	/* Low level API */
+struct osmo_conv_encoder {
+	const struct osmo_conv_code *code;
+	int i_idx;	/* Next input bit index */
+	int p_idx;	/* Current puncture index */
+	uint8_t state;	/* Current state */
+};
+
+void osmo_conv_encode_init(struct osmo_conv_encoder *encoder,
+                           const struct osmo_conv_code *code);
+int  osmo_conv_encode_raw(struct osmo_conv_encoder *encoder,
+                          const ubit_t *input, ubit_t *output, int n);
+int  osmo_conv_encode_finish(struct osmo_conv_encoder *encoder, ubit_t *output);
+
+	/* All-in-one */
+int  osmo_conv_encode(const struct osmo_conv_code *code,
+                      const ubit_t *input, ubit_t *output);
+
+
+/* Decoding */
+
+	/* Low level API */
+struct osmo_conv_decoder {
+	const struct osmo_conv_code *code;
+
+	int n_states;
+
+	int len;		/* Max o_idx (excl. termination) */
+
+	int o_idx;		/* output index */
+	int p_idx;		/* puncture index */
+
+	unsigned int *ae;	/* accumulater error */
+	unsigned int *ae_next;	/* next accumulated error (tmp in scan) */
+	uint8_t *state_history;	/* state history [len][n_states] */
+};
+
+void osmo_conv_decode_init(struct osmo_conv_decoder *decoder,
+                           const struct osmo_conv_code *code, int len);
+void osmo_conv_decode_reset(struct osmo_conv_decoder *decoder);
+void osmo_conv_decode_deinit(struct osmo_conv_decoder *decoder);
+
+int osmo_conv_decode_scan(struct osmo_conv_decoder *decoder,
+                          const sbit_t *input, int n);
+int osmo_conv_decode_finish(struct osmo_conv_decoder *decoder,
+                            const sbit_t *input);
+int osmo_conv_decode_get_output(struct osmo_conv_decoder *decoder,
+                                ubit_t *output, int has_finish);
+
+	/* All-in-one */
+int osmo_conv_decode(const struct osmo_conv_code *code,
+                     const sbit_t *input, ubit_t *output);
+
+
+#endif /* __OSMO_CONV_H__ */
diff --git a/src/Makefile.am b/src/Makefile.am
index c5c8a21d..a0d639ca 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -14,7 +14,7 @@ libosmocore_la_SOURCES = timer.c select.c signal.c msgb.c bits.c \
 			 write_queue.c utils.c \
 			 logging.c logging_syslog.c rate_ctr.c \
 			 gsmtap_util.c crc16.c panic.c backtrace.c \
-			 process.c
+			 process.c conv.c
 
 if ENABLE_PLUGIN
 libosmocore_la_SOURCES += plugin.c
diff --git a/src/conv.c b/src/conv.c
new file mode 100644
index 00000000..cf820a33
--- /dev/null
+++ b/src/conv.c
@@ -0,0 +1,493 @@
+/*
+ * conv.c
+ *
+ * Generic convolutional encoding / decoding
+ *
+ * Copyright (C) 2011  Sylvain Munaut <tnt@246tNt.com>
+ *
+ * All Rights Reserved
+ *
+ * This program is free software; you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation; either version 2 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License along
+ * with this program; if not, write to the Free Software Foundation, Inc.,
+ * 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+ */
+
+#include <alloca.h>
+#include <stdint.h>
+#include <stdlib.h>
+#include <string.h>
+
+#include <osmocom/core/bits.h>
+#include <osmocom/core/conv.h>
+
+
+/* ------------------------------------------------------------------------ */
+/* Encoding                                                                 */
+/* ------------------------------------------------------------------------ */
+
+void
+osmo_conv_encode_init(struct osmo_conv_encoder *encoder,
+                      const struct osmo_conv_code *code)
+{
+	memset(encoder, 0x00, sizeof(struct osmo_conv_encoder));
+	encoder->code = code;
+}
+
+static inline int
+_conv_encode_do_output(struct osmo_conv_encoder *encoder,
+                       uint8_t out, ubit_t *output)
+{
+	const struct osmo_conv_code *code = encoder->code;
+	int o_idx = 0;
+	int j;
+
+	if (code->puncture) {
+		for (j=0; j<code->N; j++)
+		{
+			int bit_no = code->N - j - 1;
+			int r_idx = encoder->i_idx * code->N + j;
+
+			if (code->puncture[encoder->p_idx] == r_idx)
+				encoder->p_idx++;
+			else
+				output[o_idx++] = (out >> bit_no) & 1;
+		}
+	} else {
+		for (j=0; j<code->N; j++)
+		{
+			int bit_no = code->N - j - 1;
+			output[o_idx++] = (out >> bit_no) & 1;
+		}
+	}
+
+	return o_idx;
+}
+
+int
+osmo_conv_encode_raw(struct osmo_conv_encoder *encoder,
+                     const ubit_t *input, ubit_t *output, int n)
+{
+	const struct osmo_conv_code *code = encoder->code;
+	uint8_t state;
+	int i;
+	int o_idx;
+
+	o_idx = 0;
+	state = encoder->state;
+
+	for (i=0; i<n; i++) {
+		int bit = input[i];
+		uint8_t out;
+
+		out   = code->next_output[state][bit];
+		state = code->next_state[state][bit];
+
+		o_idx += _conv_encode_do_output(encoder, out, &output[o_idx]);
+
+		encoder->i_idx++;
+	}
+
+	encoder->state = state;
+
+	return o_idx;
+}
+
+int
+osmo_conv_encode_finish(struct osmo_conv_encoder *encoder,
+                        ubit_t *output)
+{
+	const struct osmo_conv_code *code = encoder->code;
+	uint8_t state;
+	int n;
+	int i;
+	int o_idx;
+
+	n = code->K - 1;
+
+	o_idx = 0;
+	state = encoder->state;
+
+	for (i=0; i<n; i++) {
+		uint8_t out;
+
+		if (code->next_term_output) {
+			out   = code->next_term_output[state];
+			state = code->next_term_state[state];
+		} else {
+			out   = code->next_output[state][0];
+			state = code->next_state[state][0];
+		}
+
+		o_idx += _conv_encode_do_output(encoder, out, &output[o_idx]);
+
+		encoder->i_idx++;
+	}
+
+	encoder->state = state;
+
+	return o_idx;
+}
+
+int
+osmo_conv_encode(const struct osmo_conv_code *code,
+                 const ubit_t *input, ubit_t *output)
+{
+	struct osmo_conv_encoder encoder;
+	int l;
+
+	osmo_conv_encode_init(&encoder, code);
+	l  = osmo_conv_encode_raw(&encoder, input, output, code->len);
+	l += osmo_conv_encode_finish(&encoder, &output[l]);
+
+	return l;
+}
+
+
+/* ------------------------------------------------------------------------ */
+/* Decoding (viterbi)                                                       */
+/* ------------------------------------------------------------------------ */
+
+#define MAX_AE 0x00ffffff
+
+void
+osmo_conv_decode_init(struct osmo_conv_decoder *decoder,
+                      const struct osmo_conv_code *code, int len)
+{
+	int n_states;
+
+	/* Init */
+	if (len <= 0)
+		len =  code->len;
+
+	n_states = 1 << (code->K - 1);
+
+	memset(decoder, 0x00, sizeof(struct osmo_conv_decoder));
+
+	decoder->code = code;
+	decoder->n_states = n_states;
+	decoder->len = len;
+
+	/* Allocate arrays */
+	decoder->ae      = malloc(sizeof(unsigned int) * n_states);
+	decoder->ae_next = malloc(sizeof(unsigned int) * n_states);
+
+	decoder->state_history = malloc(sizeof(uint8_t) * n_states * (len + decoder->code->K - 1));
+
+	/* Classic reset */
+	osmo_conv_decode_reset(decoder);
+}
+
+void
+osmo_conv_decode_reset(struct osmo_conv_decoder *decoder)
+{
+	int i;
+
+	/* Reset indexes */
+	decoder->o_idx = 0;
+	decoder->p_idx = 0;
+
+	/* Initial error (only state 0 is valid) */
+	decoder->ae[0] = 0;
+	for (i=1; i<decoder->n_states; i++) {
+		decoder->ae[i] = MAX_AE;
+	}
+}
+
+void
+osmo_conv_decode_deinit(struct osmo_conv_decoder *decoder)
+{
+	free(decoder->ae);
+	free(decoder->ae_next);
+	free(decoder->state_history);
+
+	memset(decoder, 0x00, sizeof(struct osmo_conv_decoder));
+}
+
+int
+osmo_conv_decode_scan(struct osmo_conv_decoder *decoder,
+                      const sbit_t *input, int n)
+{
+	const struct osmo_conv_code *code = decoder->code;
+
+	int i, s, b, j;
+
+	int n_states;
+	unsigned int *ae;
+	unsigned int *ae_next;
+	uint8_t *state_history;
+	sbit_t *in_sym;
+
+	int i_idx, p_idx;
+
+	/* Prepare */
+	n_states = decoder->n_states;
+
+	ae      = decoder->ae;
+	ae_next = decoder->ae_next;
+	state_history = &decoder->state_history[n_states * decoder->o_idx];
+
+	in_sym  = alloca(sizeof(sbit_t) * code->N);
+
+	i_idx = 0;
+	p_idx = decoder->p_idx;
+
+	/* Scan the treillis */
+	for (i=0; i<n; i++)
+	{
+		/* Reset next accumulated error */
+		for (s=0; s<n_states; s++) {
+			ae_next[s] = MAX_AE;
+		}
+
+		/* Get input */
+		if (code->puncture) {
+			/* Hard way ... */
+			for (j=0; j<code->N; j++) {
+				int idx = ((decoder->o_idx + i) * code->N) + j;
+				if (idx == code->puncture[p_idx]) {
+					in_sym[j] = 0;	/* Undefined */
+					p_idx++;
+				} else {
+					in_sym[j] = input[i_idx];
+					i_idx++;
+				}
+			}
+		} else {
+			/* Easy, just copy N bits */
+			memcpy(in_sym, &input[i_idx], code->N);
+			i_idx += code->N;
+		}
+
+		/* Scan all state */
+		for (s=0; s<n_states; s++)
+		{
+			/* Scan possible input bits */
+			for (b=0; b<2; b++)
+			{
+				int nae, ov, e;
+				uint8_t m;
+
+				/* Next output and state */
+				uint8_t out   = code->next_output[s][b];
+				uint8_t state = code->next_state[s][b];
+
+				/* New error for this path */
+				nae = ae[s];			/* start from last error */
+				m = 1 << (code->N - 1);		/* mask for 'out' bit selection */
+
+				for (j=0; j<code->N; j++) {
+					ov = (out & m) ? -127 : 127; /* sbit_t value for it */
+					e = ((int)in_sym[j]) - ov;   /* raw error for this bit */
+					nae += (e * e) >> 9;         /* acc the squared/scaled value */
+					m >>= 1;                     /* next mask bit */
+				}
+
+				/* Is it survivor ? */
+				if (ae_next[state] > nae) {
+					ae_next[state] = nae;
+					state_history[(n_states * i) + state] = s;
+				}
+			}
+		}
+
+		/* Copy accumulated error */
+		memcpy(ae, ae_next, sizeof(unsigned int) * n_states);
+	}
+
+	/* Update decoder state */
+	decoder->p_idx = p_idx;
+	decoder->o_idx += n;
+
+	return i_idx;
+}
+
+int
+osmo_conv_decode_finish(struct osmo_conv_decoder *decoder,
+                        const sbit_t *input)
+{
+	const struct osmo_conv_code *code = decoder->code;
+
+	int i, s, j;
+
+	int n_states;
+	unsigned int *ae;
+	unsigned int *ae_next;
+	uint8_t *state_history;
+	sbit_t *in_sym;
+
+	int i_idx, p_idx;
+
+	/* Prepare */
+	n_states = decoder->n_states;
+
+	ae      = decoder->ae;
+	ae_next = decoder->ae_next;
+	state_history = &decoder->state_history[n_states * decoder->o_idx];
+
+	in_sym  = alloca(sizeof(sbit_t) * code->N);
+
+	i_idx = 0;
+	p_idx = decoder->p_idx;
+
+	/* Scan the treillis */
+	for (i=0; i<code->K-1; i++)
+	{
+		/* Reset next accumulated error */
+		for (s=0; s<n_states; s++) {
+			ae_next[s] = MAX_AE;
+		}
+
+		/* Get input */
+		if (code->puncture) {
+			/* Hard way ... */
+			for (j=0; j<code->N; j++) {
+				int idx = ((decoder->o_idx + i) * code->N) + j;
+				if (idx == code->puncture[p_idx]) {
+					in_sym[j] = 0;	/* Undefined */
+					p_idx++;
+				} else {
+					in_sym[j] = input[i_idx];
+					i_idx++;
+				}
+			}
+		} else {
+			/* Easy, just copy N bits */
+			memcpy(in_sym, &input[i_idx], code->N);
+			i_idx += code->N;
+		}
+
+		/* Scan all state */
+		for (s=0; s<n_states; s++)
+		{
+			int nae, ov, e;
+			uint8_t m;
+
+			/* Next output and state */
+			uint8_t out;
+			uint8_t state;
+
+			if (code->next_term_output) {
+				out   = code->next_term_output[s];
+				state = code->next_term_state[s];
+			} else {
+				out   = code->next_output[s][0];
+				state = code->next_state[s][0];
+			}
+
+			/* New error for this path */
+			nae = ae[s];			/* start from last error */
+			m = 1 << (code->N - 1);		/* mask for 'out' bit selection */
+
+			for (j=0; j<code->N; j++) {
+				ov = (out & m) ? -127 : 127; /* sbit_t value for it */
+				e = ((int)in_sym[j]) - ov;   /* raw error for this bit */
+				nae += (e * e) >> 9;         /* acc the squared/scaled value */
+				m >>= 1;                     /* next mask bit */
+			}
+
+			/* Is it survivor ? */
+			if (ae_next[state] > nae) {
+				ae_next[state] = nae;
+				state_history[(n_states * i) + state] = s;
+			}
+		}
+
+		/* Copy accumulated error */
+		memcpy(ae, ae_next, sizeof(unsigned int) * n_states);
+	}
+
+	/* Update decoder state */
+	decoder->p_idx = p_idx;
+	decoder->o_idx += code->K - 1;
+
+	return i_idx;
+}
+
+int
+osmo_conv_decode_get_output(struct osmo_conv_decoder *decoder,
+                            ubit_t *output, int has_finish)
+{
+	const struct osmo_conv_code *code = decoder->code;
+
+	int min_ae;
+	uint8_t min_state, cur_state;
+	int i, s, n;
+
+	uint8_t *sh_ptr;
+
+	/* Find state with least error */
+	min_ae = MAX_AE;
+	min_state = 0xff;
+
+	for (s=0; s<decoder->n_states; s++)
+	{
+		if (decoder->ae[s] < min_ae) {
+			min_ae = decoder->ae[s];
+			min_state = s;
+		}
+	}
+
+	if (min_state == 0xff)
+		return -1;
+
+	/* Traceback */
+	cur_state = min_state;
+
+	n = decoder->o_idx;
+
+	sh_ptr = &decoder->state_history[decoder->n_states * (n-1)];
+
+		/* No output for the K-1 termination input bits */
+	if (has_finish) {
+		for (i=0; i<code->K-1; i++) {
+			cur_state = sh_ptr[cur_state];
+			sh_ptr -= decoder->n_states;
+		}
+		n -= code->K - 1;
+	}
+
+		/* Generate output backward */
+	for (i=n-1; i>=0; i--)
+	{
+		min_state = cur_state;
+		cur_state = sh_ptr[cur_state];
+
+		sh_ptr -= decoder->n_states;
+
+		if (code->next_state[cur_state][0] == min_state)
+			output[i] = 0;
+		else
+			output[i] = 1;
+	}
+
+	return min_ae;
+}
+
+int
+osmo_conv_decode(const struct osmo_conv_code *code,
+                 const sbit_t *input, ubit_t *output)
+{
+	struct osmo_conv_decoder decoder;
+	int rv, l;
+
+	osmo_conv_decode_init(&decoder, code, 0);
+
+	l = osmo_conv_decode_scan(&decoder, input, code->len);
+	l = osmo_conv_decode_finish(&decoder, &input[l]);
+
+	rv = osmo_conv_decode_get_output(&decoder, output, 1);
+
+	osmo_conv_decode_deinit(&decoder);
+
+	return rv;
+}