前言：雷达信号处理的一般流程：ADC数据——1D-FFT——2D-FFT——CFAR检测——测距、测速、测角。本篇笔记主要记录CFAR的原理和代码实现。
一、CFAR原理
1、一维CFAR
（1）原理：仅对距离维或多普勒维做CFAR，检测单元的左右两侧均有设定好数量的保护单元和参考单元。如果是CA-CFAR则取左右两边的参考单元的均值再取平均值，然后与
相乘（计算公式如下方代码块所示），再将结果与阈值相比较。如果大于阈值则代表有目标，否则认为无目标。
（2）特殊情况
如果检测单元位于RDM（做完2D-FFT的矩阵）的边界，可以使用补0或翻转的方式。
补零法：若检测单元位于RDM矩阵的第一个点，需要额外增加RDM谱的大小，使得能够给该检测单元分配保护单元以及参考单元，补零示意图如下所示：
(2)翻转法
翻转法，是指不采用补零的方式扩展矩阵，而是采用RDM谱边界上的点相互弥补。例如，以RDM谱第一个点为被检测单元，但是左边和上边都缺少保护单元和参考单元，这个时候可以利用RDM谱矩阵最右边的点和最下边的点来弥补不足的点，从而形成一个闭合的“球状”数据体，如图13所示。
这种方式不通过补零，直接采RDM谱的原始信息，更能保证检测的稳定性，但这种方式会使得编程变得稍微复杂，随着迭代变化，翻转单元的形式也是不断变化的。
2、二维CFAR
可以分为两种形式实现：（1）先对距离维或多普勒维做1D-CFAR，然后对多普勒维或距离维再做1D-CFAR；（2）直接同时对距离维和多普勒维做CFAR，保护单元和参考单元不再是一维的了，变成上下左右形成二维的，如下图所示：
二、代码实现
转载自文章2：
%% 微信公众号：调皮的连续波
%% 知乎：调皮连续波
%% 作者：调皮哥
%% 时间：2022年5月
%% 功能：FMCW雷达发射信号、回波信号、混频、距离维FFT、速度维FFT建模仿真。
%%=========================================================================
clear all;
close all;
clc;
%% 雷达系统参数设置
maxR = 200; % 雷达最大探测目标的距离
rangeRes = 1; % 雷达的距离分率
maxV = 70; % 雷达最大检测目标的速度
fc= 77e9; % 雷达工作频率 载频
c = 3e8; % 光速
%% 用户自定义目标参数
r0 = 90; % 目标距离设置 (max = 200m)
v0 = 20; % 目标速度设置 (min =-70m/s, max=70m/s)
%% FMCW波形参数设置
B = c / (2*rangeRes); % 发射信号带宽 (y-axis) B = 150MHz
Tchirp = 5.5 * 2 * maxR/c; % 扫频时间 (x-axis), 5.5= sweep time should be at least 5 o 6 times the round trip time
slope = B / Tchirp; %调频斜率
endle_time=6.3e-6; %空闲时间
f_IFmax= (slope*2*maxR)/c ; %最高中频频率
f_IF=(slope*2*r0)/c ; %当前中频频率
Nd=128; %chirp数量
Nr=1024; %ADC采样点数
vres=(c/fc)/(2*Nd*(Tchirp+endle_time));%速度分辨率
% Nr=1024*256; %和频信号点数设置
Fs=Nr/Tchirp; %模拟信号采样频率
t=linspace(0,Nd*Tchirp,Nr*Nd); %发射信号和接收信号的采样时间，在MATLAB中的模拟信号是通过数字信号无限采样生成的。
Tx=zeros(1,length(t)); %发射信号
Rx=zeros(1,length(t)); %接收信号
Mix = zeros(1,length(t)); %差频、差拍、拍频、中频信号
r_t=zeros(1,length(t));
td=zeros(1,length(t));
%% 动目标信号生成
for i=1:length(t)
r_t(i) = r0 + v0*t(i); % 距离更新
td(i) = 2*r_t(i)/c; % 延迟时间
Tx(i) = cos(2*pi*(fc*t(i) + (slope*t(i)^2)/2)); % 发射信号 实数信号
Rx(i) = cos(2*pi*(fc*(t(i)-td(i)) + (slope*(t(i)-td(i))^2)/2)); %接收信号 实数信号
if i<=1024
freq(i)=fc+slope*i; %发射信号时频图 只取第一个chirp
freq_echo(i)=fc+slope*i;%回波信号频谱延迟
end
Mix(i) = Tx(i).*Rx(i);%差频、差拍、拍频、中频信号
end
%发射信号时域图
figure;
plot(Tx(1:1024));
xlabel('点数');
ylabel('幅度');
title('TX发射信号时域图');
% %发射信号时频图
figure;
plot(t(1:1024),freq);
xlabel('时间');
ylabel('频率');
title('TX发射信号时频图');
%接收信号时域图
figure;
plot(Rx(1:1024));
xlabel('点数');
ylabel('幅度');
title('RX接收信号时域图');
%接收信号与发射信号的时频图
figure;
plot(t(1:1024),freq);
hold on;
plot(td(1:1024)+t(1:1024),freq);
xlabel('时间');
ylabel('频率');
title('接收信号与发射信号时频图');
legend ('TX','RX');
%中频信号频谱 和频信号观察
%figure;
% plot(db(abs(fft(Mix(1:1024*256)))));%查看宽带的和频信号 将chirp的点数改为1024*256即可看到有一个门信号，但注意计算机内存。
% xlabel('频率');
% ylabel('幅度');
% title('中频信号频谱');
figure;
plot(db(abs(fft(Mix(1:1024)))));%查看宽带的和频信号 将chirp的点数改为1024*256即可看到有一个门信号，但注意计算机内存。
xlabel('频率');
ylabel('幅度');
title('中频信号频谱');
%% 低通滤波 截止频率30MHz 采样频率120MHz
% Mix=lowpass(Mix(1:1024*256),30e6,120e6);
% plot(db(abs(fft(Mix(1:1024*256)))));
% xlabel('频率');
% ylabel('幅度');
% title('中频信号低通滤波器');
%reshape the vector into Nr*Nd array. Nr and Nd here would also define the size of
%Range and Doppler FFT respectively.
signal = reshape(Mix,Nr,Nd);
figure;
mesh(signal);
xlabel('脉冲数')
ylabel('距离门数');
title('中频信号时域');
%% 距离维FFT
sig_fft = fft(signal,Nr)./Nr;
sig_fft = abs(sig_fft);
sig_fft = sig_fft(1:(Nr/2),:);
figure;
plot(sig_fft(:,1));
xlabel('距离（频率）');
ylabel('幅度')
title('第一个chirp的FTF结果')
%% 距离FFT结果谱矩阵
figure;
mesh(sig_fft);
xlabel('距离（频率）');
ylabel('chirp脉冲数')
zlabel('幅度')
title('距离维FTF结果')
%% 速度维FFT
Mix=reshape(Mix,[Nr,Nd]);
sig_fft2 = fft2(Mix,Nr,Nd);
sig_fft2 = sig_fft2(1:Nr/2,1:Nd);
sig_fft2 = fftshift (sig_fft2);
RDM = abs(sig_fft2);
RDM = 10*log10(RDM) ;
doppler_axis = linspace(-100,100,Nd);
range_axis = linspace(-200,200,Nr/2)*((Nr/2)/400);
figure;
mesh(doppler_axis,range_axis,RDM);
xlabel('多普勒通道'); ylabel('距离通道'); zlabel('幅度（dB）');
title('速度维FFT 距离多普勒谱');
%% CFAR检测算法
%通过完整的距离多普勒图滑动窗口
%在两个维度中选择参考单元的数量
Tr = 8;
Td = 4;
%选择被测单元（CUT）周围两个维度的保护单元数量，以进行准确CFAR检测
Gr = 4;
Gd = 2;
%以 dB 为单位的 SNR 值偏移阈值
snr_offset = -10*log10(0.25);
%为参考单元上的每次迭代创建一个向量来存储noise_level（噪声电平）
r_margin = 2*(Tr+Gr);
d_margin = 2*(Td+Gd);
r_grid_length = Nr/2-d_margin;
d_grid_length = Nd-r_margin;
noise_level = zeros(r_grid_length,d_grid_length);
%设计一个循环，通过在参考单元和保护单元的边缘提供边距，
%使 CUT 在距离多普勒图上滑动。 对于每次迭代，求所有参考单元中信号电平的和。
%使用 db2pow 函数将值从对数转换为线性，并平均所有使用的参考单元格的总和值。
%平均后使用 pow2db 将其转换回对数。 进一步添加偏移量以确定阈值。
%接下来，将 CUT 下的信号与此阈值进行比较。 如果 CUT > 阈值，则为其分配值 1，否则置为0。
gridSize = (2*Tr+2*Gr+1)*(2*Td+2*Gd+1);
numGcells = (2*Gr+1)*(2*Gd+1);
numTcells = gridSize - numGcells; %参考单元
% 上述过程将生成一个阈值块，
%该块小于距离多普勒图，因为 CUT 不能位于矩阵的边缘。
%因此，很少有单元不会被阈值化。 要保持MAP大小相同，请将这些值设置为 0。
sig_CFAR = zeros(size(RDM));
% 此方法计算量大，高于常规的在速度维做一次CFAR 然后在距离维再做一次CFAR。
%另外，此方法对于距离-多普勒谱矩阵的边界处存在未检测，因此对于处于边界的目标效果不友好，
%可以采用补零法，以及旋绕法检测边界目标
for i = 1:r_grid_length % 距离边界
for j = 1:d_grid_length % 多普勒边界
% 使用 db2pow 将值从对数转换为线性
% sig_pow = db2pow(RDM(i:i+d_margin,j:j+r_margin)); %从db变为线性
sig_pow = db2pow(RDM(i:i+d_margin,j:j+r_margin)); %从db变为线性
G_pow = db2pow(RDM(i+Td:i+Td+Gd*2,j+Tr:j+Tr+Gr*2)); %保护单元
% 所有参考单元内的总和信号
sig_sum = sum(sum(sig_pow))-sum(sum(G_pow));
% 对所有使用的参考单元格的总和值求平均
noise_level(i,j) = pow2db(sig_sum/numTcells);
% 添加偏移量
sig_threshold = noise_level(i,j) + snr_offset;
% 将 CUT 与阈值进行比较
if (RDM(i+d_margin/2, j+r_margin/2) > sig_threshold)
sig_CFAR(i+d_margin/2, j+r_margin/2) = 1;
end
end
end
%display the CFAR output using the Surf function like we did for Range
figure;
mesh(doppler_axis,range_axis,sig_CFAR);
colorbar;
xlabel('多普勒'); ylabel('距离');
title('CFAR 结果');
三、参考文章
1、雷达无线电系列（二）经典CFAR算法图文解析与实现（matlab） - 万雨 - 博客园 (cnblogs.com)2、干货 | FMCW雷达信号处理的二维CFAR（2D CFAR、十字CFAR）检测算法 - 知乎 (zhihu.com)