clear all
close all
clc
tic
rows=512; 
cols=1024; 
H=gen_ldpc(rows,cols);%generates a 100 x 200 H matrix( the IG parity 
                      %check matrix) 
[newh,rearranged_cols]=rearrange_cols(H); 
% newh holds the rearranged version of H and rearranged_cols holds 
%the rearrangement order of columns 
%%%%%%%%%%%%Get A and B matrices %%%%%%%%%%%%%%%% 
for i=1:rows 
    for j=1:rows 
        A(i,j)=newh(i,j); 
    end 
end 
for i=1:rows 
    for j=rows+1:cols 
        B(i,j-rows)=newh(i,j); 
    end 
end 

%%%%%%%%%%%%%%message bits%%%%%%%%%%%%%%%%%%%%% 

%generate message bits randomly 
%for j=1:10%simulate 1000 flame
for i=1:cols-rows 
    s(i)=round(rand); 
end 
s=s';
% make it a row vector 
%%%%%%% or load s from file %%%%type "help load" at matlab prompt for details 

%%%%%%%%%%%%%%get check bits( Eqn.5)%%%%%%%%%%%%%%%% 

d=mul_GF2(inv_GF2(A),B); 
c=mul_GF2(d,s); 
%%%%%%%%%%%%%%%%%form u'%%%%%%%%%%%%% 
u1=c; 
for i=rows+1:cols 
    u1(i)=s(i-rows); 
end 
%%%%%%%%%%%%%%% reorder u?to get the CW u %%%%%%%%%%%%%%%%%%% 
u=reorder_bits(u1,rearranged_cols); 



ind=find(H==1);
[r,c]=ind2sub(size(H),ind);
[rows,cols]=size(H);
h=sparse(H);                         % for use with Matlab-based LDPC Decoder
n=cols;
k=n-rows;


% Find 
% 1: maximum check degree
% 2: column indeces in each row which contain '1'
% 3: maximum variable degree
% 4: find column indeces in each row which contain '1'
[max_check_degree,check_node_ones,BIGVALUE_COLS,max_variable_degree,variable_node_ones,BIGVALUE_ROWS]=one_finder(H);



rand('seed',584);
randn('seed',843);

dB=[0:3];                           % range of SNR values in dB
SNRpbit=10.^(dB/10);                % Eb/No conversion from dB to decimal
No_uncoded=1./SNRpbit;              % since Eb=1
R=k/n;                              % code rate 
No=No_uncoded./R;



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
max_iter=5;                   %maximum number of decoder iterations
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
maximum_blockerror=30;          % maximum blockerrors per SNR point
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

blockerrors=0;
biterrors=0;
block=0;
FER=zeros(1,length(dB));        % array for Frame Error Rate
BER=zeros(1,length(dB));        % array for Channel Error Rate
block_array=zeros(1,length(dB));

