均匀背景下常用恒虚警检测推导及比较(二)：附录+代码验证

各均衡器的主体理论推导过程见《均匀背景下常用恒虚警检测推导及比较(一)：理论推导部分》

附录

(1)CA-CFAR中间推导过程

已知$f_{x_i}(x_i)={\lambda }_0{e}^{-{\lambda }_0{x}_i},z=\frac{1}{N}{\sum }_{i=1}^N{x}_i$，求$z$的概率密度函数$f_z(z)$。

$f_{x_i}(x_i)$对应的特征函数计算如下：
$${\varphi }_{x_i}(t)={\int }_0^{\infty }{e}^{it{x}_i}{f}_{x_i}(x_i)d{x}_i={\lambda }_0{\int }_0^{\infty }{e}^{(it-{\lambda }_0)x_i}d{x}_i=(1-\frac{it}{{\lambda }_0}{)}^{-1}=\frac{{\lambda }_0}{{\lambda }_0-it}\text{\hspace{1em}(1)}$$
令$z^{\prime }={\sum }_{i=1}^N{x}_i$，则$f_{z^{\prime }}(z^{\prime })$对应的特征函数为${\varphi }_{z^{\prime }}(t)={\varphi }_{x_i}^N(t)=(\frac{{\lambda }_0}{{\lambda }_0-it}{)}^N$。所以
$$f_{z^{\prime }}(z^{\prime })=\frac{1}{2\pi }{\int }_{-\infty }^{\infty }{\varphi }_{z^{\prime }}(t)e^{-it{z}^{\prime }}dt=\frac{1}{2\pi }{\int }_{-\infty }^{\infty }(\frac{{\lambda }_0}{{\lambda }_0-it}{)}^N{e}^{-it{z}^{\prime }}dt\text{\hspace{1em}(2)}$$
上面的积分式并不好求，可以通过查找积分表获得结果。另外地，可以利用一些现有结论得到想要的结果。若随机变量$X$服从Gamma分布，记为$X\backsim \mathrm{G}(n,\lambda )$，则该随机变量的概率密度函数可以表示为
$$X\backsim \mathrm{G}(n,\lambda ):f(x)=\frac{{\lambda }^n}{\Gamma (n)}{x}^{n-1}{e}^{-\lambda x}\text{\hspace{1em}(3)}$$
显然指数分布是当$n=1$时特殊Gamma分布，即若$X\backsim E(\lambda )\to X\backsim \mathrm{G}(1,\lambda )$。当$n$为整数时$\Gamma (n)=(n-1)!$。指数分布和Gamma分布具有如下关系：若$X_1,X_2,...,X_n\backsim E(\lambda )$，则$z={\sum }_{i=1}^n{x}_i\backsim \mathrm{G}(n,\lambda )$。

利用上述结论，显然有(2)中$z^{\prime }\backsim \mathrm{G}(N,{\lambda }_0)$，所以
$$z^{\prime }\backsim \mathrm{G}(N,{\lambda }_0):f_{z^{\prime }}(z^{\prime })=\frac{{\lambda }_0^N}{\Gamma (N)}(z^{\prime }{)}^{N-1}{e}^{-{\lambda }_0{z}^{\prime }}\text{\hspace{1em}(4)}$$
由于$z=\frac{1}{N}{z}^{\prime }$，所以$z\backsim \mathrm{G}(N,N{\lambda }_0)$，所以$z$的概率密度函数可以表示为
$$z\backsim \mathrm{G}(N,N{\lambda }_0):f_z(z)=\frac{(N{\lambda }_0{)}^N}{\Gamma (N)}{z}^{N-1}{e}^{-N{\lambda }_{0}z}\text{\hspace{1em}(5)}$$
(2)Gamma分布的矩母函数

若$X\backsim \mathrm{G}(n,\lambda ),(n,\lambda >0)$，则$f(x)=\frac{{\lambda }^n}{\Gamma (n)}{x}^{n-1}{e}^{-\lambda x}$，其中$\Gamma (n)={\int }_0^{\infty }{t}^{n-1}{e}^{-t}dt,n>0$，则其矩母函数可以表示为
$$\begin{array}{c}\begin{array}{r}M(t)=\mathrm{E}[e^{-tx}]={\int }_0^{\infty }{e}^{-tx}f(x)dx={\int }_0^{\infty }\frac{{\lambda }^n}{\Gamma (n)}{x}^{n-1}{e}^{-\lambda x}{e}^{-tx}dx=\frac{{\lambda }^n}{\Gamma (n)}{\int }_0^{\infty }{x}^{n-1}{e}^{(t+\lambda )x}dx\\ \underrightarrow{(t+\lambda )x=-v}\frac{1}{\Gamma (n)}(\frac{\lambda }{\lambda +t}{)}^n{\int }_0^{\infty }{v}^{n-1}{e}^{-v}dv=(\frac{\lambda }{\lambda +t}{)}^n=(1+\frac{t}{\lambda }{)}^{-n}\end{array}\end{array}\text{\hspace{1em}(6)}$$
(3)不完全Gamma函数的相关常用关系

原始Gamma、上不完全Gamma函数及下不完全Gamma函数的定义分别如下所示
$$\begin{gathered} \Gamma (n)={\int }_0^{\infty }{t}^{n-1}{e}^{-t}dt \\ \Gamma (n,x)={\int }_x^{\infty }{t}^{n-1}{e}^{-t}dt \\ \gamma (n,x)={\int }_0^x{t}^{n-1}{e}^{-t}dt\text{\hspace{1em}(1)} \end{gathered}$$
他们三者满足下述关系
$$\gamma (n,\lambda )+\Gamma (n,\lambda )=\Gamma (n)\text{\hspace{1em}(2)}$$
当$n$为整数时，$\Gamma (n,\lambda )$可以表示为$\Gamma (n,\lambda )=(n-1)!e^{-\lambda }{\sum }_{k=0}^{n-1}\frac{{\lambda }^k}{k!}$，$\Gamma (n)=(n-1)!$，所以$\gamma (n,\lambda )=(n-1)!(1-e^{-\lambda }{\sum }_{k=0}^{n-1}\frac{{\lambda }^k}{k!})$。

