【车载雷达信号处理】利用sinc函数实现扣点

针对信号处理流程中多次FFT输出的频谱结果，在特殊的场景下，可能存在针对某一特定频点的固定"虚警"，所以针对某一个特定频点进行“扣点”的操作是常有的信号处理流程需求。不仅如此，针对最大能量值的扣点也能在不适合使用cfar进行目标检测的情况下（诸如插入了大量的0从而进行的3DFFT），扣取当前FFT结果能量最大点，然后将它的能量分别与FFT扣点剩余的所有频点能量的均值乘以一个系数后进行比较，不断循环从而筛选出目标的角度信息。本文介绍了利用sinc函数对FFT结果进行扣点的操作。

卷积

信号重建

利用sinc函数实现扣点

卷积

所谓两个函数的卷积，本质上就是先将一个函数翻转，然后进行滑动叠加。在连续情况下，叠加指的是对两个函数的乘积求积分，在离散情况下就是加权求和。多次滑动得到的一系列叠加值，最终构成卷积函数。下面我在网上看到的对卷积的一种通俗理解以及相关的例子：

从“积”的过程可以看到，我们得到的叠加值，是个全局的概念。以信号分析为例，卷积的结果是不仅跟当前时刻输入信号的响应值有关，也跟过去所有时刻输入信号的响应都有关系，考虑了对过去的所有输入的效果的累积。在图像处理的中，卷积处理的结果，其实就是把每个像素周边的，甚至是整个图像的像素都考虑进来，对当前像素进行某种加权处理。所以说，“积”是全局概念，或者说是一种“混合”，把两个函数在时间或者空间上进行混合。
进行“卷”（翻转）的目的其实是施加一种约束，它指定了在“积”的时候以什么为参照。在信号分析的场景，它指定了在哪个特定时间点的前后进行“积”，在空间分析的场景，它指定了在哪个位置的周边进行累积处理。

下面对卷积这种操作从信号分析的角度举个例子：

如下图所示，输入信号是 f(t) ，是随时间变化的。系统响应函数是 g(t) ，图中的响应函数是随时间指数下降的，它的物理意义是说：如果在 t=0 的时刻有一个输入，那么随着时间的流逝，这个输入将不断衰减。换言之，到了 t=T时刻，原来在 t=0 时刻的输入f(0)的值将衰减为f(0)g(T)。

考虑到信号是连续输入的，也就是说，每个时刻都有新的信号进来，所以，最终输出的是所有之前输入信号的累积效果。如下图所示，在T=10时刻，输出结果跟图中带标记的区域整体有关。其中，f(10)因为是刚输入的，所以其输出结果应该是f(10)g(0)，而时刻t=9的输入f(9)，只经过了1个时间单位的衰减，所以产生的输出应该是 f(9)g(1)，如此类推，即图中虚线所描述的关系。这些对应点相乘然后累加，就是T=10时刻的输出信号值，这个结果也是f和g两个函数在T=10时刻的卷积值。

信号重建

将模拟信号转换为数字信号的过程，在信号处理领域叫做采样(Sampling)。我们可以形象地把这个过程理解为使用一连串宽度非常窄的脉冲和输入信号相乘。而得到的结果则是一连串时间上不连续的脉冲。我们在信号重建的时候，只需要在频谱上做一个低通滤波，把那些因为对时域信号采样引起的频域周期延展出来的频率过滤掉，得到的就是原始的信号。而根据傅立叶变换的性质，在频域上乘积，等价于在时域上的卷积。而低通滤波器，可以近似看为一个矩形函数。Sinc函数的傅立叶变换，正好近似矩形函数。

下面的matlab代码实现了利用sinc函数恢复信号的过程，读者可以调整参数体会整个插值恢复信号的过程。

%参数设定
T = 2;
f1 = 5;
f2 = 10;
fs = 40;
%可得参数
N = T * fs;
t = linspace(-T/2, T/2, 100*N);
f_ori = 0.5 * cos(2 * pi * f1 * t) - sin(2 * pi * f2 * t) + 1;
figure(1);
plot(t, f_ori);xlim([-T/10, T/10]);title('原始信号');
%采样
t2 = linspace(-T/2, T/2, N);
f_sam = 0.5 * cos(2 * pi * f1 * t2) - sin(2 * pi * f2 * t2) + 1;
figure(2);
stem(t2, f_sam);xlim([-T/10, T/10]);title('采样信号');
%恢复
y = [];
for i = 1 : length(t)x = t(i);h = sinc((x - t2).*fs);g = dot(f_sam, h);y = [y,g];
end
figure(3);
plot(t, y);xlim([-T/10, T/10]);title('恢复信号');

下面的图展示了采样信号与sinc函数卷积获得原始信号的过程，供读者对整个的恢复过程有一个直观的印象。

利用sinc函数实现扣点

扣点可以理解成插值的逆操作，我们不仅需要把我们想扣的点置0，还需要将这个点对其余点的影响消除。我们以对FFT的结果将能量最大的点进行扣选，主要的步骤如下：

