极简 FMCW 毫米波雷达 Matlab 仿真

草草不工

已于 2024-08-26 07:42:40 修改

阅读量1k

点赞数 28

分类专栏：信息与信号文章标签： matlab 信号处理汽车信息与通信

于 2024-08-26 07:41:44 首次发布

本文链接：https://blog.csdn.net/qq_44257974/article/details/141233909

版权

信息与信号专栏收录该内容

4 篇文章 0 订阅

订阅专栏

文章目录

引言
FMCW 雷达基本原理——如何测距
MATLAB仿真——测距
FMCW 雷达基本原理——如何测速
MATLAB仿真——测速
FMCW 雷达基本原理——雷达目标检测
MATLAB仿真——雷达目标检测

引言

本文旨在从信号角度仿真一个极简的 FMCW 毫米波雷达，用以理解毫米波雷达的基本原理，也为相关的数字信号处理打下基础。我们会一边看毫米波雷达实现的基本流程及原理，一边会对应着Matlab代码一点点仿真。

代码参考：https://github.com/davidscmx/radar-target-generation-and-detection

FMCW 雷达基本原理——如何测距

TI 的这篇文章讲的相当之好，以至于这一章节完全可以跳过，但此处还是略微赘述几笔，我们先只考虑距离的测量。

一个电磁波发射出去，再撞到外界物体被反弹回来，被弹回的波和发送的波之间一定有时间差，这个时间差就是我们测量距离的关键。问题是，我们应该发送一个怎样的波？FMCW（Frequency-Modulated Continuous-Wave，频率调制连续波）雷达发送的是一个频率不断变化的波形，如下图所示。

在这里插入图片描述

我们如果从频率角度来看这个波形，这是一个频率线性增长的波形，即 $f(t)=f_0+\Delta ft$ ，其中 $f_0$ 是初始频率， $\Delta f$ 是频率变化率， $t$ 是时间。可以做出下图。

在这里插入图片描述

对于接收到的波形，它会与发射波形有个时延，如下图所示。这时可以发现，在同一时刻，两个波形的频率差是相同的，这意味着在接收时，我们需要想方设法获得这个频率差，并以此判定两者的距离。如何通过频率差来判定距离呢？知道频率差，知道下图中的斜率，我们也就知道了时间差 τ，而电磁波传播的速度是光速，距离由此可得到。

在这里插入图片描述

现在的问题是，频率从多少开始，到多少结束？现在主流的车载毫米波雷达都在24GHz或是77GHz，这是频谱资源、测量精度传感器尺寸等多方面共同决定的。那如果开始是77GHz，频率应该到多少结束呢，这两者的差值也就是毫米波雷达的带宽。

带宽是上图中纵坐标的差值，其决定了雷达的距离分辨率；而雷达所能探测的最大距离，则由横轴上的 $T_{chirp}$ 决定。当然，对于chirp信号而言，带宽、周期和斜率三项中，由两项决定则波形就被确定了。

我们先来看雷达的最大作用距离：

$R_{max}≤\frac{cT_{chirp}}{2}$

当然这个仅仅是最理想条件下的，是基于信号数学意义上的决定式，真实的最大作用距离必须从功率的角度进行评估。上式理解起来也很容易，将上图中蓝色的点不断右移，直到发送和接受的两个信号在完全不同的时间段。

下面来看雷达的距离分辨率，首先，在面临多个物体放射波时，会出现下图情况。想想一下，若是两个物体相距很近，返回的两个频率自然也非常接近，在频率域上，可能就会无法分辨。

在这里插入图片描述

一个结论是：对回波信号的观察周期越长，信号处理后得到的频率分辨效果就越好。一般来说，观察时长为 $T$ 的一个观察窗可以分离相隔超过 $\frac{1}{T} Hz$ 的频率分量。

在这里插入图片描述

MATLAB仿真——测距

定义雷达参数

定义基本的雷达参数，包括：