同时下不完全Gamma函数可用级数展开进行表示
$$\gamma (n,\lambda )=\sum _{k=0}^{\infty }\frac{{\lambda }^n{e}^{-\lambda }{\lambda }^k}{s(s+1)...(s+k)}={\lambda }^n\Gamma (n)e^{-\lambda }\sum _{k=0}^{\infty }\frac{{\lambda }^k}{\Gamma (n+k+1)}\text{\hspace{1em}(3)}$$

(4)Gamma分布的分布函数

若$X\backsim \mathrm{G}(n,\lambda ),(n,\lambda >0)$，则$f(x)=\frac{{\lambda }^n}{\Gamma (n)}{x}^{n-1}{e}^{-\lambda x}$，其中$\Gamma (n)={\int }_0^{\infty }{t}^{n-1}{e}^{-t}dt,n>0$，则其分布函数可以表示为
$$F(x)=\frac{1}{\Gamma (n)}\gamma (n,\lambda x)=1-e^{-\lambda x}\sum _{k=0}^{n-1}\frac{(\lambda x{)}^k}{k!}=e^{-\lambda x}\sum _{k=n}^{\infty }\frac{(\lambda x{)}^k}{k!}\text{\hspace{1em}(1)}$$
在上面的表达式中$\gamma (n,\lambda x)$表示下不完全Gamma函数，它的表达式可以参考附录(3)。

(5)GO-CFAR中矩母函数$M_z(t)$的求解

在进行推导之前，首先给出如下表达式，后面的推导过程要利用到该结论
$${\int }_0^b{y}^m{e}^{-qy}dy=m!q^{-(m+1)}[1-\sum _{k=0}^m(yq{)}^k{e}^{-qy}/k!]\underrightarrow{b\to \infty }{\int }_0^b{y}^m{e}^{-qy}dy=m!q^{-(m+1)}\text{\hspace{1em}(1)}$$
如第四节公式(9)所示，$z$的概率密度函数可以表示为
$$f_z(z)=2\frac{(N{\lambda }_0/2{)}^{N/2}}{\Gamma (N/2)}{z}^{N/2-1}{e}^{-N{\lambda }_{0}z/2}[1-e^{-N{\lambda }_{0}z/2}\sum _{k=0}^{N/2-1}\frac{(N{\lambda }_{0}z/2{)}^k}{k!}]\text{\hspace{1em}(2)}$$
为了表述简洁，令$n=N/2$，则(1)可以表示为
$$f_z(z)=2\frac{(n{\lambda }_0{)}^n}{\Gamma (n)}{z}^{n-1}{e}^{-n{\lambda }_{0}z}[1-e^{-n{\lambda }_{0}z}\sum _{k=0}^{n-1}\frac{(n{\lambda }_{0}z{)}^k}{k!}]\text{\hspace{1em}(3)}$$
则
$$\begin{array}{c}\begin{array}{rl}{\mathrm{M}}_z(t)& ={\int }_0^{\infty }{f}_z(z)e^{-tz}dz\\ & =2\frac{(n{\lambda }_0{)}^n}{\Gamma (n)}{\int }_0^{\infty }{z}^{n-1}{e}^{-n{\lambda }_{0}z}{e}^{-tz}dz-2\frac{(n{\lambda }_0{)}^n}{\Gamma (n)}{\int }_0^{\infty }{z}^{n-1}{e}^{-2n{\lambda }_{0}z}{e}^{-tz}\sum _{k=0}^{n-1}\frac{(n{\lambda }_{0}z{)}^k}{k!}dz\end{array}\end{array}\text{\hspace{1em}(4)}$$
其中
$$\begin{array}{c}\begin{array}{rl}A& =2\frac{(n{\lambda }_0{)}^n}{\Gamma (n)}{\int }_0^{\infty }{z}^{n-1}{e}^{-n{\lambda }_{0}z}{e}^{-tz}dz=2\frac{(n{\lambda }_0{)}^n}{\Gamma (n)}{\int }_0^{\infty }{z}^{n-1}{e}^{-(n{\lambda }_0+t)z}dz\\ & \underrightarrow{m=n-1,q=n{\lambda }_0+t}\\ A& =2\frac{(n{\lambda }_0{)}^n}{\Gamma (n)}{\int }_0^{\infty }{z}^m{e}^{-qz}dz=2\frac{(n{\lambda }_0{)}^n}{\Gamma (n)}m!q^{-(m+1)}\\ & =2\frac{(n{\lambda }_0{)}^n}{\Gamma (n)}(n-1)!(n{\lambda }_0+t{)}^{-n}=2(n{\lambda }_0{)}^n(n{\lambda }_0+t{)}^{-n}\end{array}\end{array}\text{\hspace{1em}(5)}$$

