Wax 和 Kailath [6] 提出了一种统计方法，该方法基于通常用于模型选择的 AIC [7] 和 MDL [8] 准则来解决信源数量检测问题。由于其确定过程不需要主观判断，当 AIC 或 MDL 准则达到最小时，信源数量便可以被自然地确定。Zhang 等人 [9] 对 AIC 和 MDL 的统计性能进行了研究。随后，Wong 等人 [10] 利用由 AIC 和 MDL 的改进准则形成的新对数似然函数，实现了性能的提升。Wax 和 Ziskind [11] 以及 Wu 和 Fuhrmann [12] 则采用新颖的参数化方法来处理此问题，他们分别将估计的特征值和微分残差作为广义准则中的变量。Wax-Ziskind 方法是 Wax 和 Kailath [6] 所提概念的扩展，而在 [12] 中，定义了一组新的微分残差，它使用了具有不同假定信源数量的对数似然函数。最近，Wax [13] 针对检测问题提出了一种新颖的解决方案，该方案同时适用于相干和非相干信号。该解决方案基于将 Rissanen 的 MDL 原理 [8] 应用于随机信号模型。此外，Wax [14] 还提出了一种新技术，其中 MDL 被应用于空间非白噪声中的相干和非相干信号。

总的来说，AIC 和 MDL 方法，包括它们的改进版本，仍然是估计信源数量最广泛使用的方法。由于这些方法基于零均值、统计独立的髙斯随机噪声的假设，它们可被归类为依赖于模型且基于特征值的信源数量检测方法。而对于经常出现在小样本或真实数据样本中的非高斯和非白噪声模型，上述算法无法正确检测信源的数量。此外，它们没有利用特征向量中的有价值信息。

在本文中，我们提出了两种借助格什戈林半径 [1] 进行信源数量检测的方法。在第二节中，我们简要回顾了基于特征值的信源数量估计器，例如 AIC 和 MDL 方法。在第三节中，我们使用酉变换 [1] 介绍了格什戈林圆盘定理 [16]。我们建议，变换后的协方差矩阵的格什戈林半径可用于信源数量检测。在第四节中，我们探讨了使用格什戈林半径进行信源数量检测的似然法和启发式方法。似然法将现有的检测准则与格什戈林半径有效结合，以改进其对信源数量的确定。启发式方法是根据正交投影的概念发展而来的，无需信源和噪声模型的先验知识，便可以正确地检测信源数量。在第五节中，我们提供了几个仿真和真实数据示例，以评估知名估计器和我们提出的估计器的检测性能。

II. EIGENVALUE-BASED SOURCE NUMBER DETECTORs

III. UNITARY TRANSFORMED GERSCHGORIN RAD

在讨论我们提出的方法之前，我们将简要回顾格什戈林圆盘定理，该定理主要用于根据格什戈林中心（Gerschgorin centers）和格什戈林半径来估计矩阵特征值的位置。

对于一个元素表示为 $a_{ij}$ 的实数或复数 $L\times L$ 矩阵 $\mathbf{A}$，格什戈林找到了一种从 $a_{ij}$ 的值来估计特征值位置的简单方法。第 $i$ 行向量中除第 $i$ 个元素外的所有元素的模长之和，定义为
$$r_i=\sum _{\begin{array}{c}j=1\\ j\ne i\end{array}}^L|a_{ij}|\text{\hspace{1em}(12)}$$

其中 $i=1,\dots ,L$。第 $i$ 个圆盘 ${\mathbf{O}}_i$ 则被定义为复平面中与 $a_{ii}$ 的距离至多为 $r_i$ 的点的集合。换句话说，${\mathbf{O}}_i$ 代表具有以下性质的复数 $z$ 的集合
$$r_i\ge |z-a_{ii}|.\text{\hspace{1em}(13)}$$

格什戈林 [16] 已经证明，矩阵 $\mathbf{A}$ 的特征值包含在所有圆盘 ${\mathbf{O}}_i$ 的并集之内。这些圆盘 ${\mathbf{O}}_i(i=1,2,\dots ,L)$ 被称为矩阵 $\mathbf{A}$ 的格什戈林圆盘。中心 $c_i=a_{ii}$ 和半径 $r_i$ 分别被称为格什戈林中心和格什戈林半径。此外，如果 $k$ 个圆盘的集合与其他圆盘相隔离，那么在该集合中恰好存在 $k$ 个矩阵 $\mathbf{A}$ 的特征值 [16]。

在格什戈林圆盘定理的一个例子中，一个矩阵由下式给出
$$\mathbf{A}=\left(\begin{array}{ccc}10& i& -i\\ -i& 5& 0.5i\\ i& -0.5i& 4\end{array}\right)\text{\hspace{1em}(14)}$$

其格什戈林圆盘可由 (13) 式得到，分别为 O 1 = { z : 2 ≥ ∣ z − 10 ∣ } \mathbf{O}_1 = \{z : 2 \ge |z - 10|\} O1​={z:2≥∣z−10∣}， O 2 = { z : 1.5 ≥ ∣ z − 5 ∣ } \mathbf{O}_2 = \{z : 1.5 \ge |z - 5|\} O2​={z:1.5≥∣z−5∣}，以及 O 3 = { z : 1.5 ≥ ∣ z − 4 ∣ } \mathbf{O}_3 = \{z : 1.5 \ge |z - 4|\} O3​={z:1.5≥∣z−4∣}。由于 ${\mathbf{O}}_2\cup {\mathbf{O}}_3$ 与 ${\mathbf{O}}_1$ 不相交，因此矩阵 $\mathbf{A}$ 有两个特征值必须位于 ${\mathbf{O}}_2\cup {\mathbf{O}}_3$ 内，另一个则必须位于 ${\mathbf{O}}_1$ 内。(14) 式中矩阵 $\mathbf{A}$ 的三个特征值被确定为 10.35, 5.03 和 3.62，并且满足上述论点。结果如图1所示。根据上述格什戈林圆盘定理，(3) 式中协方差矩阵的每个特征值都包含在复平面中一系列圆盘的并集之内，这些圆盘以 $c_i=c_{ii}$ 为中心，半径为 $r_i={\sum }_{j\ne i}|c_{ij}|$，其中 $c_{ij}$ 是 $\mathbf{C}$ 的一个元素。由于协方差矩阵 $\mathbf{C}$ 是厄米特矩阵，其特征值和格什戈林中心都是实数。

当计算出实的格什戈林中心和 $L$ 个对应的格什戈林半径后，我们可以在实轴上定位可能的特征值。然而，由于大的格什戈林半径和相似的中心，原始协方差矩阵的格什戈林圆盘对于信源数量检测仍然没有帮助。例如，在一个仿真协方差矩阵的案例中，传感器数量为六（即 $L=6$），两个信源（即 $M=2$）是不相关的，并从 -10° 和 10° 入射（即 DOA = [-10° 10°]）。信噪比均为 2 dB（即 SNR = [2 2] dB），样本数量选择为 $N=100$。图2说明了该矩阵的格什戈林圆盘紧密重叠，其特征值分布在零到六之间的一个很大范围内。在不改变真实特征值的情况下，需要对协方差矩阵进行适当的酉变换，以便有效地利用格什戈林圆盘进行信源数量检测。

