matlab和C混合编程实现脉冲压缩

环境：matlab2015b  vs2010  cuda7.5
 
步骤一：创建ComplexMul.h文件
 
#include"cuComplex.h" //因为处理的数据是复数，需要添加这个头文件
#ifndef__COMPLEXMUL_H__
#define__COMPLEXMUL_H__
 
//声明所需要的函数。extern表示该函数在其他文件中执行
externvoid ComplexMul(cuComplex* A,  cuComplex*B, cuComplex* C, int size);//实现矩阵点乘
externvoid Scale(cuComplex* Input,  cuComplex*Output, int size);//实现归一化
 
#endif//__COMPLEXMUL_H__
 
步骤二：在ComplexMul.cu中实现声明的函数
 
#include "ComplexMul.h"
#include "mex.h"
#include "cuComplex.h"
 
__global__ void ComplexMulMask(cuComplex* A,  cuComplex* B, cuComplex* C, int size)
{
    const intnumThreads = blockDim.x * gridDim.x;
    const intthreadID = blockIdx.x * blockDim.x + threadIdx.x;
    for (int i= threadID; i < size; i += numThreads)
    {
        C[i].x= A[i].x * B[i].x - A[i].y * B[i].y;
        C[i].y= A[i].x * B[i].y + A[i].y * B[i].x;   
    }
}
__global__ void ScaleMask(cuComplex* input, floatscale, int size)
{
    const intnumThreads = blockDim.x * gridDim.x;
    const intthreadID = blockIdx.x * blockDim.x + threadIdx.x;
    for (int i= threadID; i < size; i += numThreads)
    {
       input[i].x *=scale;
        input[i].y*=scale;
    }
 
}
void ComplexMul(cuComplex* A,  cuComplex* B, cuComplex* C, int size)
{
    cuComplex*devPtrA=0, *devPtrB=0, *devPtrC=0;
    cudaMalloc(&devPtrA,sizeof(cuComplex)*size);
    cudaMalloc(&devPtrB,sizeof(cuComplex)*size);
    cudaMalloc(&devPtrC,sizeof(cuComplex)*size);
 
    cudaMemcpy(devPtrA,A,sizeof(cuComplex)*size, cudaMemcpyHostToDevice );
    cudaMemcpy(devPtrB,B,sizeof(cuComplex)*size, cudaMemcpyHostToDevice );
 
    ComplexMulMask<<<128,128>>>(devPtrA,devPtrB, devPtrC, size);
 
    cudaMemcpy(C,devPtrC,sizeof(cuComplex)*size,cudaMemcpyDeviceToHost);
 
    cudaFree(devPtrA);
    cudaFree(devPtrB);
    cudaFree(devPtrC);
 
}
 
 
 void Scale(cuComplex*Input,  cuComplex* Output, int size)
{
    cuComplex*devPtrIn=0, *devPtrOut=0;
    cudaMalloc(&devPtrIn,sizeof(cuComplex)*size);
    cudaMalloc(&devPtrOut,sizeof(cuComplex)*size);
 
    cudaMemcpy(devPtrIn,Input,sizeof(cuComplex)*size, cudaMemcpyHostToDevice );
 
    ScaleMask<<<128,128>>>(devPtrIn,1.0f/size, size);
 
    cudaMemcpy(Output,devPtrIn,sizeof(cuComplex)*size,cudaMemcpyDeviceToHost);
 
    cudaFree(devPtrIn);
    cudaFree(devPtrOut);
   
 
}
 
步骤三：使用-c选项编译CUDA代码，生成目标文件。
 
在matlab命令窗口输入：
>>system(‘nvcc–c ComplexMul.cu’)
编译成功后生成ComplexMul.obj文件
 
步骤四：创建mex函数，在matlab里面创建PulseCompression.cpp文件
 
//脉冲压缩
 
#include "mex.h"
#include <cuda_runtime.h>
#include <cufft.h>
#include "ComplexMul.h"
 
