ldpcpl1dec.c

LDPC码的译码

#include <mex.h>
#include <matrix.h>         //for Matlab mx and mex fuctions
#include <math.h>

void decode( double *vhat, double *iter,
             double *pl1, double *h1i, double *h1j,
             double rows, double cols, double max_iter, double h1num )
{
    double *qmn1;
    double *dqmn;
    double *rmn0;
    double *rmn1;
    double *qn1;

    int i, j, k;
    int iteration;

    double prod_dqmn;
    double prod_rmn0;
    double prod_rmn1;
    double const1;
    double const2;
    double const3;
    double const4;

    int done;
    int temp;

//    mexPrintf("\npl1:\n");
//    for (j=0; j<cols; j++)
//    {
//        mexPrintf("%f\t", pl1[j]);
//    }
//
//    mexPrintf("\nh1i:\n");
//    for (k=0; k<h1num; k++)
//    {
//        mexPrintf("%d\t", (int)h1i[k]);
//    }
//
//    mexPrintf("\nh1j:\n");
//    for (k=0; k<h1num; k++)
//    {
//        mexPrintf("%d\t", (int)h1j[k]);
//    }
//
//    mexPrintf("\nrows=%d\tcols=%d\tmax_iter=%d\th1num=%d\n",(int)rows,(int)cols,(int)max_iter,(int)h1num);

    qmn1    = (double *)mxMalloc((size_t) (h1num*sizeof(double)));
    dqmn    = (double *)mxMalloc((size_t) (h1num*sizeof(double)));
    rmn0    = (double *)mxMalloc((size_t) (h1num*sizeof(double)));
    rmn1    = (double *)mxMalloc((size_t) (h1num*sizeof(double)));
    qn1     = (double *)mxMalloc((size_t) (h1num*sizeof(double)));

    //初始化：对特定的信道预设信息比特的先验概率。
    for (k=0; k<h1num; k++)
    {
        qmn1[k] = pl1[(int)h1j[k]];
    }

    for (iteration=1; iteration<=(int)max_iter; iteration++)
    {   //main iteration loop
        //横向步骤：由信息节点的先验概率按置信传播算法得出各校验节点的后验概率。
        for (k=0; k<h1num; k++)
        {
            //计算概率差
            dqmn[k]=1-qmn1[k]-qmn1[k];
        }
        for (i=0; i<rows; i++)
        {
            prod_dqmn = 1;
            for (k=0; k<h1num; k++)
            {
                if (h1i[k]==i)
                {
                    prod_dqmn *= dqmn[k];
                }
            }
            for (k=0; k<h1num; k++)
            {
                if (h1i[k]==i)
                {
                    rmn0[k] = (1 + prod_dqmn/dqmn[k])/2;
                    rmn1[k] = 1 - rmn0[k];
                }
            }
        }

        //纵向步骤：由校验节点的后验概率推算出信息节点的后验概率。
        for (j=0; j<cols; j++)
        {
            prod_rmn0 = 1;
            prod_rmn1 = 1;
            for (k=0; k<h1num; k++)
            {
                if (h1j[k]==j)
                {
                    prod_rmn0 *= rmn0[k];
                    prod_rmn1 *= rmn1[k];
                }
            }
            for (k=0; k<h1num; k++)
            {
                if (h1j[k]==j)
                {
                    const3 = (1-pl1[j])*prod_rmn0;
                    const4 = pl1[j]*prod_rmn1;
                    qn1[k] = const4/(const3+const4);
                    const1 = const3/rmn0[k];
                    const2 = const4/rmn1[k];
                    qmn1[k] = const2/(const1+const2);
                    if (qn1[k]>0.5)
                    {
                        vhat[j]=1;
                    }
                    else
                    {
                        vhat[j]=0;
                    }
                }//if (h1j[k]==j)
            }//for (k=0; k<h1num; k++)
        }//for (j=0; j<cols; j++)

        //如果判决条件满足，译码结束，否则继续迭代
        done = 1;
        for (i=0; i<rows; i++)
        {
            temp = 0;
            for (k=0; k<h1num; k++)
            {
                if (h1i[k]==i)
                {
                    temp += (int)vhat[(int)h1j[k]];
                }
            }
            temp = temp%2;
            if (temp!=0)
            {
                done = 0;
                break;
            }
        }

        *iter = iteration;
        if (done)
        {
//            mexPrintf("\niter=%d\n",(int)*iter);
//            mexPrintf("vhat:\n");
//            for (j=0; j<cols; j++)
//            {
//                mexPrintf("%d\t", (int)vhat[j]);
//            }

            return;
        }

//        mexPrintf("\niter=%d\n",(int)*iter);
//        mexPrintf("vhat:\n");
//        for (j=0; j<cols; j++)
//        {
//            mexPrintf("%d\t", (int)vhat[j]);
//        }

    }//main iteration loop end
}


//  0     1            0   1   2   3    4      5
//[vhat,iter]=ldpcdec(pl1,h1i,h1j,rows,cols,max_iter)
void mexFunction( int nlhs, mxArray *plhs[],
                  int nrhs, const mxArray*prhs[] )
{
    double *pl1;    //pointer variable for input matrix pl1
    double *h1i;    //pointer variable for input matrix h1i
    double *h1j;    //pointer variable for input matrix h1j
    double rows;    //H matrix rows
    double cols;    //H matrix cols
    double max_iter;//maxim iteration

    double *vhat;   //pointer variable for output matrix vhat
    double *iter;   //iteration output

    double h1num;   //H matrix 1's number

    pl1 = mxGetPr(prhs[0]);     //pointer to pl1
    h1i = mxGetPr(prhs[1]);     //pointer to pl1
    h1j = mxGetPr(prhs[2]);     //pointer to pl1
    rows = mxGetScalar(prhs[3]);//value of rows
    cols = mxGetScalar(prhs[4]);//value of cols
    max_iter = mxGetScalar(prhs[5]);//value of max_iter

    plhs[0] = mxCreateDoubleMatrix(1, cols, mxREAL); //matrix for output
    vhat = mxGetPr(plhs[0]);    //pointer to vhat
    plhs[1] = mxCreateDoubleScalar(0);
    iter = mxGetPr(plhs[1]);

    h1num = mxGetN(prhs[1]);    //number of cols of h1i

    decode(vhat, iter, pl1, h1i, h1j, rows, cols, max_iter, h1num);

}