B. Unitary Transfomation

考虑到先前的限制，我们应该设计一种酉变换来旋转协方差矩阵，使其格什戈林圆盘可以形成两个不同的信号和噪声集合。具有较大格什戈林半径的信源集合将恰好包含 $M$ 个最大的信号特征值，而具有小格什戈林半径的噪声集合将包含剩余的噪声特征值。换句话说，我们必须选择一种能够使噪声格什戈林圆盘尽可能小，并尽可能远离信源格什戈林圆盘的酉变换，从而可以通过对圆盘大小进行分类来检测信源数量。

协方差矩阵 $\mathbf{C}$ 可以通过特征分解表示为
$$\mathbf{C}={\mathbf{U}}_d\mathbf{D}{\mathbf{U}}_d^H\text{\hspace{1em}(15)}$$

其中 ${\mathbf{U}}_d$ 是一个由 $\mathbf{C}$ 的特征向量构成的 $L\times L$ 酉矩阵，
$${\mathbf{U}}_d=[{\mathbf{q}}_1{\mathbf{q}}_2\dots {\mathbf{q}}_M\dots {\mathbf{q}}_L]\text{\hspace{1em}(16)}$$

并且 $\mathbf{D}$ 是一个对角矩阵，其对角元素由 $\mathbf{C}$ 的特征值组成，如下所示
$$\mathbf{D}=\text{diag}({\lambda }_1,{\lambda }_2,\dots ,{\lambda }_M,\dots ,{\lambda }_L).\text{\hspace{1em}(17)}$$

如果执行 ${\mathbf{U}}_d$ 的酉变换，变换后的协方差矩阵将变成一个对角矩阵，即 ${\mathbf{C}}_d={\mathbf{U}}_d^H\mathbf{C}{\mathbf{U}}_d=\mathbf{D}$。因此，变换后的协方差矩阵 ${\mathbf{C}}_d$ 的格什戈林半径全部变为零。其对角元素，即格什戈林中心，变成了特征值的精确位置。变换后的协方差矩阵提供了特征值的精确位置信息。对于真实的协方差矩阵，${\mathbf{C}}_d$ 的噪声格什戈林圆盘形成一个半径为零的单一集合。然而，对于有限数据样本的情况，${\mathbf{C}}_d$ 的 $L$ 个格什戈林圆盘是分离且清晰的。在这种情况下，应使用众所周知的基于特征值的方法 [5]–[15] 来检测信源数量。

在我们提出新的变换之前，首先将协方差矩阵分区如下
$$\mathbf{C}=\left(\begin{array}{cccc}{c}_{11}& c_{12}& \cdots & c_{1L}\\ c_{21}& c_{22}& \cdots & c_{2L}\\ ⋮& & \ddots & \\ c_{L1}& c_{L2}& \cdots & c_{LL}\end{array}\right)=\left(\begin{array}{cc}{\mathbf{C}}_1& \mathbf{c}\\ {\mathbf{c}}^H& c_{LL}\end{array}\right)\text{\hspace{1em}(18)}$$

其中 ${\mathbf{C}}_1$ 是 $\mathbf{C}$ 的一个 $(L-1)\times (L-1)$ 的前导主子矩阵，通过删除 $\mathbf{C}$ 的最后一行和最后一列得到。在物理上，这可以被看作是移除了第 $L$ 个传感器。因此，${\mathbf{C}}_1$ 成为剩余 $(L-1)$ 个传感器的降维协方差矩阵。值得注意的是，(18) 式中的 $\mathbf{c}$ 可以表示为

$$\mathbf{c}=[c_{1L}{c}_{2L}\dots c_{(L-1)L}{]}^T=[{\mathbf{b}}_1{\mathbf{b}}_2\dots {\mathbf{b}}_{L-1}{]}^T{\mathbf{C}}_s{\mathbf{b}}_L^{\ast }={\mathbf{A}}_1{\mathbf{C}}_s{\mathbf{b}}_L^{\ast }\text{\hspace{1em}(19)}$$

其中 ${\mathbf{b}}_i$ 如 (4) 式所定义，且 ${\mathbf{A}}_1=[{\mathbf{b}}_1{\mathbf{b}}_2\dots {\mathbf{b}}_{L-1}{]}^T$ 是 ${\mathbf{C}}_1$ 的方向矩阵。降维的协方差矩阵 ${\mathbf{C}}_1$ 也可以通过其特征结构分解为
$${\mathbf{C}}_1={\mathbf{U}}_1{\mathbf{D}}_1{\mathbf{U}}_1^H\text{\hspace{1em}(20)}$$

其中 ${\mathbf{U}}_1$ 是一个由 ${\mathbf{C}}_1$ 的特征向量构成的 $(L-1)\times (L-1)$ 酉矩阵，
$${\mathbf{U}}_1=[{{\mathbf{q}}^{\prime }}_1{{\mathbf{q}}^{\prime }}_2\dots {{\mathbf{q}}^{\prime }}_M\dots {{\mathbf{q}}^{\prime }}_{L-1}]\text{\hspace{1em}(21)}$$

并且 ${\mathbf{D}}_1$ 是由相应特征值构成的对角矩阵，如下所示
$${\mathbf{D}}_1=\text{diag}({\lambda }_1^{\prime },{\lambda }_2^{\prime },\dots ,{\lambda }_M^{\prime },\dots ,{\lambda }_{L-1}^{\prime }).\text{\hspace{1em}(22)}$$

与 (6) 式类似，特征值按降序排列如下
$${\lambda }_1^{\prime }\ge {\lambda }_2^{\prime }\ge \cdots \ge {\lambda }_M^{\prime }\ge {\lambda }_{M+1}^{\prime }\ge \cdots \ge {\lambda }_{L-1}^{\prime }\text{\hspace{1em}(23)}$$

由于 (22) 式中的 ${\lambda }_i^{\prime }$ 是 $\mathbf{C}$ 的前导主子矩阵（leading principal submatrix）的特征值，如 [17] 所示，它们在 (6) 和 (23) 中的特征值满足如下所示的交错特性
$${\lambda }_1\ge {\lambda }_1^{\prime }\ge {\lambda }_2\ge \cdots \ge {\lambda }_M^{\prime }\ge {\lambda }_{M+1}\ge {\lambda }_{M+1}^{\prime }\ge \cdots \ge {\lambda }_{L-1}\ge {\lambda }_{L-1}^{\prime }\ge {\lambda }_L.\text{\hspace{1em}(24)}$$

对于真实的协方差矩阵，可以证明 ${\lambda }_{M+1}^{\prime }=\cdots ={\lambda }_{L-1}^{\prime }={\sigma }_n^2$。现在，我们构造一个重要的 $L\times L$ 酉变换矩阵 $\mathbf{U}$（$\mathbf{U}{\mathbf{U}}^H=\mathbf{I}$），表示为
$$\mathbf{U}=\left(\begin{array}{cc}{\mathbf{U}}_1& 0\\ 0^T& 1\end{array}\right).\text{\hspace{1em}(25)}$$

