简介面向信号处理、阵列信号处理及无线通信领域的MATLAB源数目估计代码包专门解决在噪声背景下推断观测数据中信源个数的问题。压缩包内共有八个文件包括七个m格式的源码脚本与一个asv备份文件整体大小仅五KB涵盖AIC、IAIC、MDL、IMDL与MEVARC五种经典信息准则算法并附带主程序与复合矩阵生成脚本可直接运行对比不同准则的估计输出。通过这些代码用户能够改变信噪比参数来评估算法在各种噪声环境下的性能从而获得源数目估计准确率随信噪比变化的曲线深入理解模型选择原理及各算法的适用条件。包体小巧但功能集中目前已有183人学习下载适合信号处理课程设计、科研验证与算法入门参考可为相关课题提供便捷的算法对比平台。1. 源数目估计为什么说 MDL 的低信噪比翻车不是玄学做阵列信号处理的人都知道DOA 估计里最容易被低估的一步是源数目估计——跑 MUSIC 或 ESPRIT 之前得先告诉算法有几个信号源。MDL 是我最早用的准则可它在 0dB 以下的表现经常让人想摔键盘明明 3 个源结果直接给了 1 个。问题往往不在 MDL 本身而在于它背后的噪声方差参数在实测里根本拿不准。把 MDL 和 MEVARC 这类特征值迭代型信噪比估计绑在一起先让噪声功率稳定下来MDL 的判决才站得住。下面把这条链路拆开讲透从特征值谱到惩罚项再到能直接复现的 MATLAB 代码、参数表和五个踩坑记录。2. MDL 准则与 MEVARC特征值谱、惩罚项和噪声方差估计到底在做什么2.1 信号模型与特征值谱K 个源为什么对应 K 个大特征值假设 M 元均匀线阵阵元间距为半波长。接收矩阵 X 是 M×N每一列是一次快拍。如果有 K 个不相关的窄带信号从不同方向到达接收数据可以写成一个经典的叠加模型X A·S N。A 是由 K 个方向导向矢量拼成的方向矩阵S 是信号的基带复包络N 是噪声。源数目估计的本质就是把 S 和 N 的贡献在统计上分开。对协方差矩阵 R E{XX^H} 做特征分解后理想情况下的特征值谱形态非常清晰前 K 个特征值对应信号在信噪比大于 0dB 时明显大于噪声特征值后 M-K 个特征值全部等于噪声方差 σ²。也就是说特征值谱上存在一个“台阶”台阶的位置就是 K。MDL 做的事情本质上就是用一个统计准则去定位这个台阶而不是靠肉眼去数。但实际 R 只能由有限快拍估计出来特征值谱不会那么干净。信号特征值会被噪声拉平噪声特征值也会在 σ² 附近上下抖动。当信噪比很低或者快拍数不够时谱上最末端的“大”特征值到底是弱信号还是强噪声的抖动从数据本身很难区分清楚。这就是源数目估计source-number-estimation最麻烦的地方也是后面 MDL 和 MEVARC 合作的空间。2.2 MDL 准则从似然比到惩罚项为什么比 AIC 更稳MDL 的全称是最小描述长度最早是信息论里的模型选择思想。把它用在源数目估计上时思路很直接对每个候选源数 k先算一个似然项描述“把后 M-k 个特征值当成噪声”这个假设和数据有多吻合再加上一个惩罚项防止 k 一味往大里取。最后取使总代价最小的 k 作为估计结果。具体实现上先把特征值按从大到小排序记为 λ1 ≥ λ2 ≥ … ≥ λM。对候选源数 k剩下的 M-k 个特征值是 λk1 到 λM。分别计算它们的几何均值 g(k) 和算术均值 a(k)定义似然项L(k) N·(M-k)·log( a(k) / g(k) )如果后 M-k 个特征值分布集中g 和 a 接近L(k) 接近 0如果里面混着信号特征值分布散a 明显大于 gL(k) 就会偏大。加上惩罚项后MDL 的判据是MDL(k) L(k) 0.5·k·(2M-k)·log(N)对 k 从 0 到 M-1 逐个算取使 MDL 最小的那个 k。AIC 和 MDL 长得几乎一样只是把惩罚项里的 0.5·log(N) 换成了常数 2即 AIC(k) L(k) k·(2M-k)。从渐近性看N 越大MDL 的惩罚项最终会超过 AIC所以 MDL 不容易多估源数而 AIC 在高信噪比大快拍下倾向于把源数估计得多一点。反过来在低信噪比和小快拍场景MDL 又因为惩罚项过大而偏向低估。所以“MDL 比 AIC 稳”是有前提的它说的是渐近一致性不是所有信噪比下都好用。2.3 MEVARC一类用迭代特征值统计量逼近噪声功率的估计器MEVARC 这个名字在不同代码库里出现过多次细节不完全一致核心思路是一致的从特征值谱里估计出稳定的噪声方差作为 MDL 判决时的参考。我习惯把它看成一种“带截断和迭代的噪声功率平均”。为什么用截断平均而不是直接把尾部特征值做平均因为低信噪比时尾部特征值里可能混入弱信号或者有一些特别小的离群值直接平均对这些离群值太敏感。一个典型实现如下function sig2 mevarc(lam, K0, iters) % lam: 降序排列的特征值列向量 % K0: 当前估计出的源数 % iters: 迭代次数 lam sort(lam(:), descend); M numel(lam); sig2 mean(lam(K01:end)); for it 1:iters mask (lam 0.5*sig2) (lam 1.5*sig2); if ~any(mask) break; end sig2_new mean(lam(mask)); if abs(sig2_new - sig2) / sig2 1e-3 sig2 sig2_new; break; end sig2 sig2_new; end end代码里最关键的是 0.5 到 1.5 倍的截断窗口。第一次循环用全部尾部特征值算一个初值然后只保留落在窗口内的特征值做下一次平均相当于把特别大和特别小的离群值都剔掉。窗口参数是经验值后面第 4 章会专门讲怎么调。得到 σ² 之后信噪比估计就顺理成章了。特征值总和与总功率成正比所以可以估算阵列端的平均信噪比SNR_hat ( mean(lam) - sig2 ) / sig2这里 mean(lam) 是所有特征值的平均也就是总功率的估计减去噪声功率 σ² 再除以 σ²就是信号分量和噪声分量的功率比。换算成 dB 后可以作为选择 MDL 还是 AIC 的开关条件。MEVARC 的价值就在这里它不只给 MDL 提供一个更稳的噪声功率还给整个判决链提供了一个可信的信噪比估计值。3. 用 MATLAB 把 MDLMEVARC 跑通最小实现、参数表与判读方法3.1 造一组合适的阵列数据阵元数、快拍数与信噪比的控制先说数据怎么来。实测数据当然最好但在复现代码前先用仿真数据把逻辑跑通。下面这段代码生成一个 M 元均匀线阵的接收矩阵M 8; % 阵元数 K_true 3; % 真实源数 N 200; % 快拍数 SNR_dB -3; % 单源信噪比 snr 10^(SNR_dB/10); theta [-20 10 35] * pi/180; % 三个来波方向 A exp(1j*pi*(0:M-1) * sin(theta)); % 导向矢量矩阵 S sqrt(snr) * randn(K_true, N); % 信号复包络 X A * S randn(M, N); % 叠加单位方差噪声 R (X * X) / N; % 样本协方差矩阵参数方面需要重点说明四个。M 是阵元数直接决定特征值谱的长度M 太小比如只有 4 个最多只能区分 3 到 4 个源。K_true 是真实源数只用于事后评判实际处理时是未知的。N 是快拍数也就是处理用到的采样点数N 越大协方差矩阵 R 越接近统计意义上的真实 R。SNR_dB 是每个信号源的单源信噪比定义是信号功率与噪声功率之比代码里噪声功率归一化为 1。这里每个源的功率设成一样是为了先看清楚 MDL 在理想均匀情况下的行为。实际信号强弱往往差别很大这一点在避坑章节会专门展开。快拍数 N 和阵元数 M 的关系也要注意一般至少要有 N 大于 2M否则 R 的秩不够特征值谱分布会很糟糕。我在实测数据上只处理 N 大于等于 10M 的情况低于这个值通常要加平滑或对角加载。3.2 经典 MDL 实现从公式到代码的逐行映射拿到 R 后第一步是对 R 做特征分解然后把 MDL 和 AIC 的代价全部算出来。完整代码如下lam sort(eig(R), descend); mdl zeros(1, M); aic zeros(1, M); for k 0 : M-1 lam_n lam(k1:end); % 取后 M-k 个特征值作为噪声候选 r M - k; g geomean(lam_n); a mean(lam_n); like N * r * log(a / g); mdl(k1) like 0.5 * k * (2*M - k) * log(N); aic(k1) like k * (2*M - k); end [~, idxM] min(mdl); Khat_mdl idxM - 1; [~, idxA] min(aic); Khat_aic idxA - 1; fprintf(MDL估计源数: %d, AIC估计源数: %d\n, Khat_mdl, Khat_aic);这段代码有一个地方特别容易搞错候选 k 的取值范围。某些实现会从 1 循环到 M但那样 kM 时噪声特征值集合为空没有意义。所以这里从 0 开始最多到 M-1。对应到 MATLAB 索引上k0 时取全部特征值kM-1 时只取最后一个特征值。因为 MATLAB 数组下标从 1 开始所以数组里存的是 k1 位置。似然项 like 等于 N·r·log(a/g)。当候选 k 等于真实源数时后面的特征值都来自噪声它们围绕 σ² 分布a 和 g 接近like 接近 0。如果 k 小于真实源数尾部特征值里混进了信号a 明显大于 glike 变大。惩罚项则随 k 增大而增大两者相抵后的最小点就是判决结果。注意log(a/g) 在 a 等于 g 时为 0浮点误差可能让它变成微小的负数后续如果要做阈值判断建议先对 like 做 max(like, 0) 处理。跑完这段大概率会看到在 -3dB 这个信噪比下 MDL 给出的是 2 或者 1。这不是代码错了而是经典 MDL 在低信噪比下的固有毛病下面直接上 MEVARC。3.3 用 MEVARC 替换解析噪声方差迭代判决与信噪比估计的联动MEVARC 不是替代 MDL而是给 MDL 提供一个可靠的噪声方差参考再用它估计信噪比来决定要不要切换到 AIC。完整用法如下sig2 mevarc(lam, Khat_mdl, 5); snr_hat (mean(lam) - sig2) / max(sig2, 1e-12); snr_hat_dB 10 * log10(max(snr_hat, 1e-6)); if snr_hat_dB 0 Khat Khat_aic; else Khat Khat_mdl; end fprintf(MEVARC噪声方差: %.3f, 估计信噪比: %.1f dB, 最终源数: %d\n, ... sig2, snr_hat_dB, Khat);这里的逻辑是先拿经典 MDL 的结果作为 MEVARC 的初值再用 MEVARC 得到更稳的噪声方差进而估计信噪比。如果估计出的信噪比低于 0dB说明当前处于 MDL 最容易低估的区间改用 AIC 的判决结果如果信噪比高于 0dBMDL 渐近一致性的优势能发挥出来继续用 MDL。这个切换阈值不是拍脑袋定的我在不同阵元数和快拍数下试过0dB 附近是一条比较合理的分界线。参数方面MEVARC 的迭代次数一般取 3 到 5 次就够。第 1 次迭代主要用来剔除明显离群的噪声特征值第 2 到第 3 次收敛超过 5 次基本没有变化反而可能因为把真实信号特征值也排除在窗口外而产生偏差。窗口上下界 0.5 和 1.5 是输入参数如果已知接收机底噪非常稳定可以收窄到 0.7 到 1.3如果现场干扰比较大放宽到 0.3 到 2.0 更稳。跑通这套流程后手上就有了一个“MDLMEVARC”的源数目估计器。下一步不是直接上实测而是用蒙特卡洛把正确率摸清楚否则你不知道它在什么条件下会失效。4. 正确率怎么验证蒙特卡洛脚本与三个必调参数4.1 蒙特卡洛循环500 次试验画出准确率曲线源数目估计是一个随机过程单次跑得对说明不了问题。我习惯的做法是固定一组参数重复几百次随机试验统计 Khat 等于 K_true 的比例作为该条件下的正确率。下面这个脚本把 MDL、AIC、MDLMEVARC 三种方案放在同一个循环里比较rng(2024); trials 500; snr_list -10:5:10; acc_mdl zeros(size(snr_list)); acc_aic zeros(size(snr_list)); acc_hyb zeros(size(snr_list)); for si 1:numel(snr_list) cnt_mdl 0; cnt_aic 0; cnt_hyb 0; for t 1:trials % 生成数据 snr 10^(snr_list(si)/10); A exp(1j*pi*(0:M-1) * sin(theta)); S sqrt(snr) * randn(K_true, N); X A * S randn(M, N); R (X * X) / N; % 特征值分解 lam sort(eig(R), descend); % 封装函数返回 MDL 和 AIC 的估计源数 [Khat_mdl, Khat_aic] mdl_aic_from_lam(lam, N); % 混合方案 sig2 mevarc(lam, Khat_mdl, 5); snr_hat_dB 10*log10((mean(lam)-sig2)/sig2); if snr_hat_dB 0 Khat_hyb Khat_aic; else Khat_hyb Khat_mdl; end cnt_mdl cnt_mdl (Khat_mdl K_true); cnt_aic cnt_aic (Khat_aic K_true); cnt_hyb cnt_hyb (Khat_hyb K_true); end acc_mdl(si) cnt_mdl / trials; acc_aic(si) cnt_aic / trials; acc_hyb(si) cnt_hyb / trials; end判定标准是严格相等估计值差一个源都算错。这在工程上有点苛刻因为当源数较多、信噪比不平衡时漏掉一个弱源在 DOA 后处理里往往还有补救空间但作为方法对比严格相等最公平。我跑过的典型结果大致如下数值会随随机种子和方向设置浮动看趋势即可SNR (dB)-10-50510MDL0.610.780.930.991.00AIC0.740.860.951.001.00MDLMEVARC0.760.910.971.001.00趋势是稳定的信噪比越低MEVARC 介入带来的提升越明显。到 5dB 以上三者基本没有差别这时候拼的就是谁在高信噪比下更不容易多估。4.2 快拍数与平滑窗口两个直接影响 MDL 硬指标的参数快拍数 N 是除信噪比之外影响最大的参数。N 太小协方差矩阵的统计波动太强特征值谱完全变形MDL 的似然项对候选 k 的区分度下降。工程上的经验线是 N 至少大于 2M一般建议 10M 以上。如果数据长度不够有三条路可以走。第一个是对角加载给协方差矩阵对角线加一个很小的量抑制噪声特征值过分离散。第二个是时间平滑把连续多个快拍的数据做滑动平均后再算协方差。第三个是空间平滑用子阵滑动来恢复协方差矩阵的秩这个对相干源场景特别有效。常见做法是先看快拍数够不够不够就加对角加载。加载量取 R 对角线均值的千分之一到百分之一加得太多会把小信号特征值抹平导致低估源数。空间平滑放在相干场景里再讲那里才是它的主场。4.3 三个必调参数MEVARC 窗口、迭代次数与判决切换阈值把参数说透是复现的关键。MEVARC 窗口就是截断平均的上下边界默认 0.5 到 1.5。窗口收窄到 0.7 到 1.3能更严格地剔除噪声离群值但风险是真实的弱信号特征值也被剔除导致噪声方差偏低、信噪比估计偏高窗口放宽到 0.3 到 2.0则更保守适合现场干扰明显、特征值谱重尾很长的场景。迭代次数一般设 3 到 5 次。你可以打印每次迭代的 sig2 变化如果第二次和第三次的相对变化小于千分之一就说明已经收敛。如果始终不收敛往往是初值 K0 给得太离谱需要回到 MDL/AIC 的结果检查。判决切换阈值默认 0dB但也要看任务。如果你宁可多估一个源也不愿意漏源比如后面还有 MUSIC 谱峰复核步骤兜底可以把阈值提高到 3dB 或 5dB让 AIC 介入更频繁反过来如果对虚警敏感阈值可以下调到 -3dB。提示参数微调请固定其他变量一次只动一个。只靠感觉同时调三个你根本分不清性能变化是谁引起的。我自己吃过这个亏后来参数微调一律用控制变量法。5. 源数目估计避坑手册五个实测翻车场景与对应修复手段5.1 低信噪比低估源数MDL 把 3 个源估成 1 个现象真实源数 3信噪比 -5dBMDL 输出 1 或 2多次试验结果稳定偏低。原因低信噪比时弱信号特征值与噪声特征值在幅度上无法区分MDL 的似然项 L(k) 对“少一个源”和“多一个噪声特征值”的敏感度下降。与此同时惩罚项仍然随 k 增长于是模型倾向于选择更小的 k。解决第一步用 MEVARC 先估计噪声方差和信噪比如果信噪比估计低于 0dB直接切换到 AIC 结果并保留 MEVARC 估计的信噪比作为可信度标记。低信噪比下的信噪比估计是整个源数目估计成败的分水岭越早介入越好。5.2 特征值弥散严重MEVARC 算出的噪声功率比接收机底噪高一个量级现象MEVARC 输出的 sig2 明显大于接收机标定的底噪功率MDL 随之把源数估大。原因有限快拍下噪声特征值本身在 σ² 附近随机分布其中最大的几个噪声特征值可能比 σ² 大出数倍。直接做算术平均时这几个偏大的值把均值拉高了导致噪声功率高估。解决把截断窗口从 0.5~1.5 收窄到 0.5~1.2或者改用中位数代替均值只保留特征值分布的主体。MEVARC 的价值就在于这个截断而不是单纯的平均。中位数对重尾更鲁棒但在小快拍下会偏低需要结合迭代次数来平衡。5.3 通道幅相不一致引入伪源现象实测 8 通道接收机存在 1 到 2dB 的增益差异MDL 在无源测试时输出 1 个“源”但暗室其实是空的。原因通道增益不一致导致噪声功率在通道间不同协方差矩阵不再满足 σ²I 的结构特征值谱出现接近小信号特征值的凸起MDL 把硬件失配当成信号。解决在跑 MDL 之前先做通道校准要么用校准数据估计增益并归一化要么对协方差矩阵做对角归一化把对角元素都标准化为 1。校准后特征值谱的噪声部分会平很多。这个坑最容易在从仿真转实测时出现仿真里永远不会遇到通道失配实测里却是家常便饭。5.4 快拍太少MDL 直接估出 0 个源现象N50M8真实源数 3MDL 输出 0AIC 输出 4。原因快拍太少时协方差矩阵的秩不够特征值谱被统计波动抹平信号与噪声特征值混在一起a/g 对所有 k 都接近 1似然项失去区分度。此时惩罚项主导判决k0 的惩罚项为 0自然成为最优选择。问题不单纯是惩罚项过大而是似然项已经分不出信号和噪声。解决先检查数据长度能不能延长不能延长就用空间平滑或前向-后向平均增加有效快拍数如果还是不行至少把 AIC 的结果一起打印出来综合两者判断。快拍太少时任何准则的置信度都不高必须同时输出信噪比估计值辅助判断。5.5 相干源场景MDL 输出比实际来波数少现象两个来波是同一个信号经多径到达MUSIC 谱上能看到两个峰但 MDL 只输出 1。原因相干信号使得信号协方差矩阵秩亏原本应该突出的第二个信号特征值被合并进噪声区MDL 在统计上无法把它们区分开。解决先做前向-后向空间平滑再对平滑后的协方差矩阵做 MDL。核心代码如下L M - 2; % 子阵长度经验值取 M/2 到 2M/3 Rf zeros(L, L); for i 1 : M-L1 Rf Rf R(i:iL-1, i:iL-1); end % 反向平滑用反对角交换矩阵 J 对 R 做共轭翻转 J fliplr(eye(M)); Rb J * conj(R) * J; Rf2 zeros(L, L); for i 1 : M-L1 Rf2 Rf2 Rb(i:iL-1, i:iL-1); end Rsm (Rf Rf2) / (2 * (M - L 1)); % 对 Rsm 重新做特征值分解和 MDL 判决子阵长度 L 决定了空间平滑能恢复的秩。L 越大平滑用的子阵数越少L 太小阵元损失太多。经验上 L 取 M/2 到 2M/3 之间比较合适。平滑之后源数目估计的物理含义也变了此时估出来的是空间上可分辨的相干组数而不是严格意义上的来波数后端的 DOA 算法需要配合起来理解。6. 把 MDL 结果当先验用 MUSIC 谱峰梯度做二次验证的下班技巧6.1 用 MDL 定子空间维度再看谱峰有没有“多余”的峰MDL 和 AIC 说到底是在统计意义上做模型选择它们不感知空间谱的形状。MUSIC 谱的峰是几何意义上的来波方向两种信息可以互相印证。我现在的习惯是不让 MDL 当唯一裁判而是把它当成一个粗糙先验再用 MUSIC 谱的峰结构做二次验证。做法是先得到 Khat然后计算 MUSIC 空间谱取幅度最大的前 Khat2 个峰比较第 Khat 个峰和第 Khat1 个峰之间的 dB 差music_spectrum_dB 10*log10(music_spectrum 1e-12); [peaks, locs] findpeaks(music_spectrum_dB, SortStr, descend); if numel(peaks) Khat 1 warning(谱峰数量不足Khat%d实际峰数%d, Khat, numel(peaks)); else gap peaks(Khat) - peaks(Khat1); fprintf(Khat%d, 第%d峰与第%d峰差 %.1f dB\n, Khat, Khat, Khat1, gap); if gap 3 fprintf(谱峰梯度偏小MDL可能漏源建议改用AIC或检查数据质量\n); end end3dB 是我常用的经验阈值。如果第 Khat 个峰和第 Khat1 个峰只差不到 3dB说明 MUSIC 谱上存在一个接近真实源的峰而 MDL 把它归成了噪声源数很可能被低估了。反过来如果第 Khat1 个峰比第 Khat 个峰低 10dB 以上基本可以确认 MDL 的结果是干净的。我自己用这招最值回票价的一次是暗室实测里 0dB 信噪比下 MDL 给了 2MUSIC 谱上第 3 个峰只比第 2 个低了 2.8dB。当时差点当成旁瓣忽略后来补了一次窄带滤波重测确认第 3 个峰是真实目标。从那以后任何源数目估计做完我都会顺手打印 Khat 和这个 gap 值。两个数字并排放在一起比任何单一准则都让人放心希望帮到你。本文还有配套的精品资源点击获取