% 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);
D80