www.pudn.com > ldpc_gallager.rar > gallager1.m

%Gallager algorithm for LDPC

close all;
clear all;



randn('state',sum(100*clock));

load gallager1_data.mat H GT N0 K0 M0 rankA CodeRate




if rankA~=M0
    disp('rank of H should equal to the number rows of H.');    
    return; 
end    

n=N0;
k=K0;
m=M0;
clear N0 K0 M0


[childrenList,parentsList,children,parents]=RelationGraph(H);%Get the relation graph of the given LDPC

[graph,startGraph,endGraph,lengthGraph]=GetGraph(childrenList,parentsList,children,parents);


[rowsGraph,columnsGraph]=size(graph);


SourceCount=0;
SourceBlockLength=k;
ErrorCount=0;
TestLength=10000000;

maxIterativeNumber=100;

ErrorCount1=0;

%Intialize matrix
Alpha=zeros(n,1);
Beta=zeros(n,1);
F_Beta=zeros(n,1);
LLR=zeros(n,1);

Alpha1=zeros(n,1);
Beta1=zeros(n,1);
F_Beta1=zeros(n,1);
LLR1=zeros(n,1);

yDe_maxSNR=zeros(n,1);

NoiseVector=zeros(n,1);

A=1.0;

Eb=A*A/CodeRate;


EbN0=2; %EbN0 dB



N0=Eb*10^(-EbN0/10);
NoiseVar=N0/2;
NoiseRoot=NoiseVar^0.5;

BlockCount=0;
IgnoreCount=1;


r=zeros(n,1);

%Generate a sequence for A and noise variance estimate
s=(randn(k,1)>0); %source information sequence with length K0
t=mod(GT*s,2);

NoiseVector=randn(n,1)*NoiseRoot;
r=A*(1-2*t)+NoiseVector;

%Estimate A and Noise Variance
r1=abs(r);
A1=mean(r1);
NoiseVar1=var(r1);

while(1)
   
    %******* Coding ********
    s=(randn(k,1)>0); %source information sequence with length k
    
    % s=zeros(k,1); % this row should be removed
    
    
    t=mod(GT*s,2);
    
    %******** AWGN Channel ********
    NoiseVector=randn(n,1)*NoiseRoot;
    r=A*(1-2*t)+NoiseVector;
    
    %************* Decoding ******
    valueFac=2*A1/NoiseVar1; %Initialize LLR from channel
    
    %valueFac=2*A/NoiseVar;
           
    LLR=valueFac*r;
    
    
    Alpha=sign(LLR);
    Beta=abs(LLR);
    
    F_Beta=(exp(Beta)-1)./(exp(Beta)+1);
    indexZero=find(F_Beta<=eps);
    F_Beta(indexZero)=F_Beta(indexZero)+eps;
    F_Beta=-log(F_Beta);
    
    
    LLR1=LLR;

    
    maxSNR_LLR=0;
    minDe=m;
    %******** Iterative unit (start)****************
    for iterativeCount=1:1:maxIterativeNumber
    
        Alpha1=Alpha;
        Beta1=Beta;
        F_Beta1=F_Beta;
        LLR1=LLR;
        
                
        for u=1:1:n % u presents current the digit location
            currentLength=lengthGraph(u);
            if (0==currentLength)
                continue;
            end
            currentStart=startGraph(u);
            currentEnd=endGraph(u);
        
        
            currentRowLLR=0;
            for v=currentStart:1:currentEnd % v presents current sub graph for digit location u 
                sumF=0;
                sumCount=0;
                currentSign=1;
                for w=1:1:columnsGraph
                    currentIndex=graph(v,w);
                    if(0==currentIndex)
                        break;
                    end
                
                    sumF=sumF+F_Beta(currentIndex);
                    currentSign=currentSign*Alpha(currentIndex);
                    sumCount=sumCount+1;
                end
                if(0==sumCount) % current row are all zero elements
                    continue;
                end
            
                F_sumF=(exp(sumF)-1)/(exp(sumF)+1);
                indexZero=find(F_sumF<=eps);
                F_sumF(indexZero)=F_sumF(indexZero)+eps;
                F_sumF=-log(F_sumF);
                currentRowLLR=currentRowLLR+currentSign*F_sumF;
            
            end
            LLR1(u)=LLR(u)+currentRowLLR;
        end
        
        
        meanLLR=mean(abs(LLR1));
        
        varLLR=var(abs(LLR1));
        varLLR=varLLR+eps;
        LLR1=(2*meanLLR/varLLR)*LLR1; %Look it as a sequence polluted by 
                                       %Gaussian noise in the iterative process
                
        LLR=LLR1;
        Alpha=sign(LLR);
        Beta=abs(LLR);

        snrLLR=(meanLLR*meanLLR)/varLLR
        if(snrLLR>=80) %current LLR SNR is very high
            break;
        end
        iterativeCount
        
        
        %Tentative Decoding
        yDe=(LLR<=0);
        zDe=mod(H*yDe,2);
        valueDe=sum(zDe)
        if(valueDe==0)
            disp(iterativeCount)
            break;
        end
        
        if(isnan(snrLLR) & (iterativeCount==1)) 
            disp('snrLLR is NaN');
            break;
        end
        if(snrLLR>maxSNR_LLR)
        % if(minDe>zDe)   
            maxSNR_LLR=snrLLR;
            %minDe=zDe;
            yDe_maxSNR=yDe;
        end
        if(iterativeCount==maxIterativeNumber)
            yDe=yDe_maxSNR;    
        end
        
        
        
        
        
        F_Beta=(exp(Beta)-1)./(exp(Beta)+1);
        indexZero=find(F_Beta<=eps);
        F_Beta(indexZero)=F_Beta(indexZero)+eps;
        F_Beta=-log(F_Beta);  
        
    end

    %***************Iterative unit (end) ***********
    
    
    
    
    
    
    
       
    %Re-estimate A1 and the noise variance
    
        
    decoderOutput=yDe(1:k);

    decoderDirect=(r(1:k)<=0);
    
    if(BlockCount>IgnoreCount)
        
        errorVector=xor(s,decoderOutput);
        errornum=sum(errorVector);
        ErrorCount=ErrorCount+errornum;

        if ((errornum>0) & (zDe==0) )
            disp('Low weight code encounted...');
            break;
        end
        
        errorVector1=xor(s,decoderDirect);
        errornum1=sum(errorVector1);
        ErrorCount1=ErrorCount1+errornum1;
        
        SourceCount=SourceCount+SourceBlockLength;

        BER_DirectOutput=ErrorCount1/SourceCount;
        
        BER=ErrorCount/SourceCount;
        
        BER_DirectOutput
        ErrorCount1
        
        BER
        ErrorCount
        SourceCount
    end
    
    if(SourceCount>=TestLength)
        break;
    end
    BlockCount=BlockCount+1;
    iterativeTimes(BlockCount)=iterativeCount-1;
   
    
end