Hamming coding

hamming.cpp

/*
 * Author: Alexander Tarasikov <[email protected]>
 * Date: 2013-03-28
 *
 * Variant 5
 * Hamming (15, 11) code
 *
 */

#include <iostream>
#include <memory>
#include <cassert>

#define UNUSED(x) ((void)x)

typedef unsigned bit_t;

static bool isPowerOfTwo(unsigned x) {
	return x && (x & (x - 1)) == 0;
}

struct Matrix {
	unsigned _numRows;
	unsigned _numCols;
	bit_t *_data;

	Matrix(unsigned rows, unsigned cols)
		: _numRows(rows), _numCols(cols)
	{
		_data = new bit_t[_numRows * _numCols];
		std::fill(_data, _data + _numRows * _numCols, 0x0);
	}

	int lookup_column(Matrix *v) {
		int ret = -1;
		assert(v->_numRows == _numRows);

		for (unsigned col = 0; col < _numCols; col++) {
			bool ok = true;
			for (unsigned row = 0; row < _numRows; row++) {
				if (_data[row * _numCols + col] != v->_data[row]) {
					ok = false;
				}
			}
			if (ok) {
				ret = col;
				break;
			}
		}
		return ret;
	}

	static Matrix *transpose(Matrix *a) {
		Matrix *m = new Matrix(a->_numCols, a->_numRows);
		for (unsigned row = 0; row < m->_numRows; row++) {
			unsigned off_row = row * m->_numCols;
			for (unsigned col = 0; col < m->_numCols; col++) {
				m->_data[off_row + col] = a->_data[col * a->_numCols + row];
			}
		}
		return m;
	}

	static Matrix *multiply(Matrix *a, Matrix *b) {
		assert(a->_numCols == b->_numRows);
		Matrix *m = new Matrix(a->_numRows, b->_numCols);

		for (unsigned row = 0; row < m->_numRows; row++) {
			for (unsigned col = 0; col < m->_numCols; col++) {
				for (unsigned k = 0; k < a->_numCols; k++) {
					unsigned tmp = a->_data[row * a->_numCols + k]
						*
						b->_data[k * b->_numCols + col];
					m->_data[row * m->_numCols + col] += tmp;
					m->_data[row * m->_numCols + col] %= 2;
				}
			}
		}

		return m;
	}

	static Matrix *identity(int n) {
		Matrix *m = new Matrix(n, n);
		for (int i = 0; i < n; i++) {
			m->_data[i * m->_numRows + i] = 1;
		}
		return m;
	}

	static Matrix *vec_row(bit_t *data, unsigned length) {
		Matrix *m = new Matrix(1, length);
		std::copy(data, data + length, m->_data);
		return m;
	}
	
	static Matrix *vec_col(bit_t *data, unsigned length) {
		Matrix *m = new Matrix(length, 1);
		std::copy(data, data + length, m->_data);
		return m;
	}

	static Matrix *generator(unsigned numInfBits, unsigned numTotalBits)
	{
		unsigned _numRows = numInfBits;
		unsigned _numCols = numTotalBits;
		Matrix *m = new Matrix(numInfBits, numTotalBits);

		//first, set up identity matrix in data bits
		for (unsigned b = 0, row = 0; b < _numCols && row < _numRows; b++)
		{
			if (isPowerOfTwo(b + 1)) {
				continue;
			}
			m->_data[row * _numCols + b] = 1;
			row++;
		}
		
		//now, compute parity bits	
		//number of parity bits
		unsigned nParityBits = numTotalBits - numInfBits;

		//p is current parity bit
		for (unsigned p = 0; p < nParityBits; p++) {
			unsigned parity_col = 1 << p;
			
			for (unsigned row = 0; row < _numRows; row++) {
				//offset of current row in m->_data
				unsigned offset = row * _numCols;
				//current data bit
				for (unsigned b = 0; b < _numCols; b++) {
					if (isPowerOfTwo(b + 1)) {
						continue;
					}
					if ((b + 1) & (parity_col)) {
						m->_data[offset + parity_col - 1] ^= m->_data[offset + b];
					}
				}
			}
		}
		//now, move parity columns to end
		for (unsigned row = 0; row < _numRows; row++) {
			unsigned off_row = row * _numCols;
			for (unsigned p = 0; p < nParityBits; p++) {
				unsigned parity_col_idx = (1 << p) - 1;
				unsigned dest_col = numInfBits + p;
				
				std::swap(m->_data[off_row + parity_col_idx],
					m->_data[off_row + dest_col]);
			}
		}

		//finally, write identity matrix to start
		for (unsigned row = 0; row < _numRows; row++) {
			unsigned off_row = row * _numCols;
			for (unsigned col = 0; col < numInfBits; col++) {
				m->_data[off_row + col] = (row == col);
			}
		}
		return m;
	}
	
