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

//采用基于概率似然比的迭代譯碼算法，
//詳見《電路與系統(tǒng)學(xué)報》2004年2月第9卷第1期
//《基于低密度校驗碼的OFDM編碼調(diào)制譯碼算法》一文
//該文作者為徐志江、李式巨、官軍
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)));

    //初始化：對特定的信道預(yù)設(shè)信息比特的概率似然比。
    for (k=0; k<h1num; k++)
    {
        prmn[k] = pr[(int)h1j[k]];
    }

    for (iteration=1; iteration<=(int)max_iter; iteration++)
    {   //main iteration loop
        //橫向步驟：由信息節(jié)點的先驗概率似然比按置信傳播算法得出各校驗節(jié)點的后驗概率似然比。
        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);
                }
            }
        }

        //縱向步驟：由校驗節(jié)點的后驗概率似然比推算出信息節(jié)點的后驗概率似然比。
        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++)

        //如果判決條件滿足，譯碼結(jié)束，否則繼續(xù)迭代
        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);

}