FMCW 雷达室内多目标人员MATLAB仿真

分享一则代码，主要用于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中。这是由于高版本不兼容低版本导致 。文章来源地址https://www.toymoban.com/news/detail-481592.html

到了这里，关于FMCW 雷达室内多目标人员MATLAB仿真的文章就介绍完了。如果您还想了解更多内容，请在右上角搜索TOY模板网以前的文章或继续浏览下面的相关文章，希望大家以后多多支持TOY模板网！