★ FEATURED ARTICLE 简介本资源是鄢社锋《优化阵列信号处理》前三章关键案例的Matlab实现代码面向学习阵列信号处理、波束形成与DOA估计的高年级本科生、研究生及科研人员帮助读者把书中理论落到可运行、可调试的仿真程序上。压缩包共27个文件以26个m脚本和1个txt说明为主整体约36KB脚本覆盖波束形成、方向图绘制、三维可视化与优化求解等典型环节便于按章节对照运行。目前已有8839人学习下载热度较高。读者可借助这些代码直观观察加权与相位调整对波束指向和旁瓣的影响复现最大似然、最小均方误差及DOA估计等算法并利用surf、slice、contourf等图形函数呈现功率谱密度与空间分布从而加深对阵列增益、干扰抑制与优化技术实际作用的理解提升动手仿真与算法验证能力。1. 阵列信号处理前三章到底在解决什么问题很多人第一次翻开《优化阵列信号处理》前三章看到波束形成、协方差矩阵、MVDR 这些词第一反应是公式能看懂但落到 Matlab 里就不知道从哪下手。我当初也是这样对着推导能点头打开编辑器就卡住。前三章的核心其实就三件事把阵列接收到的信号用数学描述清楚、用最优准则设计加权向量、再把这个加权向量变成能跑的矩阵运算。难点不在推导而在“公式到代码”这一步的映射关系。这篇笔记就是把这层映射拆开给出可复现的 Matlab 实现路径适合已经学过信号处理基础、想把书上前三章案例真正跑起来的人。你不需要先读完后面章节前三章的案例自成体系跑通之后对波束形成、MVDR、LCMV 这几个核心概念会有完全不同的体感。2. 均匀线阵的接收信号模型从公式到矩阵2.1 为什么信号模型是后面所有案例的地基前三章所有案例都建立在同一个前提上阵列接收到的数据可以写成一个矩阵形式。均匀线阵ULA是最常见的起点假设有 M 个阵元间距 d远场有 K 个窄带信号源入射角度分别是 θ₁ 到 θ_K。第 m 个阵元在 t 时刻的接收信号可以写成所有源信号的叠加加上噪声。这个描述听起来简单但落到代码里关键是要把“导向矢量”这个概念变成可操作的矩阵。导向矢量 a(θ) 描述的是当一个信号从 θ 方向入射时它在各个阵元上产生的相位差。对于 ULA第 m 个阵元的相位偏移是 exp(-j·2π·d·(m-1)·sin(θ)/λ)。这里 λ 是波长d 通常取半波长。把 M 个阵元的相位偏移排成一列就是导向矢量。K 个信号源的导向矢量排成矩阵就是阵列流形矩阵 A。接收信号矩阵 X A·S N其中 S 是源信号矩阵N 是噪声矩阵。这个公式在书上是几行推导但在代码里需要仔细处理维度。我见过不少人在这里翻车原因是把阵元数和快拍数的维度搞反了。X 的维度应该是 M×NM 是阵元数N 是快拍数。A 是 M×KS 是 K×NN 是 M×N。这样乘出来才是 M×N。如果维度不对后面求协方差矩阵 R X·X/N 就会直接报错或者算出莫名其妙的结果。2.2 用 Matlab 构造导向矢量和接收信号下面这段代码构造一个 8 阵元 ULA两个信号源分别从 -10° 和 20° 入射快拍数 200信噪比 10 dB。这是前三章最基础的仿真场景后面所有案例都可以在这个基础上改。% 均匀线阵接收信号模型 % 8阵元半波长间距两个信号源200快拍 M 8; % 阵元数 d 0.5; % 阵元间距单位波长 N 200; % 快拍数 theta [-10, 20]; % 信号入射角度单位度 K length(theta); % 信号源个数 SNR 10; % 信噪比 dB % 构造导向矢量矩阵 A维度 M×K A zeros(M, K); for k 1:K A(:, k) exp(-1j * 2 * pi * d * (0:M-1) * sind(theta(k))); end % 生成源信号维度 K×N S (randn(K, N) 1j * randn(K, N)) / sqrt(2); % 生成噪声维度 M×N noise (randn(M, N) 1j * randn(M, N)) / sqrt(2); noise noise * 10^(-SNR/20); % 接收信号矩阵 X维度 M×N X A * S noise; % 样本协方差矩阵 R维度 M×M R X * X / N;这段代码里几个参数需要说明。d 取 0.5 是半波长间距这是为了避免栅瓣也是书上前三章默认的设定。theta 用角度制sind 是 Matlab 里对角度求正弦的函数不要用 sin否则会得到错误结果。SNR 控制噪声功率10^(-SNR/20) 是把 dB 转成幅度比例。R 是样本协方差矩阵用 X·X/N 估计这是最直接的做法书上前三章用的也是这个。注意sind 和 sin 的区别在这里是致命的。如果你用 sin(theta) 而 theta 是角度值算出来的导向矢量完全错误后面波束图会乱成一团。2.3 协方差矩阵的估计与对角加载样本协方差矩阵 R 在快拍数不够多的时候会接近奇异求逆的时候数值不稳定。书上前三章在 MVDR 案例里提到了对角加载但很多人第一次看的时候会忽略这个细节。对角加载就是在 R 的对角线上加一个小的量让矩阵条件数变好。% 对角加载改善条件数 loading_factor 0.01; % 加载因子通常取噪声功率的倍数 R_loaded R loading_factor * trace(R) / M * eye(M);加载因子怎么选是个经验问题。太小起不到作用太大相当于在信号方向上也加了噪声波束图会变钝。我一般先取 0.01 倍的平均对角元素然后看波束图的主瓣宽度和旁瓣电平如果旁瓣太高就适当加大如果主瓣明显变宽就减小。这个参数没有唯一正确答案取决于你的阵列和信号环境。3. 波束形成从常规到 MVDR 的代码实现3.1 常规波束形成的加权向量与波束图常规波束形成就是最简单的加权想听哪个方向就把加权向量设成那个方向的导向矢量。w a(θ₀)/M其中 θ₀ 是你想对准的方向。波束图就是计算不同角度上的响应P(θ) |w·a(θ)|²。这个计算在 Matlab 里就是一个循环但要注意角度扫描的精度。% 常规波束形成 theta_scan -90:0.5:90; % 扫描角度范围 P_conv zeros(size(theta_scan)); theta0 0; % 对准0度 w_conv exp(-1j * 2 * pi * d * (0:M-1) * sind(theta0)) / M; for i 1:length(theta_scan) a_scan exp(-1j * 2 * pi * d * (0:M-1) * sind(theta_scan(i))); P_conv(i) abs(w_conv * a_scan)^2; end % 归一化并转dB P_conv_dB 10 * log10(P_conv / max(P_conv));这段代码里 theta_scan 的步长取 0.5°这是精度和计算量的折中。如果你要画更精细的图可以取 0.1°但计算量会大很多。P_conv_dB 归一化到最大值这样波束图的峰值是 0 dB旁瓣电平可以直接读出来。常规波束形成的旁瓣大约在 -13 dB 左右主瓣宽度取决于阵元数和间距。3.2 MVDR 波束形成的完整实现MVDR 是前三章的重点案例。它的加权向量是 w R⁻¹·a(θ₀) / (a(θ₀)·R⁻¹·a(θ₀))。这个公式看起来简单但实现的时候有几个坑。第一是 R 求逆的数值稳定性第二是分母可能很小导致加权向量数值爆炸。% MVDR波束形成 theta0 0; % 期望信号方向 a0 exp(-1j * 2 * pi * d * (0:M-1) * sind(theta0)); % 使用对角加载后的协方差矩阵 R_inv inv(R_loaded); % 计算MVDR加权向量 w_mvdr R_inv * a0 / (a0 * R_inv * a0); % 计算波束图 P_mvdr zeros(size(theta_scan)); for i 1:length(theta_scan) a_scan exp(-1j * 2 * pi * d * (0:M-1) * sind(theta_scan(i))); P_mvdr(i) abs(w_mvdr * a_scan)^2; end P_mvdr_dB 10 * log10(P_mvdr / max(P_mvdr));MVDR 的波束图在干扰方向会形成零陷这是它和常规波束形成最大的区别。零陷的深度取决于协方差矩阵估计的准确度和对角加载的大小。如果零陷不够深通常是快拍数不够或者加载因子太大。如果零陷位置偏移检查你的导向矢量角度定义是否和信号模型一致。3.3 两种波束形成的对比与参数影响把常规波束形成和 MVDR 的波束图画在一起能直观看到区别。常规波束形成的主瓣宽度由阵列孔径决定MVDR 的主瓣宽度类似但在干扰方向有明显的零陷。这个零陷是 MVDR 用协方差矩阵“学”出来的不需要知道干扰方向的具体角度。% 对比绘图 figure; plot(theta_scan, P_conv_dB, b-, LineWidth, 1.5); hold on; plot(theta_scan, P_mvdr_dB, r--, LineWidth, 1.5); xlabel(角度 (度)); ylabel(归一化功率 (dB)); legend(常规波束形成, MVDR); grid on; ylim([-60, 0]);从图上能看到MVDR 在 -10° 和 20° 附近形成了零陷因为这两个方向有信号源协方差矩阵里包含了它们的信息。常规波束形成没有这个能力它只对准 0°对其他方向的信号只是被动接收。这个对比是书上前三章的核心案例之一跑通之后对自适应波束形成的理解会深很多。注意如果你发现 MVDR 的零陷没有出现在信号方向先检查协方差矩阵是不是用 X·X/N 算的再检查导向矢量的符号约定是否一致。符号约定不一致是新手最常见的错误。4. LCMV 与稳健波束形成约束条件的代码落地4.1 LCMV 的约束矩阵构造LCMV 是在 MVDR 基础上加约束。MVDR 只约束期望方向增益为 1LCMV 可以同时约束多个方向的增益。约束条件写成 C·w f其中 C 是约束矩阵f 是约束值向量。解是 w R⁻¹·C·(C·R⁻¹·C)⁻¹·f。% LCMV波束形成 % 约束0度增益为1-30度增益为0 C zeros(M, 2); C(:, 1) exp(-1j * 2 * pi * d * (0:M-1) * sind(0)); C(:, 2) exp(-1j * 2 * pi * d * (0:M-1) * sind(-30)); f [1; 0]; % LCMV加权向量 w_lcmv R_inv * C * inv(C * R_inv * C) * f; % 计算波束图 P_lcmv zeros(size(theta_scan)); for i 1:length(theta_scan) a_scan exp(-1j * 2 * pi * d * (0:M-1) * sind(theta_scan(i))); P_lcmv(i) abs(w_lcmv * a_scan)^2; end P_lcmv_dB 10 * log10(P_lcmv / max(P_lcmv));这段代码里 C 是 M×2 的矩阵两列分别对应两个约束方向。f 是 2×1 的向量第一个元素 1 表示 0° 方向增益为 1第二个元素 0 表示 -30° 方向增益为 0。C·R⁻¹·C 是 2×2 的矩阵求逆很快。如果约束条件更多这个矩阵会变大但一般不会超过几个约束否则会消耗太多自由度。4.2 对角加载对稳健性的影响对角加载在 LCMV 里同样重要。如果 R 接近奇异R⁻¹ 会数值不稳定加权向量可能非常大。加载因子的大小直接影响波束图的形状。我一般会做一组对比加载因子取 0.001、0.01、0.1看波束图的变化。% 不同加载因子对比 loading_factors [0.001, 0.01, 0.1]; figure; for idx 1:length(loading_factors) R_l R loading_factors(idx) * trace(R) / M * eye(M); w_l inv(R_l) * C * inv(C * inv(R_l) * C) * f; P_l zeros(size(theta_scan)); for i 1:length(theta_scan) a_scan exp(-1j * 2 * pi * d * (0:M-1) * sind(theta_scan(i))); P_l(i) abs(w_l * a_scan)^2; end P_l_dB 10 * log10(P_l / max(P_l)); plot(theta_scan, P_l_dB, LineWidth, 1.5); hold on; end legend(加载0.001, 加载0.01, 加载0.1); xlabel(角度 (度)); ylabel(归一化功率 (dB)); grid on; ylim([-60, 0]);从对比图能看到加载因子太小的时候波束图可能有毛刺太大则主瓣变宽、零陷变浅。0.01 左右通常是个不错的起点。这个参数需要根据你的具体场景调没有万能值。4.3 稳健波束形成的几种常见思路书上前三章还提到了稳健波束形成。实际中导向矢量总有误差可能是阵元位置不准也可能是信号方向估计有偏差。稳健波束形成的思路是让加权向量对误差不敏感。常见做法有几种一是对角加载这个前面已经用了二是基于不确定集的稳健设计把导向矢量误差建模成一个球或椭球三是用协方差矩阵拟合。对角加载是最简单也最常用的。它的本质是在协方差矩阵上加一个各向同性的噪声相当于告诉算法“我对信号方向没那么确定”。加载因子越大稳健性越好但分辨率越差。这个折中在实际系统里需要根据误差大小来定。如果阵元位置误差在 0.1 波长以内加载因子取 0.01 到 0.05 通常够用。5. 避坑与排查前三章案例实现中的常见问题5.1 波束图对称性异常现象波束图左右不对称或者峰值不在设定的对准方向。原因通常是导向矢量的符号约定不一致。有的书写 exp(-j·...)有的写 exp(j·...)如果信号模型和波束形成用了不同的符号就会出现角度偏移。解决方法是统一符号约定我一般全部用 exp(-j·2π·d·(m-1)·sin(θ)/λ)从信号模型到波束形成不改。5.2 协方差矩阵求逆报错或结果异常现象inv(R) 报奇异矩阵警告或者求出来的加权向量数值巨大。原因是快拍数不够R 接近奇异。解决方法是加对角加载或者增加快拍数。如果快拍数已经很大还是奇异检查信号源之间是否高度相关相关信号会导致协方差矩阵秩亏。这时候需要用空间平滑等去相关技术但前三章不涉及可以先增加信号源之间的角度间隔。5.3 MVDR 零陷位置偏移现象MVDR 波束图的零陷没有出现在信号入射方向而是偏了几度。原因是角度扫描的步长太大或者信号方向估计有误差。解决方法是在零陷附近加密扫描比如用 0.1° 步长。如果零陷位置和信号方向差得比较多检查信号模型里的 theta 和波束形成里的 theta_scan 是否用了同一个角度定义。5.4 对角加载后零陷变浅现象加了加载因子之后MVDR 的零陷深度明显变浅。原因是对角加载相当于在协方差矩阵里加了白噪声削弱了干扰方向的信息。解决方法是减小加载因子或者用更精细的加载策略比如只在噪声子空间加载。如果加载因子已经很小零陷还是不够深检查快拍数是否足够快拍数太少的时候协方差矩阵估计不准零陷自然浅。5.5 波束图旁瓣电平异常高现象常规波束形成的旁瓣应该在 -13 dB 左右但实际算出来高很多。原因是阵元间距不是半波长或者阵元数太少。阵元间距大于半波长会产生栅瓣旁瓣电平会升高。解决方法是把 d 设成 0.5或者用非均匀阵列。如果 d 已经是 0.5 旁瓣还是高检查加权向量是否归一化w a(θ₀)/M 里的除以 M 不能漏。6. 从跑通到用对几个让代码更接近实战的技巧跑通前三章的案例只是第一步真正要用起来还需要注意几个细节。第一个是角度扫描的精度和计算量的平衡。0.5° 步长在大多数场景够用但如果你要精确测量零陷深度0.1° 更合适。计算量增加 5 倍但现代机器完全能承受。我一般先用 0.5° 快速看整体形状确认没问题再用 0.1° 出图。第二个是协方差矩阵估计的改进。书上前三章用的是样本协方差矩阵实际中可以用前向-后向平均来改善估计质量。做法是把 X 和它的共轭反转版本拼在一起算协方差这样等效快拍数翻倍对相干信号也有一定的去相关效果。% 前向-后向平均协方差矩阵 X_fb [X, conj(flipud(X))]; R_fb X_fb * X_fb / (2 * N);这段代码里 flipud(X) 是把 X 上下翻转conj 取共轭。拼起来之后等效快拍数是 2N。这个方法在快拍数少的时候特别有用能让 MVDR 的零陷更深更稳。第三个是加权向量的归一化。MVDR 和 LCMV 的加权向量不需要像常规波束形成那样除以 M因为公式里已经包含了归一化。但如果你要比较不同算法的输出功率需要确保期望方向的增益一致。我一般会在计算完加权向量后检查 w·a(θ₀) 是否等于 1如果不是就手动归一化。% 检查期望方向增益 gain_check w_mvdr * a0; if abs(gain_check - 1) 1e-6 w_mvdr w_mvdr / gain_check; end这个检查在调试的时候很有用。如果增益不是 1说明公式实现有误或者对角加载改变了约束。LCMV 的约束是硬约束增益应该严格等于 f如果不等检查 C 和 f 的构造是否正确。最后一个技巧是关于波束图的显示。归一化到 0 dB 是常见做法但如果你要比较不同算法的绝对性能应该用相同的参考。我一般会保留原始功率值只在绘图的时候归一化。这样如果某个算法的输出功率整体偏低你能从原始值里看出来。这些技巧在书上前三章不会全部展开但实际跑代码的时候会一个个遇到。我的习惯是每跑通一个案例就把参数改一改看波束图怎么变。改得多了对参数的理解就从“书上说取 0.01”变成“我知道为什么取 0.01也知道取 0.1 会怎样”。这个体感是看公式得不到的只能靠一行行代码跑出来。希望帮到你。本文还有配套的精品资源点击获取
阅读完成 · 觉得有帮助？
