分享一则代码，主要用于FMCW雷达室内多目标MATLAB仿真，涉及到的内容和算法模块有如下：
1、目标参数设置
2、雷达参数设置
3、目标运动状态设置
4、雷达信号建模（IQ信号）
5、雷达近场收发几何位置偏差校正
6、距离维FFT
7、速度维FFT
8、2D均值滤波
9、2D峰值搜索
10、直流分量去除
11、MATLAB图像保存为视频等内容。
另外，还有很多细节部分，这里暂时不展开了，具体看代码的细节，大部分都写了注释，花点功夫应该看的明白。
静态效果如下所示：
动态效果如下所示：
室内多目标场景
代码如下：
（1）主程序
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 功能：FMCW雷达室内多目标人员仿真 %
% 时间：2023-06-11 %
% 作者：调皮连续波 %
clc;clear;close all;
%% 设置目标个数
target_num = 5;
disp(['仿真目标个数: [',num2str(target_num),']']);
%% 生成仿真数据
disp('生成仿真数据...');
[RX1_CHIRPS_I,RX1_CHIRPS_Q,RX2_CHIRPS_I,RX2_CHIRPS_Q,TargetTracks]...
= radar_simulation_core_multi_target(target_num);%存在目标
disp('仿真数据生成成功！');
%% 雷达信号处理部分
disp('雷达信号处理...');
radar_data_processing_multi_target(RX1_CHIRPS_I,RX1_CHIRPS_Q,RX2_CHIRPS_I,RX2_CHIRPS_Q,TargetTracks,target_num);
disp('雷达信号处理完成！');
（2）
[RX1_CHIRPS_I,RX1_CHIRPS_Q,RX2_CHIRPS_I,RX2_CHIRPS_Q,TargetTracks]…
= radar_simulation_core_multi_target(target_num);%存在目标
function [RX1_CHIRPS_I,RX1_CHIRPS_Q,RX2_CHIRPS_I,RX2_CHIRPS_Q,TargetTracks]= radar_simulation_core_multi_target(target_num)
SimulationTime = 5; % 雷达系统仿真时间
SampleRate = 2e6; % Chirp ADC采样率
dt = 1/SampleRate; % Chirp ADC采样时间间隔
ChirpPeriod = 1e-3; % Chirp PRT
ChirpNum = SimulationTime/ChirpPeriod; % Chirp总数
ChirpDataLen = ChirpPeriod*SampleRate; % chirp ADC数据总长度
N = SimulationTime*SampleRate; % 仿真时间内的采样点数
n = 1:N; % 仿真时间内的采样序列点索引
t = (n-1)*dt; % 仿真时间内的采样时间索引
c = 3e8; % 光速
B = 4e9; % 雷达有效带宽
S = B/ChirpPeriod; % 调频斜率
Dres = c/2/B; % 理论距离分辨率（瑞利距离分辨率）
Lamda = 0.012; % 发射信号波长
Dx = 2*(Lamda/2); % 接收阵元间距：rx1和rx2
x0 = 0; % tx 坐标 x
y0 = 0; % tx 坐标 y
z0 = 0; % tx 坐标 z
x1 = -Dx/2; % rx1 坐标 x
y1 = 0; % rx1 坐标 y
z1 = 0; % rx1 坐标 z
x2 = Dx/2; % rx2 坐标 x
y2 = 0; % rx2 坐标 y
z2 = 0; % rx2 坐标 z
%此处坐标可用于Monopluse天线阵列调整阵元位置
%设置目标在xy上的分速度
Vx = zeros(target_num,N);
Vy = zeros(target_num,N);
for k = 1: target_num %多个目标的速度是变化的，不是恒定速度，通过作图看速度的曲线近似于正弦变化
Vx(k,:) = (0.6*rand()+0.3)*cos((3*rand()+1.5)*t + 20*pi*rand())+0.01*sin(30*t + 20*pi*rand());
Vy(k,:) = (0.8*rand()+0.2)*cos((3*rand()+1.5)*t + 20*pi*rand())+0.01*sin(30*t + 20*pi*rand());
end
% plot(Vx(1,:))
X = zeros(target_num,N);
Y = zeros(target_num,N);
Z = ones(target_num,N)*1.5;%目标初始位置Z，恒定不变
for m = 1: target_num %这里的距离对应后续距离维FFT的边界划分
X(m,1) = 4*(rand()-0.5); %目标初始位置X
Y(m,1) = 4*(rand()-0.5); %目标初始位置Y
for k = 2:N
X(m,k) = X(m,k-1) + Vx(m,k-1)*dt;%目标位置X变化
Y(m,k) = Y(m,k-1) + Vy(m,k-1)*dt;%目标位置Y变化
%限定XY边界
if X(m,k)>2
X(m,k) = 2;
elseif X(m,k)<-2
X(m,k) = -2;
end
if Y(m,k)>2
Y(m,k) = 2;
elseif Y(m,k)<-2
Y(m,k) = -2;
end
end
end
%近场校正雷达tx和rx的几何位置偏差，用于提高距离精度
d0 = ((X-x0).^2 + (Y-y0).^2 + (Z-z0).^2).^0.5;
d1 = ((X-x1).^2 + (Y-y1).^2 + (Z-z1).^2).^0.5;
d2 = ((X-x2).^2 + (Y-y2).^2 + (Z-z2).^2).^0.5;
%计算中频信号频率
if1 = S*(d0+d1)/c;
if2 = S*(d0+d2)/c;
%计算中频信号角频率
w1 = 2*pi*if1;
w2 = 2*pi*if2;
%计算接收阵列初始相位
rx1_wt_phi = zeros(target_num,N);
rx2_wt_phi = zeros(target_num,N);
for m = 1: target_num
for k = 1:N
if mod(k,ChirpDataLen)==1
rx1_wt_phi(m,k) = 0;
rx2_wt_phi(m,k) = 0;
else
rx1_wt_phi(m,k) = rx1_wt_phi(m,k-1)+w1(m,k-1)*dt;
rx2_wt_phi(m,k) = rx2_wt_phi(m,k-1)+w2(m,k-1)*dt;
end
end
end
%接收阵列目标回波信号相位
rx1_phi = rx1_wt_phi + (d0+d1)*2*pi/Lamda;
rx2_phi = rx2_wt_phi + (d0+d2)*2*pi/Lamda;
%接收回波信号
rx1 = exp(1j*rx1_phi)./d0.^2./d1.^2; %归一化
rx2 = exp(1j*rx2_phi)./d0.^2./d2.^2;
% plot(real(rx1(1,:)))
%多目标信号叠加，组合为一路回波信号
rx1_sum = zeros(1,N);
rx2_sum = zeros(1,N);
for k = 1: N
rx1_sum(k) = sum(rx1(:,k));
rx2_sum(k) = sum(rx2(:,k));
end
RX1_I = real(rx1_sum);
RX1_Q = imag(rx1_sum);
RX2_I = real(rx2_sum);
RX2_Q = imag(rx2_sum);
% plot(RX1_I)
% Chirp发射信号幅度调制
RX1_CHIRPS_I = zeros(ChirpNum,ChirpDataLen);
RX1_CHIRPS_Q = zeros(ChirpNum,ChirpDataLen);
RX2_CHIRPS_I = zeros(ChirpNum,ChirpDataLen);
RX2_CHIRPS_Q = zeros(ChirpNum,ChirpDataLen);
for k = 1:ChirpNum
chirp_start = (k-1)*ChirpDataLen + 1;
chirp_end = k*ChirpDataLen;
RX1_CHIRPS_I(k,:) = RX1_I(chirp_start:chirp_end);
RX1_CHIRPS_Q(k,:) = RX1_Q(chirp_start:chirp_end);
RX2_CHIRPS_I(k,:) = RX2_I(chirp_start:chirp_end);
RX2_CHIRPS_Q(k,:) = RX2_Q(chirp_start:chirp_end);
for m=1:ChirpDataLen
RX1_CHIRPS_I(k,m) = RX1_CHIRPS_I(k,m)*(m/ChirpDataLen)^0.33;
RX1_CHIRPS_Q(k,m) = RX1_CHIRPS_Q(k,m)*(m/ChirpDataLen)^0.33;
RX2_CHIRPS_I(k,m) = RX2_CHIRPS_I(k,m)*(m/ChirpDataLen)^0.33;
RX2_CHIRPS_Q(k,m) = RX2_CHIRPS_Q(k,m)*(m/ChirpDataLen)^0.33;
end
end
% plot(RX1_CHIRPS_I(1,:))
RX1_CHIRPS_I = RX1_CHIRPS_I(:,1:1360);%信号截取前面部分
RX1_CHIRPS_Q = RX1_CHIRPS_Q(:,1:1360);
RX2_CHIRPS_I = RX2_CHIRPS_I(:,1:1360);
RX2_CHIRPS_Q = RX2_CHIRPS_Q(:,1:1360);
% 添加随机噪声
NoiseAmp = 0.001 ;
RX1_CHIRPS_I = RX1_CHIRPS_I + rand(size(RX1_CHIRPS_I))*NoiseAmp;
RX1_CHIRPS_Q = RX1_CHIRPS_Q + rand(size(RX1_CHIRPS_Q))*NoiseAmp;
RX2_CHIRPS_I = RX2_CHIRPS_I + rand(size(RX2_CHIRPS_I))*NoiseAmp;
RX2_CHIRPS_Q = RX2_CHIRPS_Q + rand(size(RX2_CHIRPS_Q))*NoiseAmp;
TargetTracks.X = X;
TargetTracks.Y = Y;
TargetTracks.Z = Z;
end
（3）radar_data_processing_multi_target(RX1_CHIRPS_I,RX1_CHIRPS_Q,RX2_CHIRPS_I,RX2_CHIRPS_Q,TargetTracks,target_num);
function radar_data_processing_multi_target(I1, Q1,I2,Q2,TargetTracks,target_num)
RANGE_MIN = 16;
RANGE_MAX = 64;
DPL_FFT_POINTS = 32;
DPL_FFT_INTVERAL = 32;
DPL_FFT_ELSIZE = 128;
[TotalChirpNum,ChirpLen] = size(I1);
%信号实部虚部组合
rx1_c = I1 + Q1*1j;
rx2_c = I2 + Q2*1j;
% Range FFT
if ~exist('RANGE_MIN','var')
RANGE_MIN = 0;
end
if ~exist('RANGE_MAX','var')
RANGE_MAX = 16;
end
rx1_range_fft = zeros(TotalChirpNum,1024);
rx2_range_fft = zeros(TotalChirpNum,1024);
tmp = zeros(1,ChirpLen);
for k=1:TotalChirpNum
for m=1:ChirpLen
tmp(m)=rx1_c(k,m)/((m+1)/ChirpLen)^0.33; %去掉幅度调制 RX1_CHIRPS_I(k,m) = RX1_CHIRPS_I(k,m)*(m/ChirpDataLen)^0.33;
end
rx1_range_fft(k,:) = fft((tmp(301:1324)-mean(tmp(301:1324))).*hamming(1024)');
for m=1:ChirpLen
tmp(m)=rx2_c(k,m)/((m+1)/ChirpLen)^0.33; %去掉幅度调制 RX2_CHIRPS_I(k,m) = RX2_CHIRPS_I(k,m)*(m/ChirpDataLen)^0.33;
end
rx2_range_fft(k,:) = fft((tmp(301:1324)-mean(tmp(301:1324))).*hamming(1024)'); %去直流后FFT
end
rx1_range_fft(:,1:RANGE_MIN) = 0;
rx2_range_fft(:,1:RANGE_MIN) = 0;
% Dopplor FFT
if ~exist('DPL_FFT_POINTS','var')
DPL_FFT_POINTS = 32;
end
if ~exist('DPL_FFT_INTVERAL','var')
DPL_FFT_INTVERAL = 32;
end
if ~exist('DPL_FFT_ELSIZE','var')
DPL_FFT_ELSIZE = 128;
end
rx1_dpl_fft = zeros(DPL_FFT_ELSIZE,RANGE_MAX);
rx2_dpl_fft = zeros(DPL_FFT_ELSIZE,RANGE_MAX);
%保存视频
SAVE_VIDEO = 1;
if SAVE_VIDEO==1
video_file_name = '室内多目标场景.avi';
video_obj = VideoWriter(video_file_name);
video_obj.FrameRate = fix(0.4*1000/DPL_FFT_INTVERAL);
open(video_obj)
end
%距离维数据提取，然后TotalChirpNum循环
for chirp_idx=1:DPL_FFT_INTVERAL:TotalChirpNum-DPL_FFT_POINTS-DPL_FFT_INTVERAL
rx1_rfft_clip = rx1_range_fft(chirp_idx:chirp_idx+DPL_FFT_POINTS-1,1:RANGE_MAX);
rx2_rfft_clip = rx2_range_fft(chirp_idx:chirp_idx+DPL_FFT_POINTS-1,1:RANGE_MAX);
rx1_rfft_clip(DPL_FFT_ELSIZE,:) = 0;
rx2_rfft_clip(DPL_FFT_ELSIZE,:) = 0;
for k=1:RANGE_MAX
rx1_dpl_fft(:,k) = fftshift(fft(rx1_rfft_clip(:,k)));
rx2_dpl_fft(:,k) = fftshift(fft(rx2_rfft_clip(:,k)));
end
eng_rx1_df = abs(rx1_dpl_fft);
eng_rx2_df = abs(rx2_dpl_fft);
eng_rx_df = eng_rx1_df.*eng_rx2_df; %两天线数据相与，得到的谱图代表真实目标
% mesh(abs(eng_rx_df));
eng_rx_df = avg_filter_2D(eng_rx_df,1);%2维均值滤波
for k=1:RANGE_MAX %幅度缩放
eng_rx_df(:,k)=eng_rx_df(:,k).*k^9;
end
for k=2:RANGE_MAX %距离维去直流
eng_rx_df(:,k-1)=eng_rx_df(:,k)-eng_rx_df(:,k-1);
end
eng_rx_df=-eng_rx_df(:,1:RANGE_MAX-1);
eng_rx_df(eng_rx_df<(1e3+max(max(eng_rx_df))*1e-1))=0;
for k=2:DPL_FFT_ELSIZE %速度维去直流
eng_rx_df(k-1,:)=eng_rx_df(k,:)-eng_rx_df(k-1,:);
end
% mesh(abs(eng_rx_df));
eng_rx_df=-eng_rx_df(1:DPL_FFT_ELSIZE-1,:);
eng_rx_df(eng_rx_df<(1e3+max(max(eng_rx_df))*1e-1))=0;
%2D峰值搜索
eng_rx_pk = find_peaks_2D(eng_rx_df);
eng_rx_pk = avg_filter_2D(eng_rx_pk,1);
eng_rx_pk = erode_2D(eng_rx_pk,1); %细化
fh = figure(1);
imagesc(-31:31,1:127,eng_rx_df);
axis xy
xlabel('多普勒索引');ylabel('距离索引');zlabel('幅度');
title('距离速度谱');
figure(2)
for k = 1:target_num
plot(TargetTracks.X(k,chirp_idx*2000),TargetTracks.Y(k,chirp_idx*2000),'bo');
xlim([-5 5]);ylim([-5 5]);
grid on
hold on
end
xlabel('X坐标(m)');ylabel('Y坐标(m)')
title(['人员统计 : ','\fontsize{12}\color{blue}',num2str(target_num)]);
hold off
if SAVE_VIDEO==1
frame = getframe(fh);
writeVideo(video_obj,frame);
end
end
if SAVE_VIDEO==1
close(video_obj);
end
end
（4）补充代码：
mat_ret = avg_filter_2D(matrix_2d,marg)
function mat_ret = avg_filter_2D(matrix_2d,marg)
[M,N] = size(matrix_2d);
matrix_fil = zeros([M+2*marg,N+2*marg]);
matrix_fil(1+marg:M+marg,1+marg:N+marg) = matrix_2d; %矩阵扩充，填零
mat_ret = matrix_fil;
%求原来矩阵的均值
for m=1+marg:M+marg
for n=1+marg:N+marg
mat_ret(m,n) = mean(matrix_fil(m,n-marg:n+marg));
end
end
matrix_fil = mat_ret;
for m=1+marg:M+marg
for n=1+marg:N+marg
mat_ret(m,n) = mean(matrix_fil(m-marg:m+marg,n));
end
end
mat_ret = mat_ret(1+marg:M+marg,1+marg:N+marg);
end
img_er = erode_2D(img_in,ep)
function img_er = erode_2D(img_in,ep)
[row,col]=size(img_in);
img_er = img_in;
%边界置零
img_er(1:2,:) = 0;
img_er(row-1:row,:) = 0;
img_er(:,1:2) = 0;
img_er(:,col-1:col) = 0;
% erode
for k=1:ep
img_tmp = img_er;
for i=3:row-2
for j=3:col-2
if img_tmp(i-1,j) && img_tmp(i+1,j) && img_tmp(i,j-1) && img_tmp(i,j+1)
% img_er(i,j) = 1;
else
img_er(i,j) = 0;
end
end
end
end
end
[img_b,num] = erode_tobin_2D(img_in)
function [img_b,num] = erode_tobin_2D(img_in)
img_b = img_in;
img_b(img_in>0) = 1;
img_b = bwmorph(img_b,'shrink',inf); %二值化
num = sum(sum(img_b));
end
find_peaks_2D(img_in)
function img_pk = find_peaks_2D(img_in)
[row,col]=size(img_in);
img_pk = img_in;
img_pk(1:2,:) = 0;
img_pk(row-1:row,:) = 0;
img_pk(:,1:2) = 0;
img_pk(:,col-1:col) = 0;
% erode
for i=3:row-2
for j=3:col-2
if img_pk(i,j)<img_pk(i-1,j) || ...
img_pk(i,j)<img_pk(i+1,j) || ...
img_pk(i,j)<img_pk(i,j-1) || ...
img_pk(i,j)<img_pk(i,j+1) || ...
img_pk(i,j)<img_pk(i-1,j-1) || ...
img_pk(i,j)<img_pk(i-1,j+1) || ...
img_pk(i,j)<img_pk(i+1,j-1) || ...
img_pk(i,j)<img_pk(i+1,j+1)
img_pk(i,j) = 0;
end
end
end
好了，本期内容就结束了，上述代码可以直接运行得到动态图，并保存视频。本代码运行环境为MATLAB2022a，如果出现乱码，可以先用记事本打开，然后复制到MATLAB中。这是由于高版本不兼容低版本导致 。
