ldpcllrdec.c

LDPC码的译码

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

void decode( double *vhat, double *iter,
             double *llr, double *h1i, double *h1j,
             double rows, double cols, double max_iter, double h1num )
{
    double *Qij;
    double *Rij;
    double *Qj;
    double *tanhQij;

    int i, j, k;
    int iteration;

    double prod_tanh;
    double sum_rmn;

    double tanh_temp;

    int done;
    int temp;

    Qij     = (double *)mxMalloc((size_t) (h1num*sizeof(double)));
    Rij     = (double *)mxMalloc((size_t) (h1num*sizeof(double)));
    Qj      = (double *)mxMalloc((size_t) (h1num*sizeof(double)));
    tanhQij = (double *)mxMalloc((size_t) (h1num*sizeof(double)));

    for (k=0; k<h1num; k++)
    {
        Qij[k] = llr[(int)h1j[k]];
    }

    for (iteration=1; iteration<=(int)max_iter; iteration++)
    {   //main iteration loop
        //横向步骤：由信息节点的先验概率对数似然比按置信传播算法得出各校验节点的后验概率对数似然比。
        for (i=0; i<rows; i++)
        {
            for (k=0; k<h1num; k++)
            {
                tanhQij[k] = tanh(Qij[k]*0.5);
            }
            prod_tanh = 1;
            for (k=0; k<h1num; k++)
            {
                if (h1i[k]==i)
                {
                    prod_tanh *= tanhQij[k];
                }
            }
            for (k=0; k<h1num; k++)
            {
                if (h1i[k]==i)
                {
                    tanh_temp = prod_tanh/tanhQij[k];
                    Rij[k] = log((1+tanh_temp)/(1-tanh_temp));
                }
            }
        }

        //纵向步骤：由校验节点的后验概率对数似然比推算出信息节点的后验概率对数似然比。
        for (j=0; j<cols; j++)
        {
            sum_rmn = 0;
            for (k=0; k<h1num; k++)
            {
                if (h1j[k]==j)
                {
                    sum_rmn += Rij[k];
                }
            }
            for (k=0; k<h1num; k++)
            {
                if (h1j[k]==j)
                {
                    Qj[k] = llr[j]+sum_rmn;
                    Qij[k] = Qj[k]-Rij[k];
                    if (Qj[k]>0)
                    {
                        vhat[j]=0;
                    }
                    else
                    {
                        vhat[j]=1;
                    }
                }//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)
        {
            return;
        }

    }//main iteration loop end
}


//  0     1               0   1   2   3    4      5
//[vhat,iter]=ldpcllrdec(llr,h1i,h1j,rows,cols,max_iter)
void mexFunction( int nlhs, mxArray *plhs[],
                  int nrhs, const mxArray*prhs[] )
{
    double *llr;     //pointer variable for input matrix llr
    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

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

}