简介：这份资源是鄢社锋《优化阵列信号处理》前三章关键案例的Matlab实现代码，面向正在学习阵列信号处理、波束形成与DOA估计的高年级本科生、研究生及工程技术人员。压缩包共27个文件，以26个m脚本和1个txt说明为主，整体约36KB，脚本按书中案例编号组织，便于对照章节逐一运行验证。内容覆盖波束形成加权、方向图与功率谱的三维可视化、以及梯度下降等优化方法在阵列处理中的应用，读者可借助这些代码直观观察信号空间分布、复现书中算法流程，并在此基础上修改参数开展自己的仿真实验。目前已有8839人学习下载，适合希望把理论推导落到代码实践、提升动手能力的学习者参考使用。
1. 阵列信号处理入门：为什么前三章案例值得逐行复现
很多人学阵列信号处理，卡在同一个地方：公式能看懂，推一遍也顺，但真到写代码就不知道从哪下手。波束形成、空间谱估计、自适应零陷这些概念，书上写得清楚，可一到工程落地就变成玄学——协方差矩阵怎么估、对角加载加多少、快拍数不够会怎样，这些细节没人替你踩坑。某高校实验室里流传的《优化阵列信号处理》前三章案例，恰好是打通这条链路的入口：它不讲空泛理论，而是用最小可运行的 Matlab 脚本，把均匀线阵的波束图、MVDR 自适应波束形成、MUSIC 空间谱估计三件事串起来。适合谁？刚接触阵列信号的研究生、做雷达或声呐预研的工程师、以及想把课本公式变成可调参数 Demo 的开发者。这一篇不抄书，只讲我照着前三章案例复现时，哪些参数必须调、哪些结果不能信、哪些坑绕不过去。
2. 均匀线阵波束形成：从阵列流形到波束图的完整链路
2.1 阵列流形矩阵的构造逻辑与参数含义
均匀线阵（ULA）是所有阵列处理案例的起点。核心就一个矩阵：阵列流形矩阵 A，尺寸为 N×M，N 是阵元数，M 是扫描角度数。每一列对应一个来波方向的方向矢量。很多新手直接调steervec就完事，但前三章案例要求手写，目的是让你看清相位差的来源。
假设阵元间距 d，波长 lambda，来波方向 theta，以第一个阵元为参考，第 n 个阵元的相位延迟是-2*pi*d*(n-1)*sin(theta)/lambda。注意符号约定：不同教材正负号相反，导致波束图左右镜像。我一般固定用“负号对应来波方向”，这样扫描角度为正时波束指向正角度，不容易搞混。
% 均匀线阵流形矩阵构造 N = 16; % 阵元数 d = 0.5; % 阵元间距，单位波长 lambda = 1; % 归一化波长 theta = -90:0.5:90; % 扫描角度，度 A = zeros(N, length(theta)); for m = 1:length(theta) for n = 1:N A(n,m) = exp(-1j*2*pi*d*(n-1)*sind(theta(m))/lambda); end end这段代码里d=0.5是半波长间距，避免栅瓣。如果你改成d=1，波束图会在 ±90 度附近出现虚假峰值，这就是空间混叠。theta步长 0.5 度足够画图，再细对结果没本质提升，但计算量翻倍。sind用的是角度制，别错写成sin把弧度塞进去，否则波束指向全错。
2.2 常规波束形成与 MVDR 的对比实现
常规波束形成就是拿流形矩阵的共轭转置去加权，权向量w = A(:,idx)/N，其中 idx 是期望方向对应的角度索引。MVDR 则要用协方差矩阵求逆：w = R^-1 * a / (a' * R^-1 * a)。案例里通常先仿真一个单信源加噪声的场景，再对比两种波束图。
% 仿真接收信号：单信源 + 高斯白噪声 K = 200; % 快拍数 theta_s = 20; % 信源方向 a_s = exp(-1j*2*pi*d*(0:N-1)'*sind(theta_s)/lambda); S = (randn(1,K) + 1j*randn(1,K))/sqrt(2); X = a_s * S + (randn(N,K) + 1j*randn(N,K))/sqrt(2); % 协方差矩阵估计 R = X * X' / K; % 常规波束形成权向量 idx = find(theta == theta_s); w_cbf = A(:,idx) / N; P_cbf = abs(w_cbf' * A).^2; % MVDR 权向量 w_mvdr = (R \ a_s) / (a_s' * (R \ a_s)); P_mvdr = abs(w_mvdr' * A).^2;关键参数是快拍数 K。K 小于阵元数 N 时，协方差矩阵秩亏，R不可逆，MVDR 直接翻车。我一般要求 K ≥ 2N，最好 5N 以上。对角加载就是给 R 加一个sigma*eye(N)，sigma 取噪声功率的 1% 到 10%，能缓解小快拍下的病态求逆。案例里没强调这点，但实际跑的时候 K=50、N=16 就会看到 MVDR 波束图畸变，主瓣偏移，这就是血泪经验。
2.3 波束图绘制与主瓣宽度、旁瓣电平的读取
画图不是目的，读参数才是。主瓣宽度通常看 -3dB 宽度，旁瓣电平看第一旁瓣比主瓣低多少 dB。均匀加权时第一旁瓣约 -13.2dB，切比雪夫加权可以压到 -30dB 但主瓣变宽。案例前三章一般只做均匀加权，但你要知道这个基准。
% 归一化并画图 P_cbf_dB = 10*log10(P_cbf / max(P_cbf)); P_mvdr_dB = 10*log10(P_mvdr / max(P_mvdr)); plot(theta, P_cbf_dB, 'b', theta, P_mvdr_dB, 'r'); xlabel('角度 (度)'); ylabel('归一化功率 (dB)'); legend('CBF','MVDR'); grid on;读图时注意：MVDR 在信源方向形成主瓣，在干扰方向形成零陷。如果仿真里只有一个信源，MVDR 和 CBF 差别不大；要体现 MVDR 优势，必须加一个干扰源，比如 theta_i = -30 度，功率比信号高 10dB。这时候你会看到 MVDR 在 -30 度处深深凹下去，而 CBF 只是略微降低。这个对比才是前三章案例想让你看到的。
3. MUSIC 空间谱估计：从子空间分解到谱峰搜索
3.1 协方差矩阵特征分解与信号子空间/噪声子空间划分
MUSIC 的核心假设是：信号子空间和噪声子空间正交。对协方差矩阵 R 做特征分解，大特征值对应信号，小特征值对应噪声。信源数估计不准，MUSIC 就废了。案例里通常假设信源数已知，但实际中你要用 AIC 或 MDL 准则去估。
% 特征分解 [V, D] = eig(R); eig_val = diag(D); [~, idx_sort] = sort(eig_val, 'descend'); V = V(:, idx_sort); eig_val = eig_val(idx_sort); % 假设信源数 P=1 P = 1; En = V(:, P+1:end); % 噪声子空间eig出来的特征值顺序是乱的，必须排序。En的列数等于 N-P。如果 P 估多了，噪声子空间里混入信号分量，谱峰变钝；P 估少了，信号子空间不完整，谱峰偏移。我一般先画特征值分布图，看有没有明显的拐点。
3.2 MUSIC 谱函数计算与谱峰搜索的代码实现
MUSIC 谱函数是P_music = 1 / (a' * En * En' * a)，对每个扫描角度算一遍。注意分母可能接近零，导致数值爆炸，实际实现时加一个极小量 eps。
P_music = zeros(size(theta)); for m = 1:length(theta) a = A(:,m); denom = a' * (En * En') * a; P_music(m) = 1 / (real(denom) + eps); end P_music_dB = 10*log10(P_music / max(P_music)); plot(theta, P_music_dB); xlabel('角度 (度)'); ylabel('MUSIC 谱 (dB)'); grid on;谱峰搜索可以用findpeaks，但要注意设置最小峰高和最小峰间距。如果两个信源靠得很近，比如间隔 5 度，MUSIC 的分辨力取决于快拍数和信噪比。K=200、SNR=10dB 时，5 度间隔能分开；K=50 时可能就糊成一个峰。这不是代码错，是物理限制。
3.3 信源数、快拍数、信噪比对 MUSIC 分辨力的影响
这三个参数是 MUSIC 的命门。信源数估错直接导致子空间划分错误；快拍数不够，协方差矩阵估计不准，特征分解不稳定；信噪比低，信号特征值和噪声特征值混在一起，分不开。案例前三章一般给的是理想条件，但你要自己扫参数。
我一般做一张表：固定 SNR=10dB，K 从 50 到 500，看两个间隔 5 度的信源能否分开。K=100 是临界点，K=200 稳定分辨。如果 SNR 降到 0dB，K 要 500 以上才稳。这些数字不是绝对，但趋势对你有用。
提示：MUSIC 谱峰高度不代表信号功率，只代表正交性程度。别拿谱峰高度去估功率，那是另一套方法。
4. 避坑与排查：前三章案例复现时最容易翻车的五个点
4.1 角度单位混用导致波束指向完全错误
现象：波束图主瓣指向和设定方向差几十度，或者左右镜像。原因：sin和sind混用，或者流形矩阵符号约定和扫描角度定义不一致。解决：统一用sind处理角度制，固定符号约定，画图前先验证单信源指向。
4.2 协方差矩阵秩亏导致 MVDR 求逆报错或结果异常
现象：R \ a_s返回 Inf 或 NaN，或者波束图完全乱掉。原因：快拍数 K 小于阵元数 N，R 不满秩。解决：增加 K 到 2N 以上，或加对角加载R = R + 0.01*trace(R)/N*eye(N)。
4.3 特征值排序遗漏导致噪声子空间取错
现象：MUSIC 谱没有明显峰值，或者峰值出现在错误角度。原因：eig输出顺序未排序，直接取后 N-P 列可能取到信号分量。解决：按特征值降序排列，再取后 N-P 列。
4.4 谱函数分母加 eps 位置不当引入虚假峰
现象：MUSIC 谱在无信源方向出现尖峰。原因：分母接近零时数值噪声被放大。解决：加eps到分母实数部分，或对谱函数做平滑。
4.5 扫描角度步长过粗导致峰值漏检
现象：信源在 20 度，但谱峰出现在 19.5 或 20.5 度。原因：theta步长 1 度太粗。解决：步长降到 0.1 度，或对谱函数做插值。
5. 进阶技巧：把前三章案例改成可复用的参数化脚本
前三章案例的代码通常是脚本式，参数写死。我习惯把它改成函数，输入 N、d、K、SNR、theta_s、theta_i，输出波束图和 MUSIC 谱。这样扫参数不用改代码。
function [P_cbf_dB, P_mvdr_dB, P_music_dB] = array_demo(N, d, K, SNR, theta_s, theta_i) lambda = 1; theta = -90:0.1:90; A = exp(-1j*2*pi*d*(0:N-1)'*sind(theta)/lambda); a_s = exp(-1j*2*pi*d*(0:N-1)'*sind(theta_s)/lambda); a_i = exp(-1j*2*pi*d*(0:N-1)'*sind(theta_i)/lambda); S = (randn(1,K)+1j*randn(1,K))/sqrt(2); I = (randn(1,K)+1j*randn(1,K))/sqrt(2) * 10^(SNR/20); X = a_s*S + a_i*I + (randn(N,K)+1j*randn(N,K))/sqrt(2); R = X*X'/K; % CBF idx = find(abs(theta-theta_s)<0.05,1); w_cbf = A(:,idx)/N; P_cbf = abs(w_cbf'*A).^2; % MVDR w_mvdr = (R\a_s)/(a_s'*(R\a_s)); P_mvdr = abs(w_mvdr'*A).^2; % MUSIC [V,D] = eig(R); [~,idx_sort] = sort(diag(D),'descend'); V = V(:,idx_sort); En = V(:,2:end); P_music = zeros(size(theta)); for m = 1:length(theta) a = A(:,m); P_music(m) = 1/(real(a'*(En*En')*a)+eps); end P_cbf_dB = 10*log10(P_cbf/max(P_cbf)); P_mvdr_dB = 10*log10(P_mvdr/max(P_mvdr)); P_music_dB = 10*log10(P_music/max(P_music)); end这个函数里 SNR 只作用于干扰源，信号源功率固定。你可以改成两个源都受 SNR 控制。调用时array_demo(16,0.5,200,10,20,-30)就能一次出三张图。我一般还会加一个输出：MVDR 在干扰方向的零陷深度，用来量化自适应能力。
验证方法：固定其他参数，只改 K，看 MVDR 零陷深度是否稳定。K=50 时零陷可能只有 -20dB，K=200 时能到 -40dB 以下。如果 K 增大零陷反而变浅，检查协方差矩阵估计是否用了X*X'/K而不是X*X'，后者会随 K 增大而幅度膨胀，导致求逆数值问题。
最后一个习惯：每次跑完案例，把关键参数和结果记一行日志，比如N=16, K=200, SNR=10, 零陷=-38dB, MUSIC 分辨=5度。攒多了你就知道边界在哪，别人问你参数怎么设，你直接翻日志，不用重新跑。希望帮到你。
本文还有配套的精品资源，点击获取
