fastldpcdec.c

LDPC码的译码

//fastldpcdec.c

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

//采用基于概率似然比的迭代译码算法，
//详见《电路与系统学报》2004年2月第9卷第1期
//《基于低密度校验码的OFDM编码调制译码算法》一文
//该文作者为徐志江、李式巨、官军


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

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

    double h1num;   //H matrix 1's number

///////////////////////////////////////////////////////////////
    double *prmn;
    double *dqmn;
    double *lrmn;
    double prn;

    int i, j, k;
    int iteration;
    int rowbase;
    int colbase;

    double prod_dqmn;
    double prod_lrmn;

    double dmn;

    int done;
    int temp;
///////////////////////////////////////////////////////////////

    pr = mxGetPr(prhs[0]);          //pointer to pr
    h1i2j = mxGetPr(prhs[1]);       //pointer to h1i2j
    h1j = mxGetPr(prhs[2]);         //pointer to h1j
    rows = mxGetScalar(prhs[3]);    //value of rows
    cols = mxGetScalar(prhs[4]);    //value of cols
    row1num = mxGetPr(prhs[5]);     //pointer to row1num
    col1num = mxGetPr(prhs[6]);     //pointer to col1num
    max_iter = mxGetScalar(prhs[7]);//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[2]);    //number of cols of h1j

///////////////////////////////////////////////////////////////
    prmn    = (double *)mxMalloc((size_t) (h1num*sizeof(double)));
    dqmn    = (double *)mxMalloc((size_t) (h1num*sizeof(double)));
    lrmn    = (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]);
        }

        rowbase = 0;
        for (i=0; i<rows; i++)
        {
            prod_dqmn = 1;
            for (k=0; k<row1num[i]; k++)
            {
                prod_dqmn *= dqmn[(int)h1i2j[k+rowbase]];
            }
            for (k=0; k<row1num[i]; k++)
            {
                dmn = prod_dqmn/dqmn[(int)h1i2j[k+rowbase]];
                lrmn[(int)h1i2j[k+rowbase]] = (1 - dmn)/(1 + dmn);
            }
            rowbase += row1num[i];
        }

        //纵向步骤：由校验节点的后验概率似然比推算出信息节点的后验概率似然比。
        colbase = 0;
        for (j=0; j<cols; j++)
        {
            prod_lrmn = 1;
            for (k=0; k<col1num[j]; k++)
            {
                prod_lrmn *= lrmn[k+colbase];
            }
            prn = pr[j]*prod_lrmn;
            if (prn>1)
            {
                vhat[j]=1;
            }
            else
            {
                vhat[j]=0;
            }

            for (k=0; k<col1num[j]; k++)
            {
                prmn[k+colbase] = prn/lrmn[k+colbase];
            }
            colbase += col1num[j];
        }

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

        *iter = iteration;
        if (done)
        {
            return;
        }

    }//main iteration loop end

}