变换后的协方差矩阵变为
$$\begin{array}{rl}\mathbf{S}={\mathbf{U}}^H\mathbf{CU}& =\left(\begin{array}{cc}{\mathbf{U}}_1^H{\mathbf{C}}_1{\mathbf{U}}_1& {\mathbf{U}}_1^H\mathbf{c}\\ {\mathbf{c}}^H{\mathbf{U}}_1& c_{LL}\end{array}\right)\\ & =\left(\begin{array}{cc}{\mathbf{D}}_1& {\mathbf{U}}_1^H\mathbf{c}\\ {\mathbf{c}}^H{\mathbf{U}}_1& c_{LL}\end{array}\right)\\ & =\left(\begin{array}{cccccc}{\lambda }_1^{\prime }& 0& 0& \cdots & 0& {\rho }_1\\ 0& {\lambda }_2^{\prime }& 0& \cdots & 0& {\rho }_2\\ 0& 0& {\lambda }_3^{\prime }& \cdots & 0& {\rho }_3\\ ⋮& & & \ddots & & ⋮\\ 0& 0& 0& \cdots & {\lambda }_{L-1}^{\prime }& {\rho }_{L-1}\\ {\rho }_1^{\ast }& {\rho }_2^{\ast }& {\rho }_3^{\ast }& \cdots & {\rho }_{L-1}^{\ast }& c_{LL}\end{array}\right)\end{array}\text{\hspace{1em}(26a)(26b)(26c)}$$

其中
$${\rho }_i={{\mathbf{q}}^{\prime }}_i^H\mathbf{c}={{\mathbf{q}}^{\prime }}_i^H{\mathbf{A}}_1{\mathbf{C}}_s{\mathbf{b}}_L^{\ast }\text{\hspace{1em}(27)}$$

其中 $i=1,2,\dots ,L-1$。

如前一节所示，格什戈林圆盘定理可用于估计 $\mathbf{S}$ 的特征值。很明显，前 $(L-1)$ 个格什戈林圆盘（即 ${\mathbf{O}}_1,{\mathbf{O}}_2,\dots ,{\mathbf{O}}_{L-1}$）拥有如下的格什戈林半径
$$r_i=|{\rho }_i|=|{{\mathbf{q}}^{\prime }}_i^H{\mathbf{A}}_1{\mathbf{C}}_s{\mathbf{b}}_L^{\ast }|=|{{\mathbf{q}}^{\prime }}_i^H\mathbf{c}|\text{\hspace{1em}(28)}$$

其中 $i=1,2,\dots ,L-1$。有必要验证，当 $i=(M+1),(M+2),\dots ,(L-1)$ 时，所有的 ${\rho }_i$ 值都等于零，这是因为噪声特征向量 ${{\mathbf{q}}^{\prime }}_i$ 与 ${\mathbf{C}}_1$ 的方向矩阵 ${\mathbf{A}}_1$ 正交。此外，因为信号特征向量 ${{\mathbf{q}}^{\prime }}_i$ 与 ${\mathbf{A}}_1$ 不正交，且 ${\mathbf{C}}_s$ 是满秩的，所以当 $i=1,2,\dots ,M$ 时，(28) 式中的 $r_i$ 是非零的。因此，变换后，真实的M信源协方差矩阵变为
$$\mathbf{S}=\left(\begin{array}{ccccccc}{\lambda }_1^{\prime }& & 0& \cdots & & 0& {\rho }_1\\ & \ddots & & & & & ⋮\\ 0& & {\lambda }_M^{\prime }& & & 0& {\rho }_M\\ ⋮& & & {\sigma }_n^2& & & 0\\ & & & & \ddots & & ⋮\\ 0& & 0& & & {\sigma }_n^2& 0\\ {\rho }_1^{\ast }& \dots & {\rho }_M^{\ast }& 0& \dots & 0& c_{LL}\end{array}\right)\text{\hspace{1em}(29)}$$

由于 $\mathbf{S}$ 是 $\mathbf{C}$ 的一个酉变换矩阵，它们将共享相同的特征值。前 $(L-1)$ 个格什戈林圆盘的集合 ${\mathbf{O}}_i$ 包含位于 $c_i={\lambda }_i^{\prime }$（即 (22) 式中 ${\mathbf{C}}_1$ 的特征值）的格什戈林中心和对应的格什戈林半径 $r_i=|{\rho }_i|,i=1,2,\dots ,(L-1)$。

• 半径为零的圆盘（即 ${\mathbf{O}}_{M+1},{\mathbf{O}}_{M+2},\dots ,{\mathbf{O}}_{L-1}$）被视为噪声格什戈林圆盘的集合。
• 剩余的包含非零半径和大圆心的圆盘（即 ${\mathbf{O}}_1,{\mathbf{O}}_2,\dots ,{\mathbf{O}}_M$）被认为是信源格什戈林圆盘。

因此，在无限样本的情况下，我们可以通过计算非零格什戈林半径的数量来确定信源数量。此外，我们也可以使用 ${\mathbf{C}}_1$ 的 $(L-1)$ 个特征值来确定信源数量。

例如，如图2中使用的相同估计协方差矩阵，通过建议的酉变换进行变换（如 (26) 所示），其格什戈林圆盘以（中心，半径）对的形式变为 {2.55, 0.69}, {2.38, 0.71}, {1.10, 0.13}, {1.02, 0.09} 和 {0.92, 0.20}。结果如图3所示。现在，显著的是格什戈林圆盘形成了两个分离的集合。信源集合包含具有大半径的 ${\mathbf{O}}_1\cup {\mathbf{O}}_2$，噪声集合则包含具有小半径的 ${\mathbf{O}}_3\cup {\mathbf{O}}_4\cup {\mathbf{O}}_5$。

V. SOURCE NUMBER ESTIMULATORS USING GERSCHGORIN RADII

