function h = ldpcenc(H)
%LDPCENC  LDPC encoder.
%   L = FEC.LDPCENC(H) constructs an LDPC encoder object L for a binary
%   systematic LDPC code with a parity-check matrix H.
%
%   H must be a sparse zero-one matrix. N and N-K are the number of columns
%   and the number of rows in H. The last N-K columns in H must be an
%   invertible matrix in GF(2).
%
%   An LDPC encoder object has the following properties. Only ParityCheckMatrix
%   is writable. All other properties are derived from it.
%
%   ParityCheckMatrix  - Parity-check matrix of the LDPC code, i.e., H.
%                        Stored as a sparse logical matrix.
%
%   BlockLength        - Total number of bits in a codeword, i.e., N.
%
%   NumInfoBits        - Number of information bits in a codeword, i.e., K.
%
%   NumParityBits      - Number of parity bits in a codeword, i.e., N-K.
%
%   EncodingAlgorithm  - Method for solving the parity-check equation to compute
%                        the parity bits using the information bits.
%
%                        'Forward Substitution' - if the last N-K columns
%                        in H are a lower triangular matrix.
%
%                        'Backward Substitution' - if the last N-K columns
%                        in H are an upper triangular matrix.
%
%                        'Matrix Inverse' - otherwise.
%
%   L = FEC.LDPCENC constructs an LDPC encoder object with a default
%   parity-check matrix (32400-by-64800) corresponding to an irregular LDPC
%   code with the following structure:
%
%     Row     No. of 1's per row        Column       No. of 1's per column
%   ------------------------------    --------------------------------------
%      1               6                  1 to 12960            8
%   2 to 32400         7              12961 to 32400            3
%
%   Columns 32401 to 64800 are a lower triangular matrix. Only the elements
%   on its main diagonal and the subdiagonal immediately below are 1's.
%   This LDPC code is used in conjunction with a BCH code in Digital Video
%   Broadcasting standard (DVB-S.2) to achieve a packet error rate below 10^-7
%   at about 0.7 dB to 1 dB from the Shannon limit.
%
%   Examples:
%
%     % Construct an LDPC encoder object
%     i = [1  3  2  4  1  2  3  3  4];   % row indices of 1's
%     j = [1  1  2  2  3  4  4  5  6];   % column indices of 1's
%     H = sparse(i,j,ones(length(i),1)); % This is just an example and is not
%     l = fec.ldpcenc(H);                % a good LDPC code.
%
%     % Construct an LDPC encoder object with the default parity-check matrix
%     l = fec.ldpcenc;
%
%   See also FEC.LDPCENC/ENCODE, FEC.LDPCDEC, FEC.LDPCDEC/DECODE, SPARSE.

%   @fec/@ldpcenc

%   Copyright 2006 The MathWorks, Inc.
%   $Revision: 1.1.6.2 $  $Date: 2006/06/23 19:24:46 $

error(nargchk(0, 1, nargin));

h = fec.ldpcenc;
if nargin == 0
    h.ParityCheckMatrix = DVBS2code(1/2);
else
    h.ParityCheckMatrix = H;
end

%%--------------------------------------------------------------------------
%% Local functions

function H = DVBS2code(R)

lenCodeWord = 64800;  % Length of codeword for DVB-S.2
NB = 360;  % Node indices parameter for DVB-S.2.

numInfoBits = lenCodeWord * R;
numParityBits = lenCodeWord - numInfoBits;

[ct1, ct2] = getchecknodetable(R);

ck1 = nodeindices(ct1, numParityBits, NB);
ck2 = nodeindices(ct2, numParityBits, NB);

d = [size(ck1,2) size(ck1,1) size(ck2,2) size(ck2,1) numParityBits-1 2 1 1];
r = [ck1(:); ck2(:); 0; reshape(ones(2,1)*(1:numParityBits-1),[],1)];
S = zeros(length(r),1);
numGroup = length(d)/2;
n = 0;
ncol = 1;
for i = 1:numGroup
    p = d(2*i-1)*d(2*i);
    S(n+1:n+p) = reshape(ones(d(2*i),1)*(ncol:ncol+d(2*i-1)-1),p,1);
    ncol = ncol + d(2*i-1);
    n = n + p;
