generateParameter.m%% 雷达参数设置functionparametergenerateParameter(targetnum)%%% 雷达仿真参数设置点目标设置 % targetnum 1-生成点目标2-生成飞机目标 %%%定义物理常数与雷达硬件指标parameter.kB1.38e-23;% 玻尔兹曼常数 parameter.T290;% 环境温度(K)parameter.Pt0.0158;% 发射功率12dBm parameter.G_ant10^(15/10);% 天线增益15dBi parameter.Loss10^(4/10);% 系统损耗4dB parameter.NoiseFigure10^(10/10);% 接收机噪声系数10dB parameter.c3e8;%光速 parameter.stratFreq77e9;%起始频率iftargetnum1%1-点目标 parameter.Tr40e-6;% 10e-6;60e-6;%【扫频时间】else% 飞机目标 parameter.Tr60e-6;40e-6;end parameter.Idle_time10e-6;% 100e-6;%空闲时间 parameter.Tcparameter.Trparameter.Idle_time;%Chirp之间的间隔 parameter.Fs25.6e6;%10e6;%【采样率】 parameter.Samplesparameter.Tr * parameter.Fs;% 采样点信号生成单个chirp采样点数扫频时间*采样率 parameter.rangeBinparameter.Samples;%2^nextpow2(parameter.Samples);%rangebin 【为了计算效率通常将 rangeBin 设置为大于 Samples 的最小的2的幂例如采样200点FFT 用256点】 parameter.Chirps512;%128;%【】chirp数在一帧数据内雷达连续发射的chirp总数 parameter.dopplerBinparameter.Chirps;%dopplerbin parameter.Slope29.982e12;%chirp斜率 parameter.Bandwidthparameter.Slope * parameter.Tr;%发射信号带宽 parameter.BandwidthValidparameter.Tr * parameter.Slope;%发射信号有效带宽 parameter.centerFreqparameter.stratFreq parameter.Bandwidth /2;%中心频率 parameter.lambdaparameter.c / parameter.centerFreq;%波长中心频率波长 parameter.numFrames50;%50;%10;% 帧数 parameter.T_frame80e-3;% 帧周期 parameter.txAntennaones(1,3);%发射天线个数 parameter.rxAntennaones(1,4);%接收天线个数 parameter.txNumlength(parameter.txAntenna);parameter.rxNumlength(parameter.rxAntenna);parameter.virtualAntennalength(parameter.txAntenna)* length(parameter.rxAntenna);parameter.angleBin180;parameter.dzparameter.lambda /2;%接收天线俯仰间距 parameter.dxparameter.lambda /2;%接收天线水平间距 parameter.doaMethod3;% 测角方法选择1:fft2:dbf3:1Dmusic4:RDMusic5:Capon parameter.clusterMethod1;% 点云聚类方法选择1-DBSCAN2-Kmeans3-X-Means %定义目标参数【注意给不同的目标添加唯一的targetID】 % 设置目标速度时需要根据雷达理论可测的最大不模糊速度保证设定速度范围在这之内iftargetnum1%1-点目标 raw_points[100-60;%target1 range speed horizontalangle203-25;%target2 range speed horizontalangle20630;%target2 range speed horizontalangle];num_ptssize(raw_points,1);% 为点目标赋予1,2,3... 的独立 Target_ID 编号 parameter.target[raw_points,(1:num_pts)]; else % 飞机目标 % % 构建飞机目标的参数配置矩阵 % 每一行代表一架飞机目标: [中心距离(m), 中心角度(°), 总速度(m/s), 侧滑角(°)] % 中心角度方位角是以 Y 轴为 0° 基准。向右偏为正x 变大向左偏为负。 % airplanes_config [ % [距离, 角度, 速度, 侧滑角] 60, 10, -8, 225; % 飞机 1: 中距离迎面快速斜切 100, -20, -4, 110; % 飞机 2: 远距离左侧斜穿 40, 35, 7, 75; % 飞机 3: 近距离大角度脱离 ]; num_airplanes size(airplanes_config, 1); all_airplane_targets []; % 用于拼接所有飞机散射点的终极矩阵 % 循环生成每架飞机的点阵并融合成一个统一的雷达观测矩阵 for airplane_idx 1:num_airplanes c_range airplanes_config(airplane_idx, 1); c_angle airplanes_config(airplane_idx, 2); v_total airplanes_config(airplane_idx, 3); h_error airplanes_config(airplane_idx, 4); % 调用你原有的单飞机点阵生成函数 single_airplane generate_airplane_target(c_range, c_angle, v_total, h_error); num_scatterers size(single_airplane, 1); id_column ones(num_scatterers, 1) * airplane_idx; single_airplane_with_id [single_airplane, id_column]; % 垂直拼接把所有飞机的散射点汇聚到一起 all_airplane_targets [all_airplane_targets; single_airplane_with_id]; end % 将拼接后的数百个散射点赋给统一的 target 变量无缝对接到你的雷达回波生成器 parameter.target all_airplane_targets; % % 可视化飞机初始位置可以解开此段注释 % current_targets parameter.target; % R current_targets(:, 1); V_radial current_targets(:, 2); Th current_targets(:, 3); % X_radar_mesh R .* cosd(Th); Y_radar_mesh R .* sind(Th); % figure(1); scatter(Y_radar_mesh, X_radar_mesh, 30, V_radial, filled, MarkerEdgeColor, k); % grid on; box on; colormap(jet); colorbar; axis equal; hold on; % plot(0, 0, rp, MarkerSize, 12, MarkerFaceColor, r);end endSamples单个chirp的采样点数 扫频时间Tr*采样率Fs。为了计算效率通常将 rangeBin 设置为大于 Samples 的最小的 2 的幂例如采样 200 点FFT 用 256 点。发射信号的有效带宽。单chirp采样点数/采样率*chirp斜率。可改进点固定发射带宽求斜率保证设置的点目标及飞机目标对比时距离分辨率是相同的。波长用中心频率计算的起始频率带宽/2不是有效带宽。我这里的带宽和有效带宽值是相同的因为没考虑到射频稳定过渡期工业中ADC并非瞬间自启动需要改进的话可将Tr减去启动延时得到有效采样时间。或者换一种看法我的稳定过渡期涵盖在了idel空闲时间范围内。为了卡尔曼滤波帧数和帧周期的设定。帧周期应远大于一帧内所有发射chirp信号的时间和剩余的帧间间歇期用于信号算法处理等。设置目标坐标和速度时注意速度需要保证处于根据设定雷达参数计算到的最大不模糊速度范围内。【或者后期算法速度解模糊】设置仿真目标时注意添加目标的TargetID多目标区分。functiontarget_matrixgenerate_airplane_target(center_range, center_angle, v_total, heading_error)%% heading_error(航向偏差角单位度)% 当 heading_error0时飞机笔直撞向雷达 % 当 heading_error5或10时【飞机产生斜切角机身点云就会显现】 %X_local[];Y_local[];% ---1. 机身边缘纺锤体机身 --- x_fuse_steps-5:0.5:5;% 飞机全长10mforxx_fuse_steps r_fuse0.4* cos(x/6);% 约束机身宽度头尾窄中间宽 X_local[X_local, x, x];Y_local[Y_local, r_fuse, -r_fuse];% 机身对称最宽时为0.4*20.8m两侧边缘最窄时0.27*20.54m % figure(11);scatter(x, r_fuse,30,r,filled);xlim([-6,6]);ylim([-4.5,4.5]);% hold on;scatter(x, -r_fuse,30,r,filled);xlim([-6,6]);ylim([-4.5,4.5]);set(gcf,Color,white);box on;end% ---2. 后掠机翼轮廓线 --- y_wing_steps[-4:0.4:-0.4,0.4:0.4:4];% 机翼展宽8m跳过(-0.4,0.4)这个区间让机翼点云从机身外壳两侧向外生长foryy_wing_steps x_front-abs(y)*0.5 0.5;% 机翼前边缘向两侧时|y|增大x坐标向后减小x轴正方向朝向机头 wing_chord1.5*(1- abs(y)/5);% 机翼宽度/弦长。机翼靠近机身的地方很宽越往翼尖走越窄 x_backx_front - wing_chord;% 机翼的后边缘坐标用前缘坐标减去当前位置的机翼宽度 X_local[X_local, x_front, x_back];Y_local[Y_local, y, y];hold on;scatter(x_front, y,30,b,filled);xlim([-6,6]);ylim([-4.5,4.5]);hold on;scatter(x_back, y,30,b,filled);xlim([-6,6]);ylim([-4.5,4.5]);end% ---3. 箭形尾翼轮廓线 --- y_tail_steps[-1.5:0.4:-0.4,0.4:0.4:1.5];% 尾翼总宽度为3m同样跳过中间区间尾翼点云接在机身尾部两侧foryy_tail_steps x_front_tail-4.5- abs(y)*0.3;% 后掠前缘线尾翼是从距离机尾最后一米-4.5 米处开始向外生长的向两侧时|y|增大x坐标向后减小 tail_chord0.8*(1- abs(y)/2);% 尾翼宽度根部最宽尾翼尖端最小 x_back_tailx_front_tail - tail_chord;% 后边缘坐标用前缘坐标减去当前位置的机翼宽度 X_local[X_local, x_front_tail, x_back_tail];Y_local[Y_local, y, y];% hold on;scatter(x_front_tail, y,30,g,filled);xlim([-6,6]);ylim([-4.5,4.5]);% hold on;scatter(x_back_tail, y,30,g,filled);xlim([-6,6]);ylim([-4.5,4.5]);end% ---4. 封闭翼尖与尾部边界 --- X_local[X_local, -1.5, -1.5];Y_local[Y_local,4, -4];X_local[X_local, -5.5];Y_local[Y_local,0];hold on;scatter(-1.5,4,30,c,filled);xlim([-6,6]);ylim([-4.5,4.5]);hold on;scatter(-1.5, -4,30,c,filled);xlim([-6,6]);ylim([-4.5,4.5]);hold on;scatter(-5.5,0,30,c,filled);xlim([-6,6]);ylim([-4.5,4.5]);%% ---5. 姿态旋转矩阵(引入了 heading_error)--- % 当 heading_error0时飞机就是绝对笔直地、自杀式地对着雷达天线冲过来。此时雷达只能看到一个正前方的“机头剪影”机翼面积被严重压缩。 % “机身点云就会显现”的原因只要 heading_error 不为0飞机就是横着翅膀、侧着身子在做斜切进场。从雷达的视线LOS看过去 % 原本被隐藏的侧面机身、巨大的后掠翼表面全都在雷达视野里暴露无遗散射点云瞬间变得极其丰富 % 根据旋转矩阵顺时针旋转角度为负逆时针为正的原理设置旋转角 rot_angle90- heading_error;R_matrix[cosd(rot_angle), -sind(rot_angle);% 二维逆时针旋转矩阵 sind(rot_angle), cosd(rot_angle)];% 把旋转后的飞机整体搬运到指定的空间位置 % 角度定义center_angle方位角是以 Y 轴为0° 基准。向右偏为正x 变大向左偏为负。 rotated_pointsR_matrix *[X_local;Y_local];% 将飞机上的所有散射点绕着自身中心旋转rot_angle角度 % figure(22);scatter(rotated_points(1,:), rotated_points(2,:),30,k,filled);%xlim([-8,8]);ylim([-8,8]);% hold on;plot(0,0,rp,MarkerSize,12,MarkerFaceColor,r);set(gcf,Color,white);box on;axis equal;X_radarrotated_points(1,:) center_range * sind(center_angle);Y_radarrotated_points(2,:) center_range * cosd(center_angle);% figure(33);scatter(X_radar, Y_radar,30,k,filled);%xlim([-8,8]);ylim([-8,8]);% hold on;plot(0,0,rp,MarkerSize,12,MarkerFaceColor,r);set(gcf,Color,white);box on;axis equal;下图是绕着中心点旋转-30°的效果设置的三个飞机目标绘制出来的效果如下对应上距离角度速度及航向%% ---6. 解算回雷达参数(计算由于斜切导致的径向速度改变)--- num_scattererslength(X_radar);target_matrixzeros(num_scatterers,3);fori1:num_scatterers r_isqrt(X_radar(i)^2 Y_radar(i)^2);% 每一个具体点如左翼尖、右翼尖、机尾【相对于雷达原点】的即时视线夹角 % 因为飞机有10米长、8 米宽在近距离下左翼尖和右翼尖的 angle_i 甚至能差出好几度 angle_iatand(X_radar(i)/ Y_radar(i));% 【多普勒径向速度投影】。核心物理修正因为飞机斜着飞各个散射点到雷达的投影速度不同 v_radial_pureabs(v_total)* cosd(angle_i - heading_error);%v_total * cosd(angle_i - flight_dir);ifv_total0% 如果顶层设置 v_total 为负靠近趋势确保投影出的速度整体呈现负值 v_radial-abs(v_radial_pure);else% 如果顶层设置 v_total 为正远离趋势确保呈现正值 v_radialabs(v_radial_pure);end target_matrix(i,:)[r_i, v_radial, angle_i];end end设计飞机点目标时飞机机头朝向与航向角Y轴为正方向一致。由于旋转矩阵是以自身为中心顺时针旋转为负逆时针为正。所以航向角应当取负号而初始飞机目标机头朝向为x轴正向所以还需要再加上90°。计算飞机上所有点目标的径向速度时需要同时考虑飞机与雷达的中心角度与航向角两者差用来计算径向速度。同时为了保持与设定的概念速度为负表示靠近雷达给计算出的径向速度符号做出相应调整。GitHub仓库https://github.com/ReddingZHL/4D-mmWave-Radar-Pipeline