% Routine to check the performance of LDPC codes

clear;

% AWGN channel: x is the transmitted vector, y is received.
% y = exp(2*n/ sigman^2); sigman^2=N0/2; 10*log(Eb/N0) = SNR

% dc = [6;1]
% dv = [3;1]
% M = 1000;
% N = 2000;

MAX_ITERATIONS = 10;

SNR = 1.5;
%SNR = 3.5:0.5:6;

%disp('Making LDPC code');
%H = makeLDPC(M,N,dc,dv,noCycleTries,maxTooLong);
%save 'Hmatrix.mat' H

%load 'H2000.mat' H;  % (N=2000, K=1000) Rate = 1/2  dc=[6;1]; dv=[3;1]; 4cycles = no
%load 'H3000.mat' H;  % (N=3000, K=2000) Rate = 1/3 dc=[9;1]; dv=[6;1]; 4cycles = no
%load 'H3000b.mat' H;  % (N=3000, K=2000) Rate = 1/3 dc=[9;1]; dv=[6;1]; 4cycles = no
%load 'H6000b.mat' H;  % (N=6000, K=3000) Rate = 1/2 dc=[6;1]; dv=[3;1]; 4cycles = no
%load 'H2000bad.mat' H; % (N=2000, K=1000) Rate = 1/2 dc=[6;1]; dv=[3;1]; 4cycles = yes
load 'H-EG2D-4.mat' H;

[M N]=size(H);
K = N-M;

% Form distributions and interleaver
numedges = 0; 
for i=1:M
    numedges = numedges + nnz(H(i,:));
    % Compute the check node distribution explicitly.
    explicitCnDist(i) = length(nonzeros(H(i,:)));
end;

% Compute variable node distribution explicitly.
for i=1:N
    explicitVnDist(i) = length(nonzeros(H(:,i)));
end;

% compute interleaver from the parity check matrix
cumCnDist = cumsum(explicitCnDist);
rowPosition = zeros(1,M);
index = 1;
for i=1:N
    nzVindex = find(H(:,i));
    if length(nzVindex)==0
        pause;
    end;
    for j=1:length(nzVindex)
        if nzVindex(j)==1
            interleaver(index)= rowPosition(nzVindex(j))+1;
        else
            interleaver(index)=cumCnDist(nzVindex(j)-1)+(1+rowPosition(nzVindex(j)));
        end;
        index = index + 1;
        rowPosition(nzVindex(j))=rowPosition(nzVindex(j))+1;
    end;
end;

fp = fopen('resultEG2D4-10.dat','w');
fprintf(fp, 'SNR\tN_ITERATIONS\tNBITS\t  BERRORS\t  FERRORS  \tBER\t\tFER\n');

% assume that all zero vector has been transmitted. 0->+1, 1->-1

if MAX_ITERATIONS==0
    N=N-M;              % uncoded
end;

c=2*zeros(1,N)-1;

Pe = zeros(1,length(SNR));
bth = 10;
maxblocks = 10000;

for s=1:length(SNR)
    disp('SNR');
    disp(SNR(s));
%     disp('Blocks...');
    
    berrors = 0;
    ferrors = 0;
    blocks = 1;
    N0 = 1/10^(SNR(s)/10); % Eb is 1;
    
%     if SNR(s)>=2.0
%         bth=50;
%     end;
%     if SNR(s)>=3.0
%         bth=25;
%     end;
%     if SNR(s)>=4.0
%         bth=10;
%     end;
    
    i=1;
    while(i<=blocks)
        %disp(blocks);
        n=sqrt(N0/2)*randn(1,N);      % should be N0/2
        
        x = c + n;
        
        y = 2*x/(N0/2);               % should be N0/2
        Lap = 0;
        if MAX_ITERATIONS==0
            diff=length(find(x>0)); % uncoded
        else
            %interleaver = randperm(length(interleaver));            % Random H matrix
            [ld, lc] = ldpcDecode(y, interleaver, explicitCnDist, explicitVnDist, MAX_ITERATIONS, Lap);
            diff = length(find(ld>0));
        end;
        
        if diff>0
            ferrors = ferrors + 1;
        end;
        berrors = berrors + diff;
        if berrors > bth
            disp('Bit Errors.');
            disp(berrors);
            break;
        else
            blocks = blocks + 1;
        end;
        if blocks > maxblocks
            break;
        end;
    i=i+1;    
    end;
    Pe(s)=berrors/(blocks*N);
    disp(Pe(s));
    fprintf(fp, '%1.1f\t%12d\t%7d%7d\t%7d\t%1.7e\t%1.7e\n\n',SNR(s),MAX_ITERATIONS,blocks*N,berrors,ferrors,berrors/(blocks*N),ferrors/blocks);
%     if Pe(s) < 4e-05
%         break;
    end;
end;

fclose(fp);