end

% Parity check matrix (sparse) for DVB-S.2
H = sparse(double(r+1), S, 1);



%--------------------------------------------------------------------------
function ck = nodeindices(ct, M, NB)
% ct: check node table (single group)
% M: number of parity bits
% NB: block size
[N, D] = size(ct);
q = (M/NB);
b = (1:NB);
bq = (b-1).'*q;
ck = zeros(D, NB*N);
for r=1:N
    ck(:, NB*(r-1)+1:NB*r) = mod(addcr(bq, ct(r,:)), M)';
end


%--------------------------------------------------------------------------
function A = addcr(c, r)
M = length(c);
N = length(r);
A = zeros(M, N);
for m = 1:M
    A(m, :) = r + c(m);
end


%--------------------------------------------------------------------------
function [ct1, ct2] = getchecknodetable(R)
% Returns check node table for LDPC coding in DVB.S2 standard.
% R = code rate.
switch R
    case 1/4
        %Table B.4: rate 1/4 (nldpc = 64 800)
        ct1 = [23606 36098 1140 28859 18148 18510 6226 540 42014 20879 23802 47088
            16419 24928 16609 17248 7693 24997 42587 16858 34921 21042 37024 20692
            1874 40094 18704 14474 14004 11519 13106 28826 38669 22363 30255 31105
            22254 40564 22645 22532 6134 9176 39998 23892 8937 15608 16854 31009
            8037 40401 13550 19526 41902 28782 13304 32796 24679 27140 45980 10021
            40540 44498 13911 22435 32701 18405 39929 25521 12497 9851 39223 34823
            15233 45333 5041 44979 45710 42150 19416 1892 23121 15860 8832 10308
            10468 44296 3611 1480 37581 32254 13817 6883 32892 40258 46538 11940
            6705 21634 28150 43757 895 6547 20970 28914 30117 25736 41734 11392
            22002 5739 27210 27828 34192 37992 10915 6998 3824 42130 4494 35739
            8515 1191 13642 30950 25943 12673 16726 34261 31828 3340 8747 39225
            18979 17058 43130 4246 4793 44030 19454 29511 47929 15174 24333 19354
            16694 8381 29642 46516 32224 26344 9405 18292 12437 27316 35466 41992
            15642 5871 46489 26723 23396 7257 8974 3156 37420 44823 35423 13541
            42858 32008 41282 38773 26570 2702 27260 46974 1469 20887 27426 38553];
        ct2 = [22152 24261 8297
            19347 9978 27802
            34991 6354 33561
            29782 30875 29523
            9278 48512 14349
            38061 4165 43878
            8548 33172 34410
            22535 28811 23950
            20439 4027 24186
            38618 8187 30947
            35538 43880 21459
            7091 45616 15063
            5505 9315 21908
            36046 32914 11836
            7304 39782 33721
            16905 29962 12980
            11171 23709 22460
            34541 9937 44500
            14035 47316 8815
            15057 45482 24461
            30518 36877 879
            7583 13364 24332
            448 27056 4682
            12083 31378 21670
            1159 18031 2221
            17028 38715 9350
            17343 24530 29574
            46128 31039 32818
            20373 36967 18345
            46685 20622 32806];
    case 1/3
        %Table B.4: rate 1/3 (nldpc = 64 800)
        ct1 = [34903 20927 32093 1052 25611 16093 16454 5520 506 37399 18518 21120
            11636 14594 22158 14763 15333 6838 22222 37856 14985 31041 18704 32910
            17449 1665 35639 16624 12867 12449 10241 11650 25622 34372 19878 26894
            29235 19780 36056 20129 20029 5457 8157 35554 21237 7943 13873 14980
            9912 7143 35911 12043 17360 37253 25588 11827 29152 21936 24125 40870
            40701 36035 39556 12366 19946 29072 16365 35495 22686 11106 8756 34863
            19165 15702 13536 40238 4465 40034 40590 37540 17162 1712 20577 14138
            31338 19342 9301 39375 3211 1316 33409 28670 12282 6118 29236 35787
            11504 30506 19558 5100 24188 24738 30397 33775 9699 6215 3397 37451
            34689 23126 7571 1058 12127 27518 23064 11265 14867 30451 28289 2966
            11660 15334 16867 15160 38343 3778 4265 39139 17293 26229 42604 13486
            31497 1365 14828 7453 26350 41346 28643 23421 8354 16255 11055 24279
            15687 12467 13906 5215 41328 23755 20800 6447 7970 2803 33262 39843
            5363 22469 38091 28457 36696 34471 23619 2404 24229 41754 1297 18563
            3673 39070 14480 30279 37483 7580 29519 30519 39831 20252 18132 20010
            34386 7252 27526 12950 6875 43020 31566 39069 18985 15541 40020 16715
            1721 37332 39953 17430 32134 29162 10490 12971 28581 29331 6489 35383
            736 7022 42349 8783 6767 11871 21675 10325 11548 25978 431 24085
            1925 10602 28585 12170 15156 34404 8351 13273 20208 5800 15367 21764
            16279 37832 34792 21250 34192 7406 41488 18346 29227 26127 25493 7048];
        ct2 = [39948 28229 24899
            17408 14274 38993
            38774 15968 28459
            41404 27249 27425
            41229 6082 43114
            13957 4979 40654
            3093 3438 34992
            34082 6172 28760
            42210 34141 41021
            14705 17783 10134
            41755 39884 22773
            14615 15593 1642
            29111 37061 39860
            9579 33552 633
            12951 21137 39608
            38244 27361 29417
            2939 10172 36479
            29094 5357 19224
            9562 24436 28637
            40177 2326 13504
            6834 21583 42516
            40651 42810 25709
            31557 32138 38142
            18624 41867 39296
            37560 14295 16245
            6821 21679 31570
            25339 25083 22081
            8047 697 35268
            9884 17073 19995
            26848 35245 8390
            18658 16134 14807
            12201 32944 5035
            25236 1216 38986
            42994 24782 8681
            28321 4932 34249
            4107 29382 32124
            22157 2624 14468
            38788 27081 7936
            4368 26148 10578
            25353 4122 39751];      
    case 2/5
        %Table B.4: rate 2/5 (nldpc = 64 800)
        ct1 = [31413 18834 28884 947 23050 14484 14809 4968 455 33659 16666 19008
            13172 19939 13354 13719 6132 20086 34040 13442 27958 16813 29619 16553
            1499 32075 14962 11578 11204 9217 10485 23062 30936 17892 24204 24885
            32490 18086 18007 4957 7285 32073 19038 7152 12486 13483 24808 21759
            32321 10839 15620 33521 23030 10646 26236 19744 21713 36784 8016 12869
            35597 11129 17948 26160 14729 31943 20416 10000 7882 31380 27858 33356
            14125 12131 36199 4058 35992 36594 33698 15475 1566 18498 12725 7067
            17406 8372 35437 2888 1184 30068 25802 11056 5507 26313 32205 37232
            15254 5365 17308 22519 35009 718 5240 16778 23131 24092 20587 33385
            27455 17602 4590 21767 22266 27357 30400 8732 5596 3060 33703 3596
            6882 873 10997 24738 20770 10067 13379 27409 25463 2673 6998 31378
            15181 13645 34501 3393 3840 35227 15562 23615 38342 12139 19471 15483
            13350 6707 23709 37204 25778 21082 7511 14588 10010 21854 28375 33591
            12514 4695 37190 21379 18723 5802 7182 2529 29936 35860 28338 10835
            34283 25610 33026 31017 21259 2165 21807 37578 1175 16710 21939 30841
            27292 33730 6836 26476 27539 35784 18245 16394 17939 23094 19216 17432
            11655 6183 38708 28408 35157 17089 13998 36029 15052 16617 5638 36464
            15693 28923 26245 9432 11675 25720 26405 5838 31851 26898 8090 37037
            24418 27583 7959 35562 37771 17784 11382 11156 37855 7073 21685 34515
            10977 13633 30969 7516 11943 18199 5231 13825 19589 23661 11150 35602
            19124 30774 6670 37344 16510 26317 23518 22957 6348 34069 8845 20175
            34985 14441 25668 4116 3019 21049 37308 24551 24727 20104 24850 12114
            38187 28527 13108 13985 1425 21477 30807 8613 26241 33368 35913 32477
            5903 34390 24641 26556 23007 27305 38247 2621 9122 32806 21554 18685];
        ct2 = [17287 27292 19033
            25796 31795 12152
            12184 35088 31226
            38263 33386 24892
            23114 37995 29796
            34336 10551 36245
            35407 175 7203
            14654 38201 22605
            28404 6595 1018
            19932 3524 29305
            31749 20247 8128
            18026 36357 26735
            7543 29767 13588
            13333 25965 8463
            14504 36796 19710
            4528 25299 7318
            35091 25550 14798
            7824 215 1248
            30848 5362 17291
            28932 30249 27073
            13062 2103 16206
            7129 32062 19612
            9512 21936 38833
            35849 33754 23450
            18705 28656 18111
            22749 27456 32187
            28229 31684 30160
            15293 8483 28002
            14880 13334 12584
            28646 2558 19687
            6259 4499 26336
            11952 28386 8405
            10609 961 7582
            10423 13191 26818
            15922 36654 21450
            10492 1532 1205
            30551 36482 22153
            5156 11330 34243
            28616 35369 13322
            8962 1485 21186
            23541 17445 35561
            33133 11593 19895
            33917 7863 33651
            20063 28331 10702
            13195 21107 21859
            4364 31137 4804
            5585 2037 4830
            30672 16927 14800];
    case 1/2
        ct1 = [54 9318 14392 27561 26909 10219 2534 8597
            55 7263 4635 2530 28130 3033 23830 3651
            56 24731 23583 26036 17299 5750 792 9169
            57 5811 26154 18653 11551 15447 13685 16264
            58 12610 11347 28768 2792 3174 29371 12997
            59 16789 16018 21449 6165 21202 15850 3186
            60 31016 21449 17618 6213 12166 8334 18212
            61 22836 14213 11327 5896 718 11727 9308
            62 2091 24941 29966 23634 9013 15587 5444
            63 22207 3983 16904 28534 21415 27524 25912
            64 25687 4501 22193 14665 14798 16158 5491
            65 4520 17094 23397 4264 22370 16941 21526
            66 10490 6182 32370 9597 30841 25954 2762
            67 22120 22865 29870 15147 13668 14955 19235
            68 6689 18408 18346 9918 25746 5443 20645
            69 29982 12529 13858 4746 30370 10023 24828
            70 1262 28032 29888 13063 24033 21951 7863
            71 6594 29642 31451 14831 9509 9335 31552
            72 1358 6454 16633 20354 24598 624 5265
            73 19529 295 18011 3080 13364 8032 15323
            74 11981 1510 7960 21462 9129 11370 25741
            75 9276 29656 4543 30699 20646 21921 28050
            76 15975 25634 5520 31119 13715 21949 19605
            77 18688 4608 31755 30165 13103 10706 29224
            78 21514 23117 12245 26035 31656 25631 30699
            79 9674 24966 31285 29908 17042 24588 31857
            80 21856 27777 29919 27000 14897 11409 7122
            81 29773 23310 263 4877 28622 20545 22092
            82 15605 5651 21864 3967 14419 22757 15896
            83 30145 1759 10139 29223 26086 10556 5098
            84 18815 16575 2936 24457 26738 6030 505