MATLAB调用C程序、调试和LDPC译码

摘要:
MATLAB是一个非常有用的工具。起初,我还注意到MATLAB可以调用C程序来加快程序运行,但利用量化代码有时调用更强大的内置函数的想法尚未付诸实践。

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代码如下

MATLAB调用C程序、调试和LDPC译码第1张MATLAB调用C程序、调试和LDPC译码第2张
#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中,可以进行逐语句的调试,如下图所示

image

   

    此时,无法操作MATLAB,可以在Visual Studio中进行操作。如果发现自动窗口中的变量取值不正确,调试无法正常进行,那多半是MATLAB数据转化过程中出现了问题,尤其是指针问题。这不仅可能导致运行结果出错,同时可能会倒是MATLAB崩溃。

4. 其他

    这是一种比较简单的调用C程序的方法,只需要对已有的C函数进行简单的修改即可。还有其他的方法,譬如调用动态链接库,可以自行查看MATLAB的帮助。

免责声明:文章转载自《MATLAB调用C程序、调试和LDPC译码》仅用于学习参考。如对内容有疑问,请及时联系本站处理。

上篇iOS 中的 HotFix 方案总结详解【文档】ORACLE RAC 修改SCAN监听端口、VIP网络配置参考文档下篇

宿迁高防,2C2G15M,22元/月;香港BGP,2C5G5M,25元/月 雨云优惠码:MjYwNzM=

相关文章

matlab练习程序(动感模糊)

  其实在matlab中调用系统函数fspecial和imfilter这两个函数就能很简单的实现动感模糊,不过我可不想就这样简单的实现,所以就自己从头写了。动感模糊最复杂的就是构造卷积矩阵了,卷积矩阵由两个参数决定,模糊半径r和模糊角度theta,通过这两个参数,就能构造不同的模板矩阵,而图像卷积时则和普通的卷积没有什么区别。   构造模板矩阵,首先是通过...

matlab绘图方法汇总

Matlab画图 强大的画图功能是Matlab的特点之中的一个,Matlab提供了一系列的画图函数,用户不须要过多的考虑画图的细节,仅仅须要给出一些基本參数就能得到所需图形,这类函数称为高层画图函数。此外,Matlab还提供了直接对图形句柄进行操作的低层画图操作。这类操作将图形的每个图形元素(如坐标轴、曲线、文字等)看做一个独立的对象,系统给每个对象分配一...

chromium浏览器开发系列第四篇:如何调试最新chromium源码

附上上几篇文章地址,方便大家查看: 下载源码 编译源码 目录结构 接二连三的事情,时间比较紧张,但是还是没有把这个系列的文章丢掉,因为这也是对自己知识的总结吧。提倡大家多写写,以后再看的时候会有种莫名的小激动。 上周写的是chromium的目录结构,好像大家不太感兴趣,在我看来这部分很重要。开头有链接地址,大家想看可以再看看。 从源码下载到编译,到目录结构...

Matlab中编译C++文件

今天在跑《Robust Object Tracking via Sparsity-based Collaborative Model》这篇文章的代码时候,发现出现如下错误: 发现错误时由于vgg_kmiter这个c类型的函数引起的,于是百度查看先关知识,原来是没有编译vgg_kmiter引起的错误,并且在此指出:由于高matlab版本不再支持mxCrea...

4.VREP常见问题1

1.VREP无法导入贴图 在vrep中无法正确识别中文路径。但是,可以正确打开位于中文路径下的场景文件。在这种情况下调试场景文件,不会出现错误。 vrep在导入贴图时,出现贴图无法导入。由于贴图所存在的路径中包含中文。 图21 导入贴图错误 解决方法:项目路径存在英文目录下。 2.Vrep软件运行错误 场景中的对象,经过复制之会取名带后缀#0 或者其他...

C++第三十七篇 -- 调试驱动程序

上一篇写的KMDF程序是通过串口进行配置的,那么我们在VS中Attach to process外,可以直接用Winbdg进行调试,winbdg.exe所在路径为C:Program Files (x86)Windows Kits10Debuggersx64。怎么使用它呢。 如果是串口方式配置的话:windbg.exe -b -k com:pipe,port=...