MATLAB調用C程序、調試和LDPC譯碼

MATLAB是一個很好用的工具。利用MATLAB腳本進行科學計算也特別方便快捷。但是代碼存在較多循環時，MATLAB運行速度極慢。如果不想放棄MATLAB中大量方便使用的庫，又希望代碼能迅速快捷的運行，可以考慮將循環較多的功能采用C編寫，MATLAB調用。本文將概述這一過程。雖然本文以LDPC譯碼算法為例，但不懂該算法不影響本文閱讀。

1. 起因

    最開始用MATLAB寫的LDPC譯碼算法中，其中一個版本是這里，里面有三重循環，運行速度極慢。后來考慮了MATLAB的向量化操作，通過算法的合理划分以及內置函數調用，成功將三重循環修改為1層，具體這一版本的代碼可見這里。通過這一手段，函數的運行速度提高了幾倍乃至幾十倍。雖然這一方法下運行速度依舊比不過MATLAB工具箱中的comm.LDPCDecoder，遠比不上利用GPU的comm.gpu.LDPCDecoder，但勝在可明確算法並具有一定擴展性。

    起初也注意到可以通過MATLAB調用C程序來加速程序運行，但向量化后的代碼湊活能用，加上有時也可調用更為強大的內置函數，這一想法一直沒有付諸實踐。這幾天想好好整理一下代碼，遂萌發了寫一個C版本譯碼算法的想法。代碼現在的狀態是“能用”，這里把相關經驗總結分析在此。

2. MATLAB調用C程序

    這一部分的內容在劉曉輝的matlab調用C程序中已經有較為詳細的介紹了，想要正確調用C程序，關鍵概括為2點。

機器上裝有MATLAB編譯器，可通過在MATLAB命令行窗口輸入mex -setup進行具體設置。

有一個正確的接口子程序mexFunction完成MATLAB和C程序之間的數據轉換和程序調用

    這里給出我寫得mexFunction（注意這個代碼寫得不好，沒有任何判斷，沒有健壯性……）

void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {
    double* llr = (double*)mxGetPr(prhs[0]);
    int* rownum = mxGetPr(prhs[1]);
    int* colnum = mxGetPr(prhs[2]);
    int* trans = mxGetPr(prhs[3]);
    double* state = mxGetPr(prhs[4]);
    plhs[0] = mxCreateDoubleMatrix(1, state[1], mxREAL);
    double* r =mxGetPr(plhs[0]);
    ldpcDec( r ,llr, rownum,colnum, trans,state);
}

    mexFunction的規范在劉曉輝的matlab調用C程序一文中已有提及，即

nlhs：輸出參數數目
plhs：指向輸出參數的指針
nrhs：輸入參數數目
prhs：指向輸入參數的指針
例如，在matlab命令行中使用
[a,b]=test(c,d,e)
調用mex函數test時，傳給test的這四個參數分別是
      2，plhs，3，prhs
其中：
prhs[0]=c
prhs[1]=d
prhs[2]=e

    由此可以解釋上述mexFunction，而命令plhs[0] = mxCreateDoubleMatrix(1, state[1], mxREAL) 則定義了一大小為1 × state[1]的矩陣，做為函數的返回值。最后調用的ldpcDec是一個C程序，運行C程序后plhs[0]指向的內存空間存儲的就是滿足要求的計算結果。ldpcDec代碼如下