$$\begin{array}{c}\begin{array}{rl}B& =2\frac{(n{\lambda }_0{)}^n}{\Gamma (n)}{\int }_0^{\infty }{z}^{n-1}{e}^{-2n{\lambda }_{0}z}{e}^{-tz}\sum _{k=0}^{n-1}\frac{(n{\lambda }_{0}z{)}^k}{k!}dz=2\frac{(n{\lambda }_0{)}^n}{\Gamma (n)}\sum _{k=0}^{n-1}[\frac{(n{\lambda }_0{)}^k}{k!}{\int }_0^{\infty }{z}^{n-1+k}{e}^{-(2n{\lambda }_0+t)z}dz]\\ & \underrightarrow{m=n-1+k,q=2n{\lambda }_0+t}\\ B& =2\frac{(n{\lambda }_0{)}^n}{\Gamma (n)}\sum _{k=0}^{n-1}[\frac{(n{\lambda }_0{)}^k}{k!}{\int }_0^{\infty }{z}^m{e}^{-qz}dz]=2\frac{(n{\lambda }_0{)}^n}{\Gamma (n)}\sum _{k=0}^{n-1}[\frac{(n{\lambda }_0{)}^k}{k!}(n-1+k)!(2n{\lambda }_0+t{)}^{-(n+k)}]\\ & =2(n{\lambda }_0{)}^n\sum _{k=0}^{n-1}\left(\begin{array}{c}n-1+k\\ k\end{array}\right)(n{\lambda }_0{)}^k(2n{\lambda }_0+t{)}^{-(n+k)}\end{array}\end{array}\text{\hspace{1em}(6)}$$

将(5)(6)代入(4)可得
$${\mathrm{M}}_z(t)=A-B=2(n{\lambda }_0{)}^n(n{\lambda }_0+t{)}^{-n}-2(n{\lambda }_0{)}^n\sum _{k=0}^{n-1}\left(\begin{array}{c}n-1+k\\ k\end{array}\right)(n{\lambda }_0{)}^k(2n{\lambda }_0+t{)}^{-(n+k)}\text{\hspace{1em}(7)}$$

(6)次序统计量的概率密度函数推导

随机变量的分布函数、概率、概率密度函数之间存在一定的关系，设随机变量$X$的pdf和cdf分别为$F(x),f(x)$，则有
$$F(x)=P(X\le x)={\int }_0^{x}f(x)dx\text{\hspace{1em}(1)}$$
这个表达式是直接根据$F(x)$和$f(x)$的定义得到的。所以
$$\begin{array}{c}\begin{array}{r}P(x\le X\le x+\Delta x)=F(x+\Delta x)-F(x)\to \\ \underset{\Delta x\to 0}{\lim }\frac{P(x\le X\le x+\Delta x)}{\Delta x}=P(X=x)=\underset{\Delta x\to 0}{\lim }\frac{F(x+\Delta x)-F(x)}{\Delta x}=f(x)dx\end{array}\end{array}\text{\hspace{1em}(2)}$$
上面的推导过程得到了在满足连续分布的随机变量等于特定值的概率，它将被应用到后续的推导过程中。

对于一个有序的序列$x_{(1)}\le x_{(2)}\le ...\le x_{(N)}$，现在要求其第$k$个元素的概率密度函数$f(y)$。显然对于事件$x_{(k)}=y$来说，其等价的命题是：序列中存在$k-1$个元素的值小于$y$，存在一个元素的值等于$y$，存在$N-k$个元素值大于$y$。上述等价命题其实包含了3小部分的内容，以下一一进行处理：

序列中存在$k-1$个元素的值小于$y$，这一事件对应的概率可以表示为$\left(\begin{array}{c}N\\ k-1\end{array}\right)[P(X\le y){]}^{k-1}=\left(\begin{array}{c}N\\ k-1\end{array}\right)[F(y){]}^{k-1}$；(从原始$N$个样本中随机选择$k-1$个，使它们的值全部小于$y$)

序列中存在一个元素的值等于$y$，这一事件对应的概率可以表示为$\left(\begin{array}{c}N-k+1\\ 1\end{array}\right)P(X=y)=\left(\begin{array}{c}N-k+1\\ 1\end{array}\right)f(y)dy$；(从剩余的$N-k+1$个样本中随机选择1个样本，让它的值等于$y$)

序列中存在$N-k$个元素的值大于$y$，这一事件对应的概率可以表示为$C_{N-k}^{N-k}[P(X>y){]}^{N-k}=\left(\begin{array}{c}N-k\\ N-k\end{array}\right)[1-P(X\le y){]}^{N-k}=\left(\begin{array}{c}N-k\\ N-k\end{array}\right)[1-F(y){]}^{N-k}$

综合上面三个事件可得
$$\begin{array}{c}\begin{array}{rl}P(x_{(k)}=y)& =f_{x_{(k)}}(y)dy=\left(\begin{array}{c}N\\ k-1\end{array}\right)[F(y){]}^{k-1}\left(\begin{array}{c}N-k+1\\ 1\end{array}\right)f(y)dy\left(\begin{array}{c}N-k\\ N-k\end{array}\right)[1-F(y){]}^{N-k}\to \\ f_{x_{(k)}}(y)& =\frac{N!}{(k-1)!(N-k)!}[F(y){]}^{k-1}[1-F(y){]}^{N-k}f(y)\end{array}\end{array}\text{\hspace{1em}(3)}$$
(7)OS-CFAR中统计量对应的矩母函数推导

