简介：面向雷达信号处理与自适应阵列研究者的Matlab算法实现包，聚焦LMS、RLS、SMI三种自适应波束形成方法，解决动态环境下信号检测与干扰抑制的工程调参问题。压缩包内共3个文件，全部为Matlab源码（.m），分别对应LMS最小均方误差、RLS递归最小二乘、SMI符号大小积分算法的主程序框架，包体仅3KB，轻量便于快速部署与二次修改。已有416人学习下载，适合通信、雷达方向的学生及工程师对照理论推导进行仿真验证。通过阅读和运行这三个脚本，可以直观比较不同算法的收敛速度、稳态误差与计算复杂度，掌握学习速率、滤波器长度等关键参数的设置方法，并结合多径传播、低信噪比等场景理解算法选型依据。源码结构简单清晰，可直接作为自适应波束形成课程设计或项目预研的起点，也为后续算法改进提供了可扩展基础。
1. 自适应波束形成为什么绕不开 LMS、RLS 与 SMI
雷达信号处理里，天线阵接收到的往往是期望信号、旁瓣干扰和噪声的混合，阵列输出质量完全取决于各阵元的复加权系数。自适应波束形成的目标就是实时估计干扰协方差、把零陷对准干扰方向，而工程中最常用的三条路线正好是 LMS、RLS 与 SMI：LMS 用梯度下降逐点修正权值，算力最低但收敛慢；RLS 用递归最小二乘换取更快的收敛；SMI 直接对采样协方差矩阵求逆，一次快拍块就能解出权值。如果你在调相控阵仿真、做雷达信号处理课程设计，或者想在 MATLAB 里快速对比三种自适应算法的波束图，这套LMS.m、RLS.m、SMI.m源码就是最直接的起点。你不用关心推导细节，照着这三个文件的逻辑去改参数，就能看到波束图的变化。
2. 阵列模型与三种自适应算法的选型逻辑
2.1 均匀线阵输出模型：先定导向矢量与快拍矩阵
自适应波束形成的第一步不是写循环，而是把阵列输出建模成矩阵形式。以均匀线阵为例，设阵元数为M，阵元间距为d，载波波长为lambda，取d = lambda / 2避免栅瓣。来自方向theta的信号在相邻阵元间产生的相位差为2 * pi * d * sin(theta) / lambda，指向theta0的导向矢量写作：
a(theta0) = [1, exp(j*2*pi*d*sin(theta0)/lambda), ..., exp(j*2*pi*(M-1)*d*sin(theta0)/lambda)]^T
在 MATLAB 里我一般这样生成导向矢量：
M = 8; % 阵元数 lambda = 0.3; % 波长，单位 m d = lambda / 2; % 阵元间距 theta0 = 0; % 期望信号方向，单位度 % 匿名函数：输入角度，输出该方向的导向矢量 a_theta = @(theta) exp(1j * 2 * pi * d / lambda * (0:M-1)' * sind(theta)); a0 = a_theta(theta0); % 期望方向导向矢量代码逻辑：(0:M-1)'把阵元序号变成M x 1列向量，sind(theta)直接以度为单位计算正弦，省去deg2rad的重复调用。生成导向矢量后，阵列在某一个快拍内的接收向量为x = a0 * s + sum(ai * ji) + n，其中s是期望信号复包络，ji是第i个干扰，n是噪声。多快拍数据拼成矩阵X = [x(1), x(2), ..., x(N)]，每个列向量是一个快拍。
在 L 波段雷达里，lambda = 0.3 m对应 1 GHz，8 元阵的孔径只有约 1.05 m，波束宽度约 12.8 度；如果换成 32 元阵，孔径变成 4.65 m，主瓣明显变窄。仿真前先想清楚要看 8 元还是 32 元，因为后续 LMS 的步长和 SMI 的快拍数需求都跟着M变。
2.2 三种算法对比：收敛速度、计算量与适用场景
LMS 本质是随机梯度下降，权重沿瞬时梯度方向更新，公式为w(n+1) = w(n) + mu * conj(e(n)) * x(n)，每个快拍只做O(M)次乘加。它的问题在于收敛步长受输入自相关矩阵特征值散布影响，特征值比大时收敛非常慢。RLS 引入代价函数中的遗忘因子，用递归方式估计自相关矩阵的逆，每步更新量达到O(M^2)，但收敛速度基本不再受特征值散布的限制。SMI 是块处理算法，直接利用N个快拍估计协方差矩阵并求逆，理论上收敛最快，代价是单次运算O(M^3)且矩阵求逆可能病态。
三者的对比可以直接看表格：
| 算法 | 单次更新复杂度 | 收敛速度 | 快拍数需求 | 典型使用场景 |
|---|---|---|---|---|
| LMS | O(M) | 慢，受特征值散布影响 | 逐快拍更新 | 实时性强、算力受限的嵌入式系统 |
| RLS | O(M^2) | 快，接近理论最优 | 几十个快拍内收敛 | 强干扰、动态变化目标 |
| SMI | O(M^3)（矩阵求逆） | 最快，块处理 | 建议 N >= 2M 到 5M | 雷达批处理、离线分析与课程设计 |
提示：表中的快拍数需求针对不含对角加载的理想情况。实际工程里快拍数不够时，SMI 的主瓣会畸变，这时一般要加对角加载，见第 4 章。
2.3 SMI 名字的常见混淆：别把两种实现搞混
搜索 SMI 时你可能会看到两种解释：一种是摘要里写的 Sign-Magnitude-Integration，另一种是经典文献中的 Sample Matrix Inversion（采样矩阵求逆）。在自适应波束形成领域，SMI 几乎都指后者，即用采样协方差矩阵求逆后计算权值。拿到源码后可以先用type SMI.m命令检查，看到inv或\运算就是采样矩阵求逆版本；如果看到sign函数，那才是符号大小积分变体。二者目标相同，但数学框架不同，做实验对比时不要混用。我在这个项目场景里习惯把SMI.m按采样矩阵求逆来实现，后面的代码也都按这个版本给出。
3. LMS 波束形成：梯度下降实现与步长坑点
3.1 迭代公式与参考信号的选取
在波束形成场景里，LMS 的参考信号d(n)不是凭空给出的目标波形，而是期望方向的导向矢量与发射信号的乘积，或者直接用一个与信号形式匹配的导引序列。这一点初学者常搞错，随手用一个随机序列当d，结果权值发散。我的习惯是：训练阶段用本地副本d(n) = s(n)（例如雷达线性调频信号），同时把阵列接收置为x(n)；两者同步后，LMS 会让输出尽量逼近参考信号，等效于在干扰方向形成零陷。
权值迭代的完整形式为：
w(n+1) = w(n) + mu * conj(e(n)) * x(n)
其中误差e(n) = d(n) - w(n)' * x(n)。mu是步长，conj取共轭是因为 MATLAB 默认变量为复数，梯度方向需要对误差取共轭才能保持维度对齐。如果信号是实基带，conj没有影响；但雷达多普勒处理后的数据基本都是复数，去掉conj会导致权值更新方向错误。
3.2 LMS.m 的最小实现
压缩包里的LMS.m建议按下面的骨架实现，保证输入输出和主脚本解耦：
function w = LMS(x, d, mu, M) % x : M x N 阵列接收矩阵，每列是一个快拍 % d : 1 x N 参考信号 % mu: 步长，建议初值 0.001~0.01 % w : M x 1 自适应权值 [N, snap] = size(x); w = zeros(M, 1); y = zeros(1, snap); for n = 1:snap x_n = x(:, n); y(n) = w' * x_n; % 当前阵列加权输出 e = d(n) - y(n); % 参考信号与实际输出的误差 w = w + mu * conj(e) * x_n; % 沿瞬时梯度方向更新 end end代码逻辑：y = w' * x_n是当前阵列加权输出；e是误差；权值更新沿瞬时梯度方向conj(e) * x_n前进。注意M作为参数传入，避免在函数里靠size(x, 1)硬编码，后续换阵元数时不用改函数体。步长mu越大更新越快，但超过稳定上限就会振荡，稳定条件为0 < mu < 2 / trace(Rxx)，Rxx是输入自相关矩阵。
y数组在这里用于后续画学习曲线：把e的平方记录下来，就能看到收敛过程。如果只想返回权值，去掉y相关行即可，但保留它对验证算法很有利。
3.3 步长选择与特征值散布问题
对于 8 元均匀线阵，输入功率归一化到 1、干噪比 30 dB 时，trace(Rxx)大概在 16 左右，mu取 0.005 是安全值。阵元数增加到 32，mu要同比例缩小到 0.001 量级。更稳的做法是使用归一化 LMS，权重更新改为：
w = w + beta * conj(e) * x_n / (delta + x_n' * x_n)
其中delta是防止除零的小常数，beta取 0.1~0.3。归一化 LMS 的步长随输入功率自动缩放，不必每次手动调整mu。一组经验参数参考如下：
| 阵元数 M | 输入功率 | 建议 mu | 收敛到 -20 dB 所需快拍 |
|---|---|---|---|
| 8 | 1 | 0.005 | 约 400~600 |
| 16 | 1 | 0.002 | 约 800~1200 |
| 32 | 1 | 0.001 | 约 1500~2000 |
表中数值来自干噪比 30 dB、期望信号 SNR 0 dB 的仿真，实际以学习曲线为准。如果拿到LMS.m源码后发现迭代不收敛，优先检查三件事：参考信号是否与期望信号对齐、mu是否超过2 / trace(Rxx)、是否忘记把输入数据转成复基带（直接传中频实数信号会导致频率偏移对消掉）。
4. RLS 与 SMI 实现：递归求逆与协方差矩阵稳定化
4.1 RLS 递推公式与 P 矩阵更新
RLS 不直接求解自相关矩阵，而是维护其逆P(n)，权值更新分为三步：计算增益向量、更新权值、更新逆矩阵。遗忘因子lambda决定了对历史数据的记忆长度，越接近 1，数据窗越长，稳态失调越小，但跟踪能力变差。雷达目标快速机动时我一般取lambda = 0.98；阵列静止且干扰缓慢变化时取lambda = 0.999。
function w = RLS(x, d, lambda, M) % lambda: 遗忘因子 0.98~0.999 [N, snap] = size(x); w = zeros(M, 1); P = eye(M) / 1e-3; % 初始逆相关矩阵，取小分母避免初始增益过小 for n = 1:snap x_n = x(:, n); k = (P * x_n) / (lambda + x_n' * P * x_n); % 增益向量 e_pri = d(n) - w' * x_n; % 先验误差 w = w + k * conj(e_pri); % 权值更新 P = (P - k * x_n' * P) / lambda; % 逆矩阵递推 end endP的初值取eye(M) / 1e-3是工程习惯：P(0)越大，初始几步的步长越大，适应速度快，但太大会让前几十个快拍剧烈抖动。e_pri是更新前的先验误差，和 LMS 里用更新后的误差不同，这是 RLS 数学推导的特点。若仿真中出现权值 NaN，多半是P矩阵递推失去正定性，常见原因包括lambda太小或输入含纯零快拍。加一行P = P + 1e-6 * eye(M)可以缓解。
4.2 SMI 块处理实现：从样本协方差到最优权
SMI 的核心思想是块处理：采集 N 个快拍，估计样本协方差矩阵，然后直接计算最优权值。以最小方差无失真响应形式给出权值：
w_smi = (R_hat \ a0) / (a0' * (R_hat \ a0))
对应代码：
function w = SMI(x, a0) % x : M x N 快拍矩阵 % a0 : M x 1 期望方向导向矢量 [M, N] = size(x); R_hat = (x * x') / N; % 样本协方差矩阵 w = (R_hat \ a0) / (a0' * (R_hat \ a0)); end代码逻辑：R_hat = (x * x') / N是最大似然意义下的协方差估计；R_hat \ a0用左除求解线性方程，比inv(R_hat) * a0数值更稳定，尤其当R_hat接近奇异时，inv会放大浮点误差，而\走 LU 或 Cholesky 路径。干扰方向完全被抑制的前提是快拍数足够多，若N < 2M，R_hat不满秩，权值会严重畸变，主瓣指向偏移。SMI 技术在实际雷达系统中通常会结合多普勒滤波后的数据使用，先做多脉冲积累再估计协方差，这样单个距离单元的快拍数会显著增加。
4.3 对角加载：让 SMI 在低快拍下不炸
低快拍时直接对R_hat求逆，特征值散布极大，噪声对应的最小特征值会被放大，波束图凹坑变浅，输出信干噪比下降。常见做法是对角加载，把R_hat替换为：
R_dl = R_hat + gamma * eye(M)
加载量gamma我一般选噪声方差的 10 倍，或R_hat最大特征值的 1/100。这样约束了最小特征值下限，让解在保持干扰零陷的同时不过度放大噪声。加入对角加载后 SMI 权值计算变为：
gamma = 0.1 * mean(eig(R_hat)); % 经验值，也可用 10 * sigma_n^2 R_dl = R_hat + gamma * eye(M); w_smi = (R_dl \ a0) / (a0' * (R_dl \ a0));mean(eig(R_hat))近似于平均输入功率，用它的 10% 做加载量在大多数仿真里都能既保留零陷深度，又压低旁瓣。加载量过大，波束会退化为常规延迟相加波束形成，自适应能力被磨平；过小则低快拍时优化效果不明显，这是 SMI 调参中最需要权衡的一组矛盾。对应场景的参数推荐如下：
| 场景 | 随手设置的效果 | 推荐参数 |
|---|---|---|
| 静止阵列、N=100 快拍 | lambda=0.99 输出平稳但跟踪慢 | lambda=0.995，P0=eye(M)/1e-3 |
| 目标在 3 秒内走完 15 度 | lambda=0.995 跟踪滞后明显 | lambda=0.97~0.98 |
| 快拍 N=8 且 M=8 | SMI 矩阵奇异，主瓣偏移 | 对角加载 gamma=0.1*mean(eig(R_hat)) |
| N=200 且 M=8 | 收敛稳定，旁瓣约 -13 dB | 无需加载，直接用 R_hat \ a0 |
4.4 三个算法放在同一份数据上的预期表现
同一组快拍数据下，LMS 大概在 500 个快拍后才进入稳态，RLS 在 50 个快拍内已把权值拉到位，SMI 则是算出多少快拍就用多少信息，快拍充足时基线最好。但这不代表 SMI 永远最优：快拍块长度增加，协方差矩阵估计越准，但也意味着对干扰变化的响应越慢，本质上是一个时间窗选择问题。RLS 的遗忘因子同样在调节时间窗，只是通过指数衰减实现。LMS 没有明确时间窗，完全靠步长控制记忆长度，所以它对非平稳环境的适应最粗糙。
排错时优先级从高到低：先把数据格式统一成M x N复矩阵，再检查a0是否归一化，最后看权值模值。权值模值突然变成NaN或接近 1e15，直接怀疑矩阵求逆失败，回到对角加载那一步。
5. 波束图验证与参数调试技巧
5.1 用一张波束图同时验证三种算法
把三种算法跑完后，统一计算波束图并叠加对比：
% 计算三种权值 w_lms = LMS(x, s, 0.005, M); w_rls = RLS(x, s, 0.999, M); w_smi = SMI(x, a0); % 扫描角度并计算方向图 theta = -90:0.1:90; A = a_theta(theta); % M x length(theta) 的导向矢量矩阵 pattern_lms = abs(w_lms' * A); pattern_rls = abs(w_rls' * A); pattern_smi = abs(w_smi' * A); % 归一化后画 dB 图 figure; plot(theta, 20*log10(pattern_lms / max(pattern_lms)), 'LineWidth', 1.2); hold on; plot(theta, 20*log10(pattern_rls / max(pattern_rls)), 'LineWidth', 1.2); plot(theta, 20*log10(pattern_smi / max(pattern_smi)), 'LineWidth', 1.2); legend('LMS', 'RLS', 'SMI'); xlabel('角度/deg'); ylabel('归一化幅度/dB'); grid on;代码逻辑：A = a_theta(theta)一次性生成所有扫描角度的导向矢量矩阵，w' * A就是阵列响应。归一化只在画图时做，不影响权值本身。三个算法期望方向都指向 0 度时，主瓣应该在 0 dB 附近重叠，区别主要体现在干扰方向 30 度附近的零陷深度。
5.2 两个最容易忽略的验收指标
第一个指标是干扰方向响应：读20*log10(pattern(30))处的值，低于 -40 dB 才算零陷有效。如果只看到 -20 dB，优先怀疑快拍数不足或 LMS 未收敛。第二个指标是cond(R_hat)：条件数超过 1e12 时，SMI 计算结果基本不可信，必须做对角加载。这两个指标比单纯看波束图形状更有说服力，评审或答辩时拿出来也是加分项。
5.3 快速定位“为什么波束没对准期望方向”
出现波束主瓣偏移时，先检查a0的生成代码是不是把sind(theta)写成了sin(theta)，这一步错得最隐蔽，因为 0 度时两者结果相同，换成 30 度立即偏差。再看权值向量的相位分布，正常时相邻阵元权值相位差应接近 0，出现随机跳变说明迭代没有收敛。做课程设计时尽量把cond(R_hat)、trace(Rxx)、学习曲线这三项打印出来，这些中间量比只看一张波束图更容易定位问题，也更能体现对自适应波束形成算法的理解深度。
本文还有配套的精品资源，点击获取