#include<stdio.h>
#include<math.h>
void ldpcDec(double*r,double* llr, int* rownum, int* colnum, int* trans, double* state){
    //列有序，trans為映射關系
    //rownum[i]-rownum[i-1],第i+1行的行重
    //colnum[i]-colnum[i-1],第i+1列的列重
    //state[0]:maxiter state[1]:llr & colnum 長度 state[2] rownum 長度,
    //state[3]:H中非零元素個數 state[4]: alpha
    double* temp;
    double* decodedtemp;
    temp = (double*)malloc(sizeof(double)*state[3]);
    decodedtemp = (double*)malloc(sizeof(double)*state[3]);
    //init
    int ii = 0;
    for (int i = 0; i<state[1]; i++){
        while (ii<colnum[i]){
            temp[ii] = llr[i];
            ii++;
        }
    }
    //iter decode
    int iter;
    for (iter = 0; iter<state[0]; iter++){
        // rowupdate;
        for (int i = 0; i<state[2]; i++){
            // temp[] trans[rownum[i-1]]~trans[rownum[i]]
            int high = rownum[i];
            int low = i>0 ? rownum[i - 1] : 0;

            double minval = fabs(temp[trans[low]]);
            double subminval = fabs(temp[trans[low + 1]]);
            for (int j = low + 1; j<high; j++){
                if (fabs(temp[trans[j]])<minval){
                    subminval = minval;
                    minval = fabs(temp[trans[j]]);
                }
                else if (fabs(temp[trans[j]])<subminval){
                    subminval = fabs(temp[trans[j]]);
                }
            }
            int mark = 1;
            for (int j = low; j < high; j++){
                if (temp[trans[j]] < 0)
                    mark = -mark;
            }
            for (int j = low; j<high; j++){
                if (fabs(temp[trans[j]]) == minval)
                    if (temp[trans[j]]<0)
                        temp[trans[j]] = -mark * state[4] * subminval;
                    else
                        temp[trans[j]] = mark*state[4] * subminval;
                else
                    if (temp[trans[j]]<0)
                        temp[trans[j]] = -mark * state[4] * minval;
                    else
                        temp[trans[j]] = mark*state[4] * minval;
            }
        }
        // colupdate;
        for (int i = 0; i<state[1]; i++){
            int high = colnum[i];
            int low = i>0 ? colnum[i - 1] : 0;
            double colsum = llr[i];
            for (int j = low; j<high; j++){
                colsum = colsum + temp[j];
            }
            if (colsum>0)
                r[i] = 0;
            else
                r[i] = 1;
            for (int j = low; j<high; j++){
                temp[j] = colsum - temp[j];
                decodedtemp[j] = r[i];
            }
        }

        // check equation
        int errflag = 0;
        for (int i = 0; i<state[2]; i++){
            int high = rownum[i];
            int low = i>0 ? rownum[i - 1] : 0;
            int sumval = 0;
            for (int j = low; j<high; j++){
                sumval = sumval + decodedtemp[trans[j]];
            }
            if (sumval % 2 != 0){
                errflag = 1;
                break;
            }
        }
        //
        if (errflag == 0)
            break;
    }
    free(temp);/*釋放指針pointer*/
    temp = NULL;
    free(decodedtemp);/*釋放指針pointer*/
    decodedtemp = NULL;
    return;
}

    上述代碼就是就是一個標准的C函數。

    如果程序無誤，使用起來是極其方便的。完整的代碼如下所示，存儲為ldpc_dec.c文件。

    在MATLAB命令行窗口輸入mex ldpc_dec.c，運行可得到文件ldpc_dec.mexw64（依平台不同可能不同）。需要使用時輸入

r = ldpc_dec(receiveSignal,rowNum,colNum,HRowNum,state);

即可。

3. MATLAB調試C程序

    一般而言，c程序可以事先調試正確，而mexFunction接口函數較為簡單，不容易出錯。然而，有時還是出現一些錯誤，此時可以通過MATLAB對C程序進行調試。以已安裝Visual Studio 和 MATLAB的電腦為例，打開MATLAB和Visual Studio。首先准備好需要調試的c代碼“ldpc_dec.c”，運行命令“mex ldpc_dec.c -g”表示后續需要對C程序進行調試（參考http://blog.csdn.net/ayw_hehe/article/details/6790147）。

    在Visual Studio中點擊“調試”-“附加到進程”，選擇MATALB，在Visual Studio中打開需要調試的C文件並設置斷點，在MATLAB中運行該程序，即輸入“ldpc_dec(receiveSignal,rowNum,colNum,HRowNum,state)”后，在設置斷點處即會中斷。此時進入Visual Studio中，可以進行逐語句的調試，如下圖所示

   

    此時，無法操作MATLAB，可以在Visual Studio中進行操作。如果發現自動窗口中的變量取值不正確，調試無法正常進行，那多半是MATLAB數據轉化過程中出現了問題，尤其是指針問題。這不僅可能導致運行結果出錯，同時可能會倒是MATLAB崩潰。

4. 其他

    這是一種比較簡單的調用C程序的方法，只需要對已有的C函數進行簡單的修改即可。還有其他的方法，譬如調用動態鏈接庫，可以自行查看MATLAB的幫助。