for z=1:length(SNRpbit), % loop for testing over range of SNR values

    biterrors=0;
    blockerrors=0;
    block=0;

    while(blockerrors<maximum_blockerror)  %while loop

       
        tx_waveform=bpsk(u,1);

        sigma=sqrt(No(z)/2);
        rx_waveform=tx_waveform + sigma*randn(1,length(tx_waveform));
        gamma_n=(4/No(z))*rx_waveform;
       [vhat,iteration]=decode_ldpc_matlab(rx_waveform,No(z),h,rows,cols,ind,r,c,max_iter);
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

        
        Errors=zeros(1,length(vhat));
        Errors(find(u'~=vhat))=1;
      
        if sum(Errors)~=0
            blockerrors=blockerrors+1;
            
            if rem(blockerrors,5)==0; %save statistics after every 5 blockerrors
                
                ber=biterrors/(block*n);
                fer=blockerrors/block;
                
                s=sprintf('%dx%d_regular_results.txt',n-k,n);
                fid = fopen(s,'a');
                fprintf(fid,'\n');
                fprintf(fid,'%s %2.1EdB\n','SNR =',dB(z));
                fprintf(fid,'%s %6.3E\n','BER =',ber);
                fprintf(fid,'%s %6.3E\n','FER =',fer);
                fprintf(fid,'%s %d\n','blocks =',block);
                fprintf(fid,'%s %d\n','blockerrors =',blockerrors);
                fprintf(fid,'%s %d\n','biterrors =',biterrors);
                fclose(fid);
            end
            
            
        end
        
        biterrors=biterrors+sum(Errors);
        block=block+1;

    end     %while loop

    block_array(z)=block; %counting blocks per SNR point
    BER(z)=biterrors/(block*n);
    FER(z)=blockerrors/block;
    fprintf(1,'\n\n Simulation finished for SNR: %d \n',dB(z))
    
%     save results.mat dB BER FER block_array
    
end      % loop for testing over range of SNR values



semilogy(dB,BER,'b-o')
title('Bit Error Rate')
ylabel('BER')
xlabel('Eb/No (dB)')
grid
hold on

rows=512; 
cols=1024; 
H=gen_ldpc(rows,cols);%generates a 100 x 200 H matrix( the IG parity 
                      %check matrix) 
[newh,rearranged_cols]=rearrange_cols(H); 
% newh holds the rearranged version of H and rearranged_cols holds 
%the rearrangement order of columns 
%%%%%%%%%%%%Get A and B matrices %%%%%%%%%%%%%%%% 
for i=1:rows 
    for j=1:rows 
        A(i,j)=newh(i,j); 
    end 
end 
for i=1:rows 
    for j=rows+1:cols 
        B(i,j-rows)=newh(i,j); 
    end 
end 

%%%%%%%%%%%%%%message bits%%%%%%%%%%%%%%%%%%%%% 

%generate message bits randomly 
%for j=1:10%simulate 1000 flame
for i=1:cols-rows 
    s(i)=round(rand); 
end 
s=s';
% make it a row vector 
%%%%%%% or load s from file %%%%type "help load" at matlab prompt for details 

%%%%%%%%%%%%%%get check bits( Eqn.5)%%%%%%%%%%%%%%%% 

d=mul_GF2(inv_GF2(A),B); 
c=mul_GF2(d,s); 
%%%%%%%%%%%%%%%%%form u'%%%%%%%%%%%%% 
u1=c; 
for i=rows+1:cols 
    u1(i)=s(i-rows); 
end 
%%%%%%%%%%%%%%% reorder u?to get the CW u %%%%%%%%%%%%%%%%%%% 
u=reorder_bits(u1,rearranged_cols); 



ind=find(H==1);
[r,c]=ind2sub(size(H),ind);
[rows,cols]=size(H);
h=sparse(H);                         % for use with Matlab-based LDPC Decoder
n=cols;
k=n-rows;


% Find 
% 1: maximum check degree
% 2: column indeces in each row which contain '1'
% 3: maximum variable degree
% 4: find column indeces in each row which contain '1'
[max_check_degree,check_node_ones,BIGVALUE_COLS,max_variable_degree,variable_node_ones,BIGVALUE_ROWS]=one_finder(H);



rand('seed',584);
randn('seed',843);

dB=[0:3];                           % range of SNR values in dB
SNRpbit=10.^(dB/10);                % Eb/No conversion from dB to decimal
No_uncoded=1./SNRpbit;              % since Eb=1
R=k/n;                              % code rate 
No=No_uncoded./R;



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
max_iter=10;                   %maximum number of decoder iterations
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
maximum_blockerror=30;          % maximum blockerrors per SNR point
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

blockerrors=0;
biterrors=0;
block=0;
FER=zeros(1,length(dB));        % array for Frame Error Rate
BER=zeros(1,length(dB));        % array for Channel Error Rate
block_array=zeros(1,length(dB));

for z=1:length(SNRpbit), % loop for testing over range of SNR values

    biterrors=0;
    blockerrors=0;
    block=0;

    while(blockerrors<maximum_blockerror)  %while loop

       
        tx_waveform=bpsk(u,1);

        sigma=sqrt(No(z)/2);
        rx_waveform=tx_waveform + sigma*randn(1,length(tx_waveform));
        gamma_n=(4/No(z))*rx_waveform;
       [vhat,iteration]=decode_ldpc_matlab(rx_waveform,No(z),h,rows,cols,ind,r,c,max_iter);
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

        
        Errors=zeros(1,length(vhat));
        Errors(find(u'~=vhat))=1;
      
        if sum(Errors)~=0
            blockerrors=blockerrors+1;
            
            if rem(blockerrors,5)==0; %save statistics after every 5 blockerrors
                
                ber=biterrors/(block*n);
                fer=blockerrors/block;
                
                s=sprintf('%dx%d_regular_results.txt',n-k,n);
                fid = fopen(s,'a');
                fprintf(fid,'\n');
                fprintf(fid,'%s %2.1EdB\n','SNR =',dB(z));
                fprintf(fid,'%s %6.3E\n','BER =',ber);
                fprintf(fid,'%s %6.3E\n','FER =',fer);
                fprintf(fid,'%s %d\n','blocks =',block);
                fprintf(fid,'%s %d\n','blockerrors =',blockerrors);
                fprintf(fid,'%s %d\n','biterrors =',biterrors);
                fclose(fid);
            end
            
            
        end
        
        biterrors=biterrors+sum(Errors);
        block=block+1;

    end     %while loop

    block_array(z)=block; %counting blocks per SNR point
    BER(z)=biterrors/(block*n);
    FER(z)=blockerrors/block;
    fprintf(1,'\n\n Simulation finished for SNR: %d \n',dB(z))
    
%     save results.mat dB BER FER block_array
    
end      % loop for testing over range of SNR values



semilogy(dB,BER,'b-o')
title('Bit Error Rate')
ylabel('BER')
xlabel('Eb/No (dB)')
grid