ldpcdec.c

此代码是ＬＤＰＣ码进行BP算法的重要参考代码

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

//采用基于概率似然比的迭代译码算法，
//详见《电路与系统学报》2004年2月第9卷第1期
//《基于低密度校验码的OFDM编码调制译码算法》一文
//该文作者为徐志江、李式巨、官军
void decode( double *vhat, double *iter,
             double *pr, double *h1i, double *h1j,
             double rows, double cols, double max_iter, double h1num )
{
    double *prmn;
    double *dqmn;
    double *lrmn;
    double *prn;

    int i, j, k;
    int iteration;

    double prod_dqmn;
    double prod_lrmn;

    double dmn;

    int done;
    int temp;

//    mexPrintf("\npr:\n");
//    for (j=0; j<cols; j++)
//    {
//        mexPrintf("%f\t", pr[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);

    prmn    = (double *)mxMalloc((size_t) (h1num*sizeof(double)));
    dqmn    = (double *)mxMalloc((size_t) (h1num*sizeof(double)));
    lrmn    = (double *)mxMalloc((size_t) (h1num*sizeof(double)));
    prn     = (double *)mxMalloc((size_t) (h1num*sizeof(double)));

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

    for (iteration=1; iteration<=(int)max_iter; iteration++)
    {   //main iteration loop
        //横向步骤：由信息节点的先验概率似然比按置信传播算法得出各校验节点的后验概率似然比。
        for (k=0; k<h1num; k++)
        {
            dqmn[k]=(1-prmn[k])/(1+prmn[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)
                {
                    dmn = prod_dqmn/dqmn[k];
                    lrmn[k] = (1 - dmn)/(1 + dmn);
                }
            }
        }

        //纵向步骤：由校验节点的后验概率似然比推算出信息节点的后验概率似然比。
        for (j=0; j<cols; j++)
        {
            prod_lrmn = 1;
            for (k=0; k<h1num; k++)
            {
                if (h1j[k]==j)
                {
                    prod_lrmn *= lrmn[k];
                }
            }
            for (k=0; k<h1num; k++)
            {
                if (h1j[k]==j)
                {
                    prn[k] = pr[j]*prod_lrmn;
                    prmn[k] = prn[k]/lrmn[k];
                    if (prn[k]>1)
                    {
                        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(pr,h1i,h1j,rows,cols,max_iter)
void mexFunction( int nlhs, mxArray *plhs[],
                  int nrhs, const mxArray*prhs[] )
{
    double *pr;     //pointer variable for input matrix pr
    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

    pr = mxGetPr(prhs[0]);      //pointer to pr
    h1i = mxGetPr(prhs[1]);     //pointer to h1i
    h1j = mxGetPr(prhs[2]);     //pointer to h1j
    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, pr, h1i, h1j, rows, cols, max_iter, h1num);

}