光速： $3\times 10^8 m/s$
最大测量范围：200m（这个参数会影响扫频周期 Tchirp）
目标速度：-20m
目标距离范围：110m
距离分辨率：1m（这个参数与观测窗口T有关，会决定带宽）

%% Radar Specifications 
%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Frequency of operation = 77GHz
% Max Range = 200m
% Range Resolution = 1 m
% Max Velocity = 100 m/s
%%%%%%%%%%%%%%%%%%%%%%%%%%%

c = 3e8;
range = 110;
vel = -20;
max_range = 200;
range_res = 1;
max_vel = 100;

下面定义单周期chirp波形的基本参数：

通过计算得到了FMCW波形的带宽B
扫频时间
斜率

这三个参数（准确的说是其中任意两个）共同决定了连续调频波的波形，chirp 信号的三角波由此决定。

%% FMCW Waveform Generation

% *%TODO* :
%Design the FMCW waveform by giving the specs of each of its parameters.
% Calculate the Bandwidth (B), Chirp Time (Tchirp) and Slope (slope) of the FMCW
% chirp using the requirements above.
B = c / (2*range_res);
Tchirp = 5.5 * 2 * (max_range/c) ; 
slope = B/Tchirp;

下面定义chirp序列的一些参数：

雷达频率：77GHz，这也是chirp信号的起始频率
Nd 代表一个序列中有多少个chirp信号
Nr 代表在每个chirp中有多少个采样点

%Operating carrier frequency of Radar 
fc= 77e9;             %carrier freq
                                                          
%The number of chirps in one sequence. Its ideal to have 2^ value for the ease of running the FFT
%for Doppler Estimation. 
Nd=128;                   % #of doppler cells OR #of sent periods % number of chirps

%The number of samples on each chirp. 
Nr=2^15;                  %for length of time OR # of range cells

下面定义了后续所需要的一些向量空间，包括：

时间轴：从0到chirp总序列完成的时间，共有 Nr*Nd 个点
Tx、Rx、Mix（混频后的信号），依据时间轴的点数进行向量创建
rt 是每个时刻的距离，由初始位置，速度与时间决定
td 是时间延迟，等于两倍的距离除以光速

% Timestamp for running the displacement scenario for every sample on each
% chirp
t=linspace(0,Nd*Tchirp,Nr*Nd); %total time for samples

%Creating the vectors for Tx, Rx and Mix based on the total samples input.
Tx=zeros(1,length(t)); %transmitted signal
Rx=zeros(1,length(t)); %received signal
Mix = zeros(1,length(t)); %beat signal

%Similar vectors for range_covered and time delay.
r_t=zeros(1,length(t));
td=zeros(1,length(t));

生成发射及接收波形

现在可以生成chrip波形了：

发射信号： $cos(2\pi(f_c\cdot t(i) + \frac{slope\cdot t(i)^2}{2} ) )$
接收信号是对于发射信号的右移： $cos(2\pi(f_c\cdot(t(i) -t_d(i)) + \frac{slope\cdot(t(i)-t_d(i))^2}{2} ) )$

%% Signal generation and Moving Target simulation
% Running the radar scenario over the time. 

for i=1:length(t)         
  % *%TODO* 
  %For each time stamp update the Range of the Target for constant velocity. 
  r_t(i) = range + (vel*t(i));
  td(i) = (2 * r_t(i)) / c;

  % *%TODO* :
  %For each time sample we need update the transmitted and
  %received signal. 

  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 ) );
    
  % *%TODO* :
  %Now by mixing the Transmit and Receive generate the beat signal
  %This is done by element wise matrix multiplication of Transmit and
  %Receiver Signal

  Mix(i) = Tx(i) .* Rx(i);
end

绘制信号：

直接在时域上分辨77GHz与77.15GHz是困难的，所以我们把chirp信号求一个FFT，以此来观察其是否符合我们的预期。
Matlab的FFT函数返回的是一个复数数组，我们需要根据采样率计算得到频率域函数的横坐标；还需要对FFT结果求模以得到频率域函数的纵坐标。
下图中可以看到我们生成出的Tx在频谱上情况。

