From: drowe67 Date: Mon, 25 Apr 2016 00:53:10 +0000 (+0000) Subject: GPL compliant replacement for golay23.c, thank you so much Tomas! X-Git-Url: http://git.whiteaudio.com/gitweb/?a=commitdiff_plain;h=ca564c8631042d90943978af1417df12bd1b46b7;p=freetel-svn-tracking.git GPL compliant replacement for golay23.c, thank you so much Tomas! git-svn-id: https://svn.code.sf.net/p/freetel/code@2794 01035d8c-6547-0410-b346-abe4f91aad63 --- diff --git a/codec2-dev/src/golay23.c b/codec2-dev/src/golay23.c index 9c31d3ab..d3821799 100644 --- a/codec2-dev/src/golay23.c +++ b/codec2-dev/src/golay23.c @@ -1,419 +1,311 @@ /*---------------------------------------------------------------------------*\ FILE........: golay23.c - AUTHOR......: Robert Morelos-Zaragoza & David Rowe + AUTHOR......: Tomas Härdin & David Rowe DATE CREATED: 3 March 2013 To test: - src$ gcc golay23.c -o golay23 -Wall -DGOLAY23_UNITTEST -DRUN_TIME_TABLES - src$ ./golay23 + src$ gcc golay23.c -o golay23 -Wall -O3 -DGOLAY23_UNITTEST && ./golay23 + src$ gcc golay23.c -o golay23 -Wall -O3 -DGOLAY23_UNITTEST -DRUN_TIME_TABLES && ./golay23 + src$ gcc golay23.c -o golay23 -Wall -O3 -DGOLAY23_UNITTEST -DNO_TABLES && ./golay23 To generate tables: - src$ gcc golay23.c -o golay23 -Wall -DGOLAY23_MAKETABLES -DRUN_TIME_TABLES + src$ gcc golay23.c -o golay23 -Wall -O3 -DGOLAY23_MAKETABLES && ./golay23 \*---------------------------------------------------------------------------*/ -/* File: golay23.c - * Title: Encoder/decoder for a binary (23,12,7) Golay code - * Author: Robert Morelos-Zaragoza (robert@spectra.eng.hawaii.edu) - * Date: August 1994 - * - * The binary (23,12,7) Golay code is an example of a perfect code, that is, - * the number of syndromes equals the number of correctable error patterns. - * The minimum distance is 7, so all error patterns of Hamming weight up to - * 3 can be corrected. The total number of these error patterns is: - * - * Number of errors Number of patterns - * ---------------- ------------------ - * 0 1 - * 1 23 - * 2 253 - * 3 1771 - * ---- - * Total number of error patterns = 2048 = 2^{11} = number of syndromes - * -- - * number of redundant bits -------^ - * - * Because of its relatively low length (23), dimension (12) and number of - * redundant bits (11), the binary (23,12,7) Golay code can be encoded and - * decoded simply by using look-up tables. The program below uses a 16K - * encoding table and an 8K decoding table. - * - * For more information, suggestions, or other ideas on implementing error - * correcting codes, please contact me at (I'm temporarily in Japan, but - * below is my U.S. address): - * - * Robert Morelos-Zaragoza - * 770 S. Post Oak Ln. #200 - * Houston, Texas 77056 - * - * email: robert@spectra.eng.hawaii.edu - * - * Homework: Add an overall parity-check bit to get the (24,12,8) - * extended Golay code. - * - * COPYRIGHT NOTICE: This computer program is free for non-commercial purposes. - * You may implement this program for any non-commercial application. You may - * also implement this program for commercial purposes, provided that you - * obtain my written permission. Any modification of this program is covered - * by this copyright. - * - * == Copyright (c) 1994 Robert Morelos-Zaragoza. All rights reserved. == - */ +/* + Copyright (C) 2016 Tomas Härdin & David Rowe + + All rights reserved. + + This program is free software; you can redistribute it and/or modify + it under the terms of the GNU Lesser General Public License version 2.1, as + published by the Free Software Foundation. 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. -#include "golay23.h" + You should have received a copy of the GNU Lesser General Public License + along with this program; if not,see . +*/ #include -#include -#include -#define X22 0x00400000 /* vector representation of X^{22} */ -#define X11 0x00000800 /* vector representation of X^{11} */ -#define MASK12 0xfffff800 /* auxiliary vector for testing */ -#define GENPOL 0x00000c75 /* generator polinomial, g(x) */ - -/* Global variables: - * - * pattern = error pattern, or information, or received vector - * encoding_table[] = encoding table - * decoding_table[] = decoding table - * data = information bits, i(x) - * codeword = code bits = x^{11}i(x) + (x^{11}i(x) mod g(x)) - * numerr = number of errors = Hamming weight of error polynomial e(x) - * position[] = error positions in the vector representation of e(x) - * recd = representation of corrupted received polynomial r(x) = c(x) + e(x) - * decerror = number of decoding errors - * a[] = auxiliary array to generate correctable error patterns - */ -static int inited = 0; +#ifdef GOLAY23_MAKETABLES +#define RUN_TIME_TABLES +#endif +#ifndef NO_TABLES #ifdef RUN_TIME_TABLES -static int encoding_table[4096], decoding_table[2048]; +int static encoding_table[4096]; +int static decoding_table[2048]; +static int inited = 0; #else +//default is to use precomputed tables #include "golayenctable.h" #include "golaydectable.h" #endif - -#ifdef GOLAY23_UNITTEST -static int position[23] = { 0x00000001, 0x00000002, 0x00000004, 0x00000008, - 0x00000010, 0x00000020, 0x00000040, 0x00000080, - 0x00000100, 0x00000200, 0x00000400, 0x00000800, - 0x00001000, 0x00002000, 0x00004000, 0x00008000, - 0x00010000, 0x00020000, 0x00040000, 0x00080000, - 0x00100000, 0x00200000, 0x00400000 }; #endif -#ifdef RUN_TIME_TABLES -static int arr2int(int a[], int r) -/* - * Convert a binary vector of Hamming weight r, and nonzero positions in - * array a[1]...a[r], to a long integer \sum_{i=1}^r 2^{a[i]-1}. - */ -{ - int i; - long mul, result = 0, temp; - - for (i=1; i<=r; i++) { - mul = 1; - temp = a[i]-1; - while (temp--) - mul = mul << 1; - result += mul; - } - return(result); +//since we want to avoid bit-reversing inside syndrome() we bit-reverse the polynomial instead +#define GOLAY_POLYNOMIAL 0xC75 //AE3 reversed + +static int syndrome(int c) { + //could probably be done slightly smarter, but works + int x; + for (x = 11; x >= 0; x--) { + if (c & ((1<<11) << x)) { + c ^= GOLAY_POLYNOMIAL << x; + } + } + return c; +} + +#ifdef __GNUC__ +#define popcount __builtin_popcount +#elif defined(__MSC_VER) +#include +#define popcount __popcnt +#else +static int popcount(unsigned int c) { + int ret = 0; + while (c) { + if (c & 1) { + ret++; + } + c >>= 1; + } + return ret; } #endif -void nextcomb(int n, int r, int a[]) -/* - * Calculate next r-combination of an n-set. - */ -{ - int i, j; - - a[r]++; - if (a[r] <= n) - return; - j = r - 1; - while (a[j] == n - r + j) - j--; - for (i = r; i >= j; i--) - a[i] = a[j] + i - j + 1; - return; +#if defined(NO_TABLES) || defined(RUN_TIME_TABLES) +static int golay23_encode_no_tables(int c) { + c <<= 11; + return syndrome(c) | c; } +#endif -int get_syndrome(int pattern) -/* - * Compute the syndrome corresponding to the given pattern, i.e., the - * remainder after dividing the pattern (when considering it as the vector - * representation of a polynomial) by the generator polynomial, GENPOL. - * In the program this pattern has several meanings: (1) pattern = infomation - * bits, when constructing the encoding table; (2) pattern = error pattern, - * when constructing the decoding table; and (3) pattern = received vector, to - * obtain its syndrome in decoding. - */ -{ - int aux = X22; - - if (pattern >= X11) - while (pattern & MASK12) { - while (!(aux & pattern)) - aux = aux >> 1; - pattern ^= (aux/X11) * GENPOL; - } - return(pattern); +#ifdef NO_TABLES +static int unrotate(unsigned int c, int x) { + return ((c << x) & 0x7FFFFF) | (c >> (23 - x)); } -/*---------------------------------------------------------------------------*\ +static int golay23_decode_no_tables(int c) { + //TODO: optimize? + int x; + c = unrotate(c, 12); - FUNCTION....: golay23_init() - AUTHOR......: David Rowe - DATE CREATED: 3 March 2013 + for (x = 0; x < 23; x++) { + int t; + int s = syndrome(c); - Call this once when you start your program to init the Golay tables. + if (popcount(s) <= 3) { + return unrotate(c ^ s, x) & 0xFFF; + } -\*---------------------------------------------------------------------------*/ + for (t = 0; t < 23; t++) { + int c2 = c ^ (1 << t); + int s = syndrome(c2); + + if (popcount(s) <= 2) { + return unrotate(c2 ^ s, x) & 0xFFF; + } + } + + //rotate + c = (c >> 1) | ((c & 1) << 22); + } + + //shouldn't reach here.. + assert("Something is wrong with golay23_decode_no_tables().."); + return c & 0xFFF; +} +#endif void golay23_init(void) { #ifdef RUN_TIME_TABLES - int i; - long temp; - int a[4]; - int pattern; - - /* - * --------------------------------------------------------------------- - * Generate ENCODING TABLE - * - * An entry to the table is an information vector, a 32-bit integer, - * whose 12 least significant positions are the information bits. The - * resulting value is a codeword in the (23,12,7) Golay code: A 32-bit - * integer whose 23 least significant bits are coded bits: Of these, the - * 12 most significant bits are information bits and the 11 least - * significant bits are redundant bits (systematic encoding). - * --------------------------------------------------------------------- - */ - for (pattern = 0; pattern < 4096; pattern++) { - temp = pattern << 11; /* multiply information by X^{11} */ - encoding_table[pattern] = temp + get_syndrome(temp);/* add redundancy */ - } + int x, y, z; + inited = 1; + for (x = 0; x < 4096; x++) { + encoding_table[x] = golay23_encode_no_tables(x); + } - /* - * --------------------------------------------------------------------- - * Generate DECODING TABLE - * - * An entry to the decoding table is a syndrome and the resulting value - * is the most likely error pattern. First an error pattern is generated. - * Then its syndrome is calculated and used as a pointer to the table - * where the error pattern value is stored. - * --------------------------------------------------------------------- - * - * (1) Error patterns of WEIGHT 1 (SINGLE ERRORS) - */ decoding_table[0] = 0; - decoding_table[1] = 1; - temp = 1; - for (i=2; i<= 23; i++) { - temp *= 2; - decoding_table[get_syndrome(temp)] = temp; + //1-bit errors + for (x = 0; x < 23; x++) { + int d = 1<= 0 && c <= 0xFFF); +#ifdef RUN_TIME_TABLES assert(inited); - assert(data <= 0xfff); +#endif - //printf("data: 0x%x\n", data); - return encoding_table[data]; +#ifdef NO_TABLES + return golay23_encode_no_tables(c); +#else + return encoding_table[c]; +#endif } -/*---------------------------------------------------------------------------*\ - - FUNCTION....: golay23_decode() - AUTHOR......: David Rowe - DATE CREATED: 3 March 2013 - - Given a 23 bit received codeword, returns the 12 bit corrected data. - -\*---------------------------------------------------------------------------*/ - -int golay23_decode(int received_codeword) { +int golay23_decode(int c) { + assert(c >= 0 && c <= 0x7FFFFF); +#ifdef RUN_TIME_TABLES assert(inited); - assert((received_codeword < (1<<23)) && (received_codeword >= 0)); - - //printf("syndrome: 0x%x\n", get_syndrome(received_codeword)); - return received_codeword ^= decoding_table[get_syndrome(received_codeword)]; -} - -int golay23_count_errors(int recd_codeword, int corrected_codeword) -{ - int errors = 0; - int diff, i; - - diff = recd_codeword ^ corrected_codeword; - for(i=0; i<23; i++) { - if (diff & 0x1) - errors++; - diff >>= 1; - } +#endif - return errors; +#ifdef NO_TABLES + //duplicate old golay23_decode()'s shift + return unrotate(golay23_decode_no_tables(c), 11); +#else + //message is shifted 11 places left in the return value + return c ^ decoding_table[syndrome(c)]; +#endif } -#ifdef GOLAY23_UNITTEST - -static int golay23_test(int error_pattern) { - int data; - int codeword; - int recd; - int pattern; - int decerror; - int i, tests; - - decerror = 0; - tests = 0; - - for (data = 0; data<(1<<12); data++) { - - codeword = golay23_encode(data); - recd = codeword ^ error_pattern; - recd = golay23_decode(recd); - pattern = (recd ^ codeword) >> 11; - for (i=0; i<12; i++) - if (pattern & position[i]) - decerror++; - if (decerror) { - printf("data: 0x%x codeword: 0x%x recd: 0x%x\n", data, codeword, recd); - printf("there were %d decoding errors\n", decerror); - exit(1); - } - tests++; - } - - return tests; +int golay23_count_errors(int recd_codeword, int corrected_codeword) { + return popcount(recd_codeword ^ corrected_codeword); } -int main(void) -{ - int i; - int tests; - int a[4]; - int error_pattern; +/** + * Table generation and testing code below + */ - golay23_init(); +#ifdef GOLAY23_MAKETABLES +#include - /* --------------------------------------------------------------------- - * Generate DATA - * --------------------------------------------------------------------- - */ +int main() { + int x; + //generate and dump + golay23_init(); - /* Test all combinations of data and 1,2 or 3 errors */ + FILE *enc = fopen("golayenctable.h", "w"); + FILE *dec = fopen("golaydectable.h", "w"); - tests = 0; - error_pattern = 1; - for (i=0; i< 23; i++) { - //printf("error_pattern: 0x%x\n", error_pattern); - tests += golay23_test(error_pattern); - error_pattern *= 2; - } - printf("%d 1 bit error tests performed OK!\n", tests); - - tests = 0; - a[1] = 1; a[2] = 2; - error_pattern = arr2int(a,2); - tests += golay23_test(error_pattern); - for (i=1; i<253; i++) { - nextcomb(23,2,a); - error_pattern = arr2int(a,2); - //printf("error_pattern: 0x%x\n", error_pattern); - tests += golay23_test(error_pattern); + fprintf(enc, "/* Generated by golay23.c -DGOLAY23_MAKETABLE */\n\ +\n\ +const int static encoding_table[]={\n"); + for (x = 0; x < 4096; x++) { + fprintf(enc, x < 4095 ? " 0x%x,\n" : " 0x%x\n", encoding_table[x]); } - printf("%d 2 bit error tests performed OK!\n", tests); - - tests = 0; - a[1] = 1; a[2] = 2; a[3] = 3; - error_pattern = arr2int(a,3); - tests += golay23_test(error_pattern); - for (i=1; i<1771; i++) { - nextcomb(23,3,a); - error_pattern = arr2int(a,3); - //printf("error_pattern: 0x%x\n", error_pattern); - tests += golay23_test(error_pattern); + fprintf(enc, "};\n"); + + fprintf(dec, "/* Generated by golay23.c -DGOLAY23_MAKETABLE */\n\ +\n\ +const int static decoding_table[]={\n"); + for (x = 0; x < 2048; x++) { + fprintf(dec, x < 2047 ? " 0x%x,\n" : " 0x%x\n", decoding_table[x]); } - printf("%d 3 bit error tests performed OK!\n", tests); + fprintf(dec, "};\n"); + + fclose(enc); + fclose(dec); return 0; } -#endif -#ifdef GOLAY23_MAKETABLES -int main(int argc, char *argv[]) { - FILE *f; - int i; +#elif defined(GOLAY23_UNITTEST) +#include +#include +#include + +int main() { + int c; golay23_init(); - f=fopen("golayenctable.h","wt"); - assert(f != NULL); + //keep track of whether every single codeword has been checked + char *checkmask = malloc(1<<23); + memset(checkmask, 0, 1<<23); - fprintf(f,"/* Generated by golay23.c -DGOLAY23_MAKETABLE */\n\n"); - fprintf(f,"const int static encoding_table[]={\n"); + //step through all possible messages + for (c = 0; c < (1<<12); c++) { + int g23 = golay23_encode(c); + int x,y,z; + checkmask[g23] = 1; + int c2 = golay23_decode(g23) >> 11; - for (i=0; i<4095; i++) - fprintf(f," 0x%x,\n", encoding_table[i]); - fprintf(f, " 0x%x\n};\n", encoding_table[i]); - fclose(f); + printf("%03x -> %06x %03x\n", c, g23, c2); - f=fopen("golaydectable.h","wt"); - assert(f != NULL); + if (c != c2) { + printf("Bad!\n"); + exit(1); + } - fprintf(f,"/* Generated by golay23.c -DGOLAY23_MAKETABLE */\n\n"); - fprintf(f,"const int static decoding_table[]={\n"); + //test the code by flipping every combination of one, two and three bits + for (x = 0; x < 23; x++) { + int flipped = g23 ^ (1<> 11; + if (c != c2) { + printf("Bad!\n"); + + exit(1); + } + } + + for (x = 0; x < 22; x++) { + for (y = x+1; y < 23; y++) { + int flipped = g23 ^ (1<> 11; + if (c != c2) { + printf("Bad!\n"); + + exit(1); + } + } + } + + for (x = 0; x < 21; x++) { + for (y = x+1; y < 22; y++) { + for (z = y+1; z < 23; z++) { + int flipped = g23 ^ (1<> 11; + if (c != c2) { + printf("Bad!\n"); + exit(1); + } + } + } + } + } - for (i=0; i<2047; i++) - fprintf(f," 0x%x,\n", decoding_table[i]); - fprintf(f, " 0x%x\n};\n", decoding_table[i]); - fclose(f); + //did we check every codeword? + for (c = 0; c < (1<<23); c++) { + if (checkmask[c] != 1) { + printf("%06x unchecked!\n", c); + exit(1); + } + } + printf("Everything checks out\n"); + free(checkmask); return 0; } - #endif - -