Low density parity check (LDPC) codes are among the best performing codes in practice, being capable of correcting errors close to the Shannon limit. Magma provides facilities for the construction, decoding, simulation and analysis of LDPC codes.
LDPC codes come in two main varieties, regular and irregular, defined by the row and column weights of the sparse parity check matrix. If all columns in the parity check matrix have some constant weight a, and all rows have some constant weight b, then the LDPC code is said to be (a, b)-regular. When either the columns or the rows have a distribution of weights, the LDPC code is said to be irregular.
Currently, there do not exist many techniques for the explicit construction of LDPC codes. More commonly, these codes are selected at random from an ensemble, and their properties determined through simulation.
Given a sparse binary matrix H, return the LDPC code which has H as its parity check matrix.
Return a random (a, b)-regular LDPC code of length n, using Gallager's original method of construction. The row weight a must dividethe length n.
Return a random code from the ensemble of (a, b)-regular binary LDPC codes.
Given (unnormalized) distributions for the variable and check weights, return length n irregular LDPC codes whose degree distributions match the given distribution. The arguments Sv and Sc are sequences of real numbers, where the i-th entry indicates what percentage of the variable (resp. check) nodes should have weight i.Note that the distributions will not be matched perfectly unless everything is in complete balance.
Return the (3,6)-regular binary LDPC code of length 2(p3 - p) using the group-based construction of Margulis.
> C1 := RegularLDPCEnsemble(10, 2, 4);
> C2 := RegularLDPCEnsemble(10, 2, 4);
> C1 eq C2;
false
> LDPCMatrix(C1):Magma;
SparseMatrix(GF(2), 5, 10, [
4, 2,1, 3,1, 4,1, 6,1,
4, 1,1, 7,1, 9,1, 10,1,
4, 1,1, 2,1, 3,1, 7,1,
4, 5,1, 6,1, 8,1, 9,1,
4, 4,1, 5,1, 8,1, 10,1
])
> H := LDPCMatrix(C1);
> C3 := LDPCCode(H);
> C3 eq C1;
true
Since a code can have many different parity check matrices, the matrix which defines a code as being LDPC must be assigned specifically. Any parity check matrix can be assigned for this purpose, and once a code is assigned an LDPC matrix it is considered by Magma to be an LDPC code (regardless of the density or other properties of the matrix). The matrix must be of sparse type (MtrxSprs).
Return true if C is an LDPC code (which is true if it has been assigned an LDPC matrix).
Given a sparse matrix H which is a parity check matrix of the code C, assign H as the LDPC matrix of C.
Given an LDPC code C, return the sparse matrix which has been assigned as its low density parity check matrix.
Given an LDPC code C, return the density of the sparse matrix which has been assigned as its low density parity check matrix.
Returns true if C is an LDPC code and has regular column and row weights. If true, the row and column weights are also returned.
For an LDPC code C, return its Tanner graph. If there are n variables and m checks, then the graph has n + m nodes, the first n of which are the variable nodes.
For an LDPC code C, return the girth of its Tanner graph.
Return the theoretical rate of LDPC codes from the ensemble described by the given inputs.
> C := RandomLinearCode(GF(2),100,50); > IsLDPC(C); false > H := SparseMatrix(ParityCheckMatrix(C)); > H; Sparse matrix with 50 rows and 100 columns over GF(2) > AssignLDPCMatrix(~C, H); > IsLDPC(C); true > LDPCDensity(C); 0.253400000000000014122036873232The density of the parity check matrices of LDPC codes is much lower than that of randomly generated codes.
> C1 := RegularLDPCEnsemble(100,3,6); > C1:Minimal; [100, 50] Linear Code over GF(2) > LDPCDensity(C1); 0.0599999999999999977795539507497
The impressive performance of LDPC codes lies in their iterative decoding algorithm. Magma provides facilities to decode using LDPC codes, as well as simulating transmission over a binary symmetric or white Gaussian noise channels.
The binary symmetric channel transmits binary values and is defined by p < 0.5. Each individual bit independently sustains a "bit-flip" error with probability p.
The Gaussian channel is analog, transmitting real-values, and is defined by a standard deviation σ. Binary values are mapped to -1 and 1 before being transmitted. Each value independently sustains an errors which are normally distributed about 0 with standard deviation σ.
Channel: MonStgElt Default: "BinarySymmetric"
p: RngReSubElt Default: 0.1
StdDev: RngReSubElt Default: 0.25
Iterations: RngIntElt Default: Dimension(C)
For an LDPC code C and a received vector v, decode v to a codeword of C using the LDPC iterative decoding algorithm.The nature of the channel from which v is received is described by the variable argument Channel, which can either be the BinarySymmetric channel or the Gaussian channel. Errors on the binary symmetric channel is described by the argument p, while on the Gaussian channel they are described by StdDev.
The vector v must be over a ring corresponding to the channel which is selected. For the binary symmetric channel v must be a binary vector over F2, while for the Gaussian channel it must be real-valued.
Since the decoding algorithm is iterative and does not necessarily terminate on its own, a maximum number of iterations needs to be specified using the argument Iterations. The default value is much larger than would normally be used in practice, giving maximum error-correcting performance (at possibly some cost to efficiency).
> n := 500; > C := RegularLDPCEnsemble(n, 4, 8); > e := 5; > Errs := ; > repeat Include(~Errs, Random(1,n)); until #Errs eq e; > v := Random(C); > ev := AmbientSpace(C)![ (i in Errs) select 1 else 0 : i in [1..n]]; > rec_vec := v + ev; > time res := LDPCDecode(C, rec_vec : Channel:="BinarySymmetric", p:=0.2); Time: 0.000 > res eq v; true
To simulate a Gaussian channel requires obtaining normally distributed errors. This can be done (discretely) by generating a multiset of possible errors.
> sigma := 0.5;
> MaxE := 3.0;
> N := 100;
> V := [ MaxE*(i/N) : i in [-N div 2..N div 2]];
> E := [ 0.5*(1+Erf(x/(sigma*Sqrt(2)))) : x in V ];
> Dist := {* V[i]^^Round(1000*(E[i]-E[i-1])) : i in [2..#V]*};
A codeword of an LDPC code needs to be mapped
into the real domain.
> n := 500; > C := RegularLDPCEnsemble(n, 4, 8); > v := Random(C); > R := RealField(); > RS := RSpace(R, n); > vR := RS ! [ IsOne(v[i]) select 1 else -1 : i in [1..n]];Normally distributed errors are then introduced, and the received vector decoded.
> for i in [1..n] do > vR[i] +:= Random(Dist); > end for; > time res := LDPCDecode(C, vR : Channel:="Gaussian", StdDev:=sigma); Time: 0.000 > res eq v; true
Channel: MonStgElt Default: "BinarySymmetric"
p: RngReSubElt Default: 0.1
StdDev: RngReSubElt Default: 0.25
Iterations: RngIntElt Default: Dimension(C)
For an LDPC code C, simulate N transmissions across the given channel and return the accumulated bit error rate and word error rate.The variable arguments are as described for the function LDPCDecode. The channel which is described controls not only the way the decoding algorithm functions, but also the nature of the errors introduced during the simulation.
> C := RegularLDPCEnsemble(200, 3, 6); > LDPCSimulate(C, 10000 : Channel := "Gaussian", StdDev := 0.7); 0.00118249999999999995739519143001 0.00619999999999999978211873141731 > LDPCSimulate(C, 10000 : Channel := "Gaussian", StdDev := 0.75); 0.00749499999999999992617016886243 0.0410000000000000017208456881690 > LDPCSimulate(C, 10000 : Channel := "Gaussian", StdDev := 0.8); 0.0337220000000000019735324485737 0.159700000000000008615330671091 > LDPCSimulate(C, 10000 : Channel := "Gaussian", StdDev := 0.85); 0.0856174999999999991606713933834 0.370800000000000018474111129763 > LDPCSimulate(C, 10000 : Channel := "Gaussian", StdDev := 0.9); 0.162790999999999991265653420669 0.640499999999999958255614274094 > LDPCSimulate(C, 10000 : Channel := "Gaussian", StdDev := 0.95); 0.237657499999999993756105709508 0.840099999999999957900342906214 > LDPCSimulate(C, 10000 : Channel := "Gaussian", StdDev := 1.0); 0.296526500000000026169288958044 0.944799999999999973177011725056
The asymptotic performance of ensembles of LDPC codes can be determined using density evolution. An ensemble of LDPC codes (either regular or irregular) is defined by a pair of degree distributions, corresponding to the degrees at the variable and check nodes of the Tanner graph.
Over a specific channel, the critical parameter which defines the asymptotic performance of a given ensemble is its threshold, which is a value of the channel parameter (i.e., the probability of error p for the binary symmetric channel, and the standard deviation σ for the Gaussian channel). When the channel parameter is less than the threshold, asymptotically a code from the ensemble will decode with an error probability tending to zero. However, at any channel parameter above the threshold there will be a non-vanishing finite error probability.
Determining the threshold of an ensemble over the binary symmetric channel is relatively trivial, however over the real-valued Gaussian channel it can involve extensive computations. The speed depends heavily on the granularity of the discretization which is used, though this also affects the accuracy of the result.
The default settings of the Magma implementation use a reasonably coarse discretization, emphasizing speed over accuracy. These (still quite accurate) approximate results can then be used to help reduce the workload of calculations over finer discretizations if more accuracy is required.
Precision: RngReSubElt Default: 0.00005
Determines the threshold of the described ensemble of LDPC codes over the binary symmetric channel, which is the critical value of the channel parameter above which there is a non-vanishing error probability (asymptotically). The ensemble can either be defined by two integers for (v, c)-regular LDPC codes, or by two density distributions Sv and Sc, which are sequences of non-negative real numbers. The density distributions do not ned to be normalized, though the first entry (corresponding to weight 1 nodes in the Tanner graph) should always be zero.The computation proceeds by establishing lower and upper bounds on the threshold, then narrowing this range by repeatedly performing density evolution on the midpoint. The argument Precision controls the precision to which the threshold is desired.
Perform density evolution on the binary symmetric channel using channel parameter p and determine the asymptotic behaviour for the given LDPC ensemble. The return value is boolean, where true indicates that p is below the threshold and the ensemble has error probability asymptotically tending to zero.
> time LDPCBinarySymmetricThreshold(3, 6); 0.0394140625000000000000000000000 Time: 0.010 > time LDPCBinarySymmetricThreshold(4, 8); 0.0473925781250000000000000000001 Time: 0.110 > time LDPCBinarySymmetricThreshold(4, 10); 0.0368359375000000000000000000000 Time: 0.090
Lower: RngReSubElt Default: 0
Upper: RngReSubElt Default: ∞
Points: RngIntElt Default: 500
MaxLLR: RngReSubElt Default: 25
MaxIterations: RngIntElt Default: ∞
QuickCheck: BoolElt Default: true
Precision: RngReSubElt Default: 0.00005
Determines the threshold of the described ensemble of LDPC codes over the Gaussian channel, which is the critical value of the standard deviation above which there is a non-vanishing error probability (asymptotically). The ensemble can either be defined by two integers for (v, c)-regular LDPC codes, or by two density distributions Sv and Sc, which are sequences of non-negative real numbers. The density distributions do not ned to be normalized, though the first entry (corresponding to weight 1 nodes in the Tanner graph) should always be zero.The computation proceeds by establishing lower and upper bounds on the threshold, then narrowing this range by repeatedly performing density evolution on the midpoint. If the threshold is approximately known then manually setting tight Lower and Upper bounds can reduce the length of the calculation.
The speed with which these evolutions are computed depends on how fine the discretization is, controlled by the variable argument Points. If the threshold is needed to high levels of accuracy then an initial computation with fewer points is recommended to get a reduced searched range. The specific meaning of eachvariable argument is described below.
Lower and Upper are real-valued bounds on the threshold, which (if tight) can help to reduce the search range and speed up the threshold determination. The validity of an input bound is verified before the search begins, and an error is returned if it is incorrect.
Points and MaxLLR define the discretized basis of log likelihood ratios on which density evolution is performed, and have integer and real values resp. Specifically, the probability mass function is defined on the range [ -MaxLLR, ..., MaxLLR] on 2 * Points + 1 discretized points.
MaxIterations allows the user to set a finite limit of iterations that a density evolution should perform in determining the asymptotic behaviour at each channel parameter. Although this may help reduce the time of a computation, it should be kept in mind that the result may not be valid.
QuickCheck defines the method by which the asymptotic behaviour at each channel parameter is identified. If set to false, then the probability density must evolve all the way to within an infinitesimal value of unity. When set to true, if the rate of change of the probability density is seen to be successively increasing then the asymptotic behaviour is assumed to go to unity. Empirically this method seems to give accurate results and so the default behaviour is true, however it has no theoretical justification.
Precision is a real-valued parameter defining the precision to which the threshold should be determined.
Setting the verbose mode Code prints out the bounds on the threshold as subsequent density evolutions narrow the search range.
Points: RngIntElt Default: 500
MaxLLR: RngReSubElt Default: 25
MaxIterations: RngIntElt Default: ∞
QuickCheck: BoolElt Default: true
Perform density evolution on the Gaussian channel using standard deviation σ and determine the asymptotic behaviour for the given LDPC ensemble. The return value is boolean, where true indicates that σ is below the threshold and the ensemble has error probability asymptotically tending to zero.See the description of LDPCGaussianThreshold for a description of the variable arguments.
Access a small database of density distributions defining good irregular LDPC ensembles. Returned is the published threshold of the ensemble over the Gaussian channel, along with the variable and check degree distributions. The input i is a non-negative integer, which indexes the database (in no particular order).
In this example a published value of the threshold of an ensemble (obtained using a different implementation) can be compared to the outputs from different levels of discretization.
> thresh, Sv, Sc := GoodLDPCEnsemble(5); > R4 := RealField(4); > [R4| x : x in Sv]; [ 0.0000, 0.3001, 0.2839, 0.0000, 0.0000, 0.0000, 0.0000, 0.4159 ] > [R4| x : x in Sc]; [ 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.2292, 0.7708 ] > thresh; 0.949700000000000000000000000000 > time approx1 := LDPCGaussianThreshold(Sv, Sc: > Points := 500, Precision := 0.001); Time: 19.400 > approx1; 0.960976562499999964472863211995 > time approx2 := LDPCGaussianThreshold(Sv, Sc: > Points := 3000, > Lower := approx1-0.02, Upper := approx1+0.02); Time: 873.560 > approx2; 0.951210937499999964472863211995This estimate can now be used to narrow the search range of accurate density evolution. For the very long calculation the verbose mode is used to keep track of the progress of the calculation.
> SetVerbose("Code", true);
> time approx3 := LDPCGaussianThreshold(Sv, Sc:
> Points := 5000,
> Lower := approx2-0.005, Upper := approx2+0.0005);
Determining the mapping matrix...
...mapping matrix obtained 19.10s
Threshold Determination for LDPC code ensemble:
c: (6):0.229190 (7):0.770810
v: (2):0.300130 (3):0.283950 (8):0.415920
will be found to precision 0.000050
Max LLR is 25.000000 distributed across 5000 points
Asymptotic behaviour determination is: fast
Beginning search with lb = 0.946211, ub = 0.951711
New Bounds: lb = 0.948961, ub = 0.951711 114.41s
New Bounds: lb = 0.950336, ub = 0.951711 367.19s
New Bounds: lb = 0.950336, ub = 0.951023 553.95s
New Bounds: lb = 0.950336, ub = 0.950680 814.35s
New Bounds: lb = 0.950508, ub = 0.950680 1261.46s
New Bounds: lb = 0.950508, ub = 0.950594 1557.46s
New Bounds: lb = 0.950508, ub = 0.950551 1891.52s
Time: 2136.150
The thresholds given in the database are published values
taken from other implementations, and so are not
guaranteed to match up exactly with the values
obtained using Magma.