% 绘制一个chrip的Tx信号的频率域结果
signal_fft = fft(Tx, Nr);
N = length(signal_fft);
Fs=Nr/Tchirp;
f = (0:N-1)*(Fs/N);
plot(f(1:N/2), abs(signal_fft(1:N/2)));
title('Single-Sided Amplitude Spectrum');
xlabel('Frequency (Hz)');
ylabel('|FFT|');

% 绘制混频后的信号的第一个脉冲
plot(t(1:1*Nr), Mix(1:1*Nr));
title('Beat Signal');
xlabel('Time');
ylabel('Amplitude');

在这里插入图片描述

接收信号的处理——测距

现在我们对接收到并混频后的信号进行fft转换，以确定基频信息，进一步通过这个频率差得到距离信息，我们在代码中打印出频谱图以及将横坐标转换成距离信息的距离图。可以看到，输出的距离与预设相符。

%% RANGE MEASUREMENT

 % *%TODO* :
%reshape the vector into Nr*Nd array. Nr and Nd here would also define the size of
%Range and Doppler FFT respectively.
Mix = reshape(Mix, [Nr, Nd]);
MixT = Mix';
dim = 2;
 
% fft运算并输出频谱图,这里只打印4行
signal_fft = fft(MixT, Nr, dim);

P2 = abs(signal_fft/Nr);
P1 = P2(:,1:Nr/2+1);
P1(:,2:end-1) = 2*P1(:,2:end-1);

% figure;
% for i=1:4
%     subplot(4,1,i)
%     plot(0:(Fs/Nr):(Fs/2-Fs/Nr),P1(i,1:Nr/2))
%     title("Row " + num2str(i) + " in the Frequency Domain")
% end

% 计算距离并直接体现在横坐标上,这里只打印4行
rowRange = (0:(Fs/Nr):(Fs/2-Fs/Nr))*c/slope/2;

figure;
for i=1:4
    subplot(4,1,i)
    plot(rowRange,P1(i,1:Nr/2))
    xlim([0,200]);
    xlabel('Range [m]');
    title("Row " + num2str(i) + " Range")
end

在这里插入图片描述

FMCW 雷达基本原理——如何测速

从刚刚的输出中，我们可以看到，在不同的时间上，由于时间间隔较短，输出的峰值位置并没有不同，但图中没有展示的是此时的相位，我们正是通过相位差来进行测速的。这里我们需要用到 2D FFT 来进行速度的求解，这个视频 很好的给出一些 2D FFT 的案例，可以帮助我们进行理解。2D FFT 在图像处理中的应用很多，其本质上是对一组信号在不同角度上进行观察。
对于我们的 chirp 信号组而言，对每个信号在频域上进行观察是一个视角，如果将视角转过来，对各组信号在同一相对时刻进行观察，就可以发现一个周期性变化的相位，这个相位经过换算就可以得到我们想要的速度。
语言表述很是乏力，建议看看上面那个链接的视频……

在这里插入图片描述

MATLAB仿真——测速

对于 2D FFT，我们首先需要 reshape 我们的混频后数组，将其组成一个由一个个时段的 chirp 信号组成的数组，这样才能从横纵坐标两个维度进行 FFT。此后这段代码进行了一系列的坐标转换，并最终生成了三维的 2D FFT 图像。

%% RANGE DOPPLER RESPONSE
% The 2D FFT implementation is already provided here. This will run a 2DFFT
% on the mixed signal (beat signal) output and generate a range doppler
% map.You will implement CFAR on the generated RDM

% Range Doppler Map Generation.
% The output of the 2D FFT is an image that has reponse in the range and
% doppler FFT bins. So, it is important to convert the axis from bin sizes
% to range and doppler based on their Max values.

Mix=reshape(Mix,[Nr,Nd]);

