Author: michael
Date: Mon Nov 24 23:08:33 2008
New Revision: 123
Log:
Low density parity check codes support.
Added:
trunk/noe/ldpc.c
trunk/noe/ldpc.h
trunk/noe/test_ldpc.c
Modified:
trunk/noe/Makefile
Modified: trunk/noe/Makefile
==============================================================================
--- trunk/noe/Makefile (original)
+++ trunk/noe/Makefile Mon Nov 24 23:08:33 2008
@@ -1,7 +1,7 @@
LIBOBJS = galois.o rs.o
OBJS = rw.o filerw.o noe.o
LIBS = libnoe_100.a libnoe_101.a libnoe_10001.a
-PUBLIC_HEADERS = rs.h galois.h
+PUBLIC_HEADERS = rs.h galois.h ldpc.h
LIBSRCS = $(LIBOBJS:.o=.c) $(ASM_LIBOBJS:.o=.s)
SRCS = $(OBJS:.o=.c)
INSTALL_BASE = /usr/local
@@ -20,9 +20,9 @@ CFLAGS = -g -Wall -O4 $(OPTFLAGS) -I. $
%_10001.o: %.c
$(CC) $(CFLAGS) -c -DSIZE=0x10001 -o $@ $<
-all: test_100 test_101 test_10001 mina
+all: test_100 test_101 test_10001 mina test_ldpc_100
-libnoe_100.a: $(LIBOBJS:.o=_100.o)
+libnoe_100.a: $(LIBOBJS:.o=_100.o) ldpc_100.o
libnoe_101.a: $(LIBOBJS:.o=_101.o) gfft_101.o
libnoe_10001.a: $(LIBOBJS:.o=_10001.o) gfft_10001.o
libnoe_100.a libnoe_101.a libnoe_10001.a:
@@ -37,6 +37,9 @@ test: all
test_%: test_%.o libnoe_%.a
$(CC) $(LDFLAGS) -o $@ $^
+test_ldpc_%: test_ldpc_%.o libnoe_%.a
+ $(CC) $(LDFLAGS) -o $@ $^
+
clean:
rm -f *.o *.a *~
Added: trunk/noe/ldpc.c
==============================================================================
--- (empty file)
+++ trunk/noe/ldpc.c Mon Nov 24 23:08:33 2008
@@ -0,0 +1,211 @@
+/*
+ * Copyright (C) 2008 Michael Niedermayer <michaelni(a)gmx.at>
+ *
+ * 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 <inttypes.h>
+#include <assert.h>
+#include <math.h>
+#include <malloc.h>
+#include <stdio.h>
+#include <string.h>
+
+#include "noe_internal.h"
+#include "galois_internal.h"
+#include "random_internal.h"
+#include "ldpc.h"
+
+#undef NDEBUG
+#include <assert.h>
+
+#define MAX_NZC 32
+
+typedef struct LDPCContext{
+ unsigned int (*parity_matrix)[2][MAX_NZC];
+ GFF4Element *inv_matrix;
+ int parity_len;
+ int data_len;
+ int nzc;
+}LDPCContext;
+
+typedef uint8_t ELEM;
+
+int inverse(ELEM *matrix, int width, int height, int solvew){
+ int i, j, k, m, p;
+
+ for(p=i=0; i<solvew; i++){
+ for(j=p; j<height; j++){
+ if(matrix[i + j*width]){
+ unsigned int logline[width];
+ unsigned int inve= EXT(log)[inv(matrix[i + j*width])];
+ for(m=i; m<width; m++){
+ ELEM t= matrix[m + j*width];
+ matrix[m + j*width]= matrix[m + p*width];
+ matrix[m + p*width]= EXT(exp)[EXT(log)[t] + inve];
+ logline[m]= EXT(log)[matrix[m + p*width]];
+ }
+ for(k=0; k<height; k++){
+ if(matrix[i + k*width] && k!=p){
+ unsigned int factor= EXT(log)[neg(matrix[i + k*width])];
+ for(m=i; m<width; m++)
+ matrix[m + k*width] = sum(matrix[m + k*width], EXT(exp)[logline[m] + factor]);
+ }
+ }
+ p++;
+ break;
+ }
+ }
+ fprintf(stderr, "%5d/%5d\r", i, width);
+ }
+ return p;
+}
+
+LDPCContext *EXT(initLDPC)(int data_len, int parity_len, uint32_t seed){
+ LDPCContext *c= malloc(sizeof(LDPCContext));
+ KISSState kiss;
+ int i, j, y;
+ int nzc= noe_log2(parity_len)+1;
+ int code_len= data_len + parity_len;
+
+ if(nzc > MAX_NZC)
+ goto fail;
+ assert(nzc <= parity_len);
+
+ init_random(&kiss, seed);
+
+ c->parity_matrix= malloc(code_len * sizeof(*c->parity_matrix));
+ c-> inv_matrix= malloc(parity_len*parity_len * sizeof(*c->inv_matrix));
+
+ for(i=0; i<code_len; i++){
+ uint8_t tab[parity_len];
+ memset(tab, 0, parity_len);
+ for(j=0; j<nzc; j++){
+ do{
+ y= get_random(&kiss)%parity_len;
+ }while(tab[y]);
+ tab[y]=1;
+ c->parity_matrix[i][0][j] = y;
+ c->parity_matrix[i][1][j] = EXT(log)[get_random(&kiss)%(SIZE-1)+1];
+ }
+ }
+#if 0
+ for(;;){
+ int maxrow=0;
+ int minrow=999999;
+ int maxrowx=0;
+ uint8_t tab[parity_len];
+ memset(tab, 0, parity_len);
+
+ for(j=0; j<parity_len; j++){
+ if(row[j] > maxrow){
+ maxrow= row[j];
+ maxrowx= j;
+ }
+ if(row[j] < minrow)
+ minrow= row[j];
+ }
+ do{
+ j= get_random(&kiss)%nzc;
+ i= get_random(&kiss)%code_len;
+ }while(c->parity_matrix[i][0][j] != maxrowx);
+ for(y=0; y<nzc; y++)
+ tab[ c->parity_matrix[i][0][y] ]=1;
+ do{
+ y= get_random(&kiss)%parity_len;
+ }while(tab[y]);
+ c->parity_matrix[i][0][j] = y;
+ row[y] ++;
+ row[maxrowx]--;
+// fprintf(stderr, "%d %d\n", maxrow, minrow);
+ if(maxrow - 3 <= minrow)
+ break;
+ }
+#endif
+
+ c-> data_len= data_len;
+ c->parity_len= parity_len;
+ c->nzc = nzc;
+
+ return c;
+fail:
+ free(c);
+ return NULL;
+}
+
+int EXT(init_matrixLDPC)(LDPCContext *c, int erasure_count, int erasure_pos[]){
+ int i, j, rank, x, y;
+ int nzc= c->nzc;
+ int parity_len= c->parity_len;
+ ELEM inv_matrix[2*parity_len*parity_len];
+
+ memset(inv_matrix, 0, sizeof(inv_matrix));
+ for(i=0; i<erasure_count; i++){
+ x= erasure_pos[i];
+ for(j=0; j<nzc; j++){
+ y= c->parity_matrix[x][0][j];
+ inv_matrix[i + y*2*parity_len]= EXT(exp)[c->parity_matrix[x][1][j]];
+ }
+ }
+ for(i=0; i<parity_len; i++){
+ inv_matrix[i + i*2*parity_len + parity_len]= 1;
+ }
+
+ rank= inverse(inv_matrix, 2*parity_len, parity_len, erasure_count);
+ for(j=0; j<erasure_count; j++)
+ for(i=0; i<parity_len; i++)
+ c->inv_matrix[i + j*parity_len]= EXT(log)[inv_matrix[i + j*2*parity_len + parity_len]];
+
+ return rank - erasure_count;
+}
+
+int EXT(decodeLDPC)(LDPCContext *c, uint8_t *code, int erasure_count, int erasure_pos[]){
+ int i,j;
+ int parity_len= c->parity_len;
+ int code_len= c->data_len + parity_len;
+ GFF4Element syndrom[c->parity_len];
+ memset(syndrom, 0, sizeof(syndrom));
+
+ for(i=0; i<erasure_count; i++)
+ code[ erasure_pos[i] ]= 0;
+
+ for(i=0; i<code_len; i++){
+ GFF4Element v= EXT(log)[ code[i] ];
+ for(j=0; j<c->nzc; j++){
+ int y= c->parity_matrix[i][0][j];
+ int p= c->parity_matrix[i][1][j];
+ syndrom[y]= sum(syndrom[y], EXT(exp)[p + v]);
+ }
+ }
+
+ for(i=0; i<parity_len; i++)
+ syndrom[i]= EXT(log)[ syndrom[i] ];
+
+ for(i=0; i<erasure_count; i++){
+ GFF4Element v=0;
+ for(j=0; j<parity_len; j++){
+ v= sum(v, EXT(exp)[syndrom[j] + c->inv_matrix[j + i*parity_len]]);
+ }
+ code[ erasure_pos[i] ]= v;
+ }
+
+ return 0;
+}
+
+void EXT(freeLDPC)(LDPCContext *c){
+ free(c->parity_matrix);
+ free(c->inv_matrix);
+ free(c);
+}
Added: trunk/noe/ldpc.h
==============================================================================
--- (empty file)
+++ trunk/noe/ldpc.h Mon Nov 24 23:08:33 2008
@@ -0,0 +1,38 @@
+/*
+ * Copyright (C) 2008 Michael Niedermayer <michaelni(a)gmx.at>
+ *
+ * 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
+ */
+
+struct LDPCContext;
+
+struct LDPCContext *EXT(initLDPC)(int data_len, int parity_len, uint32_t seed);
+
+/**
+ *
+ * @return 0 if sucessfull, <0 if not (will only partially decode)
+ * for encoding you should try another seed when this returns <0
+ */
+int EXT(init_matrixLDPC)(struct LDPCContext *c, int erasure_count, int erasure_pos[]);
+
+/**
+ *
+ * @param code the code word to encode/decode.
+ * @param erasure_count must match what has been passed to init_matrixLDPC
+ * @param erasure_pos must match what has been passed to init_matrixLDPC
+ */
+int EXT(decodeLDPC)(struct LDPCContext *c, uint8_t *code, int erasure_count, int erasure_pos[]);
+
+void EXT(freeLDPC)(struct LDPCContext *c);
Added: trunk/noe/test_ldpc.c
==============================================================================
--- (empty file)
+++ trunk/noe/test_ldpc.c Mon Nov 24 23:08:33 2008
@@ -0,0 +1,101 @@
+/*
+ * Copyright (C) 2008 Michael Niedermayer <michaelni(a)gmx.at>
+ *
+ * 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 <inttypes.h>
+#include <assert.h>
+#include <math.h>
+#include <stdio.h>
+#include <string.h>
+
+#include "galois_internal.h"
+#include "random_internal.h"
+#include "ldpc.h"
+
+#undef NDEBUG
+#include <assert.h>
+
+int main(){
+ int data_len, parity_len, i, x;
+ KISSState kiss;
+
+ init_random(&kiss, 0);
+ EXT(init)();
+
+ for(data_len=100; data_len < 2000; data_len<<=1){
+ for(parity_len=10; parity_len < data_len; parity_len<<=1){
+ const int code_len= data_len + parity_len;
+ int seed, erasure_count, try;
+ struct LDPCContext *enc= NULL;
+ struct LDPCContext *dec= NULL;
+ int parity_pos[parity_len];
+
+ for(i=0; i<parity_len; i++)
+ parity_pos[i]= data_len + i;
+
+ for(seed=0; 1; seed++){
+ enc= EXT(initLDPC)(data_len, parity_len, seed);
+ if(EXT(init_matrixLDPC)(enc, parity_len, parity_pos) >= 0)
+ break;
+ EXT(freeLDPC)(enc);
+ }
+
+ dec= EXT(initLDPC)(data_len, parity_len, seed);
+
+ for(erasure_count= parity_len; erasure_count > parity_len-3; erasure_count--){
+ int fail=0;
+ for(try=0; try<10; try++){
+ uint8_t code[code_len], code_bak[code_len];
+ int erasure_pos[erasure_count];
+
+ for(i=0; i<data_len; i++)
+ code[i]= code_bak[i]= get_random(&kiss);
+ EXT(decodeLDPC)(enc, code, parity_len, parity_pos);
+ memcpy(code_bak+data_len, code+data_len, parity_len);
+ for(i=0; i<code_len; i++)
+ if(code[i] != code_bak[i])
+ fprintf(stderr, "FATALXX error %X!= %X at %d\n", code[i], code_bak[i], i);
+
+ for(i=0; i<erasure_count; i++){
+ do{
+ x= get_random(&kiss) % code_len;
+ }while(code[x] != code_bak[x]);
+ code[x]+= get_random(&kiss)%255+1;
+ erasure_pos[i]= x;
+ }
+
+ if(EXT(init_matrixLDPC)(dec, erasure_count, erasure_pos) < 0){
+ fail++;
+ continue;
+ }
+
+ EXT(decodeLDPC)(dec, code, erasure_count, erasure_pos);
+ for(i=0; i<code_len; i++){
+ if(code[i] != code_bak[i])
+ fprintf(stderr, "FATAL error %X!= %X at %d\n", code[i], code_bak[i], i);
+ }
+ }
+ printf("data:%5d parity:%5d extra:%d fail:%d seed:%d\n", data_len, parity_len, parity_len-erasure_count, fail, seed);
+ if(!fail)
+ break;
+ }
+ EXT(freeLDPC)(dec);
+ EXT(freeLDPC)(enc);
+ }
+ }
+ return 0;
+}
