
简介面向雷达、声纳及无线通信研究者的FDA波束形成MATLAB程序包围绕多载频频率分集与自适应波束形成展开适合学习最小方差无失真响应、线性约束最小方差等经典算法以及多波束联合处理技术。压缩包共15个文件、约1.95MB其中9个.m脚本覆盖收发方向图绘制、自适应波束形成实现、信干噪比对比等核心环节3个.fig结果图便于直观核验输出3个.ini配置文件用于环境参数初始化。程序内包含类似lcmv自适应、SINR比较等功能模块可直接运行并修改频率间距、载频数等参数来观察波束指向和干扰抑制效果同时结合最小范数方向估计与多波束切换可辅助理解工程中动目标跟踪与多目标检测的算法细节也能用于课程设计或项目初期的算法验证。已有714人学习下载适合作为信号处理课程设计或雷达通信项目的基础参考。1. FDA波束形成程序是什么为什么它值得写在雷达仿真第一行如果你写过程序里只有角度权的相控阵波束形成第一次看到“FDA波束形成程序”这个标题多半会问频率分集阵列和相控阵到底差在哪一句话回答相控阵方向图只跟角度有关而FDA因为每个阵元发射频率有微小偏移方向图会随距离变化主瓣在距离-角度平面走出斜线。这个特性让FDA在近程探测、距离模糊抑制、多目标分辨上被重视代价是程序比相控阵多出一个距离维参数一错方向图就翻车。这篇文章写给要自己从零写FDA波束形成程序的雷达、通信工程师讲信号模型、MATLAB实现、参数边界和踩坑记录读完你能落一个能扫描距离-角度二维方向图的程序并知道每个参数改下去会发生什么。2. 从相控阵到FDA的信号模型频率增量如何让距离进入波束方向图2.1 相控阵的相位项与FDA的差异相控阵波束形成程序里每个阵元只加载角度相位补偿导向矢量写作a_phased(θ) exp(-j * 2π * f0 * n * d * sinθ / c)其中 n 是阵元序号0 到 N-1d 是阵元间距f0 是载频。这个矢量只和 θ 有关所以不管目标距离多少主瓣指向不变。FDA 则让第 n 个阵元的发射频率变成 f_n f0 n * Δf。频率不同信号到达或返回各阵元时叠加的相位差就不再只有角度项而是多了一项与 Δf、阵元序号和距离 r 相关的相位。也就是说导向矢量从一维变成二维a(r, θ)。工程上最容易卡住的是“FDA 方向图还随时间变化”。如果发射连续波波前相位除了载频项和 Δf 项还会带上时间 t脉冲体制下还会出现慢时间项。完全按时间表达式写程序里要多一个快拍维调试难度直线上升。为了让程序先立住业界通行做法是假设窄带、远场、Δf 远小于 f0只保留对波束形状影响最大的两项常规角度相位 f0 * n * d * sinθ / c 和距离相关相位 Δf * n * r / c。在慢时间维做平均或取某个参考时刻时变项可以用一组固定相位代替。这个假设不是偷懒而是 FDA 方向图在单个脉冲内的变化通常远小于波束宽度等你把单目标验证过了再逐步把慢时间项加回来。还有一个容易忽略的点FDA 的“距离依赖”来自 Δf * r 这一项而不是载频项的 r/c。载频项在远场近似中对所有阵元几乎相同不进入阵元间的相对相位真正让各阵元相位随距离产生差异的是频率偏移与传播时延的乘积。所以你在程序里写距离项时必须保证它乘的是 Δf而不是 f0。很多初学者在这里把公式抄串行导致方向图在距离维毫无反应。2.2 FDA波束形成程序要算的导向矢量长什么样写程序不是背公式而是把公式翻译成数组。FDA 波束形成程序最核心的量是一个关于距离和角度的导向矢量a(r, θ) exp( -j * 2π / c * ( f0 * d * sinθ Δf * r ) * n )其中 n 是列向量 (0:N-1)^T。这个式子把全部阵元的相位差写成了“一个标量因子乘阵元序号”。标量因子由角度和距离共同决定因此阵列流形矩阵 A 不再是一维阵列而是由所有扫描点和阵元序号组成的三维数组每个 (r, θ) 切片对应一个 N 维复数矢量。写程序时先把这个“标量因子”算对再广播到阵元序号上。线性模型里没有出现 n² 项这是大家最常见的简化。严格推导 FDA 发射导向矢量时阵元 n 的相位差里会出现 Δf * n² * d * sinθ / c 这样的二阶项因为第 n 个阵元相对参考阵元的频率偏移同时改变了发射频率而发射位置又差了一个 ndsinθ。在阵元数不大、Δf 在 kHz 量级时这个二阶项带来的相位变化只有几度对主瓣位置影响很小。但当你把阵元数加到 64 以上、Δf 加到 100kHz二阶项会让方向图旁瓣结构明显改变。我一般建议程序里先写一个可选开关use_second_order false; if use_second_order phase_extra -2*pi/c * delta_f * n.^2 * d * sind(Theta); % 注意维度需广播 else phase_extra 0; end这个开关让你在同一个代码框架下对比线性模型和严格模型等需要精算旁瓣时再打开它。两种模型的适用边界也值得按表格记一下模型相位项适用场景程序风险线性模型只保留 f0dsinθ 与 Δf*rN≤32Δf≤10kHz单目标验证旁瓣细节失真严格模型增加 Δfn²d*sinθ 等二阶项N32Δf50kHz旁瓣精算方向图不对称主瓣位置微偏初学阶段建议先用线性模型跑通主瓣轨迹再打开二阶项看差异。不要一上来就上严格模型否则方向图乱了你都不知道该查哪一项。2.3 程序里怎么生成距离-角度网格与方向图矩阵FDA 程序里最先要建立的不是循环而是距离-角度网格。角度范围按雷达覆盖定通常扫 ±90°距离范围按威力定例如近程雷达扫 5km 到 20km。网格宽度取决于你想要的分辨率程序初版不必太密够看出主瓣就行。生成网格和导向矢量的常见做法是theta_scan linspace(-90, 90, 361); r_scan linspace(5e3, 20e3, 301); [Theta, R] ndgrid(theta_scan, r_scan); n (0:N-1).; % 标量相位因子 phase_scalar -2*pi/c * (f0 * d * sind(Theta) delta_f * R); % 广播到每个阵元结果尺寸 [角度点数, 距离点数, N] phase phase_scalar .* reshape(n, 1, 1, N); a exp(1j * phase);为什么用 ndgrid 而不是 meshgrid雷达方向图习惯上把角度放行方向、距离放列方向ndgrid 的行、列顺序和习惯一致后面画图和找峰值都省心。reshape(n,1,1,N) 把阵元序号变成第三维phase_scalar 自动广播成三维这是 MATLAB 里沿阵元维向量化的标准姿势。这样生成的 a 是 361×301×16 的复数数组内存约 361×301×16×16 字节≈27MB完全可控。如果你把网格加密到 2000×2000这个数组会到 1GB 级别届时要分块生成后面避坑章会专门说。这段代码里最容易写错的是phase_scalar .* reshape(n,1,1,N)。如果直接用. * nMATLAB 会尝试把 361×301 矩阵和 16×1 向量广播结果变成 361×301×16 可能可以但 n 的方向容易反。写成 reshape(n,1,1,N) 后语义最清楚第三维是阵元序号。我早期在这里翻过车方向图看起来有主瓣但沿阵元维的相位顺序是反的导致波束指向镜像到另一侧。这种错误用第6章的峰值轨迹验证法能一眼看出来。另外要提醒网格步长角度网格步长建议不大于主瓣宽度的 1/5距离网格步长不大于 c/(Δf*N) 的一半。否则峰值位置可能正好落在网格缝隙里你提取的峰值距离和理论值总差几格误以为自己程序写错了。3. 用MATLAB写最小FDA多波束程序从单波束到权矩阵3.1 最小可运行程序参数表与方向图生成把第2章的模型落成一段可运行代码。参数按常见做法载频 10GHz、阵元数 16、阵元间距半波长、频率增量 10kHz、设计距离 10km、设计角度 -10°。下面是完整代码你直接复制就能跑。% min_fda_direction_plot.m % 最小FDA波束形成方向图生成程序 f0 10e9; % 载频 10GHz c 3e8; lambda c / f0; N 16; % 阵元数 d lambda / 2; % 半波长阵元间距 delta_f 10e3; % 频率增量 10kHz r0 10e3; % 设计点距离 10km theta0 -10; % 设计点角度 -10° theta_scan linspace(-90, 90, 361); r_scan linspace(5e3, 20e3, 301); [Theta, R] ndgrid(theta_scan, r_scan); n (0:N-1).; % 权矢量设计点处的导向矢量 phase_w -2*pi/c * (f0 * d * sind(theta0) delta_f * r0) * n; w exp(1j * phase_w); % 扫描区域的导向矢量三维数组 phase_a -2*pi/c * (f0 * d * sind(Theta) delta_f * R) .* reshape(n, 1, 1, N); a exp(1j * phase_a); % 波束形成权共轭点乘导向矢量沿阵元维求和 P abs(sum(conj(w) .* a, 3)).^2;运行后 P 是 361×301 矩阵。用下面这段把方向图画出来imagesc(theta_scan, r_scan/1e3, P.); xlabel(角度 (deg)); ylabel(距离 (km)); axis xy; colorbar;注意 imagesc 需要把 P 转置因为 imagesc 的第一个输入对应 x 轴第二个对应 y 轴矩阵的列才会映射到颜色。画出来的方向图会在 -10°、10km 附近出现一个主瓣并且主瓣沿距离-角度方向拉出一条斜脊这就是 FDA 区别于相控阵的视觉特征。如果你画出来是一条水平横条说明 delta_f 项没有生效回去检查 phase_a 里是否真的用了 delta_f * R。3.2 关键参数为什么这样设这段程序里的参数不是随便拍的。delta_f 10kHz 在 10GHz 载频下是载频的百万分之一属于窄带 FDA 常用区间。它决定了距离维波束的斜率delta_f 越大角度-距离脊线越陡距离可辨性越好但距离模糊周期 c/delta_f 越小10kHz 对应 30km我们的扫描范围 5km 到 20km 不会跨周期所有距离峰都是单值的。N16 是入门级线阵角度波束宽度约 6°够看出两个角度目标分离计算量也小。r0 和 theta0 是你要瞄准的空域位置多波束场景下就是多个这样的点。一个细节这里权矢量没有加窗。矩形权会让旁瓣在 -13dB 左右FDA 的距离维旁瓣会更高一些因为频率增量本身会带来“栅瓣式”的距离副峰。如果要压低旁瓣可以对权矢量加海明窗但要付出主瓣展宽的代价。程序里加窗很简单win hamming(N); w w .* win;加窗后主瓣会变宽峰值位置会稍偏所以一般先不加窗调试调准主瓣轨迹后再加。3.3 从单波束到多波束权矩阵一次算K个波束多波束实现最省事的方式是构造 K 个权矢量每个权矢量指向各自的 (r_k, θ_k)然后把 K 个权堆成矩阵 W。下面以 3 个波束为例K 3; r_des [8e3, 1e4, 1.2e4]; % 三个设计距离 theta_des [-30, -10, 10]; % 三个设计角度 % 生成 KxN 权矩阵每行是一个匹配权矢量 phase_W -2*pi/c * (f0 * d * sind(theta_des(:)) delta_f * r_des(:)) * n.; W exp(1j * phase_W); % 把扫描导向矢量 a 从 361x301xN 转成 N x (361*301) A reshape(a, N, []); P_multi abs(conj(W) * A).^2; % K x (361*301) P_multi reshape(P_multi, K, size(Theta,1), size(Theta,2));逻辑说明theta_des(:) 强制列向量(f0*d*sind(theta_des(:)) delta_f*r_des(:))是 K×1乘以 n. 后得到 K×N 的相位矩阵W 的每一行对应一个波束的权。A 把三维导向矢量展成二维行是阵元序号列是扫描点conj(W) * A 同时完成 K 个匹配滤波。P_multi 的每个切片是一个波束的方向图用squeeze(P_multi(k,:,:))就能单独画。这种写法的优点是增加波束数只动 WA 只算一次。如果目标数量变化频繁把 W 换成多目标导向矩阵即可后端检测可以直接复用这个矩阵乘积的输出。另外如果只想同时形成多个距离波束但角度相同可以把 r_des 设成多个距离、theta_des 设为同一个角度如果只想多角度反过来即可。注意波束设计点不要离得太近否则会出现第5章说的互相串扰。3.4 发射方向图和接收波束形成程序差在哪很多读者把方向图计算理解成发射波束其实上面代码算的是接收波束权矢量对空间各方向的目标回波做匹配滤波。如果要算发射方向图需要把每个阵元的发射相位写到观察点然后对各阵元信号求和权矢量变成发射权。程序上差别很小发射方向图输出是abs(sum(w .* a, 3)).^2接收方向图是abs(sum(conj(w) .* a, 3)).^2一个共轭号之差。但物理上差别很大发射方向图是能量合成接收方向图是信号处理增益。FDA 多波束通常关心接收因为接收时可以做多目标匹配发射时只能叠加一个波束除非用正交波形。建议把两种都写在同一个函数里用参数决定。多波束还有一种做法是多频偏给每个波束分配不同的 Δf在接收端用滤波器组分离。这种方式抗干扰更强但程序里要同时维护多个频率偏移复杂度高不少。初学者先用“多权矢量、同一Δf”方案把波束间相关性调明白再考虑多频偏。4. 频率增量、阵元数与多波束权矢量调出稳定方向图的参数边界4.1 频率增量Δf怎么选距离周期、主瓣斜率和旁瓣Δf 是 FDA 波束形成程序里第一敏感参数。它直接出现在导向矢量的距离项上。Δf 越大距离维相位变化越快方向图在距离方向的斜率越大目标距离不同时主瓣角度偏移越明显。但副作用是距离周期 c/Δf 变小。例如 Δf10kHz 时周期 30kmΔf30kHz 时周期 10km。如果你的雷达要探测 20km 外的目标距离周期至少要大于最大作用距离否则同一个角度上会出现多个距离模糊峰。反过来Δf 太小比如 1kHz距离周期 300km看起来没有模糊但距离维波束太钝几乎退化成相控阵。一个合理的选参流程是先定最大作用距离 Rmax要求 Δf ≤ c/Rmax保证不模糊再用距离分辨率需求给下界。FDA 的距离维等效波束宽度大约由 Δf*N 决定实际工程里可以直接做一个 Δf 扫描程序一次性看多组参数下的方向图delta_f_list [5e3, 10e3, 20e3, 50e3]; for idf 1:length(delta_f_list) delta_f_tmp delta_f_list(idf); % 按最小程序生成P画图 end观察距离维主瓣宽度和是否出现周期重复。这个扫描脚本值得保留因为换雷达频段时只需重跑一遍。注意Δf 增大后旁瓣也会升高尤其是靠近 0° 方向。原因是频率偏移导致有效孔径出现幅度锥削。如果旁瓣超标首先要查的并不是窗函数而是 Δf 是否已经到“距离周期 2*Rmax”的边界。这时候减小 Δf 比加窗更有效。4.2 阵元数N与阵型多波束的自由度方向图角度维主瓣宽度和相控阵一样约 λ/(N*d) 弧度。N 增大波束更锐旁瓣更低但 N 也决定你可用的多波束自由度。用匹配权矩阵做 K 个波束时K 个权矢量之间的相关性取决于设计点间隔。理论上 N 个阵元可以提供 N 个自由度的波束成形但实际受波束宽度限制K 远小于 N。比如 N16角度波束宽度约 6°在 ±90° 范围里最多放 30 个波束但要让波束间隔大于 3dB 宽度K 实际只能在 10 上下。阵型也有讲究。均匀线阵只能把距离-角度耦合方向图投影到一维角度矩形平面阵多一个俯仰角同一个距离-方位单元可以再做俯仰分辨。程序改动量很小把代码里的阵元序号 n 换成阵元位置坐标矩阵相位项里的 ndsinθ 换成位置矢量与方向余弦的点乘。例如平面阵的第 m 个阵元位置 (xm, ym)则角度项是 xmcosφcosθ ymcosφsinθ具体取决于坐标系定义。初学建议先锁死均匀线阵把距离维跑顺再扩展不然调试时多一个维度的方位混叠很难分清是模型错还是参数错。还有一个和窗函数相关的经验先调 Δf再调 N最后加窗。三者同时动的话方向图变化叠加在一起你没法判断是哪个参数导致的旁瓣抬高。我一般会固定 N 和窗只扫 Δf确定 Δf 后再固定 Δf 和窗扫 N最后才加窗。4.3 设计点怎么摆避免波束间高相关多波束程序最容易出现的“看着能跑、输出一团糟”九成原因在设计点布局。两个波束设计点如果距离很近匹配权矢量高度相关任何一个目标都会被多个波束同时响应形成伪目标。经验做法是让相邻设计点的角度差大于单波束的 3dB 宽度距离差大于距离维主瓣的 3dB 宽度如果空间受限至少保证权矩阵 W 的 Gram 矩阵非对角元素低于 0.5。检验方法corr_W abs(conj(W) * W.) / N; % KxK 波束间相关如果第 i 行第 j 列大于 0.5就把这两个设计点往远拉。这个量纲很简单却能让多波束程序少走很多弯路。如果确实需要在很近的位置放两个波束可以对 W 做 Gram-Schmidt 正交化正交化后的权矢量旁瓣会抬高方向图峰值位置不变但副瓣水平需要重新评估。这里还要提醒一个常见误用不要用“增大 Δf”来强行增加多波束的区分度。Δf 增大改变的是距离维斜率不是波束间隔离度。两个波束在距离维重叠时Δf 只会让整条脊线变陡并不会让两个脊线分开。真正要拉开的是设计点的 r_k 和 θ_k。记住这句经验能让你少走很多路。5. FDA波束形成程序常见问题排查让距离维波束翻车的5个真实原因5.1 主瓣位置整体偏移设计点对不上现象程序设了 (r010km, θ0-10°)可 imagesc 画出来最大点跑到 (-8°, 12km) 一类的位置。原因最常见是 delta_f 单位搞混比如把 10e3 写成了 10也可能是 sind 和 sin 混用theta0 已经是弧度值却仍然套了 sind。另一个可能是在生成权矢量时用了 theta0但在生成扫描导向矢量时用了 Theta 却忘了度转弧度。解决先把程序里所有角度相关调用统一成 sind再用单点数值验证。在命令行打印phase_w(8)和phase_a(1,1,8)手算这两个复数相位差是否构成共轭关系。如果设计点是 -10°峰值却偏到正角度几乎可以肯定是阵列相位符号写反-2*pi/c的负号被漏掉或是乘了 n 后方向反了。用第6章的峰值轨迹验证法把所有角度下的峰值距离提出来和理论直线对比偏差超过一个网格步长就说明相位项有问题。5.2 距离维波束展得很宽看不出距离分辨现象方向图在角度方向有主瓣沿距离方向却几乎平铺好像没加 Δf。原因delta_f 太小或者 r_scan 网格太稀。尤其在 N16 时如果 delta_f 只有几百 Hz距离相位变化在扫描范围内只有几度方向图当然看不出距离选择。另一个原因是 phase_a 中距离项误乘了 f0 而不是 delta_f等于把 Δf 项写成了载频项距离分辨率完全丢失。解决把 delta_f 提高到 10kHz 以上同时把 r_scan 步长设为 c/(delta_f*N) 量级不要让每个单元的距离跨度超过半个主瓣宽度。在程序中可以加一条断言assert(delta_f c / max(r_scan), Δf过大距离周期小于扫描范围);这句话能挡掉一大半参数错误。同时检查距离项代码是否写成了f0 * R而不是delta_f * R。5.3 多波束输出互相串扰出现伪目标现象三个波束输出里每个波束的旁瓣都在其他波束主瓣位置留下凸起甚至出现比主瓣还高的峰。原因权矢量没归一化而且设计点间隔不够多个权矢量相关性强。归一化只影响输出幅度不会改变主瓣位置真正串扰的根源是相关性。解决先对 W 每行做W(k,:) W(k,:) / norm(W(k,:));归一化再用corr_W abs(conj(W) * W.) / N;看相关性。如果相关系数超过 0.5把对应设计点往远拉。如果设计点已经受系统约束不能动对 W 做正交化[Q, ~] qr(conj(W)., 0); W conj(Q).;这个操作让 K 个波束在无噪声时绝对正交但旁瓣会抬高需要重新测一遍方向图。注意正交化只适用于权矩阵不适用于发射方向图发射端的多波束不能这么做。5.4 程序跑得慢到没法调参现象角度网格 361、距离网格 301程序要跑几十秒把网格加密到 1000×1000 后几乎不动。原因代码里还留着三重循环对角度、对距离、对阵元逐个累加。MATLAB 的 for 循环在 10 万次以上时开销显著尤其是循环里有 exp 和复数乘法。解决把阵元维向量化用 ndgrid 生成网格按第3章的矩阵乘法一次算完。如果 A 太大分块处理比如距离维分 10 块每块生成 N×块内点数 的 A做完波束形成就释放。分块后内存可控速度反而比全矩阵快因为 exp 计算时缓存更友好。实际经验是N 通常只有几十瓶颈不在 N而在重复 for 遍历每个网格点所以矩阵化优先于一切。5.5 方向图在距离维出现周期性重复峰现象设计点是 10km但 20km、30km、40km 处出现相似形状的峰距离间隔固定。原因这就是 FDA 的距离周期。周期等于 c/Δf当 Δf 较大时周期可能只有几公里扫描范围内出现多个相同形状的峰。这不是数值 bug是 FDA 本身的特性。解决如果不需要距离模糊就把 Δf 减小让周期大于扫描范围。如果需要用大 Δf 保证距离分辨率那么系统设计就要加距离窗或者用多频点解模糊。程序层面可以在 5.2 的 assert 里检查delta_f c / max(r_scan)如果触发提示用户减小 Δf。注意有时候周期性重复峰并不是周期而是旁瓣需要比较峰的幅度如果重复峰的幅度和主瓣相近就是模糊如果低 13dB 以上是旁瓣。6. 验证与进阶用峰值轨迹确认程序正确再做多波束加速6.1 峰值轨迹验证法FDA 程序写完后先别急着接检测算法用一条直线验证模型正确性。根据线性模型设计点 (r0, θ0) 处的峰值轨迹满足 f0dsinθ Δfr f0dsinθ0 Δfr0。提取程序输出的峰值距离和理论直线对比[~, idx_r_peak] max(P, [], 1); r_peak r_scan(idx_r_peak); r_theory (f0*d*sind(theta0) delta_f*r0 - f0*d*sind(theta_scan)) / delta_f; plot(theta_scan, r_peak/1e3, o, theta_scan, r_theory/1e3, -);误差在一个网格步长内就说明导向矢量和权矢量写对了。6.2 多波束程序加速与内存控制矩阵化之后A 是 N×M 复数矩阵。N16、M108k 时 A 约 27MB可接受网格到 2000×2000 后 A 超过 1GB必须分块生成。分块时注意 exp 计算比例是最大的开销避免重复计算 phase。多波束的 K 个权矩阵乘积本身很快不需要额外优化。6.3 从方向图走向检测最后是一个职业习惯多波束程序不要只画方向图把 P_multi 保存成三维数组在距离维做归一化和 CFAR 门限每个波束切片就是一路检测输入。我早期写 FDA 程序时总想一步到位加高阶项结果方向图一团乱后来先按线性模型验证主瓣轨迹再逐项加 n²、加窗才稳定下来。波束形成程序不怕简单怕的是连一条直线都验证不过去就堆复杂度。希望帮到你。本文还有配套的精品资源点击获取