test_mac.m

速度很快的低密度奇偶校验码的Matlab仿真程序

% test for 802.16e's LDPC codes
clear all;

%%%%%%%%%%%%%%%%%%%%for 802.16e LDPC code construction%%%%%%%%%%%%%%%%%%%%55
% % code rate option( 1:1/3 2:1/2 3:2/3 4:3/4)
% % generate the HBMatrix
% option=2;
% HBMatrix= LDPC_HBMatrix_std(option);
% %informaiton bit length.
% Bit_Length=288;                       
% % Call LDPC maker, the format of HMatrix refer to LDPC_maker_16e.c
% [a_tmp, b_tmp, HMatrix] = LDPC_maker_16e(HBMatrix,Bit_Length);
% HMatrix=double(HMatrix);
% % find H information
% M=HMatrix(1,1);N=HMatrix(2,1);
% % construct Ha
% Ha=zeros(M,N);
% for i=2:(M+1)
%     num=HMatrix(1,i);
%     tt=HMatrix(2:num+1,i);
%     Ha(i-1,tt+1)=1;
% end;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%for H matrix downloaded from Mackay %%%%%%%%%%%%%%%%%
% load PEGirReg504x1008Direct.mat
%load PEGirReg252x504Direct.mat
load test.mat
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[M,N]=size(Ha);
% [K,N]=size(Ga);
%最大度数
maxVarNodeDegree=max(sum(Ha));
maxChkNodeDegree=max(sum(Ha'));
%度分布情况列向量
VarNodeDegree=sum(Ha)';
CheckNodeDegree=sum(Ha')';
%H矩阵向交织器转化
alpha=H2alpha(Ha);
alpha=alpha-1;
%变量结点与校验结点所连结点信息
CheckNodeConnectionVarNode=zeros(N,maxVarNodeDegree);
for i=1:N
    tt=find(Ha(:,i));
    tt=[tt;zeros(maxVarNodeDegree-length(tt),1)];
    CheckNodeConnectionVarNode(i,:)=tt';
end;
VarNodeConnectionCheckNode=zeros(M,maxChkNodeDegree);
for i=1:M
   tt=find(Ha(i,:));
   tt=[tt,zeros(1,maxChkNodeDegree-length(tt))];
   VarNodeConnectionCheckNode(i,:)=tt;
end;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[K,N]=size(Ga);
Rate=K/N;

% crc_length=8;
%仿真参数
Lmax=20.0;
IterMax=100;
% beta=0.5;

% EbN0=[1:0.25:2,2.2:0.1:2.8];
EbN0=0;
frames=1e2;
codetype=6;

fprintf('simulation conditions:\n');
disp('code length N=');disp(N);
disp('information length K=');disp(K);
disp('Rate=');disp(Rate);
disp('IterMax=');disp(IterMax);
disp('frames=');disp(frames);

Pb=zeros(1,length(EbN0));fer=zeros(1,length(EbN0));

fr=zeros(1,length(EbN0));
for i=1:length(EbN0)
    sigma=1.0/sqrt(2*Rate*10.^(EbN0(i)/10));
    err_sum=0;fer_sum=0;err_sum1=0;fer_sum1=0;
    for j=1:frames       
        %generate bit;
        x=floor(2*rand(K,1));
        %encode
        y=ldpc_encode(x,Ga,2);
        %bpsk
        z=2*y'-1;
        %awgn
        z=z+sigma*randn(N,1);
        %demap
        inf_L=2./sigma^2*z;
        
        %for ldpc_decoder_spa
        tic;
        [inf_bit,success]=ldpc_decoder_spa(inf_L,VarNodeDegree,CheckNodeDegree,alpha,IterMax,Lmax);
        toc
        d=floor((sign(inf_bit)+1)/2);
        sum_d=sum(abs(d(N-K+1:end)-x));
        err_sum=err_sum+sum_d;
        if (sum_d>0)
            fer_sum=fer_sum+1;
        end;
          if (mod(j,50)==0)
              [i j fer_sum err_sum]
          end;
        if (fer_sum>100 || (fer_sum>50 && j>1e7))
            break;
        end;
    end;     
    Pb(i)=err_sum/j/K;      fer(i)=fer_sum/j;
    fr(i)=j;
    if (err_sum==0)
        break;
    end;
end;
semilogy(EbN0,Pb,'r-');
hold on;
semilogy(EbN0,fer,'r-');
title('LDPC decoder');
xlabel('EbN0(dB)');
ylabel('BER');

save ldpc_crc_h