A. Likelihood Approach-Gerschgorin Likelihood Estimator (GLE）

从酉变换后的数据快照 $\mathbf{z}(k)=\mathbf{Ux}(k)$ 出发，我们可以根据变换后的协方差矩阵 $\mathbf{S}=E[\mathbf{z}(k)\mathbf{z}(k{)}^H]$ 来构建用于信源数量检测的似然函数。使用由 Wax 和 Kailath [6] 建立的类似程序，我们得到对数似然函数 ${\Psi }_L(k)$ 如下
$${\Psi }_L(k)=-N\log [\det ({\mathbf{S}}^{(k)})]-\text{tr}[({\mathbf{S}}^{(k)}{)}^{-1}\hat{\mathbf{S}}]\text{\hspace{1em}(30)}$$
其中 $\hat{\mathbf{S}}$ 和 ${\mathbf{S}}^{(k)}$ 是变换后的样本协方差矩阵，分别定义为 $\hat{\mathbf{S}}={\mathbf{U}}^H\hat{\mathbf{C}}\mathbf{U}$ 和 ${\mathbf{S}}^{(k)}={\mathbf{U}}^H{\mathbf{C}}^{(k)}\mathbf{U}$。通过观察 (29)，${\mathbf{S}}^{(k)}$ 可以通过设定以下条件来近似
$${\lambda }_{k+1}^{\prime }={\lambda }_{k+2}^{\prime }\cdots ={\lambda }_{L-1}^{\prime }={\sigma }_n^2={\beta }_{L-1}^{\prime }(k)=\frac{1}{L-k-1}\sum _{i=k+1}^{L-1}{\lambda }_i^{\prime }\text{\hspace{1em}(31)}$$
并且
$${\rho }_{k+1}={\rho }_{k+2}=\cdots ={\rho }_{L-1}=0.\text{\hspace{1em}(32)}$$

因此，变换后的协方差矩阵可以被分为四个子矩阵
$${\mathbf{S}}^{(k)}=\left(\begin{array}{cc}{\mathbf{S}}_{11}^{(k)}& {\mathbf{S}}_{12}^{(k)}\\ {\mathbf{S}}_{21}^{(k)}& {\mathbf{S}}_{22}^{(k)}\end{array}\right)\text{\hspace{1em}(33)}$$
其中 ${\mathbf{S}}_{11}^{(k)}=\text{diag}({\lambda }_1^{\prime },{\lambda }_2^{\prime },\dots ,{\lambda }_k^{\prime },{\sigma }_n^2,\dots ,{\sigma }_n^2)$，${\mathbf{S}}_{12}^{(k)}={\mathbf{S}}_{21}^{(k)H}=[{\rho }_1,{\rho }_2,\dots ,{\rho }_k,0,\dots ,0{]}^T$，并且 ${\mathbf{S}}_{22}^{(k)}=c_{LL}$。使用 [18] 中的恒等式
$$\det (\mathbf{S})=\det \left[\begin{array}{cc}{\mathbf{S}}_{11}^{(k)}& {\mathbf{S}}_{12}^{(k)}\\ {\mathbf{S}}_{21}^{(k)}& {\mathbf{S}}_{22}^{(k)}\end{array}\right]=\det ({\mathbf{S}}_{11}^{(k)})\det ({\mathbf{S}}_{22}^{(k)}-{\mathbf{S}}_{21}^{(k)}({\mathbf{S}}_{11}^{(k)}{)}^{-1}{\mathbf{S}}_{12}^{(k)})\text{\hspace{1em}(34)}$$
可以很容易地验证 $\det ({\mathbf{S}}^{(k)})$ 可表示为
$$\det ({\mathbf{S}}^{(k)})\approx \left(\prod _{i=1}^k{\lambda }_i^{\prime }\prod _{i=k+1}^{L-1}{\sigma }_n^2\right)\left(c_{LL}-\sum _{i=1}^k\frac{r_i^2}{{\lambda }_i^{\prime }}\right).\text{\hspace{1em}(35)}$$
由于 ${\mathbf{S}}^{(k)}={\mathbf{U}}^H{\mathbf{C}}^{(k)}\mathbf{U}$ 且 $\hat{\mathbf{S}}={\mathbf{U}}^H\hat{\mathbf{C}}\mathbf{U}$，我们也可以证明
$$\text{tr}[({\mathbf{S}}^{(k)}{)}^{-1}\hat{\mathbf{S}}]=\text{tr}[{\mathbf{U}}^H({\mathbf{C}}^{(k)}{)}^{-1}\hat{\mathbf{C}}\mathbf{U}]=\text{tr}[({\mathbf{C}}^{(k)}{)}^{-1}\hat{\mathbf{C}}]\approx L\text{\hspace{1em}(36)}$$
通过应用 [6] 中使用的相同论证。根据 (35) 和 (36) 并忽略与 k 无关的常数项，新的对数似然函数最终变为
$$\begin{array}{rl}{\Psi }_{L,G}(k)& \approx -N\log \left(\left(\prod _{i=1}^k{\lambda }_i^{\prime }\prod _{i=k+1}^{L-1}{\sigma }_n^2\right)\left(c_{LL}-\sum _{i=1}^k\frac{r_i^2}{{\lambda }_i^{\prime }}\right)\right)\\ & \approx N(L-1-k)\log \left(\frac{{\alpha }_{L-1}^{\prime }(k)}{{\beta }_{L-1}^{\prime }(k)}\right)-N\log \left(c_{LL}-\sum _{i=1}^k\frac{r_i^2}{{\lambda }_i^{\prime }}\right)\end{array}\text{\hspace{1em}(37)}$$

其中 ${\alpha }_{L-1}^{\prime }(k)=({\prod }_{i=k+1}^{L-1}{\lambda }_i^{\prime }{)}^{1/(L-k-1)}$。我们可以看到 (37) 式右侧的第一项是在 (8) 式中陈述的 $L-1$ 传感器阵列的对数似然函数。第一项描述了噪声元素算术平均值与其几何平均值之比的单调函数。除了这个仅与特征值相关的对数似然函数外，(37) 式右侧的第二项是由格什戈林半径通过其对应特征值加权所贡献的。这个格什戈林项是信号子空间中元素之间距离的一种度量。很明显，对于有限数据样本，第二项也是一个递减函数。当 $\hat{M}=M$ 时，格什戈林似然函数的值会急剧下降，从而提高了检测性能。应用惩罚函数后，格什戈林似然估计器 (GLE) 现在变为
$$\text{GLE}(k)={\Psi }_{L,G}(k)+P(N,k,L)\text{\hspace{1em}(38)}$$

用于信源数量检测。由于所提出的似然函数使用了格什戈林半径和 ${\mathbf{S}}_1$ 的信号子空间分量的信息，可以很容易地验证 (37) 式中信号子空间的自由参数数量为 $k^2$ 个分量 [13] 以及格什戈林半径 $(r_1,\dots ,r_k)$ 的 $k$ 个自由参数。因此，如果我们为 (38) 式选择 $P(N,k,L)=k^2+k$，则 GLE 准则就变成了格什戈林 AIC (GAIC)。如果 $P(N,k,L)=\frac{1}{2}(k^2+k)\log N$，则 GLE 准则就变成了格什戈林 MDL (GMDL)。信源数量由使 GAIC 或 GMDL 准则最小化的值确定。

通过 [19] 中建议的图形证明，我们可以表明这两种准则都是一致的。在图4中，我们发现当 $\hat{M}=M$ 时，GAIC 从一个正值下降到一个负值。因此，GAIC 是强一致的，而 AIC 则倾向于高估信源数量，因为格什戈林半径信息在 $\hat{M}=M$ 时增加了似然函数的深度。我们也可以证明 GMDL 准则也是强一致的。

值得注意的是，所有的似然准则都是在高斯过程和独立数据观测的假设下推导出来的。对于非高斯和非白噪声的情况，需要一种启发式的方法。

B. Heuristic Approach-Gerschgorin Disk Estimator (GDE）

通过观察 (27)，我们得知 ${\rho }_i={{\mathbf{q}}^{\prime }}_i^H\mathbf{c}={{\mathbf{q}}^{\prime }}_i^H{\mathbf{A}}_1{\mathbf{C}}_s{\mathbf{b}}_L^{\ast }$。格什戈林半径可以解释为从 ${{\mathbf{q}}^{\prime }}_i$ 到 $\mathbf{c}$ 的投影的幅度。根据柯西-施瓦茨不等式，我们知道
$$r_i=|{\rho }_i|=|{{\mathbf{q}}^{\prime }}_i^H{\mathbf{A}}_1{\mathbf{C}}_s{\mathbf{b}}_L^{\ast }|\le |{{\mathbf{q}}^{\prime }}_i^H{\mathbf{A}}_1|\cdot |{\mathbf{C}}_s{\mathbf{b}}_L^{\ast }|=\lambda |{{\mathbf{q}}^{\prime }}_i^H{\mathbf{A}}_1|\text{\hspace{1em}(39)}$$

其中 $\lambda =|{\mathbf{C}}_s{\mathbf{b}}_L^{\ast }|$ 与 $i$ 无关。第 $i$ 个格什戈林半径 $r_i$ 实际上是由 ${{\mathbf{q}}^{\prime }}_i^H{\mathbf{A}}_1$ 的幅度决定的。如果 ${{\mathbf{q}}^{\prime }}_i$ 是噪声特征向量，那么作为噪声格什戈林半径的第 $i$ 个格什戈林半径就会显著减小。如果 ${{\mathbf{q}}^{\prime }}_i$ 是信号特征向量，那么作为信号格什戈林半径的第 $i$ 个格什戈林半径将不为零。${{\mathbf{q}}^{\prime }}_i^H{\mathbf{A}}_1$ 的幅度几乎与噪声模型以及信号和噪声功率无关。因此，格什戈林半径仅取决于 ${{\mathbf{q}}^{\prime }}_i$ 的准确估计。当快照数量较大时，估计的协方差矩阵更接近真实的协方差矩阵，从而导致更精确的特征向量估计。

为了研究快照数量和格什戈林半径之间的影响，我们使用了一个八传感器均匀线性阵列来接收两个等强度信号（SNR = 10 dB），其DOA分别为10°和20°（DOA = [10° 20°]）。图5显示了在不同快照数量下，变换后协方差矩阵的噪声（-）和信源（+）格什戈林半径。图中显示，随着快照数量的增加，变换后协方差矩阵的噪声和信源格什戈林半径变得更容易区分。当 $N$ 趋于无穷大时，它们是完全分离的。因此，只要选择了合适的阈值，就可以仅通过格什戈林半径来检测信源数量。

从图5可以看出，阈值必须根据不同的快照数量进行调整。因此，我们定义了一个启发式的决策规则，如 [1] 所示

$$\text{GDE}(k)=r_k-\frac{D(N)}{L-1}\sum _{i=1}^{L-1}{r}_i\text{\hspace{1em}(40)}$$

其中 $k$ 是闭区间 $[1,L-2]$ 内的整数。可调因子 $D(N)$ 可以是一个在 $N$ 增加时值介于1和0之间的非增函数。随着 $N$ 变大，由于噪声格什戈林半径会随着快照数量的增加而减小，阈值也会变小。如果从 $k=1$ 开始计算 GDE(k)，当达到第一个非正值的 GDE(k) 时，信源数量就被确定为 $k-1$（即 $M=k-1$）。这是因为低于可调阈值的半径值将被视为噪声集合。当传感器数量和信源的信噪比非常低时，噪声格什戈林半径的值会接近大多数信源格什戈林半径。因此，上述 GDE 规则可能会产生低估的问题。

两种提出的方法都只涉及 $(L-1)$ 维的特征分解，并且比需要 $L$ 维特征分解的著名检测器计算复杂度更低。表I显示了当通过QR迭代执行特征分解时，AIC、GAIC和GDE方法的乘法次数。

C. Improved Decision Rules-Rotation and Averaging

因为格什戈林半径 $r_i=|{\rho }_i|$ 可以被视为 $\mathbf{C}$ 的第 $L$ 列在降维协方差矩阵的第 $i$ 个特征向量上的投影，即 $r_i=|{\rho }_i|=|q_i^{H}c|$，所以变换后的半径应该取决于 $q_i$ 的精确估计以及第 $i$ 个传感器和第 $L$ 个传感器 $c$ 之间的互相关函数。降维协方差矩阵的计算可以作用于另一组行列而非最后一组。因此，存在 $L$ 个不同的酉变换。例如，我们可以按如下方式旋转协方差矩阵：

$${\mathbf{C}}^{(m+1)}={\mathbf{EC}}^{(m)}{\mathbf{E}}^{-1},\text{for }m=0,1,\cdots ,(L-1)\text{\hspace{1em}(41)}$$

其中 ${\mathbf{C}}^{(0)}=\mathbf{C}$ 且 $\mathbf{E}$ 是一个选择如下的基本交换矩阵：

$$\mathbf{E}=\left(\begin{array}{cccccc}0& 1& 0& \cdots & 0& \\ 0& 0& 1& \cdots & 0& \\ 0& 0& 0& 1& \cdots & 0\\ ⋮& ⋮& ⋮& \ddots & ⋮& \\ 0& 0& 0& 0& \cdots & 1\\ 1& 0& 0& 0& \cdots & 0\end{array}\right)\text{\hspace{1em}(42)}$$

因为 ${\mathbf{E}}^T={\mathbf{E}}^{-1}$，所以该交换矩阵也是酉矩阵。我们可以将这组新的酉变换矩阵视为：
$${\mathbf{U}}^{(m+1)}={\mathbf{E}}^H{\mathbf{U}}^{(m)},\text{for }m=0,1,\cdots ,(L-1)\text{\hspace{1em}(43)}$$

其中 ${\mathbf{U}}^{(0)}=\mathbf{U}$，如公式 (25) 所述。多个酉矩阵的乘积仍然是酉矩阵。一旦获得了新的变换后的特征结构和相应的格什戈林半径，就可以重新计算上述两个检测准则。所有旋转和变换后的协方差矩阵的格什戈林半径都可以用来检测信源数量。由于 ${\mathbf{U}}^{(m)}$ 是直接通过估计 ${\mathbf{C}}^{(m)}$ 得到的，因此该检测准则与第一个准则具有相同的期望，但统计特性不同。为了提高检测性能，改进的格什戈林估计器按以下步骤执行：

1. 令 $m=0$。通过使用建议的协方差矩阵 ${\mathbf{C}}^{(0)}=\mathbf{C}$ 的酉变换，计算格什戈林半径 $r_i^{(0)}=|{\rho }_i^{(0)}|$，$i=1,2,\cdots ,L$。

2. 计算检测准则
   $${\text{GLE}}^{(m)}(k)=-N(L-1-k)\log \left[\frac{{\alpha }_{L-1}^{\prime }(k)}{{\beta }_{L-1}^{\prime }(k)}\right]+N\log \left(c_{LL}-\sum _{i=1}^k\frac{r_i^2}{{\lambda }_i}\right)+P(N,k,L)\text{\hspace{1em}(44)}$$
   或
   $${\text{GDE}}^{(m)}(k)=r_k-\frac{D(N)}{L-1}\sum _{i=1}^{L-1}{r}_i\text{\hspace{1em}(45)}$$

3. 通过 ${\mathbf{C}}^{(m+1)}={\mathbf{EC}}^{(m)}{\mathbf{E}}^{-1}$ 旋转当前的协方差矩阵 $\mathbf{C}$，然后令 $m=m+1$。

4. 重复步骤 2，以找到对于 $m=1,2,\cdots ,L-1$ 的 ${\text{GLE}}^{(m)}(k)$ 或 ${\text{GDE}}^{(m)}(k)$。

5. 计算平均准则
   $$\text{MGLE}(k)=\frac{1}{L}\sum _{m=0}^{L-1}{\text{GLE}}^{(m)}(k)\text{\hspace{1em}(46)}$$
   对于似然法
   或
   $$\text{MGDE}(k)=\frac{1}{L}\sum _{m=0}^{L-1}{\text{GDE}}^{(m)}(k)\text{\hspace{1em}(47)}$$
   对于启发式方法

6. 当准则函数 $\text{MGLE}(k)$ 在 k ∈ { 0 , L − 1 } k \in \{0, L-1\} k∈{0,L−1} 范围内取最小值时，将信源数量确定为 $k$；或者当 $\text{GDE}(k)$ 在 k ∈ { 1 , L − 1 } k \in \{1, L-1\} k∈{1,L−1} 范围内达到第一个非正值时，将信源数量确定为 $k-1$。

MGLE 准则在 $P(N,k,L)=k^2+k$ 时被称为修正的 GAIC (MGAIC)，在 $P(N,k,L)=\frac{1}{2}(k^2+k)\log N$ 时被称为修正的 GMDL (MGMDL)。公式 (46) 或 (47) 中的修正检测准则，通过对 $L$ 次旋转结果进行平均，有助于提高准则的检测性能。然而，MGLE 准则所需的计算复杂度大约是 GLE 的 $L$ 倍。

V. SIMULATION RESULT

我们现在提供仿真结果，以展示所提出方法的性能。在第一个实验中，我们使用一个由八个各向同性传感器组成的均匀线性阵列，传感器间距为半个波长，并带有加性不相关白噪声。所提出的准则、AIC 和 MDL 被用于检测两个不相关的信源，其信噪比（SNR）均为 0 dB，分别从 0° 和 10° 入射。经过 200 次蒙特卡洛运行后，我们使用不同数量的快照计算了它们的相对错误检测频率。

这是我实际仿真出来的结果，代码我已经放到文章的最后了。可以看到和文章中的结果还是有不小差距的。MDL算法依然还是性能最好的算法。

图6描绘了 AIC, MDL, GAIC, GMDL, MGAIC, 和 MGMDL 在概率方面的错误检测性能。GAIC, GMDL, MGAIC, MGMDL 的性能优于 AIC 和 MDL。可以看出，格什戈林 AIC (GAIC) 准则能够对信源数量产生一致的估计，而格什戈林 MDL (GMDL) 准则在小或中等数据样本下沒有低估信源数量的倾向。

此外，通过对 $L$ 组格什戈林半径进行平均，我们能够提高 GAIC 和 GMDL 的检测性能。例如，在仿真数据的情况下，MGAIC 以 $L$ 倍的复杂度为代价，实现了比 GAIC 更准确的信源数量检测。通过在 MATLAB 中评估浮点运算（flop）次数，AIC, GAIC, 和 MGAIC 方法分别需要 45222, 29107 和 240416 次浮点运算。因此，我们可以得出结论，AIC 准则比 GAIC 准则的计算成本更高。MGAIC 在所有方法中需要最高的计算负载。类似地，MDL, GMDL, 和 MGMDL 也有相同的结论。

第二个实验使用了 48 组真实数据，这些数据由明尼苏达大学 Kaveh 实验室提供的经过校准的实验传感器阵列（$L=8$）采集 [15]。其中包括 30 组数据，用于两个在 -4° 和 -6° 的等功率信源（具有相同SNR）；以及 18 组数据，用于两个在 1° 和 3° 的等功率信源（具有不同SNR）。每组数据包含 64 个数据样本。通过对这 48 组数据的检测结果进行平均，我们可以计算出这四种方法在使用 8-64 范围内的不同快照数量时的错误概率，如图7所示。

我们发现，AIC 和 MDL 方法确定的检测信源数量接近七个，而当使用超过 26 个样本的数据序列时，GDE 估计器能以 100% 的准确率确定确切的信源数量。此外，使用 MGDE 准则，我们在使用超过 9 个数据样本的情况下，对所有 48 组数据都能正确估计信源数量。该实验的结果表明，在实测实验数据的情况下，所提出的 GDE 和 MGDE 方法的性能优于 AIC, MDL, 和 GAIC。

最后，第三个实验的条件与第一个实验相同，除了 SNR = [1 3] dB, DOA = [20° 30°]，并且使用了非白噪声。在图8中，通过使用自回归系数为 0.9 的 AR(1) 色高斯噪声，我们发现 AIC, MDL, GAIC, 和 GMDL 估计器都无法检测出信源数量。当快照数量 $N$ 增加时，GDE 和 MGDE 方法的检测性能变得更加准确。然而，增加数据样本并不能改善 AIC, MDL, GAIC, 和 GMDL 的性能。所提出的 GDE 和 MGDE 规则，在使用变换后协方差矩阵的格什戈林半径时，能够正确检测 AR 模型的信源数量。这个实验证明了即使在缺少噪声模型信息的情况下，GDE 和 MGDE 依然可以使用。

VI. Simulation Code

为了方便大家更好的理解，本人编写了一套仿真代码，各位可以参考。

% MATLAB Code for comparing GAIC, GMDL, AIC, and MDL algorithms
% Based on the paper: "Source number estimators using transformed Gerschgorin radii"
% by H.T. Wu, J.F. Yang, and F.K. Chen, IEEE T-SP, 1995.

clear; clc; close all;

%% 1. 仿真参数设置 (Simulation Parameters)
% --- 固定参数 (与论文图6场景类似) ---
L = 8;                  % 传感器数量 (Number of sensors)
M = 3;                  % 真实信源数量 (True number of sources)
DOA = [0, 10, 30];          % 信号入射角度 (Directions of Arrival in degrees)
SNR_dB = [5, 7, 10];        % 信噪比 (Signal-to-Noise Ratios in dB)
d = 0.5;                % 阵元间距 (半波长) (Sensor spacing in half-wavelengths)
num_runs = 200;         % 每个N值的蒙特卡洛仿真次数

% --- 变化参数 (Varying Parameter) ---
N_values = 10:10:200;   % 快照数 N 的变化范围
prob_correct_gaic = zeros(size(N_values));
prob_correct_gmdl = zeros(size(N_values));
prob_correct_aic = zeros(size(N_values));
prob_correct_mdl = zeros(size(N_values));

%% --- 主仿真循环 ---
for idx = 1:length(N_values)
    N = N_values(idx);
    correct_gaic = 0;
    correct_gmdl = 0;
    correct_aic = 0;
    correct_mdl = 0;
    
    % --- 蒙特卡洛循环 ---
    for run = 1:num_runs
        %% 2. 生成阵列接收信号 (Generate Array Signals)
        A = exp(-1j * 2 * pi * d * (0:L-1)' * sind(DOA));
        s_powers = 10.^(SNR_dB / 10);
        s = sqrt(s_powers') .* (randn(M, N) + 1j * randn(M, N)) / sqrt(2);
        noise_w = (randn(L, N) + 1j * randn(L, N)) / sqrt(2); % 高斯白噪声
        X = A * s + noise_w;

        %% 3. 实现所有算法 (Implement All Algorithms)
        C_hat = (1/N) * (X * X');

        % --- 首先计算传统AIC/MDL所需的全尺寸特征值 ---
        lambda_full = sort(eig(C_hat), 'descend');
        
        % --- GAIC/GMDL所需的酉变换和Gerschgorin半径 ---
        C1_hat = C_hat(1:L-1, 1:L-1);
        c_hat = C_hat(1:L-1, L);
        c_LL = C_hat(L, L);
        [U1_hat, D1_hat] = eig(C1_hat);
        [lambda_prime, sorted_idx] = sort(diag(D1_hat), 'descend');
        U1_hat_sorted = U1_hat(:, sorted_idx);
        r = zeros(L-1, 1);
        for i = 1:(L-1)
            r(i) = abs(U1_hat_sorted(:, i)' * c_hat);
        end
        
        % --- 为每个k计算所有四种判据 ---
        GAIC_values = zeros(1, L-1);
        GMDL_values = zeros(1, L-1);
        AIC_values = zeros(1, L);
        MDL_values = zeros(1, L);

        for k = 0:(L-1) % k是假设的信源数量
            
            % --- 计算传统AIC和MDL (k from 0 to L-1) ---
            lambda_noise_full = lambda_full(k+1 : L);
            alpha_full = geo_mean(lambda_noise_full);
            beta_full = mean(lambda_noise_full);
            
            if alpha_full < 1e-10 || beta_full < 1e-10
                log_likelihood_full = Inf;
            else
                log_likelihood_full = N * (L - k) * log(alpha_full / beta_full);
            end
            
            penalty_aic_trad = k * (2*L - k);
            penalty_mdl_trad = 0.5 * k * (2*L - k) * log(N);
            
            AIC_values(k+1) = -log_likelihood_full + penalty_aic_trad;
            MDL_values(k+1) = -log_likelihood_full + penalty_mdl_trad;

            % --- 计算GAIC和GMDL (k from 0 to L-2) ---
            if k <= L-2
                lambda_noise_prime = lambda_prime(k+1 : L-1);
                alpha_prime = geo_mean(lambda_noise_prime);
                beta_prime = mean(lambda_noise_prime);
                
                if alpha_prime < 1e-10 || beta_prime < 1e-10
                    log_term1 = Inf;
                else
                    log_term1 = N * (L - 1 - k) * log(alpha_prime / beta_prime);
                end

                sum_term = 0;
                if k > 0
                    sum_term = sum(r(1:k).^2 ./ lambda_prime(1:k));
                end
                log_arg = c_LL - sum_term;

                if log_arg <= 0
                    log_term2 = -Inf;
                else
                    log_term2 = -N * log(log_arg);
                end
                
                psi_lg_k = log_term1 + log_term2;
                
                penalty_gaic = k^2 + k;
                penalty_gmdl = 0.5 * (k^2 + k) * log(N);
                
                GAIC_values(k+1) = -psi_lg_k + penalty_gaic;
                GMDL_values(k+1) = -psi_lg_k + penalty_gmdl;
            else % 当k=L-1时, GAIC/GMDL无定义, 设为无穷大
                GAIC_values(k+1) = Inf;
                GMDL_values(k+1) = Inf;
            end
        end
        
        % 寻找最小值对应的k
        [~, M_hat_gaic_idx] = min(GAIC_values);
        M_hat_gaic = M_hat_gaic_idx - 1;
        
        [~, M_hat_gmdl_idx] = min(GMDL_values);
        M_hat_gmdl = M_hat_gmdl_idx - 1;
        
        [~, M_hat_aic_idx] = min(AIC_values);
        M_hat_aic = M_hat_aic_idx - 1;
        
        [~, M_hat_mdl_idx] = min(MDL_values);
        M_hat_mdl = M_hat_mdl_idx - 1;

        % 检查估计是否正确
        if M_hat_gaic == M, correct_gaic = correct_gaic + 1; end
        if M_hat_gmdl == M, correct_gmdl = correct_gmdl + 1; end
        if M_hat_aic == M, correct_aic = correct_aic + 1; end
        if M_hat_mdl == M, correct_mdl = correct_mdl + 1; end
    end

    % 计算并存储正确检测概率
    prob_correct_gaic(idx) = correct_gaic / num_runs;
    prob_correct_gmdl(idx) = correct_gmdl / num_runs;
    prob_correct_aic(idx) = correct_aic / num_runs;
    prob_correct_mdl(idx) = correct_mdl / num_runs;
    
    fprintf('N = %d, P_correct(GAIC)=%.2f, P_correct(GMDL)=%.2f, P_correct(AIC)=%.2f, P_correct(MDL)=%.2f\n', ...
            N, prob_correct_gaic(idx), prob_correct_gmdl(idx), prob_correct_aic(idx), prob_correct_mdl(idx));
end

%% 4. 绘图 (Plotting)
figure;
% 绘制错误检测概率 (Probability of Detection Error = 1 - P_correct)
plot(N_values, 1 - prob_correct_gaic, '-s', 'LineWidth', 1.5, 'MarkerSize', 6);
hold on;
plot(N_values, 1 - prob_correct_gmdl, '-o', 'LineWidth', 1.5, 'MarkerSize', 6);
plot(N_values, 1 - prob_correct_aic, '--^', 'LineWidth', 1.5, 'MarkerSize', 6);
plot(N_values, 1 - prob_correct_mdl, '--x', 'LineWidth', 1.5, 'MarkerSize', 6);
grid on;
title('多种算法检测性能对比 vs. 快照数 (N)', 'FontSize', 14);
xlabel('快照数 (Number of Snapshots, N)', 'FontSize', 12);
ylabel('错误检测概率 (Probability of Detection Error)', 'FontSize', 12);
legend('GAIC', 'GMDL', 'AIC', 'MDL', 'Location', 'northeast');
set(gca, 'FontSize', 11);

Fig.6 的仿真代码

% MATLAB Code for comparing MGAIC, MGMDL, GAIC, GMDL, AIC, and MDL algorithms
% Based on the paper: "Source number estimators using transformed Gerschgorin radii"
% by H.T. Wu, J.F. Yang, and F.K. Chen, IEEE T-SP, 1995.

clear; clc; close all;

%% 1. 仿真参数设置 (Simulation Parameters)
% --- 固定参数 (与论文图6场景类似) ---
L = 8;                  % 传感器数量 (Number of sensors)
M = 2;                  % 真实信源数量 (True number of sources)
DOA = [0, 10];          % 信号入射角度 (Directions of Arrival in degrees)
SNR_dB = [0, 0];        % 信噪比 (Signal-to-Noise Ratios in dB)
d = 0.5;                % 阵元间距 (半波长) (Sensor spacing in half-wavelengths)
num_runs = 200;         % 每个N值的蒙特卡洛仿真次数

% --- 变化参数 (Varying Parameter) ---
N_values = round(logspace(1, 3, 15)); % 从 10^1 到 10^3 取15个点
prob_correct_gaic = zeros(size(N_values));
prob_correct_gmdl = zeros(size(N_values));
prob_correct_aic = zeros(size(N_values));
prob_correct_mdl = zeros(size(N_values));
prob_correct_mgaic = zeros(size(N_values));
prob_correct_mgmdl = zeros(size(N_values));

%% --- 主仿真循环 ---
for idx = 1:length(N_values)
    N = N_values(idx);
    correct_gaic = 0; correct_gmdl = 0;
    correct_aic = 0; correct_mdl = 0;
    correct_mgaic = 0; correct_mgmdl = 0;
    
    % --- 蒙特卡洛循环 ---
    for run = 1:num_runs
        %% 2. 生成阵列接收信号 (Generate Array Signals)
        A = exp(-1j * 2 * pi * d * (0:L-1)' * sind(DOA));
        s_powers = 10.^(SNR_dB / 10);
        s = sqrt(s_powers') .* (randn(M, N) + 1j * randn(M, N)) / sqrt(2);
        noise_w = (randn(L, N) + 1j * randn(L, N)) / sqrt(2); % 高斯白噪声
        X = A * s + noise_w;
        C_hat = (1/N) * (X * X');

        %% 3. 实现所有算法 (Implement All Algorithms)
        % --- 计算传统AIC/MDL ---
        lambda_full = sort(eig(C_hat), 'descend');
        AIC_values = zeros(1, L); MDL_values = zeros(1, L);
        for k = 0:(L-1)
            lambda_noise_full = lambda_full(k+1 : L);
            alpha_full = geo_mean(lambda_noise_full);
            beta_full = mean(lambda_noise_full);
            if alpha_full < 1e-10, log_likelihood_full = Inf; else
                log_likelihood_full = N * (L - k) * log(alpha_full / beta_full);
            end
            AIC_values(k+1) = -log_likelihood_full + k * (2*L - k);
            MDL_values(k+1) = -log_likelihood_full + 0.5 * k * (2*L - k) * log(N);
        end
        [~, M_hat_aic_idx] = min(AIC_values); M_hat_aic = M_hat_aic_idx - 1;
        [~, M_hat_mdl_idx] = min(MDL_values); M_hat_mdl = M_hat_mdl_idx - 1;
        
        % --- 实现GAIC/GMDL 和 MGAIC/MGMDL ---
        GAIC_all_rotations = zeros(L, L-1);
        GMDL_all_rotations = zeros(L, L-1);
        C_current = C_hat;

        for m = 1:L % L次旋转
            % a. 酉变换
            C1_hat = C_current(1:L-1, 1:L-1);
            c_hat = C_current(1:L-1, L);
            c_LL = C_current(L, L);
            [U1_hat, D1_hat] = eig(C1_hat);
            [lambda_prime, sorted_idx] = sort(diag(D1_hat), 'descend');
            U1_hat_sorted = U1_hat(:, sorted_idx);

            % b. 计算Gerschgorin半径
            r = zeros(L-1, 1);
            for i = 1:(L-1)
                r(i) = abs(U1_hat_sorted(:, i)' * c_hat);
            end

            % c. 计算当次旋转的GLE判据
            for k = 0:(L-2)
                lambda_noise_prime = lambda_prime(k+1 : L-1);
                alpha_prime = geo_mean(lambda_noise_prime);
                beta_prime = mean(lambda_noise_prime);
                
                if alpha_prime < 1e-10, log_term1 = Inf; else
                    log_term1 = N * (L-1-k) * log(alpha_prime / beta_prime);
                end

                sum_term = 0;
                if k > 0, sum_term = sum(r(1:k).^2 ./ lambda_prime(1:k)); end
                
                log_arg = c_LL - sum_term;
                if log_arg <= 0, log_term2 = -Inf; else
                    log_term2 = -N * log(log_arg);
                end
                
                psi_lg_k = log_term1 + log_term2;
                GAIC_all_rotations(m, k+1) = -psi_lg_k + (k^2 + k);
                GMDL_all_rotations(m, k+1) = -psi_lg_k + 0.5 * (k^2 + k) * log(N);
            end
            
            % 旋转协方差矩阵用于下一次迭代
            perm_idx = [2:L, 1];
            C_current = C_current(perm_idx, perm_idx);
        end

        % 对L次旋转结果求平均
        MGAIC_values = mean(GAIC_all_rotations, 1);
        MGMDL_values = mean(GMDL_all_rotations, 1);
        
        % 第一次旋转的结果即为GAIC/GMDL的结果
        GAIC_values = GAIC_all_rotations(1, :);
        GMDL_values = GMDL_all_rotations(1, :);

        % 寻找最小值对应的k
        [~, M_hat_gaic_idx] = min(GAIC_values); M_hat_gaic = M_hat_gaic_idx - 1;
        [~, M_hat_gmdl_idx] = min(GMDL_values); M_hat_gmdl = M_hat_gmdl_idx - 1;
        [~, M_hat_mgaic_idx] = min(MGAIC_values); M_hat_mgaic = M_hat_mgaic_idx - 1;
        [~, M_hat_mgmdl_idx] = min(MGMDL_values); M_hat_mgmdl = M_hat_mgmdl_idx - 1;

        % 检查估计是否正确
        if M_hat_gaic == M, correct_gaic = correct_gaic + 1; end
        if M_hat_gmdl == M, correct_gmdl = correct_gmdl + 1; end
        if M_hat_aic == M, correct_aic = correct_aic + 1; end
        if M_hat_mdl == M, correct_mdl = correct_mdl + 1; end
        if M_hat_mgaic == M, correct_mgaic = correct_mgaic + 1; end
        if M_hat_mgmdl == M, correct_mgmdl = correct_mgmdl + 1; end
    end

    % 计算并存储正确检测概率
    prob_correct_gaic(idx) = correct_gaic / num_runs;
    prob_correct_gmdl(idx) = correct_gmdl / num_runs;
    prob_correct_aic(idx) = correct_aic / num_runs;
    prob_correct_mdl(idx) = correct_mdl / num_runs;
    prob_correct_mgaic(idx) = correct_mgaic / num_runs;
    prob_correct_mgmdl(idx) = correct_mgmdl / num_runs;
    
    fprintf('N = %d, P_error(MGAIC)=%.2f, P_error(MGMDL)=%.2f\n', ...
            N, 1-prob_correct_mgaic(idx), 1-prob_correct_mgmdl(idx));
end

%% 4. 绘图 (Plotting)
figure;
% 绘制错误检测概率 (Probability of Detection Error)
semilogx(N_values, 1 - prob_correct_aic, ':^', 'LineWidth', 1.5, 'MarkerSize', 5);
hold on;
semilogx(N_values, 1 - prob_correct_gaic, '-s', 'LineWidth', 1, 'MarkerSize', 6);
semilogx(N_values, 1 - prob_correct_mgaic, '-d', 'LineWidth', 1, 'MarkerSize', 7);
semilogx(N_values, 1 - prob_correct_mdl, ':x', 'LineWidth', 1.5, 'MarkerSize', 7);
semilogx(N_values, 1 - prob_correct_gmdl, '-o', 'LineWidth', 1, 'MarkerSize', 6);
semilogx(N_values, 1 - prob_correct_mgmdl, '-*', 'LineWidth', 1, 'MarkerSize', 7);

grid on;
title('多种算法检测性能对比 vs. 快照数 (N)', 'FontSize', 14);
xlabel('快照数 (Number of Snapshots, N)', 'FontSize', 12);
ylabel('错误检测概率 (Probability of Detection Error)', 'FontSize', 12);
legend('AIC', 'GAIC', 'MGAIC', 'MDL', 'GMDL', 'MGMDL', 'Location', 'northeast');
set(gca, 'FontSize', 11);