% 2D FFT using the FFT size for both dimensions.
signal_fft2 = fft2(Mix,Nr,Nd);

% Taking just one side of signal from Range dimension.
signal_fft2 = signal_fft2(1:Nr/2,1:Nd);
signal_fft2 = fftshift (signal_fft2);

RDM = abs(signal_fft2);
RDM = 10*log10(RDM) ; % 转换为分贝单位

%use the surf function to plot the output of 2DFFT and to show axis in both
%dimensions

doppler_axis = linspace(-100,100,Nd);
range_axis = linspace(-200,200,Nr/2)*((Nr/2)/400);

figure,surf(doppler_axis,range_axis,RDM,'EdgeColor', 'none');
title('Amplitude and Range From FFT2');
xlabel('Speed');
ylabel('Range');
zlabel('Amplitude');

在这里插入图片描述

FMCW 雷达基本原理——雷达目标检测

刚刚我们得到了 2D FFT 的结果，但还并没有完成目标检测的工作。雷达在判决过程中，可能会出现两类错误。第一类是在没有目标时判断为有目标，这类错误称为虚警。另一类是在有目标时判断为没有目标，这类错误称为漏警。
目标检测方法的核心是阈值法，其实说来简单——如果雷达回波大于阈值，则显示检测到目标，否则视为噪声。

在这里插入图片描述
恒虚警率（CFAR）算法的基本要素包括以下几个关键部分：

参考窗口（Reference Window）：CFAR算法需要一个参考窗口来估计背景噪声或杂波的统计特性。这个窗口通常选在认为没有目标出现的区域。
测试单元（Cell Under Test, CUT）：这是算法正在评估是否包含目标的单元。CFAR算法会将测试单元的响应与根据参考窗口估计的噪声水平进行比较。
噪声估计：在参考窗口内，算法会估计噪声的功率或强度。这可以通过平均、中值或其他统计方法来实现。
阈值计算：根据噪声估计，CFAR算法会计算一个阈值。如果测试单元的响应超过这个阈值，则可能检测到目标。
阈值偏移（Threshold Offset）：为了考虑信噪比和所需的检测概率，阈值可能会根据特定的系统性能要求进行调整。
比较和决策：测试单元的响应与阈值进行比较。如果响应高于阈值，则可能检测到目标；如果低于阈值，则可能是噪声或杂波。
保护区域（Guard Cells）：在测试单元周围可能会设置保护区域，以防止测试单元与参考窗口或相邻单元的干扰。
训练窗口（Training Window）：这是参考窗口内用于估计噪声的一部分，通常选在远离测试单元的位置。
边缘处理：在处理矩阵边缘附近的测试单元时，可能需要特别处理，因为这些位置可能没有足够的参考窗口。
输出结果：CFAR算法的输出是一个二值图像或矩阵，表示是否检测到目标（例如，1表示检测到目标，0表示未检测到）。
算法类型：有多种CFAR算法，如细胞平均CFAR（CA-CFAR）、最大值保持CFAR（MCFAR）、排序统计CFAR（OS-CFAR）等，每种算法有其特定的噪声估计和阈值计算方法。

MATLAB仿真——雷达目标检测

把这段代码贴上来，其过程和 CFAR 的基本要素很接近，让 GPT 帮忙读了一下：

