前言：雷达信号处理的一般流程：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：
%% 微信公众号：调皮的连续波%% 知乎：调皮连续波%% 作者：调皮哥%% 时间：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 = 150MHzTchirp = 5.5 * 2 * maxR/c; % 扫频时间 (x-axis), 5.5= sweep time should be at least 5 o 6 times the round trip timeslope = 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<=1024freq(i)=fc+slope*i; %发射信号时频图 只取第一个chirpfreq_echo(i)=fc+slope*i;%回波信号频谱延迟endMix(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('中频信号时域');%% 距离维FFTsig_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结果')%% 速度维FFTMix=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 endend%display the CFAR output using the Surf function like we did for Rangefigure;mesh(doppler_axis,range_axis,sig_CFAR);colorbar;xlabel('多普勒'); ylabel('距离');title('CFAR 结果');
三、参考文章
1、雷达无线电系列（二）经典CFAR算法图文解析与实现（matlab） - 万雨 - 博客园 ()2、干货 | FMCW雷达信号处理的二维CFAR（2D CFAR、十字CFAR）检测算法 - 知乎 ()
