A C implementation of LDPC over GF(q)

Seishi Takamura @ Stanford Univ./NTT Corporation

Since June 9, 2006.

Introduction

This page is on LDPC over GF(q>=2) software, absolutely with no warranty. The source codes are written in plain C. You do not need any support applications other than a C compiler. I've tested with VC/gcc on Cygwin and gcc on Linux.

With use of several tricks such as it speeded up by 30 times for GF(32) and more than 1,000 times for GF(256) compared to straightforward implementation.

This program can be used not only for multi-level signals but for (concatenated) binary sources i.e., for m-tuple signal simply GF(2m) does that.

About the speed, GF2_LDPC is fastest, even though taking into account that GFq_LDPC_NTT decodes multiple bits (i.e., GF(8) LDPC is more than three (=log28) times slower than binary LDPC. The advantage of GFq_LDPC_NTT is its potential of better error correction performance (you can easily find a lot of articles on this topic).
Comments are very much welcome. Please give any comments to .

Download

The package contains C sources and some LDPC matrix files.
LDPC-v0.1.zip (649kB) Initial version (June 9, 2006).
LDPC-v0.2.zip (649kB) Fixed a bug in GF2_LDPC.c for irregular matrix case (such as 3584x444). Special thanks: Mr. David Elkouss (Feb. 5, 2008)

Compilation and demonstration

For binary decoding demo, compile and run 
% gcc GF2_LDPC.c -lm
% ./a.out 0.08 4986.93xb.329.alist
This performs binary LDPC decoding of rate (m/n=)0.5 over BSC channel of 8% crossover probability.

For non-binary, macro-variable Log2Q is set to 5 (i.e., LDPC over GF(32)). Please override the definition with compiler's command line option.

For GF(8) decoding demo, compile and run 
% gcc -DLog2Q=3 GFq_LDPC_NTT.c -lm
% ./a.out bsc 0.147 q8.sp.6000.4000.3000.1
For GF(4) demo, compile and run 
% gcc -DLog2Q=2 GFq_LDPC_NTT.c -lm
% ./a.out bsc 0.145 q4.sp.9000.6000.4500.1
I do not have LDPC matrices for other than GF order of 4 and 8. (your contribution appreciated)

If you wish, you can specify maximum sum-product iteration using "-iter num" option as the first argument e.g., 
% ./a.out -iter 50 bsc 0.145 q4.sp.9000.6000.4500.1
Internally, the demo does
  1. read m, n values and m*n LDPC matrix H from an alist-format file
  2. create a q-ary random vector x of length n (source signal)
  3. obtain s = Hx in GF(q) of length m (syndrome signal)
  4. add BSC(for q>=2) or Laplacian(for q>2) noise onto x to yield y
  5. do sum-product loop using y and s to obtain x^ (recovered signal)
  6. if Hx^ equals s, then stop, else continue from step 5
  7. do next experiment from step 4 (x is unchanged, only additive noise is changed)
Execution snapshot looks as follows, where HamDist stands for Hamming distance and synd(x) means Hx. 
% ./GFq_LDPC_NTT bsc 0.1 ../LDPC/q8.sp.6000.4000.3000.1
Warning! GF order(32) does not match with matrix file(8)
I'll use it anyways.

GF(32) experiment 1:
m/n=0.666667, BSC channel capacity(rate) = 0.468996 (bits)
start: HamDist(x,y)=3000 HamDist(synd(x),synd(y))=7199
 1: HamDist(x,x^)=2535 HamDist(s,synd(x^))=4020
 2: HamDist(x,x^)=1756 HamDist(s,synd(x^))=2619
 3: HamDist(x,x^)=1075 HamDist(s,synd(x^))=1608
 4: HamDist(x,x^)=514 HamDist(s,synd(x^))=794
 5: HamDist(x,x^)=171 HamDist(s,synd(x^))=288
 6: HamDist(x,x^)=36 HamDist(s,synd(x^))=53
 7: HamDist(x,x^)=8 HamDist(s,synd(x^))=11
 8: HamDist(x,x^)=0 HamDist(s,synd(x^))=0
converged.

GF(32) experiment 2:
m/n=0.666667, BSC channel capacity(rate) = 0.468996 (bits)
start: HamDist(x,y)=3000 HamDist(synd(x),synd(y))=7137
 1: HamDist(x,x^)=2606 HamDist(s,synd(x^))=3887
 2: HamDist(x,x^)=1818 HamDist(s,synd(x^))=2539
 3: HamDist(x,x^)=1136 HamDist(s,synd(x^))=1556
 4: HamDist(x,x^)=627 HamDist(s,synd(x^))=917
 5: HamDist(x,x^)=247 HamDist(s,synd(x^))=379
 6: HamDist(x,x^)=75 HamDist(s,synd(x^))=130
 7: HamDist(x,x^)=2 HamDist(s,synd(x^))=6
 8: HamDist(x,x^)=0 HamDist(s,synd(x^))=0
converged.

GF(32) experiment 3:
m/n=0.666667, BSC channel capacity(rate) = 0.468996 (bits)
start: HamDist(x,y)=3000 HamDist(synd(x),synd(y))=7191
 1: HamDist(x,x^)=2661 HamDist(s,synd(x^))=4008
 2: HamDist(x,x^)=1850 HamDist(s,synd(x^))=2649
 3: HamDist(x,x^)=1104 HamDist(s,synd(x^))=1611
 4: HamDist(x,x^)=569 HamDist(s,synd(x^))=861
 5: HamDist(x,x^)=223 HamDist(s,synd(x^))=420
 6: HamDist(x,x^)=59 HamDist(s,synd(x^))=103
 7: HamDist(x,x^)=7 HamDist(s,synd(x^))=11
 8: HamDist(x,x^)=2 HamDist(s,synd(x^))=6
 9: HamDist(x,x^)=0 HamDist(s,synd(x^))=0
converged.

Notes

GF2_LDPC, a binary (GF(2)) LDPC decoder, is based on the IEEE paper of A.D. Liveris et al.: "Compression of binary sources with side information at the decoder using LDPC codes".

Low density parity matrix is in 'alist' format, see Prof. MacKay's description for details.

GFq_LDPC_NTT, a non-binary (GF(2m)) LDPC decoder, is based on Dr. Davey's dissertation, "Error-correction using Low-Density Parity-Check Codes". Please be noted that there's a typo in the first line of page 26 of the dissertation:
Ramn is the a'th coordinate
should be
Ramn is the (zm - Hmna)'th coordinate
which is already incorporated in my source code.

The original alist format is aimed for binary LDPC, however, its extension to GF(q>2) is found in GFQmatrices directory. Also the matrix examples (in fact these are the only GF(q>2) matrices I know of) can be found in the same directory.

API

Binary LDPC (GF2_LDPC)
void initdec(char *s)
s: alist file name
-- Reads in a low density parity check matrix H in alist format and sets up necessary variables such as m, n etc.
void bsc(int x[], int y[], double p, int q0[])
-- lets y[] be noisy version of x[], with crossover probability of p, and sets the values of q0[] (see description below)
int dec(int q0[], int s[], int loop_max, int x[])
q0[0<=i<n]: log likelihood ratio, equals log(Prob(x[i]=0|y[i]) / Prob(x[i]=1|y[i]))
s[m]: syndrome
loop_max: maximum loop iteration
x[n]: original input signal
return value: 0=converged, -1=not converged
-- sum-product decoder. Please be noted that information of y[] is already incorporated in q0[]
Non-binary LDPC (GFq_LDPC_NTT)
void initdec(char *s)
s: matrix file name
-- Reads in a low density parity check matrix H in extended alist format and sets up necessary variables such as m, n etc.
void bsc(int x[], int y[], double p, int **logfna)
-- lets each bit of y[] be noisy version of corresponding bit of x[], with crossover probability of p, and sets the values of logfna[][] (see below)
void lap(int x[], int y[], double stddev, int **logfna)
-- lets y[] be noisy version of x[], with Laplacian noise with standard deviation of stddev, and sets the values of logfna[][] (see below)
int dec(int **logfna, int z[], int loop_max, int x[])
logfna[0<=i<n][0<=a<q]: prior, equals log(Prob(x[i]=a|y[i])
z[m]: syndrome
loop_max: maximum loop iteration
x[n]: original input signal (not used for sum-product calculation, just for reference)
-- sum-product decoder

Tweaks and modifications

You might perhaps want to tweak the following fixed-point-number-related definitions in the source code. 
#define INT	6              // int part
#define DECI	14             // fraction part
Only up to GF(256) is supported so far. To raise this limit, just add new logq[] and expq[] arrays appropriately (your contribution appreciated).

The case for GF(prime_number) can be simply implemented. Just take modulo of the prime number after normal addition/subtraction/multiplication operations.

Bugs

You can compile GFq_LDPC_NTT.c as a binary LDPC decoder by compiling 
% gcc -DLog2Q=1 GFq_LDPC_NTT.c -lm
but this does not execute well. Please use GF2_LDPC.c for this purpose.

LDPC decoding convergence (i.e., successful/unsuccessful decoding) is now automatically determined by watching the transition of Hamming distance between given syndrome and temporarily decoded signal's syndrome. Though this practically works, it is not very intelligent.