统计量$z$的概率密度函数为
$$f_z(z)={\lambda }_0\frac{N!}{(k-1)!(N-k)!}{e}^{-{\lambda }_0(N-k+1)z}[1-e^{-{\lambda }_{0}z}{]}^{k-1}\text{\hspace{1em}(1)}$$
其矩母函数可以表示为
$${\mathrm{M}}_z(t)={\int }_0^{\infty }{f}_z(z)e^{-tz}dz={\lambda }_0\frac{N!}{(k-1)!(N-k)!}{\int }_0^{\infty }{e}^{-{\lambda }_0(N-k+1)z}[1-e^{-{\lambda }_{0}z}{]}^{k-1}{e}^{-tz}dz\text{\hspace{1em}(2)}$$
利用二项展开式$(a+b{)}^n={\sum }_{r=0}^n\left(\begin{array}{c}n\\ r\end{array}\right)a^{n-r}{b}^r$，有$[1-e^{-{\lambda }_{0}z}{]}^{k-1}={\sum }_{r=0}^{k-1}(-1{)}^r\left(\begin{array}{c}k-1\\ r\end{array}\right)e^{-{\lambda }_{0}rz}$，将其代入(2)中并进行整理可得
$$\begin{array}{c}\begin{array}{rl}{\mathrm{M}}_z(t)& ={\lambda }_0\frac{N!}{(k-1)!(N-k)!}\sum _{r=0}^{k-1}[(-1{)}^r\left(\begin{array}{c}k-1\\ r\end{array}\right){\int }_0^{\infty }{e}^{-[{\lambda }_0(N-k+1+r)+t]z}dz]\\ & ={\lambda }_0\frac{N!}{(k-1)!(N-k)!}\sum _{r=0}^{k-1}\frac{(-1{)}^r\left(\begin{array}{c}k-1\\ r\end{array}\right)}{{\lambda }_0(N-k+1+r)+t}={\lambda }_0(N-k+1)\sum _{r=0}^{k-1}\frac{(-1{)}^r\left(\begin{array}{c}k-1\\ r\end{array}\right)\left(\begin{array}{c}N\\ k-1\end{array}\right)}{{\lambda }_0(N-k+1+r)+t}\end{array}\end{array}\text{\hspace{1em}(3)}$$
上面的表达式感觉还没有化简到最简形式，但是目前还不知道怎么往下继续化简。

通过查阅参考文献[10]，发现如下积分公式
$${\int }_0^{\infty }(1-e^{-x/\beta }{)}^{v-1}{e}^{-\mu x}dx=\beta \mathrm{B}(\beta \mu ,v)\text{\hspace{1em}(4)}$$
其中$\mathrm{B}(x,y)$表示beta函数，且$\mathrm{B}(x,y)={\int }_0^1{t}^{x-1}(1-t{)}^{y-1}dt$，当$x,y$均为正整数时，$\mathrm{B}(x,y)=\frac{(x-1)!(y-1)!}{(x+y-1)!}=\frac{\Gamma (x)\Gamma (y)}{\Gamma (x+y)}$，利用(4)所示积分，对(2)进行重新整理
$$\begin{array}{c}\begin{array}{rl}{\mathrm{M}}_z(t)& ={\lambda }_0\frac{N!}{(k-1)!(N-k)!}{\int }_0^{\infty }(1-e^{-{\lambda }_{0}z}{)}^{k-1}{e}^{-[{\lambda }_0(N-k+1)+t]z}dz\\ & \underrightarrow{\mu ={\lambda }_0(N-k+1)+t,v=k,\beta =1/{\lambda }_0}\\ {\mathrm{M}}_z(t)& ={\lambda }_0\frac{N!}{(k-1)!(N-k)!}\beta \mathrm{B}(\beta \mu ,v)=\frac{N!}{(k-1)!(N-k)!}\mathrm{B}(N-k+1+t/{\lambda }_0,k)\\ & =k\left(\begin{array}{c}N\\ k\end{array}\right)\mathrm{B}(N-k+1+t/{\lambda }_0,k)\end{array}\end{array}\text{\hspace{1em}(5)}$$

(8)利用matlab生成指定均值和方差的正态分布随机变量

在上面的推导过程中，我们均假设背景噪声服从高斯分布，因此平方率检波之后噪声背景服从指数分布，关于高斯分布、指数分布、瑞利分布之间的关系，可以参考之前写的博客随机变量概率密度函数和概率分布函数相关总结。在CFAR检测器的仿真研究过程中，产生特定功率的背景噪声信号，以及特定信噪比的目标信号是很重要的(需要特别注意的是，这边所说的信噪比，是指原始信号中的目标信噪比，而不是平方率检波后的目标的信噪比，所以仿真信号的产生最好直接从高斯分布的原始信号开始，再由这个原始信号的模值平方生成平方率检波后的信号)。以下就如何按要求产生特定背景噪声及目标信号进行说明。在说明之前，先看一些结论：
$$\begin{array}{c}\begin{array}{rl}& X\backsim N(\mu ,{\sigma }^2)\to Y=aX+b\backsim N(a\mu +b,a^2{\sigma }^2)\\ & X\backsim N(0,{\sigma }^2),Y\backsim N(0,{\sigma }^2)\to Z=X^2+Y^2\backsim \mathrm{Exp}(\lambda )=\frac{1}{\lambda }{e}^{-z/\lambda }=\frac{1}{2{\sigma }^2}{e}^{-z/2{\sigma }^2}\\ & X\backsim \mathrm{Exp}(\lambda )\to \mathrm{E}(x)=\lambda ,D(x)={\lambda }^2\end{array}\end{array}\text{\hspace{1em}(1)}$$
设要求的背景噪声功率和目标信噪比分别为$nP$ dB和$SNR$ dB，则仿真代码如下：

nP=10;    %背景噪声功率(dB)
SNR=10;   %目标信噪比(dB)
N=10000;  %信号长度