将FFT的结果中对应的abs值中最大点的位置（下标）和他的频率值（复数）保存。
根据获取的位置移动sinc函数数组，然后乘以它的频率值，获得与原FFT结果一样点数的参考信号数组。
将原FFT结果与获得的参考信号数组做差，即可实现扣点操作。

下面是实现扣点的一种C实现，代码中中文注释帮助您更好的理解整个扣点的过程，代码能够实现对同一FFT结果进行多次的最大值扣点，可以看出代码将FFT结果剩余的所有频点的能量均值乘以一个系数作为阈值，筛选了当前所选取的最大点是否能代表当前目标的某一维度（一般为角度）信息：

void simpleDPK(STRUCT_DPKRST *Sresult,STRUCT_DPKCFG *dpkCfg, float *din,float *dinabs)
{float maxVal[2]={0.0};uint16_t maxIdx = 0;uint16_t maxIdx_left0 = 0;uint16_t maxIdx_right0 = 0;uint8_t Offset = 0;uint8_t dpkTime=dpkCfg->dpkTime;uint8_t thres=dpkCfg->thres;uint8_t sincLen=dpkCfg->sincLen;uint8_t sNumANT=dpkCfg->numAnt;float scalePow=1.0/sincLen;float sRes[12] = {0.0}; //默认扣最大12次Sresult->targetNum=0;arm_cmplx_mag_f32(din,dinabs, sincLen); //先做1次absfor(uint8_t k=0;k<dpkTime*3;k=k+3)  //扣点次数{	for(uint32_t j=1;j<sincLen;j++)			//拿到dinabs最大值及索引{if(dinabs[maxIdx]<dinabs[j])//拿到最大值，index是对应索引{maxIdx = j;}					}		maxVal[0] = din[maxIdx*2];		//拿最大实部maxVal[1] = din[maxIdx*2+1];	//拿最大虚部maxVal[0] = din[maxIdx*2];		//拿最大实部maxVal[1] = din[maxIdx*2+1];	//拿最大虚部//取左右两边索引值if (maxIdx == 0){maxIdx_left0 = sincLen-1;maxIdx_right0 = maxIdx+1;}else if (maxIdx == sincLen-1){maxIdx_left0 = maxIdx-1;maxIdx_right0 = 0;					}else {maxIdx_left0 = maxIdx-1;maxIdx_right0 = maxIdx+1;	}	for(uint8_t im=0;im<sincLen;im++) {Offset=im+sincLen-maxIdx;	//拿移位索引if(Offset>=sincLen){Offset=Offset-sincLen; }			SincSeq[im][0]= SincBuf[Offset][0]*maxVal[0]-SincBuf[Offset][1]*maxVal[1];SincSeq[im][1] = SincBuf[Offset][1]*maxVal[0]+SincBuf[Offset][0]*maxVal[1];din[im*2]=din[im*2]-SincSeq[im][0];din[im*2+1]=din[im*2+1]-SincSeq[im][1];}Sresult->idx[k]=maxIdx;   //先存最大目标索引Sresult->pow[k]=dinabs[maxIdx];  //先存最大目标能量值		Sresult->idx[k+1]=maxIdx_left0;   //先存最大目标索引Sresult->pow[k+1]=dinabs[maxIdx_left0];  //先存最大目标能量值				Sresult->idx[k+2]=maxIdx_right0;   //先存最大目标索引Sresult->pow[k+2]=dinabs[maxIdx_right0];  //先存最大目标能量值		arm_cmplx_mag_f32(din,dinabs, sincLen); //拿到ABS值sRes[k] = 0.0;for(uint8_t im=0;im<sincLen;im++) {sRes[k]=sRes[k]+dinabs[im];}sRes[k+1]=sRes[k];sRes[k+2]=sRes[k];
//		printf("k=%d, idx = %d, pow = %f, sRes = %f\r\n",k,Sresult->idx[k],Sresult->pow[k],sRes[k]);}for(uint8_t k=0;k<dpkTime*3;k=k+3){if(Sresult->pow[k] > sRes[0]/sincLen/(sNumANT-dpkTime) * sNumANT * thres){if (Sresult->pow[k] > 1500000) {for (uint8_t kk=k;kk<(k+3);kk++){Sresult->pow[Sresult->targetNum] = Sresult->pow[kk] /(sRes[kk]/sincLen/(sNumANT-dpkTime) * sNumANT );//		printf("case1: kk = %d, itargetNum=%d, snr = %f\r\n",kk, Sresult->targetNum,Sresult->pow[Sresult->targetNum]);	Sresult->targetNum=Sresult->targetNum+1;}}else{Sresult->pow[k] = Sresult->pow[k] /(sRes[dpkTime-1]/sincLen/(sNumANT-dpkTime) * sNumANT );//		printf("case2: kk = 1, itargetNum=%d, snr = %f\r\n", Sresult->targetNum,Sresult->pow[Sresult->targetNum]);	Sresult->targetNum=Sresult->targetNum+1;}}}
}