本文编辑：@调皮连续波，保持关注调皮哥，获得更多学习内容和建议！
感谢大家前来捧场，我是调皮哥！
今天继续来分享干货，我认为一个东西，只有把它研究彻底，了解它的全部，才能够真正理解并掌握它。今天调皮哥采用一种新的方法来处理雷达信号处理中的静态杂波，这个方法是圆拟合，粉丝问我是否可以，感谢粉丝提出的问题，那么我今天就分享我验证的结果。
首先我们需要了解什么是圆拟合？
1.圆拟合
简单的说，就是给你一堆离散的点，然后根据这些离散点的特点，拟合出一个圆，思想和最小二乘法线性拟合相似，模型如下图所示。
2.为什么可以采用圆拟合？
为什么用圆拟合算法就能够消除静态杂波干扰（直流分量）呢？这个问题是本文的核心问题，只有理解了这个才能明白圆拟合算法的本质。
我们知道，雷达回波信号经过正交采样后得到IQ复信号，时域IQ信号在复平面上的分布是随着采样点的幅度和相位散落在复平面上，如果不存在直流分量，那么离散点的是以零点为圆心，以幅度为半径，以随时间变化的相位为角度，呈圆弧分布。
如果受到静态物体产生的杂波干扰，由于信号存在直流分量，所以估计的相位随时间演变的幅值往往小于真实值，且杂波信号强度越大，偏差值越大，体现在图形中，就如下图所示，A1就是直流分量的幅度，可见A0因为受到直流分量的干扰，圆心偏离零点很远。
因此从信号处理的角度解决实际应用中静态杂波的干扰问题，因此可以利用信号后处理方法去补偿、校正相位时间序列的估计偏差，即消除A1的幅度，使得A0的圆心重回零点。
关于这部分内容的深刻理解，考虑到文章篇幅和时间关系，我暂时不做深入分析，下一篇文章会详细从底层解释，建议读者先反复深入理解傅里叶变换、复信号平面、直流和交流信号、向量叠加等内容。
3.雷达信号处理圆拟合算法步骤
如果你了解了上面的原理后，其实算法程序编写起来很简单，步骤如下：
（1）找到当前离散点在复平面上的分布
（2）对这些点进行圆拟合操作
（3）根据拟合后的圆，补偿偏移
（4）根据补偿后的相位和幅度，得到新的离散点。
4.圆拟合数据说明与处理结果
MATLAB版本：2021a
操作系统：Win10
数据格式：脉冲数*采样点数=128*128的时域信号，是两个目标在3米附近相向而行。调皮哥经过调试，得到的结果如下图所示，圆拟合采用最小二乘法实现。
5.处理效果对比与分析
说实话，其实效果一般，比起均值相消算法，圆拟合其实效果并不那么好，首先在距离维上的直流分量去除还有少量部分存在，速度维上也有少量直流分量（零速通道）没有消除，如下图所示。
总体看起来效果还能接受，不过就是计算量很大，不利于工程师实现。
6.程序分析
（1）雷达参数设置部分
主要是根据采集目标信号时所设置的参数，各位读者可以根据自己的参数修改。建议初学者不要将采样点数和脉冲数设置为一样的，便于调试。
%% 雷达参数
Tx_Number = 2; %发射天线
Rx_Number = 4; %接收天线
Range_Number = 128; %距离点数（每个脉冲128个点）
Doppler_Number = 128; %多普勒通道数(总共128个重复脉冲数)
global Params;
Params.NChirp = Doppler_Number; %1帧数据的chirp个数
Params.NChan = Rx_Number; %RxAn数,ADC通道数
Params.NSample = Range_Number; %每个chirp ADC采样数
Params.Fs = 2.5e6; %采样频率
Params.c = 3.0e8; %光速
Params.startFreq = 77e9; %起始频率
Params.freqSlope = 60e12; %chirp的斜率
Params.bandwidth = 3.072e9; %真实带宽
Params.lambda=Params.c/Params.startFreq; %雷达信号波长
Params.Tc = 144e-6; %chirp周期
global FFT2_mag;（2）坐标计算
这部分就是根据公式计算绘制图形的坐标，挺简单的，根据雷达测距测速的原理公式就能够明白了。
%% 坐标计算
[X,Y] = meshgrid(Params.c*(0:Params.NSample-1)*Params.Fs/2/Params.freqSlope/Params.NSample, ...
(-Params.NChirp/2:Params.NChirp/2 - 1)*Params.lambda/Params.Tc/Params.NChirp/2);（3）距离维直流分量去除
首先load加载原始数据，赋值给fft1d_before变量，然后建立两个会用到的空数组。
load ReIm_Data_All.mat ;
fft1d_before=ReIm_Data_All;
ai=zeros(Range_Number,Doppler_Number);
data=zeros(Range_Number,Doppler_Number,Tx_Number*Rx_Number);距离维圆拟合程序如下，其实是同时对8个天线进行拟合，因此解算的速度会慢一些，其中程序有绘图的部分我给注释了，读者可以在自己调试时打开并观察结果。拟合的算法采用的是最小二乘法，其原理请各位读者自行研究。
for antenna=1:Tx_Number*Rx_Number
for Range=1:Range_Number
%1.估计每个扫频周期时间内的基带复信号中频信号的幅值
AmR(Range,:)=fft1d_before(Range,:,antenna);
%2.幅值时间序列和已知的初始相位时间序列得到复平面上的离散点
% figure(2);
% plot(fft1d(:,doppler,1),'o');
% title([num2str(doppler)]);
%3.一个距离门上的所有多普勒点进行圆拟合
xdataR=real(AmR(Range,:));
ydataR=imag(AmR(Range,:));
%最小二乘法拟合
k0 = ones(1,3);
F = @(k)(xdataR-k(1)).^2+(ydataR-k(2)).^2-k(3)^2;
[k,resnorm] = lsqnonlin(F,k0);
%k(1)是圆心的x坐标
%k(2)是圆心的y坐标
%k(3)的绝对值是圆的半径
% r0 = [k(1),k(2)];
% R = abs(k(3));
% xx = k(1)-R:0.01*R:k(1)+R;
% y1_h = sqrt(R.^2 - (xx - r0(1)).^2) + r0(2);
% y2_h = -sqrt(R.^2 - (xx - r0(1)).^2) + r0(2);
% figure(1);
% plot(xx,y1_h,'b')
% hold on
% plot(xx,y2_h','b')
% plot(xdata,ydata,'*r')
% title('距离维圆拟合');
% xlabel('实部');
% ylabel('虚部');
% axis equal %axis square
%4.修正补偿
%获取拟合圆的圆心
x=k(1);
y=k(2);
%将圆心移到零点（0，0）
xdataR=xdataR-x;
ydataR=ydataR-y;
%5.得到新的点的时间序列相位
dataR(Range,:,antenna)=complex(xdataR,ydataR);
hold off;
end
end（4）速度维直流分量去除
速度维和距离维的算法一样，其实就是采用两次圆拟合运算，这里为了方便理解，我分开做计算，读者可以优化这个程序步骤，使得算法运算更加有效率。
%% 速度维 圆拟合
AmV=zeros(Range_Number,Doppler_Number);
dataV=zeros(Range_Number,Doppler_Number,Tx_Number*Rx_Number);
for antenna=1:Tx_Number*Rx_Number
for doppler=1:Doppler_Number
%1.估计每个扫频周期时间内的基带复信号中频信号的幅值
AmV(:,doppler)=dataR(:,doppler,antenna);
%2.幅值时间序列和已知的初始相位时间序列得到复平面上的离散点
% figure(2);
% plot(fft1d(:,doppler,1),'o');
% title([num2str(doppler)]);
%3.一个距离门上的所有多普勒点进行圆拟合
xdataV=real(AmV(:,doppler));
ydataV=imag(AmV(:,doppler));
%最小二乘法拟合
k0 = ones(1,3);
F = @(k)(xdataV-k(1)).^2+(ydataV-k(2)).^2-k(3)^2;
[k,resnorm] = lsqnonlin(F,k0);
%k(1)是圆心的x坐标
%k(2)是圆心的y坐标
%k(3)的绝对值是圆的半径
% r0 = [k(1),k(2)];
% R = abs(k(3));
% xx = k(1)-R:0.01*R:k(1)+R;
% y1_h = sqrt(R.^2 - (xx - r0(1)).^2) + r0(2);
% y2_h = -sqrt(R.^2 - (xx - r0(1)).^2) + r0(2);
% figure(1);
% plot(xx,y1_h,'b')
% hold on
% plot(xx,y2_h','b')
% plot(xdata,ydata,'*r')
% axis equal %axis square
%4.修正补偿
%获取拟合圆的圆心
x=k(1);
y=k(2);
%将圆心移到零点（0，0）
xdataV=xdataV-x;
ydataV=ydataV-y;
%5.得到新的点的时间序列相位
dataV(:,doppler,antenna)=complex(xdataV,ydataV);
hold off;
end
end（5）距离维和速度维FFT
在这里我们其实清楚地知道，两次圆拟合算法都是在FFT运算之前做的。本程序还有很多图形暂时注释，读者可以自行画出分析，比如拟合前、拟合后的图形以及信号波形等，我个人觉得读者应该有较高的编程能力，这些都不是困难的事情，因此不再繁琐列出。
%% 1D FFT
fft1d= zeros(Range_Number,Doppler_Number,Tx_Number*Rx_Number);
for antenna =1:Tx_Number*Rx_Number
for Range=1:Range_Number
fft1d(Range,:,antenna) = fft((dataV(Range,:,antenna)));
end
end
FFT1_mag=abs(fft1d(:,:,1));
figure();
mesh(FFT1_mag);
xlabel('采样点数');ylabel('脉冲数');zlabel('幅度');
title('圆拟合 1D-FFT结果');
%% 2D-FFT
fft2d= zeros(Range_Number,Doppler_Number,Tx_Number*Rx_Number);
for antenna=1:Tx_Number*Rx_Number
for doppler=1:Doppler_Number
fft2d(:,doppler,antenna) =fftshift( fft((fft1d(:,doppler,antenna))));
end
end
FFT2_mag=(abs(fft2d(:,:,1)));
figure();
mesh(X,Y,FFT2_mag);
xlabel('距离维(m)');ylabel('速度维(m/s)');zlabel('幅度');
title('圆拟合 2D-FFT结果');感谢阅读，觉得还不错，可以给调皮哥一个打赏，鼓励调皮哥继续创作，喜欢本文也可以点赞、收藏、转发，让更多的人一起学雷达（加粗的谢谢）！
参考资料：
[1]一种FMCW雷达静态杂波干扰消除的信号处理方法 彭志科
[2]Dynamic Digital Signal Processing Algorithm for Vital Signs Extraction in Continuous-Wave Radars
上述这两篇文献强烈建议大家阅读，其中就是用到了本文所论述的这种方法。
程序和数据下载方式：后台回复：“0607”获得下载链接。
