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的帮助。