#defineNX 16384 //脉压点数，必须是2的整数次幂
#defineBATCH 1
 
 
voidmexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[] )
{
  
    if( nrhs != 4 )
        mexErrMsgTxt("Invaid number ofinput arguments");//检查输入参数个数是否为2
   
    double* signal_i =(double*)mxGetData(prhs[0]); //读取数据
    double* signal_q = (double*)mxGetData(prhs[1]);
    double* filter_i =(double*)mxGetData(prhs[2]);
    double* filter_q =(double*)mxGetData(prhs[3]);
   
    cufftComplex*signal_H=(cufftComplex*)malloc(NX*BATCH*sizeof(cufftComplex));
    cufftComplex* filter_H=(cufftComplex*)malloc(NX*BATCH*sizeof(cufftComplex));
   // printf("输入信号为：\n");
    for(int i=0;i<NX;i++)
    {
        signal_H[i].x=signal_i[i];
        signal_H[i].y=signal_q[i];
      // if(i<10) printf("%.4f+%.4fi\n",signal_i[i],signal_q[i] );
    }
     //printf("输入滤波函数为：\n");
    for(int i=0;i<NX;i++)
    {
        filter_H[i].x=filter_i[i];
        filter_H[i].y=filter_q[i];
      // if(i<10) printf("%.4f+%.4fi\n",filter_i[i],filter_q[i] );
    }
    
    cufftComplex *signal_D; //存储信号
    cudaMalloc( &signal_D,  sizeof(cufftComplex)*NX*BATCH);  // malloc memory
    cudaMemcpy( signal_D, signal_H,sizeof(cufftComplex)*NX*BATCH, cudaMemcpyHostToDevice);   // copy memory from host to device
    cufftComplex *filter_D; //存储滤波器系数
    cudaMalloc( &filter_D,  sizeof(cufftComplex)*NX*BATCH);  // malloc memory
    cudaMemcpy( filter_D, filter_H,sizeof(cufftComplex)*NX*BATCH, cudaMemcpyHostToDevice);   // copy memory from host to device
   
    cufftComplex* signal_deviceOut; //存储信号fft的结果
    cudaMalloc( &signal_deviceOut,sizeof(cufftComplex)*NX*BATCH );
    cufftComplex* filter_deviceOut; //存储滤波器系数  fft的结果
    cudaMalloc( &filter_deviceOut,sizeof(cufftComplex)*NX*BATCH );
   
    //记录时间
    cudaEvent_t start, stop;
         cudaEventCreate(&start);
         cudaEventCreate(&stop);
         cudaEventRecord(start,0);//计时开始
   
    cufftHandle plan_S; //FFT的参数都相同，只需要创建一个plan
    cufftPlan1d(&plan_S, NX*BATCH,CUFFT_C2C, BATCH);
    cufftExecC2C(plan_S, signal_D,signal_deviceOut, CUFFT_FORWARD );
    cufftExecC2C(plan_S, filter_D,filter_deviceOut, CUFFT_FORWARD );
   
    cufftComplex* MatchS_device; //存储点乘的结果
    cudaMalloc( &MatchS_device,sizeof(cufftComplex)*NX*BATCH );
   ComplexMul(signal_deviceOut,filter_deviceOut, MatchS_device,NX*BATCH);//调用函数计算复数点乘
   
    cufftComplex* MatchSifft_device; //存储ifft之后的结果
    cudaMalloc( &MatchSifft_device,sizeof(cufftComplex)*NX*BATCH );
   
    cufftExecC2C(plan_S, MatchS_device,MatchSifft_device, CUFFT_INVERSE );
   
    cufftComplex* MatchS_scale; //存储归一化之后的结果
    cudaMalloc( &MatchS_scale,sizeof(cufftComplex)*NX*BATCH );
    Scale(MatchSifft_device, MatchS_scale,NX);//调用函数对ifft结果进行归一
   
    //记录结束时间
         cudaEventRecord(stop,0);
         cudaEventSynchronize(stop);//计时结束
         float elapsedTime;
         cudaEventElapsedTime(&elapsedTime,start,stop);
   //printf("%3.3fms\t",elapsedTime);//输出计时结果
         cudaEventDestroy(start);
         cudaEventDestroy(stop);
    
    float* signal_out = (float*)mxMalloc(sizeof(cufftComplex)*NX*BATCH );
    cudaMemcpy(signal_out,MatchS_scale,sizeof(cufftComplex)*NX*BATCH, cudaMemcpyDeviceToHost );
 
     plhs[0] = mxCreateNumericMatrix( BATCH,NX, mxSINGLE_CLASS, mxCOMPLEX);
     plhs[1] = mxCreateDoubleMatrix(1,1,mxREAL);   
     double *T;
     T=mxGetPr(plhs[1]);
     *T=elapsedTime;
    
    float* S_real = (float*)mxGetPr(plhs[0]);
    float* S_imag = (float*)mxGetPi(plhs[0]);
    float* S_complex = signal_out;
        for ( int c = 0; c < BATCH; ++c )
        {
            for (int r = 0; r < NX; ++r )
            {
                *S_real++ = *S_complex++;
                *S_imag++ = *S_complex++;
            }
        }
    mxFree(signal_out);
    cufftDestroy(plan_S);
    cudaFree(signal_D);
   
}
 
步骤五：编译mex，在matlab里输入以下命令：
 
>>mex PulseCompression.cpp complexMul.obj -lcudart -lcufft -L"C:\ProgramFiles\NVIDIA GPU Computing Toolkit\CUDA\v7.5\lib\x64" -v-I"C:\Program Files\NVIDIA GPU Computing Toolkit\CUDA\v7.5\include"
l  -lcudart：告知mex正在使用cuda运行时库。（l是小写的L）
l  -lcufft：表明子啊调用CUFFT库。
l  -Ldir：dir是CUDA和CUFFT库所在的目录。
l  -Idir：dir是CUDA和CUFFT头文件所在的目录。（I是大写的i）
编译成功会生成PulseCompression.mexw64文件。
 
步骤六：在matlab里面调用PulseCompression函数。编写PulseCom.m文件计算设备端脉压并与matlab脉压结果做对比。
 
clc;
clearall;
MON=1000;
fori=1:MON
%随机生成的数据
data_AI=rand(1,FFT_N)*100;
data_AQ=rand(1,FFT_N)*100;
PC_filter_I=rand(1,FFT_N)*100;
PC_filter_Q=rand(1,FFT_N)*100;
 
%设备端脉压
tic
[MatchSo_GPU,Tfft(i)]=PulseCompression(data_AI,data_AQ,PC_filter_I,PC_filter_Q);
TD(i)=toc;
%MATALB做脉压
tic
dataA=data_AI+j*data_AQ;
dataA_fft=fft(dataA,FFT_N);
PC_filter=PC_filter_I+j*PC_filter_Q;
Match_fft=fft(PC_filter,FFT_N);
MatchS=dataA_fft.*Match_fft;
MatchSo=ifft(MatchS);%脉压输出
TH(i)=toc;
end
T_fft=sum(Tfft)/MON  
Tdev=sum(TD)/MON
THos=sum(TH)/MON
 
figure();
plot(abs(MatchSo));
xlabel('距离单元')
ylabel('脉压输出')
title('Matlab端脉冲压缩');grid;
figure();
plot(abs(MatchSo_GPU));
xlabel('距离单元')
ylabel('脉压输出')
title('GPU端脉冲压缩');grid;
 
结果：
 
 
 
 
测试数据结果：

随机数据结果：