	static Matrix *parity(Matrix *gen) {
		unsigned numTotalBits = gen->_numCols;
		unsigned numInfBits = gen->_numRows;
		unsigned numParityBits = numTotalBits - numInfBits;
		Matrix *m = new Matrix(numParityBits, numTotalBits);

		//setup parity matrix as the transpose of generator parity 
		for (unsigned row = 0; row < m->_numRows; row++) {
			unsigned off_row = row * m->_numCols;
			for (unsigned col = 0; col < numInfBits; col++) {
				unsigned g_off = col * gen->_numCols + numInfBits;
				
				m->_data[off_row + col] = gen->_data[g_off + row];
			}
		}

		//setup identity matrix at the end
		for (unsigned row = 0; row < m->_numRows; row++) {
			unsigned off_row = row * m->_numCols + numInfBits;
			m->_data[off_row + row] = 1;
		}

		return m;
	}

	~Matrix() {
		delete[] _data;
	}

	friend std::ostream &operator<<(std::ostream &stream, Matrix const &m) {
		for (unsigned i = 0; i < m._numRows; i++) {
			stream << "|";
			for (unsigned j = 0; j < m._numCols; j++) {
				stream << m._data[i * m._numCols + j] << " ";
			}
			stream << "|" << std::endl;
		}
		
		return stream;
	}
};

#if 0
	static const unsigned DATA_BITS = 4;
	static const unsigned TOTAL_BITS = 7;
#else
	static const unsigned DATA_BITS = 11;
	static const unsigned TOTAL_BITS = 15;
#endif

int main(int argc, char **argv) {
	UNUSED(argc);
	UNUSED(argv);
	std::cout << "Generator(" << DATA_BITS << "," << TOTAL_BITS << ")"
		<< std::endl;
	std::auto_ptr<Matrix> g(Matrix::generator(DATA_BITS, TOTAL_BITS));
	std::cout << *g << std::endl;

	std::cout << "Parity(" << TOTAL_BITS - DATA_BITS << "," << TOTAL_BITS
		<< ")" << std::endl;
	std::auto_ptr<Matrix> h(Matrix::parity(&*g));
	std::cout << *h << std::endl;

	std::cout << "Generator transposed G'" << std::endl;
	std::auto_ptr<Matrix> gt(Matrix::transpose(&*g));
	std::cout << *gt << std::endl;

	std::cout << "H * G' test" << std::endl;
	std::auto_ptr<Matrix> gh_test(Matrix::multiply(&*h, &*gt));
	std::cout << *gh_test << std::endl;
	for (unsigned row = 0; row < gh_test->_numRows; row++) {
		for (unsigned col = 0; col < gh_test->_numCols; col++) {
			assert(gh_test->_data[row * gh_test->_numCols + col] == 0);
		}
	}

	std::cout << "Input vector v" << std::endl;
	bit_t input[DATA_BITS] = {1, 0, 1, 1};
	std::auto_ptr<Matrix> vin(Matrix::vec_col(input, DATA_BITS));
	std::cout << *vin << std::endl;
	
	std::cout << "x = G' * v" << std::endl;
	std::auto_ptr<Matrix> x(Matrix::multiply(&*gt, &*vin));
	std::cout << *x << std::endl;

	std::cout << "syndrome = H * x" << std::endl;
	std::auto_ptr<Matrix> s(Matrix::multiply(&*h, &*x));
	std::cout << *s << std::endl;

	std::cout << "Now, let's flip the last bit of code" << std::endl;
	x->_data[x->_numCols * x->_numRows - 1] ^= 1;
	std::cout << "Vector X with error " << std::endl;
	std::cout << *x << std::endl;

	std::cout << "new syndrome = H * x" << std::endl;
	std::auto_ptr<Matrix> s2(Matrix::multiply(&*h, &*x));
	std::cout << *s2 << std::endl;

	std::cout << "looking up error bit number" << std::endl;
	int bit = h->lookup_column(&*s2);
	if (bit < 0) {
		std::cout << "error bit not found, no error?" << std::endl;
	}
	else {
		std::cout << "error at bit " << bit + 1 << " c++ index "
			<< bit << std::endl;
	}

	return 0;
}