sigma=sqrt(10^(nP/10))/sqrt(2);  %单通道噪声标准差
noise=sigma*(randn(N,1)+1j*randn(N,1));  %生成特定功率的噪声信号
trgt=10^((SNR+nP)/20;   %目标幅值

(9)各CFAR算法门限乘积因子、虚警概率、检测概率的计算代码

% CA-CFAR虚警概率计算函数(直接用公式计算)
function [PFA] = func_CaPfa(alpha,N)
    PFA=(1+alpha/N)^(-N);
end

% OS-CFAR虚警概率计算函数(直接用公式计算)
function [PFA] = func_OsPfa(alpha,N,k)
   PFA=k*nchoosek(N,k)*gamma(N-k+1+alpha)*gamma(k)/gamma(N+1+alpha);
end

% GO-CFAR虚警概率计算函数(直接用公式计算)
function [PFA] = func_GoPfa(alpha,N)
  n=N/2;
  A=2*(1+alpha/n)^(-n);
  B=0;
  for k=0:n-1
      B=B+2*nchoosek(n-1+k,k)*(2+alpha/n)^(-(n+k));
  end
  PFA=A-B;
end

% SO-CFAR虚警概率计算函数(直接用公式计算)
function [PFA] = func_SoPfa(alpha,N)
  n=N/2;
  PFA=0;
  for k=0:n-1
      PFA=PFA+2*nchoosek(n-1+k,k)*(2+alpha/n)^(-(n+k));
  end
end

% CA-CFAR检测概率计算函数(直接用公式计算)
function [Pd] = func_CaPd(alpha,N,snr)
  snr=10.^(snr/10); %将snr转化成线性
  Pd=(1+alpha./(N*(1+snr))).^(-N);
end

% OS-CFAR检测概率计算函数(直接用公式计算)
function [Pd] = func_OsPd(alpha,N,k,snr)
  snr=10.^(snr/10); %将snr转化成线性
  Pd=k*nchoosek(N,k)*gamma(N-k+1+alpha./(1+snr))*gamma(k)./gamma(N+1+alpha./(1+snr));
end

% GO-CFAR检测概率计算函数(直接用公式计算)
function [Pd] = func_GoPd(alpha,N,snr)
  snr=10.^(snr/10); %将snr转化成线性
  n=N/2;
  A=2*(1+alpha./(n*(1+snr))).^(-n);
  B=0;
  for k=0:n-1
      B=B+2*nchoosek(n-1+k,k)*(2+alpha./(n*(1+snr))).^(-(n+k));
  end
  Pd=A-B;
end

% SO-CFAR检测概率计算函数(直接用公式计算)
function [Pd] = func_SoPd(alpha,N,snr)
  snr=10.^(snr/10); %将snr转化成线性
  n=N/2;
  Pd=0;
  for k=0:n-1
      Pd=Pd+2*nchoosek(n-1+k,k)*(2+alpha./(n*(1+snr))).^(-(n+k));
  end
end

% CA-CFAR门限乘积因子计算函数(直接用公式计算)
function [alpha] = func_CaAlpha(PFA,N)
    alpha=N*(PFA^(-1/N)-1);
end

% OS-CFAR门限乘积因子计算函数(构造多项式计算)
function [alpha] = func_OsAlpha(PFA,N,k)
    syms T;
    x=vpasolve(PFA-k*nchoosek(N,k)*gamma(N-k+1+T)*gamma(k)/gamma(N+T+1)==0,T);
    alpha=x(x==abs(x));
end

% GO-CFAR门限乘积因子计算函数(用二分法逼近)
function [alpha] = func_GoAlpha(PFA,N)
   scope=[0,500];    %预设门限乘积因子的范围
   res=1;            %预设门限乘积因子的计算误差
   while 1
     alpha=mean(scope);
     PFA_=func_GoPfa(alpha,N);   %根据alpha值计算虚警率
     d=1/PFA_-1/PFA;             %理想PFA与计算得到的PFA之前的差异
     if abs(d)<res || abs(scope(1)-scope(2))<0.0001
         return;
     elseif d<0
         scope(1)=alpha;
     else
         scope(2)=alpha;
     end
  end
end

% SO-CFAR门限乘积因子计算函数(用二分法逼近)
function [alpha] = func_SoAlpha(PFA,N)
   scope=[0,500];    %预设门限乘积因子的范围
   res=1;            %预设门限乘积因子的计算误差
   while 1
     alpha=mean(scope);
     PFA_=func_SoPfa(alpha,N);   %根据alpha值计算虚警率
     d=1/PFA_-1/PFA;             %理想PFA与计算得到的PFA之前的差异
     if abs(d)<res || abs(scope(1)-scope(2))<0.0001
         return;
     elseif d<0
         scope(1)=alpha;
     else
         scope(2)=alpha;
     end
   end
end

上述代码已与matlab自带的CFAR算法结果进行比对，证明了代码的正确性。

(10)各种CFAR算法的实现代码

% 函数说明:实现CA-CFAR
% 输入参数:
%    (1)x:原始样本，暂时只支持一维向量
%    (2)cutI:待检单元索引
%    (3)refN:参考窗长度
%    (4)prctN:保护窗长度
%    (5)PFA:虚警概率
% 输出参数:
%    (1)flag:该检测单元上是否存在目标，0表示不存在，1表示存在
%    (2)alpha:门限乘积因子
%    (3)threshold:检测门限
% 版本历史:
% 20210802:程序创建，只支持一维输入向量，且不考虑边缘情况
function [flag,alpha,threshold] = func_CaCFAR(x,cutI,refN,prctN,PFA)
    x=x(:);
    if cutI-refN/2-prctN/2<=0 | cutI+refN/2+prctN/2>length(x)
        error('出现边缘情况，请检查函数输入参数!');
    end

    leading=x(cutI-refN/2-prctN/2:cutI-prctN/2-1);  %前参考窗
    lagging=x(cutI+prctN/2+1:cutI+prctN/2+refN/2);  %后参考窗
    cut=x(cutI);                                    %待检单元

    alpha=func_CaAlpha(PFA,refN);            %计算门限乘积因子
    noiseP=mean([leading;lagging]);          %计算噪声功率
    threshold=alpha*noiseP;                  %计算门限
    flag=(cut>=threshold);                   %待检单元与检测门限对比获得检测结果
end

% 函数说明:实现GO-CFAR
% 输入参数:
%    (1)x:原始样本，暂时只支持一维向量
%    (2)cutI:待检单元索引
%    (3)refN:参考窗长度
%    (4)prctN:保护窗长度
%    (5)PFA:虚警概率
% 输出参数:
%    (1)flag:该检测单元上是否存在目标，0表示不存在，1表示存在
%    (2)alpha:门限乘积因子
%    (3)threshold:检测门限
% 版本历史:
% 20210802:程序创建，只支持一维输入向量，且不考虑边缘情况
function [flag,alpha,threshold] = func_GoCFAR(x,cutI,refN,prctN,PFA)
    x=x(:);
    if cutI-refN/2-prctN/2<=0 | cutI+refN/2+prctN/2>length(x)
        error('出现边缘情况，请检查函数输入参数!');
    end

    leading=x(cutI-refN/2-prctN/2:cutI-prctN/2-1);  %前参考窗
    lagging=x(cutI+prctN/2+1:cutI+prctN/2+refN/2);  %后参考窗
    cut=x(cutI);                                    %待检单元

    alpha=func_GoAlpha(PFA,refN);            %计算门限乘积因子
    noiseP=max(mean(leading),mean(lagging)); %计算噪声功率
    threshold=alpha*noiseP;                  %计算门限
    flag=(cut>=threshold);                   %待检单元与检测门限对比获得检测结果
end

% 函数说明:实现SO-CFAR
% 输入参数:
%    (1)x:原始样本，暂时只支持一维向量
%    (2)cutI:待检单元索引
%    (3)refN:参考窗长度
%    (4)prctN:保护窗长度
%    (5)PFA:虚警概率
% 输出参数:
%    (1)flag:该检测单元上是否存在目标，0表示不存在，1表示存在
%    (2)alpha:门限乘积因子
%    (3)threshold:检测门限
% 版本历史:
% 20210802:程序创建，只支持一维输入向量，且不考虑边缘情况
function [flag,alpha,threshold] = func_SoCFAR(x,cutI,refN,prctN,PFA)
    x=x(:);
    if cutI-refN/2-prctN/2<=0 | cutI+refN/2+prctN/2>length(x)
        error('出现边缘情况，请检查函数输入参数!');
    end
    
    leading=x(cutI-refN/2-prctN/2:cutI-prctN/2-1);  %前参考窗
    lagging=x(cutI+prctN/2+1:cutI+prctN/2+refN/2);  %后参考窗
    cut=x(cutI);                                    %待检单元

    alpha=func_SoAlpha(PFA,refN);            %计算门限乘积因子
    noiseP=min(mean(leading),mean(lagging)); %计算噪声功率
    threshold=alpha*noiseP;                  %计算门限
    flag=(cut>=threshold);                   %待检单元与检测门限对比获得检测结果
end

% 函数说明:实现OS-CFAR
% 输入参数:
%    (1)x:原始样本，暂时只支持一维向量
%    (2)cutI:待检单元索引
%    (3)refN:参考窗长度
%    (4)prctN:保护窗长度
%    (5)k:噪声估计单元索引
%    (5)PFA:虚警概率
% 输出参数:
%    (1)flag:该检测单元上是否存在目标，0表示不存在，1表示存在
%    (2)alpha:门限乘积因子
%    (3)threshold:检测门限
% 版本历史:
% 20210802:程序创建，只支持一维输入向量，且不考虑边缘情况
function [flag,alpha,threshold] = func_OsCFAR(x,cutI,refN,prctN,k,PFA)
    x=x(:);
    if cutI-refN/2-prctN/2<=0 | cutI+refN/2+prctN/2>length(x)
        error('出现边缘情况，请检查函数输入参数!');
    end

    leading=x(cutI-refN/2-prctN/2:cutI-prctN/2-1);  %前参考窗
    lagging=x(cutI+prctN/2+1:cutI+prctN/2+refN/2);  %后参考窗
    cut=x(cutI);                                    %待检单元

    alpha=func_OsAlpha(PFA,refN,k);          %计算门限乘积因子
    tmp=sort([leading;lagging],'ascend');
    noiseP=tmp(k);                           %计算噪声功率
    threshold=alpha*noiseP;                  %计算门限
    flag=(cut>=threshold);                   %待检单元与检测门限对比获得检测结果
end

上述代码已与matlab自带的CFAR代码进行了对比，证明了上述代码的正确性。

(11)理论性能对比代码

snr=5:30;       %信噪比(dB)
N=[16,32];      %参考窗长度
k=round(N*3/4); 
PFA=1e-6;      %虚警概率

PD_Optimal=PFA.^(1./(1+10.^(snr/10))); 
PD_CA=zeros(length(N),length(snr));
PD_SO=zeros(length(N),length(snr));
PD_GO=zeros(length(N),length(snr));
PD_OS=zeros(length(N),length(snr));
for ii=1:length(N)
    alpha_CA=func_CaAlpha(PFA,N(ii));
    alpha_OS=func_OsAlpha(PFA,N(ii),k(ii));
    alpha_GO=func_GoAlpha(PFA,N(ii));
    alpha_SO=func_SoAlpha(PFA,N(ii));
     
    PD_CA(ii,:)=func_CaPd(alpha_CA,N(ii),snr);
    PD_OS(ii,:)=func_OsPd(alpha_OS,N(ii),k(ii),snr);
    PD_SO(ii,:)=func_SoPd(alpha_SO,N(ii),snr);
    PD_GO(ii,:)=func_GoPd(alpha_GO,N(ii),snr);
end

figure;hold on;box on;
plot(PD_Optimal,'r');
% plot(PD_CA(1,:),'b.-');
plot(PD_CA(2,:),'b*-');
% plot(PD_GO(1,:),'k.-');
plot(PD_GO(2,:),'k*-');
% plot(PD_SO(1,:),'g.-');
plot(PD_SO(2,:),'g*-');
% plot(PD_OS(1,:),'y.-');
plot(PD_OS(2,:),'y*-');
legend('optimal','CA32','GO32','SO32','OS32');
axis tight;

(12)单目标情况下各CFAR算法检测门限

PFA=1e-5;          %虚警概率
N=32;              %参考窗长度
k=round(N*3/4);   
L=100;             %原始数据长度(这边暂时不考虑边缘效应，所以实际生成的数据要比这个长度稍长)

SNR=15;            %目标信噪比(dB)
noiseP=20;         %噪声背景平均功率(dB)
trgtLoc=50;        %目标所在距离单元

% 计算每种算法对应的门限乘积因子
alpha_CA=func_CaAlpha(PFA,N);
alpha_OS=func_OsAlpha(PFA,N,k);
alpha_GO=func_GoAlpha(PFA,N);
alpha_SO=func_SoAlpha(PFA,N);

% 每种算法对应的门限值(dB)
threshold_CA=zeros(L+N,1);
threshold_OS=zeros(L+N,1);
threshold_GO=zeros(L+N,1);
threshold_SO=zeros(L+N,1);

% 生成原始平方率检波数据
sigma0=sqrt(10^(noiseP/10))/sqrt(2);    %单路噪声标准差
sigma1=sqrt(10^((noiseP+SNR)/10))/sqrt(2);%单路目标标准差
x=sigma0*(randn(L+N,1)+1j*randn(L+N,1));%生成指定功率的噪声信号
cut=sigma1*(randn(1,1)+1j*randn(1,1));%待检单元回波
x(trgtLoc)=cut;                       %待检单元回波
x=abs(x).^2;                           %模值平方

for ii=N/2+1:L+N/2
    left=ii-N/2:ii-1;   %前参考窗索引
    right=ii+1:ii+N/2;  %后参考窗索引
    
    threshold_CA(ii)=10*log10(alpha_CA*mean(x([left right])));
    threshold_GO(ii)=10*log10(alpha_GO*max(mean(x(left)),mean(x(right))));
    threshold_SO(ii)=10*log10(alpha_SO*min(mean(x(left)),mean(x(right))));
    
    tmp=sort(x([left right]),'ascend');
    threshold_OS(ii)=10*log10(alpha_OS*tmp(k));
end

figure;hold on;box on;
plot(10*log10(x),'b');
plot(threshold_CA,'r');
plot(threshold_GO,'k');
plot(threshold_SO,'y');
plot(threshold_OS,'g');
xlabel('采样点');ylabel('幅度(dB)');axis tight;xlim([N/2+1 N/2+1+L]);
legend('原始采样点','CA-CFAR','GO-CFAR','SO-CFAR','OS-CFAR');

(13)蒙特卡洛仿真实验代码

PFA=1e-5;          %虚警概率
mcTimes=round(1/PFA);%蒙特卡洛仿真次数
N=32;              %参考窗长度
k=round(N*3/4);   

noiseP=10;         %背景噪声功率(dB)
SNR=5:30;          %目标信噪比(dB)

% 计算各种算法的门限乘积因子
alpha_CA=func_CaAlpha(PFA,N);
alpha_GO=func_GoAlpha(PFA,N);
alpha_SO=func_SoAlpha(PFA,N);
alpha_OS=func_OsAlpha(PFA,N,k);

% 计算各种算法理论的检测概率
CaPd_Theo=zeros(length(SNR),1);
OsPd_Theo=zeros(length(SNR),1);
SoPd_Theo=zeros(length(SNR),1);
GoPd_Theo=zeros(length(SNR),1);

for ii=1:length(SNR)
    CaPd_Theo(ii)=func_CaPd(alpha_CA,N,SNR(ii));
    OsPd_Theo(ii)=func_OsPd(alpha_OS,N,k,SNR(ii));
    SoPd_Theo(ii)=func_SoPd(alpha_SO,N,SNR(ii));
    GoPd_Theo(ii)=func_GoPd(alpha_GO,N,SNR(ii));
end

% 计算各种算法蒙特卡洛仿真的结果
rng(1000);
CaPd_Mc=zeros(length(SNR),1);
OsPd_Mc=zeros(length(SNR),1);
SoPd_Mc=zeros(length(SNR),1);
GoPd_Mc=zeros(length(SNR),1);

sigma=sqrt(10^(noiseP/10))/sqrt(2);    %单路噪声标准差
for isnr=1:length(SNR)
    trgt=10^((noiseP+SNR(isnr))/20);                   %目标幅值
    x=sigma*(randn(N+1,mcTimes)+1j*randn(N+1,mcTimes));%生成指定功率的噪声信号
    x(N/2+1,:)=trgt;
    x=abs(x).^2;

    CaDetector=phased.CFARDetector('Method','CA','NumGuardCells',0,'NumTrainingCells',N,'ProbabilityFalseAlarm',PFA);
    CaResult=CaDetector(x,N/2+1);
    CaPd_Mc(isnr)=sum(CaResult)/mcTimes;
    
    GoDetector=phased.CFARDetector('Method','GOCA','NumGuardCells',0,'NumTrainingCells',N,'ProbabilityFalseAlarm',PFA);
    GoResult=GoDetector(x,N/2+1);
    GoPd_Mc(isnr)=sum(GoResult)/mcTimes;
    
    SoDetector=phased.CFARDetector('Method','SOCA','NumGuardCells',0,'NumTrainingCells',N,'ProbabilityFalseAlarm',PFA);
    SoResult=SoDetector(x,N/2+1);
    SoPd_Mc(isnr)=sum(SoResult)/mcTimes;
    
        OsDetector=phased.CFARDetector('Method','OS','Rank',k,'NumGuardCells',0,'NumTrainingCells',N,'ProbabilityFalseAlarm',PFA);
    OsResult=OsDetector(x,N/2+1);
    OsPd_Mc(isnr)=sum(OsResult)/mcTimes;
    
    disp(strcat('第',num2str(isnr),'个信噪比计算完成!'));
end

figure;hold on;box on;
plot(SNR,CaPd_Theo,'r');
plot(SNR,CaPd_Mc,'r.-');
plot(SNR,OsPd_Theo,'b');
plot(SNR,OsPd_Mc,'b.-');
plot(SNR,GoPd_Theo,'k');
plot(SNR,GoPd_Mc,'k.-');
plot(SNR,SoPd_Theo,'g');
plot(SNR,SoPd_Mc,'g.-');
xlabel('信噪比(dB)');ylabel('P_d');axis tight;
legend('CA','CA-MC','OS','OS-MC','GO','GO-MC','SO','SO-MC');

由于在上面的附录里面已经验证了自己写的代码和matlab自带的CFAR代码结果是相吻合的，所以以后在处理的时候可以直接用matlab的库函数来实现，显得更加专业和灵活。

(14)修改目标回波生成方式后代码

PFA=1e-5;          %虚警概率
mcTimes=round(1/PFA);%蒙特卡洛仿真次数
N=32;              %参考窗长度
k=round(N*3/4);   

noiseP=10;         %背景噪声功率(dB)
SNR=5:30;          %目标信噪比(dB)

% 计算各种算法的门限乘积因子
alpha_CA=func_CaAlpha(PFA,N);
alpha_GO=func_GoAlpha(PFA,N);
alpha_SO=func_SoAlpha(PFA,N);
alpha_OS=func_OsAlpha(PFA,N,k);

% 计算各种算法理论的检测概率
CaPd_Theo=zeros(length(SNR),1);
OsPd_Theo=zeros(length(SNR),1);
SoPd_Theo=zeros(length(SNR),1);
GoPd_Theo=zeros(length(SNR),1);

for ii=1:length(SNR)
    CaPd_Theo(ii)=func_CaPd(alpha_CA,N,SNR(ii));
    OsPd_Theo(ii)=func_OsPd(alpha_OS,N,k,SNR(ii));
    SoPd_Theo(ii)=func_SoPd(alpha_SO,N,SNR(ii));
    GoPd_Theo(ii)=func_GoPd(alpha_GO,N,SNR(ii));
end

% 计算各种算法蒙特卡洛仿真的结果
rng(1000);
CaPd_Mc=zeros(length(SNR),1);
OsPd_Mc=zeros(length(SNR),1);
SoPd_Mc=zeros(length(SNR),1);
GoPd_Mc=zeros(length(SNR),1);

for isnr=1:length(SNR)
    % 回波生成方式1:高斯分布的模值平方
    %sigma0=sqrt(db2pow(noiseP)/2);                     %单路噪声标准差
    %sigma1=sqrt(db2pow(noiseP+SNR(isnr))/2);           %单路目标回波的标准差
    %x=sigma0*(randn(N+1,mcTimes)+1j*randn(N+1,mcTimes));%生成指定功率的噪声信号
    %trgt=sigma1*(randn(1,mcTimes)+1j*randn(1,mcTimes)); %生成指定信噪比的目标回波信号
    %x(N/2+1,:)=trgt;
    %x=abs(x).^2;
    
    % 回波生成方式2:直接用指数分布生成
    lambda0=db2pow(noiseP);                           %噪声指数分布参数
    lambda1=db2pow(noiseP+SNR(isnr));                 %目标回波指数分布
    x=exprnd(lambda0,N+1,mcTimes);                    %生成指定功率的噪声信号
    trgt=exprnd(lambda1,1,mcTimes);                   %生成指定信噪比的目标回波信号
    x(N/2+1,:)=trgt;

    CaDetector=phased.CFARDetector('Method','CA','NumGuardCells',0,'NumTrainingCells',N,'ProbabilityFalseAlarm',PFA);
    CaResult=CaDetector(x,N/2+1);
    CaPd_Mc(isnr)=sum(CaResult)/mcTimes;
    
    GoDetector=phased.CFARDetector('Method','GOCA','NumGuardCells',0,'NumTrainingCells',N,'ProbabilityFalseAlarm',PFA);
    GoResult=GoDetector(x,N/2+1);
    GoPd_Mc(isnr)=sum(GoResult)/mcTimes;
    
    SoDetector=phased.CFARDetector('Method','SOCA','NumGuardCells',0,'NumTrainingCells',N,'ProbabilityFalseAlarm',PFA);
    SoResult=SoDetector(x,N/2+1);
    SoPd_Mc(isnr)=sum(SoResult)/mcTimes;
    
        OsDetector=phased.CFARDetector('Method','OS','Rank',k,'NumGuardCells',0,'NumTrainingCells',N,'ProbabilityFalseAlarm',PFA);
    OsResult=OsDetector(x,N/2+1);
    OsPd_Mc(isnr)=sum(OsResult)/mcTimes;
    
    disp(strcat('第',num2str(isnr),'个信噪比计算完成!'));
end

figure;hold on;box on;
plot(SNR,CaPd_Theo,'r');
plot(SNR,CaPd_Mc,'r.-');
plot(SNR,OsPd_Theo,'b');
plot(SNR,OsPd_Mc,'b.-');
plot(SNR,GoPd_Theo,'k');
plot(SNR,GoPd_Mc,'k.-');
plot(SNR,SoPd_Theo,'g');
plot(SNR,SoPd_Mc,'g.-');
xlabel('信噪比(dB)');ylabel('P_d');axis tight;
legend('CA','CA-MC','OS','OS-MC','GO','GO-MC','SO','SO-MC');