CUDA脉冲压缩

1 前言

目前已经在cuda代码中完成RD算法，其中距离向压缩和方位向压缩fft使用了cufft库，距离徙动矫正和矩阵转置采用传统c的做法，还没有针对gpu进行优化。

对于最后的成像结果，还没有进行定量分析，后续还需要添加输出各项参数的函数。

2 CUDA脉冲压缩

2.1 C中生成回波数据及Gnuplot绘图

为了方便调试，采用c语言直接生成回波数据。首先通过以下代码生成回波信号CompData，距离向采样点数LENGTH为1024，方位向采样点数tm为800。

for(j=0;j<tm;j++)
{
	for (i = 0; i < LENGTH; i++)
	{
		Rtm = sqrt(pow(Xc,2) + pow((-v*(rm/v)/2+v*j/fa-Yc),2) + pow(H,2));	
		t = (i-LENGTH/2)/fs-2*Rtm/c;	
		CompData[j*LENGTH+i].x = cos(PI * Kr * pow(t,2))*(t>(-Tp/2-2*Rtm/c))*(t<(Tp/2-2*Rtm/c));	
		CompData[j*LENGTH+i].y = sin(PI * Kr * pow(t,2))*(t>(-Tp/2-2*Rtm/c))*(t<(Tp/2-2*Rtm/c));	
	}
}

通过Gnuplot工具，可以实现matlab中数据绘图的功能，回波数据波形图如下：

同样生成距离向参考函数，参考函数波形图如下：

2.2 脉冲压缩

根据脉冲压缩的原理，需要对回波数据和方位向参考函数分别作傅里叶变换，将结果相乘后做傅里叶逆变换。

这里的fft和ifft均采用cuFFT实现，因为回波数据矩阵形式，所以需要对每行数据分别作fft，这里直接采用for循环进行处理，后期可通过流处理加速。向量乘法这里采用编写核函数进行gpu加速，核函数代码如下：

__global__ void conj(cufftComplex *src, cufftComplex *data, cufftComplex *dest)

{
	size_t i = blockIdx.x * blockDim.x + threadIdx.x;
	if(i<LENGTH)
	{
		dest[i].x = src[i].x * data[i].x + src[i].y * data[i].y;	
		dest[i].y = src[i].y * data[i].x - src[i].x * data[i].y;	
	}
}

最后通过ifft得到脉冲压缩结果：

在nano中测得的运算时间为140ms左右，有很大的有化空间。

3 RDA成像结果

将成像结果和Matlab进行对比，左图为cuda实现的点目标成像，右图为Matlab实现的点目标成像。

4 Gunplot的使用

    在cuda编程中，需要将计算的结果进行可视化，比如回波数据、脉冲压缩的数据、最后成像结果等，因为c不像matlab本身集成了plot工具，所以需要第三方工具进行画图。

4.1 gunplot 命令行用法

-V, —version
-h, —help
-p —persist 保持窗口不关闭
-d —default-settings 不执行初始化
-s —slow 等待加载
-e “command1; command2; …” 执行命令
-c scriptfile ARG1 ARG2 … 调用脚本

为了实验该工具的用法，编写了以下程序：

环境：ubuntu20.04, gcc

#include <stdio.h>
#include <math.h>
#define pi 3.1415926
#define wLow 2 * pi *fLow / fs
#define wHight 2 * pi *fHight / fs

int fs = 4096;
int fLow = 50;
int fHight = 1000;

short buff[1024], x[16], y[1024] = {0};

const short h[16] = {
	-2095, -730, -130, 1036, 2611, 4295, 5714, 6524, 6524,
	5714, 4295, 2611, 1036, -130, -730, -2095};

short addNewData(short newData)
{
	int y0;
	short i;
	static short index = 0;
	index = index % 16;
	x[index] = newData;
	for (i = 0; i < 16; i++)
	{
		y0 += x[(i + 1 + index) % 16] * h[i];
	}
	index++;
	return (int)(y0 >> 15);
}
int main(void)
{
	short i = 0;
	//	模拟生成输入信号
	for (i = 0; i < 1024; i++)
	{
		buff[i] = 1024 * sin(wLow * i) + 1024 * sin(wHight * i);
	}
	// 每次输入一个数
	for (i = 0; i < 1024; i++)
	{
		y[i] = addNewData(buff[i]);
		printf("%d %d\n", i, y[i]);
	}

	return 0;
}

该段程序为一个循环滤波器，将程序编译后运行下面命令得到滤波后的波形图。其中u 表示使用第几列数据，l 表示用线表示。

gnuplot -p -e "plot [0:1023] '<fir2.exe' u 1:2 w l"