初始化参数： n_train_cells 和 n_train_bands 确定了用于噪声估计的训练单元在距离和多普勒维度上的数量。n_guard_cells 和 n_guard_bands 确定了测试单元（CUT）周围的保护单元数量，以避免边缘效应。offset 是加在噪声估计上的一个偏移量，用于调整最终的检测阈值。noise_level 是一个向量，用于存储每次迭代中训练单元的噪声水平。
归一化RDM： RDM 被归一化，可能是为了在后续处理中避免数值过大或过小。
滑动窗口： 双重 for 循环在RDM上滑动CUT，同时在边缘留出训练和保护单元的余量。
噪声估计： 内层双重 for 循环遍历训练单元，如果单元不在CUT的保护区域内，则累加其噪声水平（从分贝转换为线性单位）。
计算阈值： 计算所有训练单元的平均噪声水平，然后转换回分贝单位，并加上偏移量以确定阈值。
CUT的检测： 对于每个CUT，将其信号水平与阈值进行比较。如果CUT的信号水平大于阈值，则标记为1（检测到目标），否则标记为0（未检测到目标）。
处理边缘情况： 由于CUT不能位于矩阵边缘，因此处理后的RDM会比原始的RDM小。将未被阈值化处理的单元格设置为0。
显示CFAR输出： 使用 surf 函数显示经过CA-CFAR滤波后的RDM，并添加了图形的标题、轴标签和视角设置。

%% CFAR implementation

%Slide Window through the complete Range Doppler Map

% *%TODO* :
%Select the number of Training Cells in both the dimensions.
n_train_cells = 10;
n_train_bands = 8;

% *%TODO* :
%Select the number of Guard Cells in both dimensions around the Cell under 
%test (CUT) for accurate estimation
n_guard_cells = 4;
n_guard_bands = 4;

% *%TODO* :
% offset the threshold by SNR value in dB
offset = 1.4;

% *%TODO* :
%Create a vector to store noise_level for each iteration on training cells
noise_level = zeros(1,1);


% *%TODO* :
%design a loop such that it slides the CUT across range doppler map by
%giving margins at the edges for Training and Guard Cells.
%For every iteration sum the signal level within all the training
%cells. To sum convert the value from logarithmic to linear using db2pow
%function. Average the summed values for all of the training%cells used. After averaging convert it back to logarithimic using pow2db.
%Further add the offset to it to determine the threshold. Next, compare the
%signal under CUT with this threshold. If the CUT level > threshold assign
%it a value of 1, else equate it to 0.


   % Use RDM[x,y] as the matrix from the output of 2D FFT for implementing
   % CFAR

RDM = RDM / max(RDM(:));


for row0 = n_train_cells + n_guard_cells + 1 : (Nr/2) - (n_train_cells + n_guard_cells)
  for col0 = n_train_bands + n_guard_bands + 1 : (Nd) - (n_train_bands + n_guard_bands)
    %Create a vector to store noise_level for each iteration on training cells
    noise_level = zeros(1, 1);

    for row1 = row0 - (n_train_cells + n_guard_cells) : row0 + (n_train_cells + n_guard_cells)
      for col1 = col0 - (n_train_bands + n_guard_bands) : col0 + (n_train_bands + n_guard_bands)
        if (abs(row0 - row1) > n_guard_cells || abs(col0 - col1) > n_guard_bands)
          noise_level = noise_level + db2pow(RDM(row1, col1));
        end
      end
    end

    % Calculate threshold from noise average then add the offset
    thresh = pow2db(noise_level / (2 * (n_train_bands + n_guard_bands + 1) * 2 * (n_train_cells + n_guard_cells + 1) - (n_guard_cells * n_guard_bands) - 1));
    thresh = thresh + offset;

    CUT = RDM(row1,col1);

    if (CUT < thresh)
      RDM(row0, col0) = 0;
    else
      RDM(row0, col0) = 1;
    end

  end
end


% *%TODO* :
% The process above will generate a thresholded block, which is smaller 
%than the Range Doppler Map as the CUT cannot be located at the edges of
%matrix. Hence,few cells will not be thresholded. To keep the map size same
% set those values to 0. 

RDM(RDM~=0 & RDM~=1) = 0;

% *%TODO* :
%display the CFAR output using the Surf function like we did for Range
%Doppler Response output.
figure('Name', 'CA-CFAR Filtered RDM')
surf(doppler_axis, range_axis, RDM);
%colorbar;
title( 'CA-CFAR Filtered RDM surface plot');
xlabel('Speed');
ylabel('Range');
zlabel('Normalized Amplitude');
view(315, 45);