简介针对稀疏圆阵的相干信源波达方向估计问题这份代码包实现了解相干求根多信号分类算法。算法在传统波束变换基础上引入相位校正与波束域误差补偿构造共轭对称导向矢量并结合前后向平均与求根多信号分类完成高精度估计适合阵列信号处理、稀疏圆阵波达方向估计方向的学习者与研究人员参考使用。压缩包共二十二个文件、约13KB内容精简主体为十个m脚本文件和八个备份脚本覆盖宽带信号波达方向估计、仿真实验等模块备份文件保留修改过程便于对照理解。目前已有两百二十八人学习资源虽小但代码结构清晰包含mydoa宽带信号示例与仿真脚本可直接加载运行也可根据自身阵列参数对算法做进一步改进与扩展。1. 圆阵相干快拍做 Root-MUSIC直接抄线阵代码必翻车一个很常见的需求手里有一套均匀圆阵的 16 通道快拍数据信号里有相干源算法点名要 Root-MUSIC 出 DOA。把线阵代码的导向矢量换掉就跑结果通常是求根多项式病态、根满天飞。原因有两层圆阵导向矢量是 cos(θ-φ_n) 的指数求和没有 ULA 的 Vandermonde 结构Root-MUSIC 赖以求根的等比数列前提不成立相干源又把协方差矩阵的秩打掉子空间分裂失效。工程上的标准做法是把链路拆成三步相位模式激励把圆阵变成虚拟 ULA前后向空间平滑恢复秩再在模式空间做多项式求根。快拍数少到平滑撑不住时换稀疏重构把角度估计写成网格字典上的稀疏向量恢复。下面按这条链路把每一步的原理、参数边界和可复现代码讲清楚。2. 圆阵 Root-MUSIC 的模式空间预变换把 UCA 拉直成虚拟 ULA2.1 圆阵导向矢量没有 Vandermonde 结构多项式求根无从谈起Root-MUSIC 把谱搜索变成多项式求根依赖的是均匀线阵ULA导向矢量 a(θ) [1, z, z², …, z^{N-1}]ᵀ 这种等比数列结构。把 z e^{-j2πd sinθ/λ} 代进 MUSIC 谱峰表达式分母整理后是 z 的多项式单位圆附近的根就是入射角估计。它的价值是连续输出角度不做网格搜索精度只受快拍数和信噪比限制。均匀圆阵UCA的情况完全不同。N 个阵元均匀分布在半径 r 的圆周上第 n 个阵元方位角 φ_n 2πn/N一个方向为 θ 的窄带信号在阵元上产生的相移是a_n(θ) e^{j2πr cos(θ-φ_n)}角度出现在 cos 的自变量里整体写不成某个变量 z 的幂次。硬凑多项式也能凑但阶数高、系数病态低信噪比下根会飘满整个 z 平面。所以不要在元素空间硬做求根通用做法是先用相位模式激励把圆阵映射成虚拟 ULA再跑标准 Root-MUSIC。这一步是圆阵 DOA 估计所有子空间类算法的公共前置。2.2 相位模式激励用波束形成权把圆阵响应拆成谐波这一步的数学基础是 Jacobi-Anger 展开e^{j2πr cos(θ-φ)} Σ_{m-∞}^{∞} j^m J_m(2πr) e^{jm(θ-φ)}它把每个阵元的响应拆成无数个角谐波 e^{jmθ} 的叠加第 m 个谐波的幅度是 j^m J_m(2πr)。如果能把第 m 个谐波单独取出来就相当于得到一个响应为 e^{jmθ} 的虚拟阵元——这正好是虚拟 ULA 第 m 个阵元的导向响应。取出谐波的动作就是一组波束形成每个阵元乘 e^{-jmφ_n} 再求和。因为 Σ_{n0}^{N-1} e^{j(m-m)φ_n} 在 m≠m 时为 0阵元分布的均匀性会把其它谐波抵消掉只剩目标模式。实际能用的模式只有 2M1 个。M 有两个硬上限离散阵元求和最多支持到第 N/2 阶模式再高出现角度混叠Bessel 函数 J_m(2πr) 在 m 超过 2πr 后迅速衰减到接近 0强行取高阶模式等于放大噪声。经验取值是 M min(⌊2πr⌋, ⌊N/2⌋−1)。2.3 MATLAB 实现从阵元坐标到模式空间的变换矩阵% 圆阵几何半径 r 按载波波长归一化 N 16; % 阵元数 r 0.8; % 半径 0.8 个波长 phi 2*pi*(0:N-1)/N; % 阵元方位角 % 模式阶数选择见 2.2 的边界讨论 M min(floor(2*pi*r), floor(N/2)-1); % 本例取 5 m (-M:M); % 模式索引 -5..5共 11 个模式 % 相位模式束形成矩阵 F(2M1) x N F zeros(2*M1, N); for k 1:2*M1 Jk besselj(abs(m(k)), 2*pi*r); % 该模式 Bessel 幅度 F(k,:) exp(-1j*m(k)*phi.) / (1j^abs(m(k)) * Jk) / sqrt(N); end % 快拍 X 是 N x T 复矩阵变换到模式空间 Y F * X; % 尺寸 (2M1) x T等效虚拟 ULA 输出F 每行做三件事exp(-1jm(k)phi) 是相位匹配把第 m 个角谐波从各阵元输出里选出来除以 1j^abs(m(k))Jk 是补上展开系数让虚拟阵列对 e^{jmθ} 的响应是等幅的1/sqrt(N) 让噪声功率尺度基本不变。省略 Bessel 幅度补偿的话虚拟阵列各阵元增益呈 J_m 递减序列Root-MUSIC 也能跑但高阶模式噪声被放大低 SNR 下估计偏差明显所以补偿不能省。矩阵化写法是 F (1/sqrt(N)) * diag(1 ./ (1j.^abs(m) .besselj(abs(m), 2pir))) * exp(1jmphi.)循环写法便于逐行核对。2.4 模式数 M 直接决定可分辨源数别往大里堆M虚拟阵元数 2M1非相干源可分辨上限典型半径范围254约 0.3λ0.35λ498约 0.6λ0.75λ51110约 0.8λ0.95λ71514约 1.1λ1.3λ可分辨源数上限是虚拟阵元数减一而虚拟阵元数等于 2M1。加了相干平滑之后还要再扣掉子阵造成的孔径损失所以表格右侧上限在相干场景基本达不到。另一个约束来自圆阵半径小圆阵r 0.3λ能用的模式太少虚拟孔径撑不起来r 过大时阵元互耦和模式混叠同时恶化。实践中圆阵半径设计一般落在 0.5λ1λM 按 ⌊2πr⌋ 取不要为了追求虚拟阵元数强行加大 M。3. 相干快拍下的协方差秩亏空间平滑次数与子阵长度怎么定3.1 相干信号让信号子空间分裂失效两个完全相干的信源例如同一发射信号的直达径和反射径快拍模型是 s₂(t) α·s₁(t)。此时信源协方差 R_s E[s sᴴ] 的秩只有 1接收协方差 R A R_s Aᴴ σ²I 的特征分解里信号子空间只剩一个大特征值。MUSIC 和 Root-MUSIC 都依赖信号子空间与噪声子空间的正交性秩亏直接导致两个源的特征向量混叠成一个谱峰变成一个求根也只剩一个有效根。这不是实现层面的 bug是数据模型层面的信息缺失。解法不是换谱函数而是先对协方差做空间平滑把秩恢复过来再进求根。平滑处理的时机有讲究一定要在模式空间变换之后做因为虚拟 ULA 是等间距线性结构滑窗子阵才有意义直接对圆阵元素空间做平滑子阵之间没有均匀相位递进效果很差。3.2 平滑次数 L、子阵长度 P 与可处理相干源数的关系前向空间平滑把虚拟 ULA 的 P₀ 2M1 个阵元划分成 L 个长度为 P P₀−L1 的滑窗子阵。同一个角度的信号在每个子阵上有不同的相位偏移子阵协方差平均之后相干源之间的固定相位关系被打散R_s 的秩恢复。两个硬条件L 不小于最大相干组的信源数且 P 不小于总信源数加一。后向平滑把数据反向共轭再平均一次等效把平滑次数翻倍L 可以放宽到相干源数的一半工程上保守起见仍按 L ≥ 相干源数取。平滑的代价是有效阵列长度从 P₀ 缩到 PMUSIC 分辨力随之下降。P 太小两个相邻角度的根在低 SNR 下分不开L 太大孔径损失又反噬精度。这几个参数是同一个 P₀ 下的跷跷板只能按目标相干组规模取最小满足值。虚拟阵元 P₀子阵长 P平滑次数 L可处理相干组平滑后最多可辨源数1184271510639158857实际系统里判断 L 够不够看特征值谱比看谱图更直接平滑后信号特征值仍然只比噪声大一两个量级说明 L 不够需要加大平滑次数。注意这里说的判断前提是信源数估计准确否则特征值个数本身就是错的。3.3 模式空间 FBSS 落地代码与求根多项式构造function [Rfb, P] fbss_virtual(Y, L) % Y : 模式空间快拍 (2M1) x T来自相位模式变换 % L : 前向平滑子阵数 P0 size(Y, 1); P P0 - L 1; % 子阵长度 Rf zeros(P, P); for l 1:L Yl Y(l:lP-1, :); % 第 l 个子阵快拍 Rf Rf (Yl * Yl) / size(Y, 2); end Rf Rf / L; J flip(eye(P)); % 反对角交换矩阵 Rfb 0.5 * (Rf J * conj(Rf) * J); % 叠加后向平滑 end代码逻辑上Y 的模式索引按 −M..M 顺序排列等同于一个阵元顺序固定的虚拟 ULAY(l:lP-1,:) 就是对这个虚拟阵列做滑窗切片后向平滑用 J 实现J * conj(Rf) * J 等价于把阵列反向、快拍共轭后的协方差与 Rf 平均后平滑次数翻倍。平滑完成后再做 Root-MUSIC 求根[V, D] eig(Rfb); [~, idx] sort(diag(D), descend); Ks 2; % 信源数估计值 En V(:, idx(Ks1:end)); % 噪声子空间 P x (P-Ks) % 构造 Root-MUSIC 多项式系数噪声特征向量自相关求和 c zeros(1, 2*P-1); for k 1:size(En, 2) c c conv(En(:,k), conj(flip(En(:,k)))).; end rt roots(c); % 2(P-1) 个根 [~, ii] sort(abs(abs(rt)-1), ascend); % 按离单位圆距离排序 z_sel rt(ii(1:Ks)); doa mod(-angle(z_sel) * 180/pi, 360); % 虚拟 ULA 的空间频率就是方位角多项式系数的构造依据是 p(z) z^{P-1} aᵀ(z^{-1}) EₙEₙᴴ a(z)展开后每个系数恰好对应某个噪声特征向量自相关在相应时延上的叠加所以用 conv 做自相关再求和即可。求出的 2(P-1) 个根关于单位圆成共轭倒数对离单位圆最近的那个半组就是信号根。角度换算用 doa −arg(z)因为虚拟 ULA 导向是 e^{-jmθ}复数变量 z e^{-jθ}。不要按 ULA 的 asin 公式换算虚拟阵列的空间频率是方位角本身不是 sinθ。注意L 每加 1P 就减 1。平滑次数覆盖相干组数之后不是越大越好低 SNR 下分辨力随孔径收缩先升后降。调参时盯着特征值谱看比盯谱峰图更可靠。4. 稀疏视角的 DOA 估计网格字典、分组稀疏与 l1-SVD4.1 相干和少快拍是稀疏模型的天然主场空间平滑救回秩代价是快拍数和孔径的双重需求。快拍降到个位数或者信噪比很低时协方差估计本身就不稳平滑救不回来。这时换一种建模方式不在协方差域做子空间分解直接对快拍本身做稀疏重构。把角度平面离散成 G 个网格点DOA 估计变成在过完备字典 Φ 上找一个稀疏向量 b使每个快拍 x(t) ≈ Φ·b(t)。D 个入射源对应 b 只有 D 个非零元素这是典型的稀疏向量恢复问题。它的优点在于模型里根本不出现协方差矩阵信源之间相不相关无所谓单快拍也能算。Φ 是 (2M1)×G 的复矩阵也叫冗余字典本质是一张有结构的稀疏矩阵每一列是候选角度的虚拟 ULA 导向矢量相邻网格列之间相关性经常超过 0.99。列相关性高意味着 OMP 这类贪心算法容易在相邻网格间反复横跳一般直接用 l1 范数正则或贝叶斯类方法前者实现简单后者免调正则系数但收敛慢。4.2 字典构造与 l1-SVD 的 CVX 求解多快拍不能逐个快拍恢复再叠加单快拍恢复的毛刺会在叠加时放大。常见做法是 l1-SVD对快拍矩阵做 SVD只保留前 Ks 个右奇异向量对应的投影把噪声维压掉再对降维结果做分组稀疏恢复。每个网格点在 Ks 个有效快拍上要么同时激活要么同时为零这个约束用行 2-范数求和实现。% 模式空间字典尺寸 (2M1) x G比元素空间字典小一个量级 G 360; % 网格数1° 步进 grid (0:G-1) * 2*pi / G; % 方位角网格 Phi exp(-1j * m * grid); % m 是模式索引列向量 % 快拍降维SVD 取前 Ks 个有效奇异方向 [~, ~, Vv] svd(Y, econ); Y_sv Y * Vv(:, 1:Ks); % (2M1) x Ks % 分组稀疏恢复 cvx_begin quiet variable B(G, Ks) complex minimize( sum(norms(B, 2, 2)) ) % 每行跨 Ks 列的 2-范数求和 subject to norm(Phi * B - Y_sv, fro) eta; % 噪声界约束 cvx_end pw sum(abs(B).^2, 2); % 每个网格角度的能量 [~, pk] sort(pw, descend); doa_sparse grid(pk(1:Ks)) * 180/pi;minimize 里的 sum(norms(B,2,2)) 是分组稀疏group sparsity的复数域写法先对每行求跨 Ks 个有效快拍的 2-范数再把所有行的范数求和。这个目标让解倾向于整行激活或整行为零避免同一角度在不同快拍上时有时无。约束里的 eta 是噪声界经验取值在 eta ≈ σ̂·√(2M1 Ks·√(2M1)) 附近σ̂ 用协方差矩阵最小特征值的均方根估计。eta 调太大解会从一个稀疏向量变成稠密向量全部网格都有能量调太小则过拟合噪声。这个稀疏与稠密的平衡是稀疏求解唯一的旋钮实际调试时以残差范数接近已知噪声能量为准。4.3 和 Root-MUSIC 怎么选精度、稀疏算力成本与网格失配对比项模式空间 Root-MUSIC FBSS稀疏字典重构相干信号处理必须平滑损失孔径天然免疫不依赖协方差秩最小快拍数数十个起协方差才稳定单快拍可运行角度精度无网格连续输出精度高受网格步进限制存在量化误差计算成本毫秒级稀疏算力开销小CVX 秒级G 和 Ks 增大后明显信源数估计错误根选错但可排查能量分布形态直接暴露问题选型判据很直接快拍多于 100、信噪比高于 5dB、只要方位角用 Root-MUSIC单快拍、强相干、低信噪比用稀疏重构。近两年 SubspaceNet 这类把学习型网络接进 DOA 任务的方案也有不少评估基准它们主要解决特征分解在阵列模型失配下的鲁棒性和求根、稀疏两条经典线路是互补关系替换的是子空间估计环节不是整个链路。稀疏方法最大的坑是网格失配真实角度落在两个网格点之间时稀疏解会把能量摊到相邻网格上出现双峰。工程上的修复是两轮细化第一轮用 1° 或 0.5° 粗网格定位峰值第二轮在峰值周围 ±1° 重建 0.01° 步进的小字典再解一次。第二轮字典只有几百列稀疏算力开销可接受精度能逼近 Root-MUSIC。小字典同样可以套用分组稀疏目标函数代码只需把 grid 换成细化区间。5. 快拍数据端到端验证从相干源生成到 RMSE 蒙特卡罗统计5.1 一条直接能跑的完整链路rng(1); N16; r0.8; phi2*pi*(0:N-1)/N; T200; th_true[40 130]*pi/180; % 两个真实入射角 Aexp(1j*2*pi*r*cos(th_true-phi)); % 圆阵流形 s(randn(1,T)1j*randn(1,T))/sqrt(2); S[s; s*exp(1j*pi/4)]; % s2 是 s1 的相干副本相位差 45° XA*S0.1*(randn(N,T)1j*randn(N,T)); % 噪声标准差 0.1 Mmin(floor(2*pi*r), floor(N/2)-1); m(-M:M); Fzeros(2*M1,N); for k1:2*M1 F(k,:)exp(-1j*m(k)*phi.)/(1j^abs(m(k))*besselj(abs(m(k)),2*pi*r))/sqrt(N); end YF*X; % 模式空间快拍 L3; P2*M1-L1; % 平滑次数 3子阵长 9 Rfzeros(P,P); for l1:L YlY(l:lP-1,:); RfRfYl*Yl/T; end RfRf/L; Jflip(eye(P)); Rfb0.5*(RfJ*conj(Rf)*J); [V,D]eig(Rfb); [~,idx]sort(diag(D),descend); EnV(:,idx(3:end)); % 已知两个源取第 3 个特征值之后 czeros(1,2*P-1); for k1:size(En,2) ccconv(En(:,k),conj(flip(En(:,k)))).; end rtroots(c); [~,ii]sort(abs(abs(rt)-1)); doamod(-angle(rt(ii(1:2)))*180/pi,360);这段代码跑完doa 应该在 [39.x, 130.x] 附近。相干副本的相位差 π/4 不影响秩亏性质只改变恢复后特征向量的相位分布这是检查实现是否正确的好抓手改相位差、改相干幅度 α估计角度不应有明显漂移。5.2 三个最容易翻车的参数平滑次数 L。L 小于相干组源数时特征值谱里第二个信号特征值会沉进噪声doa 只剩一个角。定位方法把 Rfb 的特征值打印出来看第 Ks1 个特征值和第 1 个噪声特征值是否在同一量级是则逐步加大 L。信源数 Ks。相干低 SNR 下 MDL/AIC 经常过估特征值谱没有明显膝部。经验做法是取谱膝部之前的个数然后 Root-MUSIC 多选 12 个候选根把落在雷达覆盖范围外或幅度异常的根剔除剩下的按离单位圆距离排序。半径归一化 r。代码里 r 是相对波长的值实际由物理半径除以波长得到。r 取错直接导致 M 和 Bessel 补偿全错表现为 F 条件数异常大、Y 能量比 X 小几个量级。上板前先用单已知源验证一遍变换矩阵的增益这一步能省掉大半调参时间。5.3 RMSE 校验与角度卷绕的细节for mc 1:200 X A*S 0.1*(randn(N,T)1j*randn(N,T)); % 每次只重采噪声 % ... 复制 5.1 的完整链路得到 doa ... err abs(doa - th_true*180/pi); err min(err, 360 - err); % 角度卷绕 errs(mc,:) err; end rms sqrt(mean(errs.^2, 1));角度卷绕这行不能省40° 和 400° 是同一个方向abs 直接减会在跨越 0°/360° 边界时把误差算成接近 360°RMSE 直接失真。评估这套链路的设计保证能力不要看单次谱峰图要看蒙特卡罗下的分辨概率和 RMSE 曲线横轴是信噪比从 0dB 到 20dB每条曲线固定快拍数。调试时如果 RMSE 在高 SNR 下不降反升优先怀疑根筛选逻辑其次才是平滑参数。把稀疏重构的粗估计角度作为 Root-MUSIC 求根前的先验约束只保留距粗估计 2° 以内的候选根圆阵相干场景的整套 DOA 链路才算真正闭环。本文还有配套的精品资源点击获取