clc; 
clear all;
% LDPC matrix size, rate must be 1/2
% Warning: encoding - decoding can be very long for large LDPC matrix!
M = 10000;
N = 20000;   
% Method for creating LDPC matrix (0 = Evencol; 1 = Evenboth)
method = 1;

% Eliminate length-4 cycle
noCycle = 1;

% Number of 1s per column for LDPC matrix
onePerCol = 3;

% LDPC matrix reorder strategy (0 = First; 1 = Mincol; 2 = Minprod)
strategy = 2;

% Number of iteration;
iter = 5;
% Make the LDPC matrix
H = makeLdpc(M, N, 1, 1, onePerCol);
% Make random data (0/1)
dSource = round(rand(M, 1));
% Encoding message
[c, newH] = makeParityChk(dSource(:, 1), H, strategy);
u = [c; dSource(:, 1)];
ldpc_code=u;
   
%Lsnr, Usnr : Lower and Upper SNR bound in dB
%nsymb : Number of symbols to transmit
%modtype : 'psk', 'qask' types of modulation one can use.
%Mary : M-ary modulation
%stbc : 
%       G2  k=2, p=2, R=1       (Two Tx antennas; Alamouti's code)
%Rate is defined as: R=k/p (Tx k symbols / p time slots)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%Define all User Parameters here(STBC)
Lsnr=1;          %Lower SNR level
Usnr=10;          %Upper SRN level
nsymb=N;         %number of symbols to transmit
Mary = 4;        %Constellation size, Mary = 2^b where b is # of bits per symbol
modtype = 'psk'; %'psk' or 'qask'
stbc = 'G2';     
nRx = 2;         %Number of receive antennas to use
dobiterr = 'yes';%'yes' or 'no'; work with bit or symbol error respectively
%%%%%%%%
tic
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Save file name

fname = sprintf('stbc_%s_%d%s_%dRx',stbc,Mary,modtype,nRx);
fid = fopen([fname '.txt'],'w');
fprintf(fid,'Starting receiver simulation : \n');
fclose(fid);

sigStr=['code' stbc];
err1=ones(Usnr-Lsnr+1,1); %Define error vector
err2=ones(Usnr-Lsnr+1,1);
index=0;
for snr=Lsnr:Usnr,
    %Write to file to show routines' progress
    fid = fopen([fname '.txt'],'a');
    fprintf(fid,'Lower SNR: %d, Current SNR: %d, Upper SNR: %d, Elapsed Time: %6.0f \n',Lsnr,snr,Usnr,toc);
    fclose(fid);
   
    index=index+1;
    
    switch dobiterr
    case 'no'
        abscissa = 'SNR / (dB)';
        ordinate = 'Symbol Error Probability';
    case 'yes', 
        K=log2(Mary);
        snr=snr + 10*log10(K);
        abscissa = 'SNR / (dB)';
        ordinate = 'Bit Error Probability';
    end
    
    switch sigStr
    
    case 'codeG2'
        %Call routine for code G2 2-by-nRx
        nsymb=nsymb - mod(nsymb,2);
       [out11,in11]=STBCcodeG2(index,snr,nsymb,Mary,modtype,nRx,err2,ldpc_code);
        Rate = 1;
    
    end%switch
    % Decoding 
    tx=out11;
   
    vhat2 = decodeLDPC(tx, newH, iter);
    
    err1(index)=symerr(u,out11); 
    err2(index)=symerr(u,vhat2'); 
end%snr

save(fname)

fid = fopen([fname '.txt'],'a');
fprintf(fid,'FINISHED!!!   Total Elapsed Time: %6.0f\n',toc)    
fclose(fid);
%Plot the Results
figure(4)
speff = Rate*log2(Mary);
str = ['%d bps/Hz, %d receive antenna, Code ' stbc ' , %d ' modtype ' modulation'];
tle = [sprintf(str,speff,nRx,Mary)];
%semilogy(Lsnr:Usnr,err/nsymb,'ko-');grid
semilogy(Lsnr:Usnr,err1/nsymb, 'o-');
hold;
semilogy(Lsnr:Usnr,err2/nsymb, 'o--');
grid on;
hold off;
xlabel(abscissa);ylabel(ordinate);
title(tle)
legend('without LDPC','having LDPC',0); 
fprintf('Finished.\n')