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;
}

View Code

    上述代碼就是就是一個標準的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的幫助。