杂波边缘及多目标场景下常用恒虚警检测算法的性能推导
在常用CFAR(如CA-CFAR,GO-CFAR,SO-CFAR,OS-CFAR等)的推导过程中，均隐含了均匀背景的假设。即参考窗中所有单元是独立同分布的，且其分布与待检单元中的噪声分布是一致的。当满足这个前提假设时，从参考单元得到的估计噪声，能够很好地反映待检单元中的噪声情况，此时CFAR算法的性能较好。均匀背景下的CFAR算法的详细推导和分析过程可以参考文献[1]。
在实际应用中，上述均匀背景的假设往往难以得到满足，比如在参考窗中同时存在多个干扰目标或参考窗中包含杂波过渡区。此时CFAR将产生部分性能损失，具体地，举几种常见的情况：
(1) 当参考窗中包含多个点干扰目标时，CA-CFAR的噪声估计值偏高，造成检测概率下降；GO-CFAR选择前后参考窗中估计噪声较大的作为噪声估计值，因此也会得到偏高的估计值，造成检测概率下降；SO-CFAR选择前后参考窗中估计噪声较小的作为噪声估计值，因此若这些点目标集中分布于前、后参考窗中的一个时，SO-CFAR的性能受影响较小，当前后参考窗中同时含有干扰目标时，SO-CFAR也会得到偏高的噪声估计值，造成检测性能下降；OS-CFAR若选择的k值不在干扰目标单元上，则性能影响不大。
(2)当参考窗跨越杂波边缘时，当待检单元落在杂波区时，CA-CFAR噪声估计值偏低，容易造成虚警，GO-CFAR噪声估计值影响较小，SO-CFAR噪声估计值偏小，虚警增加，OS-CFAR噪声估计值与k值选取有关；当待检单元落在非杂波区时，CA-CFAR噪声估计值偏高，检测概率下降，GO-CFAR噪声估计值偏小，虚警增加，SO-CFAR噪声估计值影响较小，OS-CFAR噪声估计值与k值选取有关。
对于常用CFAR算法来说，上述两种非理想情况，可以统一表示多干扰目标的情况，为简化分析过程，下面仅考虑所有干扰的强度都保持一致的情况。对各CFAR算法在非理想背景条件下的一般性能进行分析。
设参考窗长为\(N\)，\(n=N/2\)表示半窗长度。设干扰目标数目为\(M\)，且\(M \leq N\)。设噪声功率为\(\sigma_0^2\)，有效目标信噪比为\({\rm SNR}\)，干噪比为\({\rm CNR}\)，则有效目标和干扰的功率分别为\(\sigma_t^2=\sigma_0^2(1+{\rm SNR})\)，\(\sigma_c^2=\sigma_0^2(1+{\rm CNR}\))。经平方率检波后，噪声单元、有效目标单元和干扰单元均服从指数分布，他们的概率密度函数分别表示为:
\[\begin{equation}
\begin{aligned}
f_n(x)=\lambda_0 e^{-\lambda_0 x},\lambda_0=1/\sigma_0^2\\
f_t(x)=\lambda_t e^{-\lambda_t x},\lambda_t=1/\sigma_t^2\\
f_c(x)=\lambda_c e^{-\lambda_c x},\lambda_c=1/\sigma_c^2
\end{aligned}
\end{equation}\tag{1}
\]
(1) CA-CFAR性能分析：
CA-CFAR将参考单元平均值作为噪声估计值，因此其估计噪声可以表示为(需要注意的是，为了推导过程更加简洁，我们将取平均的系数归到了门限乘积因子中)：
\[Z=\sum_{i=1}^r{X_i}+\sum_{i=r+1}^N{X_i}=Z_1+Z_2\tag{2}
\]
参考文献[1]易知：\(Z_1 \backsim {\rm G}(r,\lambda_c),Z_2 \backsim {\rm G}(N-r,\lambda_0)\)，且\({\rm M}_{z_1}(t)=(1+\frac{t}{\lambda_c})^{-r},{\rm M}_{z_2}(t)=(1+\frac{t}{\lambda_0})^{-(N-r)}\),\({\rm M}_z(t)={\rm M}_{z_1}(t){\rm M}_{z_2}(t)\)，所以此时CA-CFAR的虚警率和检测概率可以分别表示为
\[\begin{equation}
\begin{aligned}
P_{FA}&={\rm M}_z(\alpha \lambda_0)={\rm M}_{z_1}(\alpha \lambda_0){\rm M}_{z_2}(\alpha \lambda_0)\\
&=(1+\frac{\alpha \lambda_0}{\lambda_c})^{-r}(1+\alpha)^{-(N-r)}=[1+\alpha(1+{\rm CNR})]^{-r}(1+\alpha)^{-(N-r)}\\
P_D&={\rm M}_z(\alpha \lambda_t)={\rm M}_{z_1}(\alpha \lambda_t){\rm M}_{z_2}(\alpha \lambda_t)\\
&=(1+\frac{\alpha \lambda_t}{\lambda_c})^{-r}(1+\frac{\alpha \lambda_t}{\lambda_0})^{-(N-r)}=[1+\frac{\alpha(1+{\rm CNR})}{(1+{\rm SNR})}]^{-r}[1+\frac{\alpha}{(1+{\rm SNR})}]^{-(N-r)}
\end{aligned}
\end{equation}\tag{3}
\]
显然，当干扰目标的信噪比\({\rm CNR}=0\)时，对应均匀噪声背景的情况，此时\(P_{FA}=(1+\alpha)^{-N},P_D=[1+\alpha/(1+{\rm SNR})]^{-N}\)，与均匀噪声背景情况吻合，证明了(2)的正确性。
由于CA-CFAR是对整个参考单元取平均，因此多干扰目标情况下，其性能与目标在参考单元中的具体分布无关。
上图给出了多干扰目标情况下CA-CFAR理论性能和蒙特卡洛仿真的结果对比，其中第一个图为虚警概率与干扰数目和\({\rm CNR}\)之间的关系，可以看出随着\({\rm CNR}\)的提高，虚警概率逐渐降低，这是因为\({\rm CNR}\)越高，估计的噪声功率越大，CA-CFAR的门限值也就越高。此外，第一个图中理论结果和蒙特卡洛仿真结果好像不是完全吻合，这是因为本身虚警概率就很小(\(10^{-5}\)级别)，而蒙特卡洛仿真次数设置为\(10^6\)，理论上只会有几次虚警，因此结果起伏比较大。第二幅图中两种结果是完全吻合的，证明了(3)式理论推导结果的正确性，同时由于仿真过程中将待检单元目标信噪比设置成和干扰目标信噪比相同，因此才会出现第二幅图中检测概率随\({\rm CNR}\)增大而提高这种"怪现象"。
(2) GO-CFAR性能分析：
GO-CFAR将前、后参考窗中平均功率大的作为估计噪声，显然GO-CFAR的性能与干扰目标的分布有关，设前参考窗中有\(r_1\)个干扰，后参考窗中有\(r_2\)个干扰，同时为了书写简便，将求平均的系数统一到门限乘积因子中去，因此有
\[\begin{equation}
\begin{aligned}
Z_1&=\sum_{i=1}^{r_1}{X_i}+\sum_{i=r_1+1}^n{X_i}=Z_{1-1}+Z_{1-2}\\
Z_2&=\sum_{i=n+1}^{n+r_2}{X_i}+\sum_{i=n+r_2+1}^{2n}{X_i}=Z_{2-1}+Z_{2-2}\\
Z&={\rm max}(Z_1,Z_2)
\end{aligned}
\end{equation}\tag{4}
\]
由此可以得到以下关系：
\[\begin{equation}
\begin{aligned}
Z_{1-1} \backsim {\rm G}(r_1,\lambda_c)&:f_{z_{1-1}}(z)=\frac{\lambda_c^{r_1}}{\Gamma(r_1)}z^{r_1-1}e^{-\lambda_c z}\\
Z_{1-2} \backsim {\rm G}(n-r_1,\lambda_0)&:f_{z_{1-2}}(z)=\frac{\lambda_0^{n-r_1}}{\Gamma(n-r_1)}z^{n-r_1-1}e^{-\lambda_0 z}\\
Z_{2-1} \backsim {\rm G}(r_2,\lambda_c)&:f_{z_{2-1}}(z)=\frac{\lambda_c^{r_2}}{\Gamma(r_2)}z^{r_2-1}e^{-\lambda_c z}\\
Z_{2-2} \backsim {\rm G}(n-r_2,\lambda_0)&:f_{z_{2-2}}(z)=\frac{\lambda_0^{n-r_2}}{\Gamma(n-r_2)}z^{n-r_2-1}e^{-\lambda_0 z}\\
\end{aligned}
\end{equation}\tag{5}
\]
需要注意的是，(5)中为了保证Gamma函数是有意义的，需要满足\(1 \leq r_1 <n,1 \leq r_2 <n\)，因此如果在(5)的基础上进行推导的话，是不能涵盖所有干扰的形式的，因此需要分情况进行推导，最终给出来的性能表达式也应该是一个分段函数的形式，具体地，可以分成下述几个情况进行讨论：
① \(r_1=r_2=0\)或\(r_1=r_2=n\)：此时对应参考窗中无干扰和全干扰的情况，均可均匀背景进行处理，相关结果可参考 《均匀背景背景下常用恒虚警检测推导与比较》；
② \(1 \leq r_1 \leq n-1,1\leq r_2 \leq n-1\)：此时两个参考半窗中均存在干扰目标，且干扰目标数均小于半窗长度；
③ \(r_1=0,1\leq r_2 \leq n-1\)或\(r_1=n,1\leq r_2\leq n-1\)：此时一个半窗部分干扰，另一个半窗无干扰或全干扰
④ \(r_1=0,r_2=n\)：此时只有一个半窗存在干扰，且干扰数目恰好等于半窗长度。
以下分别进行探讨：
① 均匀背景情况
其矩母函数可在《均匀背景背景下常用恒虚警检测推导与比较》中找到，具体表示如下：
\[{\rm M}_z^{(1)}(t)=2(1+t/\lambda)^{-n}-2\sum_{k=0}^{n-1}\begin{pmatrix}n-1+k\\k\end{pmatrix}(2+t/\lambda)^{-(n+k)}\tag{6}
\]
需要注意的是，无干扰情况下\(\lambda=\lambda_0\)，全干扰情况下\(\lambda=\lambda_c\)。
② \(1 \leq r_1 \leq n-1,1\leq r_2 \leq n-1\)，两个参考半窗中均存在部分干扰
参考《两个服从Gamma分布的随机变量的和的pdf和cdf》可得\(Z_1,Z_2\)的概率密度函数和分布函数可以表示为：
\[\begin{equation}
\begin{aligned}
f_{z_1}(z)&=A_0\sum_{i=0}^{r_1-1}A_1(i)[z^i e^{-\lambda_c z}-\sum_{j=0}^{n-i-2}\frac{(\lambda_0-\lambda_c)^j}{j!}z^{i+j}e^{-\lambda_0 z}]\\
f_{z_2}(z)&=\hat{A}_0\sum_{i=0}^{r_2-1}\hat{A}_1(i)[z^i e^{-\lambda_c z}-\sum_{j=0}^{n-i-2}\frac{(\lambda_0-\lambda_c)^j}{j!}z^{i+j}e^{-\lambda_0 z}]\\
F_{z_1}(z)&=A_2-A_0\sum_{m=0}^{r_1-1}A_1(m)[m!\sum_{k=0}^m \frac{\lambda_c^{-(m-k+1)}}{k!}z^k e^{-\lambda_c z}-\sum_{p=0}^{n-m-2}\sum_{q=0}^{m+p}\frac{(m+p)!(\lambda_0-\lambda_c)^p\lambda_0^{-(m+p-q+1)}}{p!q!}z^q e^{-\lambda_0 z}]\\
F_{z_2}(z)&=\hat{A}_2-\hat{A}_0\sum_{m=0}^{r_2-1}\hat{A}_1(m)[m!\sum_{k=0}^m \frac{\lambda_c^{-(m-k+1)}}{k!}z^k e^{-\lambda_c z}-\sum_{p=0}^{n-m-2}\sum_{q=0}^{m+p}\frac{(m+p)!(\lambda_0-\lambda_c)^p\lambda_0^{-(m+p-q+1)}}{p!q!}z^q e^{-\lambda_0 z}]\\
\end{aligned}
\end{equation}\tag{7}
\]
其中
\[\begin{equation}
\begin{aligned}
A_0(\lambda_0,\lambda_c,n,r)&=\frac{\lambda_0^{n-r}\lambda_c^r}{\Gamma(n-r)\Gamma(r)}\\
A_1(\lambda_0,\lambda_c,n,r,i)&=(-1)^{r-1-i}\begin{pmatrix} r-1\\i \end{pmatrix}(n-i-2)!(\lambda_0-\lambda_c)^{-(n-i-1)}\\
A_2(\lambda_0,\lambda_c,n,r)&=A_0(\lambda_0,\lambda_c,n,r)\sum_{i=0}^{r-1}A_1(\lambda_0,\lambda_c,n,r,i)[i!\lambda_c^{-(i+1)}-\sum_{j=0}^{n-i-2}\frac{(i+j)!}{j!}(\lambda_0-\lambda_c)^j \lambda_0^{-(i+j+1)}]\\
A_0(\lambda_0,\lambda_c,n,r_1)&\overset{\Delta}{=}A_0, \quad A_0(\lambda_0,\lambda_c,n,r_2)\overset{\Delta}{=}\hat{A}_0\\
A_1(\lambda_0,\lambda_c,n,r_1,i)&\overset{\Delta}{=}A_1(i),\quad A_1(\lambda_0,\lambda_c,n,r_2,i)\overset{\Delta}{=}\hat{A}_1(i)\\
A_2(\lambda_0,\lambda_c,n,r_1)&\overset{\Delta}{=}A_2,\quad A_2(\lambda_0,\lambda_c,n,r_2)\overset{\Delta}{=}\hat{A}_2
\end{aligned}
\end{equation}\tag{8}
\]
需要注意的是，为了保证(5)中Gamma函数是有意义的，需要满足\(1 \leq r_1 \leq n-1,1 \leq r_2 \leq n-1\)，且\(\lambda_0 \neq \lambda_c\)，也就是说(6)推导得到的模型是不能统一均匀背景的情况的。此外，从《两个服从Gamma分布的随机变量的和的pdf和cdf》中我们知道，利用数值仿真可以发现\(A_2(\lambda_0,\lambda_c,n,r)\equiv 1\)但是还不清楚如何进行严密的证明，不过后续推导仍然将其当成一个已知条件来用。
所以
\[f_z(z)=f_{z_1}(z)F_{z_2}(z)+f_{z_2}(z)F_{z_1}(z)\tag{9}
\]
所以
\[\begin{equation}
\begin{aligned}
{\rm M}_z^{(2)}(t)&=\int_0^\infty f_z(z)e^{-tz}dz\\
&=\int_0^\infty f_{z_1}(z)F_{z_2}(z)e^{-tz}dz+\int_0^\infty f_{z_2}(z)F_{z_1}(z)e^{-tz}dz=I_0+I_1
\end{aligned}
\end{equation}\tag{10}
\]
其中
\[\begin{equation}
\begin{aligned}
I_0&=\int_0^\infty f_{z_1}(z)F_{z_2}(z)e^{-tz}dz\\
&=A_0 \hat{A}_2 \sum_{i=0}^{r_1-1}A_1(i)i!(\lambda_c+t)^{-(i+1)}\\
&-A_0\hat{A}_2\sum_{i=0}^{r_1-1}A_1(i)\sum_{j=0}^{n-i-2}A_3(j)(i+j)!(\lambda_0+t)^{-(i+j+1)}\\
&-A_0\hat{A}_0\sum_{i=0}^{r_1-1}A_1(i)\sum_{m=0}^{r_2-1}\hat{A}_1(m)\sum_{k=0}^m A_4(m,k)(i+k)!(2\lambda_c+t)^{-(i+k+1)}\\
&+A_0\hat{A}_0\sum_{i=0}^{r_1-1}A_1(i)\sum_{j=0}^{n-i-2}A_3(j)\sum_{m=0}^{r_2-1}\hat{A}_1(m)\sum_{k=0}^m A_4(m,k)(i+j+k)!(\lambda_0+\lambda_c+t)^{-(i+j+k+1)}\\
&+A_0\hat{A}_0\sum_{i=0}^{r_1-1}A_1(i)\sum_{m=0}^{r_2-1}\hat{A}_1(m)\sum_{p=0}^{n-m-2}\sum_{q=0}^{m+p}A_5(p,q,m)(i+q)!(\lambda_0+\lambda_c+t)^{-(i+q+1)}\\
&-A_0\hat{A}_0\sum_{i=0}^{r_1-1}A_1(i)\sum_{j=0}^{n-i-2}A_3(j)\sum_{m=0}^{r_2-1}\hat{A}_1(m)\sum_{p=0}^{n-m-2}\sum_{q=0}^{m+p}A_5(p,q,m)(i+j+q)!(2\lambda_0+t)^{-(i+j+q+1)}
\end{aligned}
\end{equation}\tag{11}
\]
\[\begin{equation}
\begin{aligned}
I_1&=\int_0^\infty f_{z_2}(z)F_{z_1}(z)e^{-tz}dz\\
&=\hat{A}_0 A_2 \sum_{i=0}^{r_2-1}\hat{A}_1(i)i!(\lambda_c+t)^{-(i+1)}\\
&-\hat{A}_0 A_2\sum_{i=0}^{r_2-1}\hat{A}_1(i)\sum_{j=0}^{n-i-2}A_3(j)(i+j)!(\lambda_0+t)^{-(i+j+1)}\\
&-\hat{A}_0 A_0\sum_{i=0}^{r_2-1}\hat{A}_1(i)\sum_{m=0}^{r_1-1}A_1(m)\sum_{k=0}^m A_4(m,k)(i+k)!(2\lambda_c+t)^{-(i+k+1)}\\
&+\hat{A}_0 A_0\sum_{i=0}^{r_2-1}\hat{A}_1(i)\sum_{j=0}^{n-i-2}A_3(j)\sum_{m=0}^{r_1-1}A_1(m)\sum_{k=0}^m A_4(m,k)(i+j+k)!(\lambda_0+\lambda_c+t)^{-(i+j+k+1)}\\
&+\hat{A}_0 A_0\sum_{i=0}^{r_2-1}\hat{A}_1(i)\sum_{m=0}^{r_1-1}A_1(m)\sum_{p=0}^{n-m-2}\sum_{q=0}^{m+p}A_5(p,q,m)(i+q)!(\lambda_0+\lambda_c+t)^{-(i+q+1)}\\
&-\hat{A}_0 A_0\sum_{i=0}^{r_2-1}\hat{A}_1(i)\sum_{j=0}^{n-i-2}A_3(j)\sum_{m=0}^{r_1-1}A_1(m)\sum_{p=0}^{n-m-2}\sum_{q=0}^{m+p}A_5(p,q,m)(i+j+q)!(2\lambda_0+t)^{-(i+j+q+1)}
\end{aligned}
\end{equation}\tag{12}
\]
在上面的表达式中
\[\begin{equation}
\begin{aligned}
A_3(\lambda_0,\lambda_c,j)&=\frac{(\lambda_0-\lambda_c)^j}{j!}\overset{\Delta}{=}A_3(j)\\
A_4(\lambda_c,m,k)&=\frac{m!\lambda_c^{-(m-k+1)}}{k!}\overset{\Delta}{=}A_4(m,k)\\
A_5(\lambda_0,\lambda_c,p,q,m)&=\frac{(m+p)!(\lambda_0-\lambda_c)^p \lambda_0^{-(m+p-q+1)}}{p!q!}\overset{\Delta}{=}A_5(p,q,m)
\end{aligned}
\end{equation}\tag{13}
\]
③ \(r_1=0,1\leq r_2 \leq n-1\)或\(r_1=n,1\leq r_2\leq n-1\)：一个半窗部分干扰，另一个半窗无干扰或全干扰
设前参考窗中没有干扰目标，后参考窗中有\(r\)个干扰，则此时前后参考窗的概率密度函数和分布函数可以分别表示为
\[\begin{equation}
\begin{aligned}
f_{z_1}(z)&=\frac{\lambda^n}{\Gamma(n)}z^{n-1}e^{-\lambda z}\\
f_{z_2}(z)&=A_0\sum_{i=0}^{r-1}A_1(i)[z^i e^{-\lambda_c z}-\sum_{j=0}^{n-i-2}\frac{(\lambda_0-\lambda_c)^j}{j!}z^{i+j}e^{-\lambda_0 z}]\\
F_{z_1}(z)&=1-\sum_{k=0}^{n-1}\frac{\lambda^k}{k!}z^k e^{-\lambda z}\\
F_{z_2}(z)&=A_2-A_0\sum_{m=0}^{r-1}A_1(m)[m!\sum_{k=0}^m \frac{\lambda_c^{-(m-k+1)}}{k!}z^k e^{-\lambda_c z}-\sum_{p=0}^{n-m-2}\sum_{q=0}^{m+p}\frac{(m+p)!(\lambda_0-\lambda_c)^p\lambda_0^{-(m+p-q+1)}}{p!q!}z^q e^{-\lambda_0 z}]\\
\end{aligned}
\end{equation}\tag{14}
\]
所以，整个参考单元的概率密度函数可以表示为
\[\begin{equation}
\begin{aligned}
f_z(z)=f_{z_1}(z)F_{z_2}(z)+f_{z_2}(z)F_{z_1}(z)
\end{aligned}
\end{equation}\tag{15}
\]
所以，其矩母函数可以表示为
\[\begin{equation}
\begin{aligned}
{\rm M}_z^{(3)}(t)=\int_0^\infty f_z(z)dz=I_0+I_1
\end{aligned}
\end{equation}\tag{16}
\]
其中
\[\begin{equation}
\begin{aligned}
I_0&=A_2(1+t/\lambda)^{-n}\\
&-A_0 \lambda^n\sum_{m=0}^{r-1}A_1(m)m!\sum_{k=0}^m \begin{pmatrix}n+k-1\\k\end{pmatrix}\lambda_c^{-(m-k+1)}(\lambda+\lambda_c+t)^{-(n+k)}\\
&+A_0 \lambda^n\sum_{m=0}^{r-1}A_1(m)\sum_{p=0}^{n-m-2}\sum_{q=0}^{m+p}\frac{(m+p)!}{p!}\begin{pmatrix}n+q-1\\q\end{pmatrix}(\lambda_0-\lambda_c)^p\lambda_0^{-(m+p-q+1)}(\lambda+\lambda_0+t)^{-(n+q)}\\
I_1&=A_0\sum_{i=0}^{r-1}A_1(i)i!(\lambda_c+t)^{-(i+1)}\\
&-A_0\sum_{i=0}^{r-1}A_1(i)\sum_{j=0}^{n-i-2}\frac{(i+j)!}{j!}(\lambda_0-\lambda_c)^j(\lambda_0+t)^{-(i+j+1)}\\
&-A_0\sum_{i=0}^{r-1}A_1(i)\sum_{k=0}^{n-1}\frac{(k+i)!}{k!}\lambda^k(\lambda+\lambda_c+t)^{-(k+i+1)}\\
&+A_0\sum_{i=0}^{r-1}A_1(i)\sum_{j=0}^{n-i-2}\frac{(\lambda_0-\lambda_c)^j}{j!}\sum_{k=0}^{n-1}\frac{(i+j+k)!}{k!}\lambda^k(\lambda+\lambda_0+t)^{-(i+j+k+1)}
\end{aligned}
\end{equation}\tag{17}
\]
需要注意的是，当\(r_1=0\)时，\(\lambda=\lambda_0\)，当\(r_1=n\)时，\(\lambda=\lambda_c\)。
④ \(r_1=0,r_2=n\)：此时只有一个半窗受干扰且干扰数目等于半窗长度
此时两个半窗的概率密度函数和分布函数可以分别表示为
\[\begin{equation}
\begin{aligned}
f_{z_1}(z)&=\frac{\lambda_0^n}{\Gamma(n)}z^{n-1}e^{-\lambda_0 z}\\
f_{z_2}(z)&=\frac{\lambda_c^n}{\Gamma(n)}z^{n-1}e^{-\lambda_c z}\\
F_{z_1}(z)&=1-\sum_{k=0}^{n-1}\frac{\lambda_0^k}{k!}z^k e^{-\lambda_0 z}\\
F_{z_2}(z)&=1-\sum_{k=0}^{n-1}\frac{\lambda_c^k}{k!}z^k e^{-\lambda_c z}\\
\end{aligned}
\end{equation}\tag{18}
\]
所以统计量的概率密度函数可以表示为
\[\begin{equation}
\begin{aligned}
f_z(z)&=f_{z_1}(z)F_{z_2}(z)+f_{z_2}(z)F_{z_1}(z)\\
&=\frac{\lambda_0^n}{\Gamma(n)}z^{n-1}e^{-\lambda_0 z}+\frac{\lambda_c^n}{\Gamma(n)}z^{n-1}e^{-\lambda_c z}\\
&-\frac{\lambda_0^n}{\Gamma(n)}\sum_{k=0}^{n-1}\frac{\lambda_c^k}{k!}z^{n+k-1}e^{-(\lambda_0+\lambda_c)z}-\frac{\lambda_c^n}{\Gamma(n)}\sum_{k=0}^{n-1}\frac{\lambda_0^k}{k!}z^{n+k-1}e^{-(\lambda_0+\lambda_c)z}
\end{aligned}
\end{equation}\tag{19}
\]
所以其矩母函数可以表示为
\[\begin{equation}
\begin{aligned}
{\rm M}_z^{(4)}(t)&=(1+t/\lambda_0)^{(-n)}+(1+t/\lambda_c)^{(-n)}\\
&-\lambda_0^n\sum_{k=0}^{n-1}\begin{pmatrix}n+k-1\\k\end{pmatrix}\lambda_c^k(\lambda_0+\lambda_c+t)^{-(n+k)}-\lambda_c^n\sum_{k=0}^{n-1}\begin{pmatrix}n+k-1\\k\end{pmatrix}\lambda_0^k(\lambda_0+\lambda_c+t)^{-(n+k)}
\end{aligned}
\end{equation}\tag{20}
\]
至此，所有情况都已推导完毕，总结如下：
\[\begin{equation}
\begin{aligned}
{\rm M}_z(t)=
\begin{cases}
{\rm M}_z^{(1)}(z),&\text{$r_1=r_2=0||r_1=r_2=n$}\\
{\rm M}_z^{(2)}(z),&\text{$1\leq r_1,r_2\leq n-1$}\\
{\rm M}_z^{(3)}(z),&\text{$r_1=0,1\leq r_2\leq n-1$||$r_1=n,1\leq r_2\leq n-1$}\\
{\rm M}_z^{(4)}(z),&\text{$r_1=0,r_2=n$}
\end{cases}
\end{aligned}
\end{equation}\tag{21}
\]
针对上面的情况进行理论计算和蒙特卡洛仿真，仿真代码如附录(3)所示，仿真结果展示如下：
从上面的结果可以看出，理论计算结果和蒙特卡洛仿真结果吻合得很好，证明了上述推导结果的正确性。
(3)SO-CFAR性能分析：
SO-CFAR对比前后参考窗均值大小，并将其中较小的作为噪声估计结果。显然，它跟GO-CFAR一样，性能与干扰目标在参考窗中的具体分布有关。它的整个建模过程与(2)中GO-CFAR的建模过程基本一致，只是其概率密度函数表达式为：
\[\begin{equation}
\begin{aligned}
f_z(z)&=f_{z_1}(z)+f_{z_2}(z)-[f_{z_1}(z)F_{z_2}(z)+f_{z_2}(z)F_{z_1}(z)]\\
&=f_{z_1}(z)[1-F_{z_2}(z)]+f_{z_2}(z)[1-F_{z_1}(z)]
\end{aligned}
\end{equation}\tag{22}
\]
由于SO-CFAR的很多推导过程和推导结果可以直接从GO-CFAR中获取，因此这边就省去中间的繁琐的推导过程，直接给出推导结果：
\[\begin{equation}
\begin{aligned}
{\rm M}_z^{(1)}(z)&=2\sum_{k=0}^{n-1}\lambda^{n+k}\begin{pmatrix}n+k-1\\k\end{pmatrix}(2\lambda+t)^{-(n+k)}\\
{\rm M}_z^{(2)}(z)&=A_0\hat{A}_0\sum_{i=0}^{r_1-1}A_1(i)\sum_{m=0}^{r_2-1}\hat{A}_1(m)\sum_{k=0}^m A_4(m,k)(i+k)!(2\lambda_c+t)^{-(i+k+1)}\\
&-A_0\hat{A}_0\sum_{i=0}^{r_1-1}A_1(i)\sum_{j=0}^{n-i-2}A_3(j)\sum_{m=0}^{r_2-1}\hat{A}_1(m)\sum_{k=0}^m A_4(m,k)(i+j+k)!(\lambda_0+\lambda_c+t)^{-(i+j+k+1)}\\
&-A_0\hat{A}_0\sum_{i=0}^{r_1-1}A_1(i)\sum_{m=0}^{r_2-1}\hat{A}_1(m)\sum_{p=0}^{n-m-2}\sum_{q=0}^{m+p}A_5(p,q,m)(i+q)!(\lambda_0+\lambda_c+t)^{-(i+q+1)}\\
&+A_0\hat{A}_0\sum_{i=0}^{r_1-1}A_1(i)\sum_{j=0}^{n-i-2}A_3(j)\sum_{m=0}^{r_2-1}\hat{A}_1(m)\sum_{p=0}^{n-m-2}\sum_{q=0}^{m+p}A_5(p,q,m)(i+j+q)!(2\lambda_0+t)^{-(i+j+q+1)}\\
&+\hat{A}_0 A_0\sum_{i=0}^{r_2-1}\hat{A}_1(i)\sum_{m=0}^{r_1-1}A_1(m)\sum_{k=0}^m A_4(m,k)(i+k)!(2\lambda_c+t)^{-(i+k+1)}\\
&-\hat{A}_0 A_0\sum_{i=0}^{r_2-1}\hat{A}_1(i)\sum_{j=0}^{n-i-2}A_3(j)\sum_{m=0}^{r_1-1}A_1(m)\sum_{k=0}^m A_4(m,k)(i+j+k)!(\lambda_0+\lambda_c+t)^{-(i+j+k+1)}\\
&-\hat{A}_0 A_0\sum_{i=0}^{r_2-1}\hat{A}_1(i)\sum_{m=0}^{r_1-1}A_1(m)\sum_{p=0}^{n-m-2}\sum_{q=0}^{m+p}A_5(p,q,m)(i+q)!(\lambda_0+\lambda_c+t)^{-(i+q+1)}\\
&+\hat{A}_0 A_0\sum_{i=0}^{r_2-1}\hat{A}_1(i)\sum_{j=0}^{n-i-2}A_3(j)\sum_{m=0}^{r_1-1}A_1(m)\sum_{p=0}^{n-m-2}\sum_{q=0}^{m+p}A_5(p,q,m)(i+j+q)!(2\lambda_0+t)^{-(i+j+q+1)}\\
M_z^{(3)}(z)&=A_0 \lambda^n\sum_{m=0}^{r-1}A_1(m)m!\sum_{k=0}^m \begin{pmatrix}n+k-1\\k\end{pmatrix}\lambda_c^{-(m-k+1)}(\lambda+\lambda_c+t)^{-(n+k)}\\
&-A_0 \lambda^n\sum_{m=0}^{r-1}A_1(m)\sum_{p=0}^{n-m-2}\sum_{q=0}^{m+p}\frac{(m+p)!}{p!}\begin{pmatrix}n+q-1\\q\end{pmatrix}(\lambda_0-\lambda_c)^p\lambda_0^{-(m+p-q+1)}(\lambda+\lambda_0+t)^{-(n+q)}\\
&+A_0\sum_{i=0}^{r-1}A_1(i)\sum_{k=0}^{n-1}\frac{(k+i)!}{k!}\lambda^k(\lambda+\lambda_c+t)^{-(k+i+1)}\\
&-A_0\sum_{i=0}^{r-1}A_1(i)\sum_{j=0}^{n-i-2}\frac{(\lambda_0-\lambda_c)^j}{j!}\sum_{k=0}^{n-1}\frac{(i+j+k)!}{k!}\lambda^k(\lambda+\lambda_0+t)^{-(i+j+k+1)}\\
{\rm M}_z^{(4)}(z)&=\lambda_0^n\sum_{k=0}^{n-1}\begin{pmatrix}n+k-1\\k\end{pmatrix}\lambda_c^k(\lambda_0+\lambda_c+t)^{-(n+k)}+\lambda_c^n\sum_{k=0}^{n-1}\begin{pmatrix}n+k-1\\k\end{pmatrix}\lambda_0^k(\lambda_0+\lambda_c+t)^{-(n+k)}
\end{aligned}
\end{equation}\tag{23}
\]
将所有结果按照(21)进行综合，可以得到最终的表达式。由于它与GO-CFAR公式中有很多相似的地方，所以这边不给出仿真代码，只展示结果
(4) OS-CFAR性能分析：
OS-CFAR将参考窗内的值进行排序，然后选择某个单元的值作为估计噪声值。其中排序操作是针对整个参考窗进行的，因此干扰目标在参考窗中的位置并不影响最终的结果。若取排序后第\(k\)个参考单元的值进行噪声功率估计(即OS-CFAR对应的统计量为\(Z=X_{(k)}\))。则参考《次序统计量的概率密度函数》中的推导思路，这边采用基于分布函数的推导思路，具体推导过程如下：
根据分布函数的定义有\(F_{X_{(k)}}(x)=P(X_{(k)} \leq x)\)，则参考附录(2)有：
\[\begin{equation}
\begin{aligned}
F_{X_{(k)}}(x)&=P(X_{(k)} \leq x)=\sum_{i=k}^N{P[有i个样本值不大于x]}\\
&=\sum_{i=k}^N{\sum_j{P[有j个不受干扰的样本值不大于x，有i-j个干扰不大于x]}}
\end{aligned}
\end{equation}\tag{24}
\]
基于(7)，现在需要对\(j\)的取值范围进行分析，显然\(j \geq 0,j \leq N-M,i-j \leq M,i-j \geq 0\)，所以可以得到\({\rm max}(0,i-M) \leq j \leq {\rm min}(i,N-M)\)，所以(7)可以写成
\[\begin{equation}
\begin{aligned}
F_{X_{(k)}}(x)=\sum_{i=k}^N{\sum_{{\rm max}(0,i-M)}^{{{\rm min}(i,N-M)}}{\begin{pmatrix} N-M \\ j \end{pmatrix} \begin{pmatrix} M \\ i-j \end{pmatrix}[F_n(x)]^j [1-F_n(x)]^{N-M-j} [F_c(x)]^{i-j} [1-F_c(x)]^{M-i+j}}}
\end{aligned}
\end{equation}\tag{25}
\]
在上面的表达式中\(F_n(x)\)和\(F_c(x)\)分别表示噪声背景和干扰目标的分布函数，其中\(F_n(x)=1-e^{-\lambda_0 x}\)，\(F_c(x)=1-e^{-\lambda_c x}\)，将其代入(8)可以得到
\[\begin{equation}
\begin{aligned}
F_{X_{(k)}}(x)=\sum_{i=k}^N \sum_{{\rm max}(0,i-M)}^{{\rm min}(i,N-M)}\begin{pmatrix} N-M \\ j \end{pmatrix} \begin{pmatrix} M \\ i-j \end{pmatrix} (1-e^{-\lambda_0 x})^j e^{-\lambda_0(N-M-j)x} (1-e^{-\lambda_c x})^{i-j} e^{-\lambda_c(M-i+j)x}
\end{aligned}
\end{equation}\tag{26}
\]
进一步地，利用二项展开式有：
\[\begin{equation}
\begin{aligned}
(1-e^{-\lambda_0 x})^j&=\sum_{l=0}^j{(-1)^l \begin{pmatrix} j \\ l \end{pmatrix} e^{-\lambda_0 l x}}\\
(1-e^{-\lambda_c x})^{i-j}&=\sum_{p=0}^{i-j}{(-1)^p \begin{pmatrix} i-j \\ p \end{pmatrix} e^{-\lambda_c p x}}
\end{aligned}
\end{equation}\tag{27}
\]
将(10)代入(9)化简得到：
\[\begin{equation}
\begin{aligned}
F_{X_{(k)}}(x)=\sum_{i=k}^N \sum_{{\rm max}(0,i-M)}^{{\rm min}(i,N-M)}\begin{pmatrix} N-M \\ j \end{pmatrix} \begin{pmatrix} M \\ i-j \end{pmatrix} [\sum_{l=0}^j \sum_{p=0}^{i-j} (-1)^{p+l} \begin{pmatrix} j \\ l \end{pmatrix} \begin{pmatrix} i-j \\ p \end{pmatrix} e^{-[\lambda_0(N-M-j+l)+\lambda_c(M-i+j+p)]x}]
\end{aligned}
\end{equation}\tag{28}
\]
则\(f_{X_{(k)}}(x)\)可以通到对\(F_{X_{(k)}}(x)\)求导得到：
\[\begin{equation}
\begin{aligned}
f_{X_{(k)}}(x)&=\sum_{i=k}^N \sum_{{\rm max}(0,i-M)}^{{\rm min}(i,N-M)}\begin{pmatrix} N-M \\ j \end{pmatrix} \begin{pmatrix} M \\ i-j \end{pmatrix} [\sum_{l=0}^j \sum_{p=0}^{i-j} (-1)^{p+l+1}\\
&[\lambda_0(N-M-j+l)+\lambda_c(M-i+j+p)] \begin{pmatrix} j \\ l \end{pmatrix} \begin{pmatrix} i-j \\ p \end{pmatrix} e^{-[\lambda_0(N-M-j+l)+\lambda_c(M-i+j+p)]x}]
\end{aligned}
\end{equation}\tag{29}
\]
所以\(f_{X_{(x)}}(x)\)对应的矩母函数可以表示为(首先将\(X_{(k)}\)用符号\(Z\)表示)：
\[\begin{equation}
\begin{aligned}
{\rm M}_z(t)&=\int_0^\infty{f_z(z) e^{-tz}dz}\\
&=...\int_0^\infty{e^{-[\lambda_0(N-M-j+l)+\lambda_c(M-i+j+p)+t]z}dz}\\
&=...[\lambda_0(N-M-j+l)+\lambda_c(M-i+j+p)+t]^{-1}
\end{aligned}
\end{equation}\tag{30}
\]
所以此时其虚警概率和检测概率可以分别表示为
\[\begin{equation}
\begin{aligned}
P_{FA}&={\rm M}_z(\alpha \lambda_0)=...[\lambda_0(N-M-j+l+\alpha)+\lambda_c(M-i+j+p)]^{-1}\\
P_D&={\rm M}_z(\alpha \lambda_t)=...[\lambda_0(N-M-j+l+\alpha)+\lambda_c(M-i+j+p)+\alpha \lambda_t]^{-1}
\end{aligned}
\end{equation}\tag{31}
\]
值得注意的是，上述推导OS-CFAR性能的过程是完全按照之前的推导思路，即得到统计量的概率密度函数之后，再求其对应的矩母函数，然后由矩母函数和虚警概率、检测概率之间的关系得到OS-CFAR的性能。虽然这种做法也是可以的，但是由于在上述求解过程中首先得到的是统计量的分布函数，所以是否可以考虑省略掉由分布函数计算概率密度函数的步骤，而直接利用分布函数求虚警概率和检测概率，显然应该是可以的，具体操作过程如下：
由之前的分析可以，虚警概率和矩母函数之间存在如下关系：
\[\begin{equation}
\begin{aligned}
P_{FA}&={\rm M}_z(\alpha \lambda_0)=\int_0^\infty{f_z(z)e^{-\alpha \lambda_0 z}dz}\\
&=\int_0^\infty{F_z'(z) e^{-\alpha \lambda_0 z}dz}=\int_0^\infty{e^{-\alpha \lambda_0 z}dF_z(z)}\\
&=F_z(z)e^{-\alpha \lambda_0 z}|_0^\infty-\int_0^\infty{F_z(z)de^{-\alpha \lambda_0 z}}\\
&=\alpha \lambda_0 \int_0^\infty{F_z(z) e^{-\alpha \lambda_0 z}dz}
\end{aligned}
\end{equation}\tag{32}
\]
同理可得\(P_D=\alpha \lambda_t \int_0^\infty{F_z(z)e^{-\alpha \lambda_t z}dz}\)，所以将\(F_z(z)\)表达式代入可以求得：
\[\begin{equation}
\begin{aligned}
P_{FA}&=\alpha \lambda_0 \sum_{i=k}^N \sum_{j={\rm max}(0,i-M)}^{{\rm min}(i,N-M)}\begin{pmatrix} N-M \\ j \end{pmatrix} \begin{pmatrix} M \\ i-j \end{pmatrix} [\sum_{l=0}^j \sum_{p=0}^{i-j} (-1)^{p+l} \begin{pmatrix} j \\ l \end{pmatrix} \begin{pmatrix} i-j \\ p \end{pmatrix} \int_0^\infty{e^{-[\lambda_0(N-M-j+l+\alpha)+\lambda_c(M-i+j+p)]z}dz}]\\
&=\alpha \lambda_0 \sum_{i=k}^N \sum_{j={\rm max}(0,i-M)}^{{\rm min}(i,N-M)}\begin{pmatrix} N-M \\ j \end{pmatrix} \begin{pmatrix} M \\ i-j \end{pmatrix} \lbrace \sum_{l=0}^j \sum_{p=0}^{i-j} (-1)^{p+l} \begin{pmatrix} j \\ l \end{pmatrix} \begin{pmatrix} i-j \\ p \end{pmatrix} \\
&[\lambda_0(N-M-j+l+\alpha)+\lambda_c(M-i+j+p)]^{-1}\rbrace
\end{aligned}
\end{equation}\tag{33}
\]
\[\begin{equation}
\begin{aligned}
P_D&=\alpha \lambda_t \sum_{i=k}^N \sum_{j={\rm max}(0,i-M)}^{{\rm min}(i,N-M)}\begin{pmatrix} N-M \\ j \end{pmatrix} \begin{pmatrix} M \\ i-j \end{pmatrix} [\sum_{l=0}^j \sum_{p=0}^{i-j} (-1)^{p+l} \begin{pmatrix} j \\ l \end{pmatrix} \begin{pmatrix} i-j \\ p \end{pmatrix} \int_0^\infty{e^{-[\lambda_0(N-M-j+l)+\lambda_c(M-i+j+p)+\alpha \lambda_t]z}dz}]\\
&=\alpha \lambda_t \sum_{i=k}^N \sum_{j={\rm max}(0,i-M)}^{{\rm min}(i,N-M)}\begin{pmatrix} N-M \\ j \end{pmatrix} \begin{pmatrix} M \\ i-j \end{pmatrix} \lbrace \sum_{l=0}^j \sum_{p=0}^{i-j} (-1)^{p+l} \begin{pmatrix} j \\ l \end{pmatrix} \begin{pmatrix} i-j \\ p \end{pmatrix} \\
&[\lambda_0(N-M-j+l)+\lambda_c(M-i+j+p)+\alpha \lambda_t]^{-1}\rbrace
\end{aligned}
\end{equation}\tag{34}
\]
(16)(17)得到的结果与参考文献[1]、[2]中的结果相吻合。
上述理论分析结果与蒙特卡洛仿真结果的对比代码参考附录(2)。
[2] Wilson S L. Two CFAR algorithms for interfering targets and nonhomogeneous clutter[J]. IEEE transactions on aerospace and electronic systems, 1993, 29(1): 57-72.
[3] Gandhi P P, Kassam S A. Analysis of CFAR processors in nonhomogeneous background[J]. IEEE Transactions on Aerospace and Electronic systems, 1988, 24(4): 427-445.
(1)多干扰目标条件下CA-CFAR算法性能代码
N=32; %参考窗长
PFA=1e-5; %预设虚警概率
mcTime=1e6; %蒙特卡洛仿真次数
M=[4,8,12]; %干扰数目
noiseP=0; %背景噪声功率
CNR=0:2:30; %杂噪比(dB)
SNR=CNR; %目标信噪比(dB)
sigma_0=db2pow(noiseP); %噪声功率
sigma_c=db2pow(noiseP+CNR); %干扰功率
sigma_t=db2pow(noiseP+SNR); %目标功率
lambda_0=1./sigma_0; %噪声指数分布参数
lambda_c=1./(1+sigma_c); %干扰指数分布参数
lambda_t=1./(1+sigma_t); %目标指数分布参数
alpha=func_CaAlpha(PFA,N); %计算门限乘积因子
theoPfa=zeros(length(M),length(CNR)); %理论计算得到的虚警概率
mcPfa=zeros(length(M),length(CNR)); %蒙特卡洛仿真得到的虚警概率
theoPd=zeros(length(M),length(CNR)); %理论计算得到的检测概率
mcPd=zeros(length(M),length(CNR)); %蒙特卡洛仿真得到的检测概率
% 计算理论性能
for m=1:length(M)
Mz1=(1+alpha*lambda_0./(N*lambda_c)).^(-M(m));
Mz2=(1+alpha/N)^(-(N-M(m)));
theoPfa(m,:)=Mz1.*Mz2;
Mz1=(1+alpha*lambda_t./(N*lambda_c)).^(-M(m));
Mz2=(1+alpha*lambda_t./(N*lambda_0)).^(-(N-M(m)));
theoPd(m,:)=Mz1.*Mz2;
end
% 蒙特卡洛仿真性能
for m=1:length(M)
for icnr=1:length(CNR)
faCnt=0;
pdCnt=0;
for jj=1:mcTime
noise=exprnd(sigma_0,N-M(m),1); %生成噪声背景
inter=exprnd(sigma_c(icnr),M(m),1); %生成干扰
noTrgtCut=exprnd(sigma_0,1,1); %生成不含目标待检单元
trgtCut=exprnd(sigma_t(icnr),1,1); %生成含目标待检单元
refWin=[noise;inter]; %生成参考窗
theshold=alpha*mean(refWin); %计算门限
faCnt=faCnt+(noTrgtCut>=theshold);
pdCnt=pdCnt+(trgtCut>=theshold);
end
mcPfa(m,icnr)=faCnt/mcTime;
mcPd(m,icnr)=pdCnt/mcTime;
end
end
figure;hold on;box on;
plot(CNR,theoPfa(1,:),'r');
plot(CNR,mcPfa(1,:),'r.-');
plot(CNR,theoPfa(2,:),'b');
plot(CNR,mcPfa(2,:),'b.-');
plot(CNR,theoPfa(3,:),'k');
plot(CNR,mcPfa(3,:),'k.-');
xlabel('CNR(dB)');ylabel('P_{FA}');
legend('theo:M=4','MC:M=4','theo:M=8','MC:M=8','theo:M=12','MC:M=12');
figure;hold on;box on;
plot(CNR,theoPd(1,:),'r');
plot(CNR,mcPd(1,:),'r.-');
plot(CNR,theoPd(2,:),'b');
plot(CNR,mcPd(2,:),'b.-');
plot(CNR,theoPd(3,:),'k');
plot(CNR,mcPd(3,:),'k.-');
xlabel('CNR(dB)');ylabel('P_{D}');
legend('theo:M=4','MC:M=4','theo:M=8','MC:M=8','theo:M=12','MC:M=12');
(2)多干扰目标条件下SO-CFAR算法性能代码
N=32; %参考窗长
k=round(N*3/4); %噪声估计单元索引
PFA=1e-5; %预设虚警概率
mcTime=1e6; %蒙特卡洛仿真次数
M=[4,12]; %干扰数目
noiseP=0; %背景噪声功率
CNR=0:1:20; %杂噪比(dB)
SNR=CNR; %目标信噪比(dB)
sigma_0=db2pow(noiseP); %噪声功率
sigma_c=db2pow(noiseP+CNR); %干扰功率
sigma_t=db2pow(noiseP+SNR); %目标功率
lambda_0=1./sigma_0; %噪声指数分布参数
lambda_c=1./(1+sigma_c); %干扰指数分布参数
lambda_t=1./(1+sigma_t); %目标指数分布参数
alpha=func_OsAlpha(PFA,N,k); %计算门限乘积因子
theoPfa=zeros(length(M),length(CNR)); %理论计算得到的虚警概率
mcPfa=zeros(length(M),length(CNR)); %蒙特卡洛仿真得到的虚警概率
theoPd=zeros(length(M),length(CNR)); %理论计算得到的检测概率
mcPd=zeros(length(M),length(CNR)); %蒙特卡洛仿真得到的检测概率
% 计算理论性能
for m=1:length(M)
for icnr=1:length(CNR)
tmp1=0;
tmp2=0;
for i=k:N
for j=max(0,i-M(m)):min(i,N-M(m))
a1=nchoosek(N-M(m),j);
a2=nchoosek(M(m),i-j);
c1=0;
c2=0;
for l=0:j
for p=0:i-j
b1=nchoosek(j,l);
b2=nchoosek(i-j,p);
b3=(lambda_0*(N-M(m)-j+l+alpha)+lambda_c(icnr)*(M(m)-i+j+p))^(-1);
b4=(lambda_0*(N-M(m)-j+l)+lambda_c(icnr)*(M(m)-i+j+p)+alpha*lambda_t(icnr))^(-1);
c1=c1+(-1)^(p+l)*b1*b2*b3;
c2=c2+(-1)^(p+l)*b1*b2*b4;
end
end
tmp1=tmp1+a1*a2*c1;
tmp2=tmp2+a1*a2*c2;
end
end
theoPfa(m,icnr)=alpha*lambda_0*tmp1;
theoPd(m,icnr)=alpha*lambda_t(icnr)*tmp2;
end
end
% 蒙特卡洛仿真性能
for m=1:length(M)
for icnr=1:length(CNR)
faCnt=0;
pdCnt=0;
for jj=1:mcTime
noise=exprnd(sigma_0,N-M(m),1); %生成噪声背景
inter=exprnd(sigma_c(icnr),M(m),1); %生成干扰
noTrgtCut=exprnd(sigma_0,1,1); %生成不含目标待检单元
trgtCut=exprnd(sigma_t(icnr),1,1); %生成含目标待检单元
refWin=[noise;inter]; %生成参考窗
tmp=sort(refWin,'ascend');
theshold=alpha*tmp(k); %计算门限
faCnt=faCnt+(noTrgtCut>=theshold);
pdCnt=pdCnt+(trgtCut>=theshold);
end
mcPfa(m,icnr)=faCnt/mcTime;
mcPd(m,icnr)=pdCnt/mcTime;
end
end
figure;hold on;box on;
plot(CNR,theoPfa(1,:),'r');
plot(CNR,mcPfa(1,:),'r.-');
plot(CNR,theoPfa(2,:),'b');
plot(CNR,mcPfa(2,:),'b.-');
xlabel('CNR(dB)');ylabel('P_{FA}');
legend('theo:M=4','MC:M=4','theo:M=12','MC:M=12');
figure;hold on;box on;
plot(CNR,theoPd(1,:),'r');
plot(CNR,mcPd(1,:),'r.-');
plot(CNR,theoPd(2,:),'b');
plot(CNR,mcPd(2,:),'b.-');
xlabel('CNR(dB)');ylabel('P_{D}');
legend('theo:M=4','MC:M=4','theo:M=12','MC:M=12');
(3)多干扰目标条件下GO-CFAR算法性能代码
function main
N=32; %窗长
n=N/2; %半窗长度
PFA=1e-6; %预设虚警
mcTimes=1e5; %蒙特卡洛仿真次数
r1=2; %前参考窗干扰目标数目
r2=2; %后参考窗干扰数目
noiseP=0; %噪声功率(dB)
SNR=20; %目标信噪比(dB)
CNR=5:2:29; %干扰目标信噪比(dB)
sigma_0=db2pow(noiseP);
sigma_t=db2pow(noiseP+SNR);
sigma_c=db2pow(noiseP+CNR);
lambda_0=1./sigma_0;
lambda_c=1./(1+sigma_c);
lambda_t=1./(1+sigma_t);
alpha=func_GoAlpha(PFA,N); %计算门限乘积因子
theoPd=zeros(length(CNR),1);%理论检测概率
mcPd=zeros(length(CNR),1); %蒙特卡洛计算得到的检测概率
% 计算理论检测概率
for icnr=1:length(CNR)
theoPd(icnr)=func_GoPd(lambda_0,lambda_c(icnr),lambda_t,n,r1,r2,alpha);
end
% 进行蒙特卡洛仿真
for icnr=1:length(CNR)
cnt=0;
for imc=1:mcTimes
noise=exprnd(sigma_0,N-r1-r2,1); %生成背景噪声
leadClutter=exprnd(sigma_c(icnr),r1,1); %前参考窗中的干扰目标
lagClutter=exprnd(sigma_c(icnr),r2,1); %后参考窗中的干扰目标
cut=exprnd(sigma_t,1,1); %待检单元
leadWin=[noise(1:n-r1);leadClutter];
lagWin=[noise(n-r1+1:end);lagClutter];
goT=alpha*max(mean(leadWin),mean(lagWin));
cnt=cnt+(cut>=goT);
end
mcPd(icnr)=cnt/mcTimes;
end
% 绘图比较
figure;hold on;box on;
plot(CNR,theoPd,'r.-');
plot(CNR,mcPd,'bo-');
xlabel('CNR(dB)');ylabel('检测概率');axis tight;ylim([0 1]);
legend('理论','蒙特卡洛仿真');
end
function A0=A_0(lambda_0,lambda_c,n,r)
A0=(lambda_0^(n-r)*lambda_c^r)/(gamma(n-r)*gamma(r));
end
function A1=A_1(lambda_0,lambda_c,n,r,i)
A1=(-1)^(r-1-i)*nchoosek(r-1,i)*factorial(n-i-2)*(lambda_0-lambda_c)^(-(n-i-1));
end
function A2=A_2(lambda_0,lambda_c,n,r)
A0=A_0(lambda_0,lambda_c,n,r);
A2=0;
for i=0:r-1
A1=A_1(lambda_0,lambda_c,n,r,i);
tmp=0;
for j=0:n-i-2
tmp=tmp+factorial(i+j)/factorial(j)*(lambda_0-lambda_c)^j*lambda_0^(-(i+j+1));
end
A2=A2+A1*(factorial(i)*lambda_c^(-(i+1))-tmp);
end
A2=A0*A2;
end
function A3=A_3(lambda_0,lambda_c,j)
A3=(lambda_0-lambda_c)^j/factorial(j);
end
function A4=A_4(lambda_c,m,k)
A4=factorial(m)*lambda_c^(-(m-k+1))/factorial(k);
end
function A5=A_5(lambda_0,lambda_c,p,q,m)
A5=factorial(m+p)*(lambda_0-lambda_c)^p*lambda_0^(-(m+p-q+1))/factorial(p)/factorial(q);
end
function I0=I_0(lambda_0,lambda_c,n,r1,r2,t)
A0=A_0(lambda_0,lambda_c,n,r1);
A0_=A_0(lambda_0,lambda_c,n,r2);
A2_=A_2(lambda_0,lambda_c,n,r2);
part1=0;
for i=0:r1-1
A1=A_1(lambda_0,lambda_c,n,r1,i);
tmp1=factorial(i)*(lambda_c+t)^(-(i+1));
tmp2=0;
for j=0:n-i-2
A3=A_3(lambda_0,lambda_c,j);
tmp2=tmp2+A3*factorial(i+j)*(lambda_0+t)^(-(i+j+1));
end
part1=part1+A1*(tmp1-tmp2);
end
part1=A0*A2_*part1;
part2=0;
for i=0:r1-1
A1=A_1(lambda_0,lambda_c,n,r1,i);
tmp1=0;
for m=0:r2-1
A1_=A_1(lambda_0,lambda_c,n,r2,m);
tmp2=0;
for k=0:m
A4=A_4(lambda_c,m,k);
tmp2=tmp2+A4*factorial(i+k)*(2*lambda_c+t)^(-(i+k+1));
end
tmp3=0;
for p=0:n-m-2
for q=0:m+p
A5=A_5(lambda_0,lambda_c,p,q,m);
tmp3=tmp3+A5*factorial(i+q)*(lambda_0+lambda_c+t)^(-(i+q+1));
end
end
tmp1=tmp1+A1_*(tmp3-tmp2);
end
part2=part2+A1*tmp1;
end
part2=A0*A0_*part2;
part3=0;
for i=0:r1-1
A1=A_1(lambda_0,lambda_c,n,r1,i);
tmp1=0;
for j=0:n-i-2
A3=A_3(lambda_0,lambda_c,j);
tmp2=0;
for m=0:r2-1
A1_=A_1(lambda_0,lambda_c,n,r2,m);
tmp3=0;
for k=0:m
A4=A_4(lambda_c,m,k);
tmp3=tmp3+A4*factorial(i+j+k)*(lambda_0+lambda_c+t)^(-(i+j+k+1));
end
tmp4=0;
for p=0:n-m-2
for q=0:m+p
A5=A_5(lambda_0,lambda_c,p,q,m);
tmp4=tmp4+A5*factorial(i+j+q)*(2*lambda_0+t)^(-(i+j+q+1));
end
end
tmp2=tmp2+A1_*(tmp3-tmp4);
end
tmp1=tmp1+A3*tmp2;
end
part3=part3+A1*tmp1;
end
part3=A0*A0_*part3;
I0=part1+part2+part3;
end
function I1=I_1(lambda_0,lambda_c,n,r1,r2,t)
A0=A_0(lambda_0,lambda_c,n,r1);
A0_=A_0(lambda_0,lambda_c,n,r2);
A2=A_2(lambda_0,lambda_c,n,r1);
part1=0;
for i=0:r2-1
A1_=A_1(lambda_0,lambda_c,n,r2,i);
tmp1=factorial(i)*(lambda_c+t)^(-(i+1));
tmp2=0;
for j=0:n-i-2
A3=A_3(lambda_0,lambda_c,j);
tmp2=tmp2+A3*factorial(i+j)*(lambda_0+t)^(-(i+j+1));
end
part1=part1+A1_*(tmp1-tmp2);
end
part1=A0_*A2*part1;
part2=0;
for i=0:r2-1
A1_=A_1(lambda_0,lambda_c,n,r2,i);
tmp1=0;
for m=0:r1-1
A1=A_1(lambda_0,lambda_c,n,r1,m);
tmp2=0;
for k=0:m
A4=A_4(lambda_c,m,k);
tmp2=tmp2+A4*factorial(i+k)*(2*lambda_c+t)^(-(i+k+1));
end
tmp3=0;
for p=0:n-m-2
for q=0:m+p
A5=A_5(lambda_0,lambda_c,p,q,m);
tmp3=tmp3+A5*factorial(i+q)*(lambda_0+lambda_c+t)^(-(i+q+1));
end
end
tmp1=tmp1+A1*(tmp3-tmp2);
end
part2=part2+A1_*tmp1;
end
part2=A0_*A0*part2;
part3=0;
for i=0:r2-1
A1_=A_1(lambda_0,lambda_c,n,r2,i);
tmp1=0;
for j=0:n-i-2
A3=A_3(lambda_0,lambda_c,j);
tmp2=0;
for m=0:r1-1
A1=A_1(lambda_0,lambda_c,n,r1,m);
tmp3=0;
for k=0:m
A4=A_4(lambda_c,m,k);
tmp3=tmp3+A4*factorial(i+j+k)*(lambda_0+lambda_c+t)^(-(i+j+k+1));
end
tmp4=0;
for p=0:n-m-2
for q=0:m+p
A5=A_5(lambda_0,lambda_c,p,q,m);
tmp4=tmp4+A5*factorial(i+j+q)*(2*lambda_0+t)^(-(i+j+q+1));
end
end
tmp2=tmp2+A1*(tmp3-tmp4);
end
tmp1=tmp1+A3*tmp2;
end
part3=part3+A1_*tmp1;
end
part3=A0_*A0*part3;
I1=part1+part2+part3;
end
% 对应r1=r2=0或r1=r2=n的情况
function Pd=func_GoPd_r1(lambda,lambda_t,n,alpha)
t=alpha*lambda_t;
A=2*(1+t/lambda)^(-n);
B=0;
for k=0:n-1
B=B+nchoosek(n+k-1,k)*(2+t/lambda)^(-(n+k));
end
Pd=A-2*B;
end
% 对应1<=r1<=n-1,1<=r2<=n-1的情况
function Pd=func_GoPd_r2(lambda_0,lambda_c,lambda_t,n,r1,r2,alpha)
Pd=I_0(lambda_0,lambda_c,n,r1,r2,alpha*lambda_t)+I_1(lambda_0,lambda_c,n,r1,r2,alpha*lambda_t);
end
% 对应r1=0,1<=r2<=n-1的情况
function Pd=func_GoPd_r3(lambda_0,lambda_c,lambda_t,n,r1,r2,alpha)
t=alpha*lambda_t;
A0=A_0(lambda_0,lambda_c,n,r2);
A2=A_2(lambda_0,lambda_c,n,r2);
if r1==0
lambda=lambda_0;
else
lambda=lambda_c;
end
I0=0;
for m=0:r2-1
A1=A_1(lambda_0,lambda_c,n,r2,m);
tmp1=0;
for k=0:m
tmp1=tmp1+nchoosek(n+k-1,k)*lambda_c^(-(m-k+1))*(lambda+lambda_c+t)^(-(n+k));
end
tmp1=tmp1*factorial(m);
tmp2=0;
for p=0:n-m-2
for q=0:m+p
tmp2=tmp2+factorial(m+p)/factorial(p)*nchoosek(n+q-1,q)*(lambda_0-lambda_c)^p*lambda_0^(-(m+p-q+1))*(lambda+lambda_0+t)^(-(n+q));
end
end
I0=I0+A1*(tmp2-tmp1);
end
I0=A0*lambda^n*I0+A2*(1+t/lambda)^(-n);
I1=0;
for i=0:r2-1
A1=A_1(lambda_0,lambda_c,n,r2,i);
tmp1=factorial(i)*(lambda_c+t)^(-(i+1));
tmp2=0;
for j=0:n-i-2
tmp2=tmp2+factorial(i+j)/factorial(j)*(lambda_0-lambda_c)^j*(lambda_0+t)^(-(i+j+1));
end
tmp3=0;
for k=0:n-1
tmp3=tmp3+factorial(k+i)/factorial(k)*lambda^k*(lambda+lambda_c+t)^(-(k+i+1));
end
tmp4=0;
for j=0:n-i-2
tmp5=0;
for k=0:n-1
tmp5=tmp5+factorial(i+j+k)/factorial(k)*lambda^k*(lambda+lambda_0+t)^(-(i+j+k+1));
end
tmp4=tmp4+(lambda_0-lambda_c)^j/factorial(j)*tmp5;
end
I1=I1+A0*A1*(tmp1-tmp2-tmp3+tmp4);
end
Pd=I0+I1;
end
% 对应r1=0,r2=n的情况
function Pd=func_GoPd_r4(lambda_0,lambda_c,lambda_t,n,alpha)
t=alpha*lambda_t;
tmp1=(1+t/lambda_0)^(-n);
tmp2=(1+t/lambda_c)^(-n);
tmp3=0;
for k=0:n-1
tmp3=tmp3+nchoosek(n+k-1,k)*lambda_c^k*(lambda_0+lambda_c+t)^(-(n+k));
end
tmp3=tmp3*lambda_0^n;
tmp4=0;
for k=0:n-1
tmp4=tmp4+nchoosek(n+k-1,k)*lambda_0^k*(lambda_0+lambda_c+t)^(-(n+k));
end
tmp4=tmp4*lambda_c^n;
Pd=tmp1+tmp2-tmp3-tmp4;
end
% 根据输入参数情况调用对应情况下的函数
function Pd=func_GoPd(lambda_0,lambda_c,lambda_t,n,r1,r2,alpha)
r_min=min(r1,r2);
r_max=max(r1,r2);
if r1==0 & r2==0
Pd=func_GoPd_r1(n*lambda_0,lambda_t,n,alpha);
elseif r1==n & r2==n
Pd=func_GoPd_r1(n*lambda_c,lambda_t,n,alpha);
elseif r_min==0 & r_max==n
Pd=func_GoPd_r4(n*lambda_0,n*lambda_c,lambda_t,n,alpha);
elseif (r_min==0 & r_max~=n)|(r_max==n & r_min~=0)
Pd=func_GoPd_r3(n*lambda_0,n*lambda_c,lambda_t,n,r1,r2,alpha);
else
Pd=func_GoPd_r2(n*lambda_0,n*lambda_c,lambda_t,n,r1,r2,alpha);
end
end
