常用恒虚警检测基本原理推导

常用恒虚警检测算法的推导过程

检测问题其实就是一个二元假设检验的问题,即判断待检测数据中是仅含噪声(干扰)还是同时含有目标。通常信号是在统计意义上描述的,因此涉及到上述两种假设检验条件下的概率密度函数\(f(y|H_0)\)\(f(y|H_1)\),其中\(y\)表示样本值。通常情况下,若是基于\(N\)个样本进行检测,则需要考虑这些样本组成的样本向量\(\boldsymbol{y}=[y_1,...,y_N]^{\rm T}\)的联合概率密度函数\(f(\boldsymbol{y}|H_0)\)\(f(\boldsymbol{y}|H_1)\)。目前常用的检测方法均是以似然比检测(Likehood Ratio Test, LRT)为基础进行后续的推导的,似然比检测的决策准则如下:

\[\frac{f(\boldsymbol{y}|H_1)}{f(\boldsymbol{y}|H_0)} \gtrless T \tag{1} \]

若根据(1)对一个样本进行检测,需要分别计算\(H_0\)\(H_1\)条件下的联合概率密度函数,然后将二者的比值与对应的门限值进行比较,然后在实际应用中往往难以得到样本值准确的概率密度函数表达式,所以严格按照(1)式进行处理可能会很难办。不同下面通过几个简单的例子我们会发现,利用(1)推导出来的等效决策准则并不需要知道准确的概率密度函数。

例子1:在零均值,方差为\(\sigma^2\)的高斯白噪声中,检测是否存在一个常数\(m\),则\(H_0\)\(H_1\)对应的样本应该满足:

\[\begin{equation} \begin{aligned} &H_0:\boldsymbol{y} \backsim N(0,\sigma^2 \boldsymbol{I}_N)\\ &H_1:\boldsymbol{y} \backsim N(m \boldsymbol{I}_N,\sigma^2 \boldsymbol{I}_N) \end{aligned} \end{equation}\tag{2} \]

所以

\[\begin{equation} \begin{aligned} f(\boldsymbol{y}|H_0)&=\prod_{n=0}^{N-1}{\frac{1}{\sqrt{2\pi \sigma^2}} {\rm exp}[-\frac{1}{2} (\frac{y_n}{\sigma})^2]}\\ f(\boldsymbol{y}|H_1)&=\prod_{n=0}^{N-1}{\frac{1}{\sqrt{2\pi \sigma^2}} {\rm exp}[-\frac{1}{2} (\frac{y_n-m}{\sigma})^2]} \end{aligned} \end{equation}\tag{3} \]

所以

\[\Lambda(\boldsymbol{y})=\frac{f(\boldsymbol{y}|H_1)}{f(\boldsymbol{y}|H_0)}=\frac{\prod_{n=0}^{N-1}{{\rm exp}[-\frac{1}{2} (\frac{y_n-m}{\sigma})^2]}}{\prod_{n=0}^{N-1}{{\rm exp}[-\frac{1}{2} (\frac{y_n}{\sigma})^2]}}\tag{4} \]

对上式左右两端同时取对数可以得到

\[{\rm ln}\Lambda(\boldsymbol{y})=\sum_{n=0}^{N-1}{[-\frac{1}{2} (\frac{y_n-m}{\sigma})^2+\frac{1}{2}(\frac{y_n}{\sigma})^2]}=\frac{1}{\sigma^2}\sum_{n=0}^{N-1}{m y_n}-\frac{1}{2\sigma^2}\sum_{n=0}^{N-1}{m^2} \tag{5} \]

将(5)代入(1)中化简得到

\[\sum_{n=0}^{N-1}{y_n} \gtrless \frac{\sigma^2}{m} {\rm ln}(T)+Nm/2=T' \tag{6} \]

从(6)可以看出,要想得到(1)的似然比检测结果,并不需要知道完整的概率密度函数,而仅仅需要利用样本值的函数(本例中\(g(\boldsymbol{y})=\sum_{n=0}^{N-1}{y_n}\)),我们称之为样本值的充分统计量与门限值进行比较即可,门限值的设置可能与样本数、虚警概率、检测概率等因素有关。

例子2:高斯白噪声背景下非起伏目标的平方律检波过程推导

设目标样本\(m=\hat{m}{\rm exp}(i\theta)\),则两种假设检验条件下的采样样本分别可以表示为

\[H_0:y_n=w_n \\ H_1:y_n=m+w_n \tag{7} \]

\(z_n=|y_n|\),则在\(H_0\)\(H_1\)条件下,\(z_n\)分别服从瑞利分布和莱斯分布,两者对应的概率密度函数可以表示为:

\[\begin{equation} \begin{aligned} f_{z_n}(z_n|H_0)&= \begin{cases} \frac{2 z_n}{\sigma^2}e^{-z_n^2/\sigma^2},&\text{$z_n \geq 0$}\\ 0,&\text{$z_n < 0$} \end{cases}\\ f_{z_n}(z_n|H_1)&= \begin{cases} \frac{2 z_n}{\sigma^2}e^{-(z_n^2+\hat{m}^2)/\sigma^2}I_0(\frac{2\hat{m}z_n} {\sigma^2}),&\text{$z_n \geq 0$}\\ 0,&\text{$z_n < 0$} \end{cases} \end{aligned} \end{equation}\tag{8} \]

需要注意的是在(8)中I,Q通道中的噪声功率相等,均为\(\sigma^2/2\),而原始采样点\(y_n\)是单通道采样点(I或Q通道)。其中\(I_0(z)\)为第一类修正贝塞尔函数。此时\(N\)个采样的联合概率密度函数可以表示为

\[\begin{equation} \begin{aligned} f_{\boldsymbol{z}}(\boldsymbol{z}|H_0)&=\prod_{n=0}^{N-1}{\frac{2 z_n}{\sigma^2}e^{-z_n^2/\sigma^2}} \\ f_{\boldsymbol{z}}(\boldsymbol{z}|H_1)&=\prod_{n=0}^{N-1}{\frac{2 z_n}{\sigma^2}e^{-(z_n^2+\hat{m}^2)/\sigma^2}I_0(\frac{2\hat{m}z_n}{\sigma^2})} \end{aligned} \end{equation}\tag{9} \]

所以其似然比和对数似然比可以分别表示为

\[\begin{equation} \begin{aligned} \Lambda(\boldsymbol{z})&=\frac{f_{\boldsymbol{z}}(\boldsymbol{z}|H_1)}{f_{\boldsymbol{z}}(\boldsymbol{z}|H_0)}=e^{-\hat{m}/\sigma^2}\prod_{n=0}^{N-1}{I_0(\frac{2\hat{m}z_n}{\sigma^2})} \gtrless T \\ {\rm ln}\Lambda(\boldsymbol{z})&=-\frac{\hat{m}^2}{\sigma^2}+\sum_{n=0}^{N-1}{{\rm ln}[I_0(\frac{2\hat{m}z_n}{\sigma^2})]} \gtrless {\rm ln}T \end{aligned} \end{equation}\tag{10} \]

所以

\[\sum_{n=0}^{N-1}{{\rm ln}[I_0(\frac{2\hat{m}z_n}{\sigma^2})]} \gtrless {\rm ln}T+\frac{\hat{m}^2}{\sigma^2}=T' \tag{11} \]

在上面的表达式中,\(g(\boldsymbol{z})=\sum_{n=0}^{N-1}{{\rm ln}[I_0(\frac{2\hat{m}z_n}{\sigma^2})]}\)为样本值的充分统计量,可以看出此时只需要知道目标信噪比和原始门限值即可实现似然比检测。

对例子2进行进一步分析,根据贝塞尔函数和对数函数展开性质:

\[I_0(z)=1+x^2/4+x^4/64+... \\ {\rm ln}(1+z)=z-z^2+z^3/3+... \tag{12} \]

对(11)进行化简可以得到

\[\sum_{n=0}^{N-1}{\frac{\hat{m}^2 z_n^2}{\sigma^4}} \gtrless T' \tag{13} \]

所以

\[\sum_{n=0}^{N-1}{z_n^2} \gtrless \frac{\sigma^4}{\hat{m}^2}T' \tag{13} \]

所以当给定一组采样值之后,我们可以通过(13)所示的式子来实现似然比检测。

上面通过给定的两个例子说明了在实际应用中往往并不需要通过确切的样本概率密度函数,而是通过样本的充分统计以及其他的样本特征参数即可实现似然比检测。不过似然比检测的另外一个重要因素就是门限值的设置,而在上面的描述中还未涉及到相关内容。通常情况下,在雷达系统中,往往需要在给定虚警概率(或检测概率)的条件下来设置对应的门限值,以下给出虚警概率、检测概率与门限值之间的关系

\[P_{FA}=\int_T^\infty{f(g(\boldsymbol{y})|H_0)d\boldsymbol{y}} \\ P_D=\int_T^\infty{f(g(\boldsymbol{y})|H_1)d\boldsymbol{y}} \tag{14} \]

以下推导平方率检波器的检测门限和虚警率之间的关系式。平方率检波器的实现过程大致如(13)所示,设有\(N\)个采样点,则将这些采样点幅度的平方进行求和,并与特定门限进行比较。在(13)中令\(r_n=z_n^2\),则有所有采样点对应的充分统计量可以表示为\(g(\boldsymbol{r})=\sum_{n=0}^{N-1}{r_n}\),则根据(14)可知,需要首先求出充分统计量在两种假设检验条件下的概率密度函数。由于噪声服从高斯分布,因此\(r_n\)\(H_0\)\(H_1\)条件下分别服从指数分布和莱斯分布。两者的概率密度函数可以分别表示为

\[\begin{equation} \begin{aligned} H_0:f_{r_n}(r_n|H_0)&= \begin{cases} \frac{1}{\sigma^2}e^{-r_n/\sigma^2},&\text{$r_n \geq 0$}\\ 0,&\text{$r_n < 0$} \end{cases}\\ H_1:f_{r_n}(r_n|H_1)&= \begin{cases} \frac{1}{\sigma^2}e^{-\frac{r_n+m^2}{\sigma^2}}I_0(2\frac{m \sqrt{r_n}} {\sigma^2}),&\text{$r_n \geq 0$}\\ 0,&\text{$r_n < 0$} \end{cases} \end{aligned} \end{equation}\tag{15} \]

为了后续推导的方便,令\(z_n'=z_n/\sigma\),即将样本点对噪声功率进行归一化,则此时重新令上面的\(x_n=z_n'^2\),则\(x_n=r_n/\sigma^2\),则归一化后\(x_n\)的概率密度函数可以表示为(具体推导过程参考附录1)

\[\begin{equation} \begin{aligned} H_0:f_{x_n}(x_n|H_0)&= \begin{cases} e^{-x_n},&\text{$x_n \geq 0$}\\ 0,&\text{$x_n < 0$} \end{cases}\\ H_1:f_{x_n}(x_n|H_1)&= \begin{cases} e^{-(x_n+\chi)}I_0(2\sqrt{y\chi}),&\text{$x_n \geq 0$}\\ 0,&\text{$x_n < 0$} \end{cases} \end{aligned} \end{equation}\tag{16} \]

\(N\)个采样点对应的功率和的概率密度函数为单个采样点功率的概率密度函数的\(N\)维折叠卷积(两个随便变量的和的概率密度函数为各自概率密度函数的卷积),直接从\(N\)维折叠卷积的定义出发计算这个概率密度函数并不容易,因此下面采用特征函数的方法进行求解。特征函数的定义及基本性质参考附录2。

\(x_n\)\(H_0\)情况下的特征函数可以表示为

\[\phi(t)=\int_0^{\infty}{e^{-x_n}e^{i t x_n}dx_n}=\frac{1}{1-it} \tag{17} \]

\(x=\sum{x_n}\)的特征函数为\(\phi(t)^N=(\frac{1}{1-it})^N\),所以其概率密度函数可以表示为

\[f_x(x|H_0)=\frac{1}{2\pi}\int_{-\infty}^{\infty}{(\frac{1}{1-it})^N e^{-i t x}dt}\tag{18}= \begin{cases} \frac{x^{N-1}}{(N-1)!}e^{-x},&\text{$x \geq 0$}\\ 0,&\text{$x < 0$} \end{cases} \]

显然,当\(N=1\)时,上式退化为指数分布。虚警概率只与目标不存在时的概率密度函数有关,因此可以计算得到虚警概率为

\[P_{FA}=\int_T^{\infty}{f_x(x|H_0)dx}=1-I(T/\sqrt{N},N-1)\tag{19} \]

其中\(I(\mu,M)=\int_0^{\mu \sqrt{M+1}}{\frac{e^{-\tau}\tau^M}{M!}d\tau}\)是不完全伽马函数的皮尔逊形式。对于\(N=1\)的单采样情况,上式退化为\(P_{FA}=e^{-T}\),对未归一化的情况应该有\(P_{FA}=e^{-T/\sigma^2}\),所以\(T=-\sigma^2 {\rm ln}P_{FA}\)

其实对于平方率检波来说,在单样本情况下,其虚警概率、检测概率与门限之间的关系可以通过下面的方法进行简单计算:

虚警概率仅与目标不存在(\(H_0\))情况下的概率密度函数有关,而此时检测单元的概率密度函数可以表示为

\[f_y(y|H_0)=\frac{1}{\sigma^2}e^{-y/\sigma^2} \tag{20} \]

所以虚警概率可计算如下

\[P_{FA}=\int_T^{\infty}{f_y(y|H_0)dy}=e^{-T/\sigma^2}\tag{21} \]

检测概率与目标存在(\(H_1\))条件下的概率密度函数有关,此时若目标的信噪比为\(\chi\),则检测单元的功率为\(\sigma^2(1+\chi)\),所以此时检测单元的概率密度函数可以表示为

\[f_y(y|H_1)=\frac{1}{\sigma^2(1+\chi)}e^{-y/[\sigma^2(1+\chi)]}\tag{22} \]

所以检测概率可计算如下

\[P_D=\int_T^{\infty}{f_y(y|H_1)dy}=e^{-T/[\sigma^2(1+\chi)]}\tag{23} \]

需要注意的是,上面推导的是平稳噪声环境的情况,只需要根据噪声情况设置一个固定的门限即可,在雷达的整个工作过程中并不用对检测门限值进行调整,即没有用到恒虚警的思想。上面的推导结果可以作为最优检测器用于与下述的各种CFAR检测器进行性能对比。

(1) CA-CFAR推导

从上面的推导过程中我们知道,对于平方率检波来说,单样本情况下检测门限与虚警概率满足\(T=-\sigma^2 {\rm ln}P_{FA}\)。需要注意的是,这个关系是在噪声(或杂波)功率精确已知的情况下得到的,但是在实际应用中\(\sigma^2\)并不知道,需要从样本中进行估计,显然估计值\(\hat{\sigma}^2\)将服从一定的分布,而\(P_{FA}\)可以表示为\(\hat{\sigma}^2\)的函数,因此它也满足一定的分布,所以实际应用条件下可以将\(T=-\sigma^2 {\rm ln}P_{FA}\)改写为\(T=\alpha \hat{\sigma}^2\)。其中\(\alpha\)是与虚警概率有关的参数。现在的任务有两个,一个是利用实测数据估计噪声功率\(\hat{\sigma}^2\),另一个是求\(\alpha\)

噪声功率通常是利用待检单元两端的\(N\)个单元进行估计的,因此一般需要满足以下两个条件:

1)临近单元所含噪声(杂波)的统计特性与待检单元一致(称为均匀干扰);

2)临近单元不包含任何目标,其仅仅存在干扰噪声。

设第\(i\)个待检单元为\(x_i\)(功率值),它服从指数分布,即

\[f_{x_i}(x_i)=\frac{1}{\sigma^2}e^{-x_i/\sigma^2} \tag{24} \]

\(N\)个样本的联合概率密度可以表示为

\[f_{\boldsymbol{x}}(\boldsymbol{x})=\prod_{i=0}^{N-1}{f_{x_i}(x_i)}=\frac{1}{\sigma^{2N}}\prod_{i=0}^{N-1}{e^{-x_i/\sigma^2}}=\frac{1}{\sigma^{2N}}{\rm exp}[-(\sum_{i=0}^{N-1}{x_i})/\sigma^2] \tag{25} \]

所以

\[{\rm ln}f_{\boldsymbol{x}}(\boldsymbol{x})=-N{\rm ln}\sigma^2-\frac{1}{\sigma^2}(\sum_{i=0}^{N-1}{x_i}) \tag{26} \]

将(22)对\(\sigma^2\)求导,并让导数等于零,可以求得\(\sigma^2\)的最大似然估计解为

\[\hat{\sigma}^2=\frac{1}{N}\sum_{i=0}^{N-1}{x_i} \tag{27} \]

所以\(\hat{T}=\alpha\hat{\sigma}^2=\frac{\alpha}{N}\sum_{i=0}^{N-1}{x_i}\),设\(z_i=\frac{\alpha}{N}x_i\),则\(z_i\)的概率密度函数可以参考附录1进行求解,而\(\hat{T}=\sum{z_i}\)率密度函数则可以利用特征函数进行求解,并且得到

\[f_{\hat{T}}(\hat{T})=(\frac{N}{\alpha \sigma^2})^N \frac{\hat{T}^{N-1}}{(N-1)!}e^{-N\hat{T}/\alpha \sigma^2} \tag{28} \]

所以\(P_{FA}=e^{-\hat{T}/\hat{\sigma}^2}\),所以可以求得\(P_{FA}\)的期望为

\[{\rm E}(P_{FA})=\int_{-\infty}^{\infty}{e^{-\hat{T}/\hat{\sigma}^2}f_{\hat{T}}(\hat{T})d\hat{T}} \tag{29} \]

可以求解得到

\[{\rm E}(P_{FA})=(1+\frac{\alpha}{N})^{-N} \tag{30} \]

所以

\[\alpha=N({\rm E}(P_{FA})^{-1/N}-1) \tag{31} \]

(2)OS-CFAR推导

OS-CFAR和CA-CFAR的差别在于估算噪声功率的方法不一样。CA-CFAR将待检单元周围单元的平均功率作为噪声功率估计值,而OS-CFAR则是将待检单元周围单元进行排序后选择其中某个单元的功率作为噪声功率估计值。需要注意的是,虽然OS-CFAR最终只用到一个单元的功率值,但是在排序过程中综合应用了所有单元的信息。推导OS-CFAR相关性能可遵循类似CA-CFAR的步骤进行,即首先需要求得待检单元的概率密度函数,然后由\(T=\alpha \hat{\sigma}^2\)求得检测门限\(T\)的概率密度函数,最后由\(P_{FA}=e^{-T/\sigma^2}\)可求得虚警概率期望值,进而求得门限乘积因子\(\alpha\)的表达式。以下具体进行求解。

设排序之后选择第\(k\)个样本的值\(x_k\)作为噪声功率估计值,则\(x_k\)的概率密度函数可以表示为

\[f_{x_k}(x)=k C_N^k f_x(x)F_x^{k-1}(x)[1-F_x(x)]^{N-k} \tag{32} \]

(目前还不清楚上面的表达式是如何推导出来的)在上面的表达式中\(f_x(x)=\frac{1}{\sigma^2}e^{-x/\sigma^2}\)表示任意一个单元的概率密度函数,\(F_x(x)=1-e^{-x/\sigma^2}\)表示任意一个单元的分布函数,将两者代入上式可以得到

\[f_{x_k}(x)=k C_N^k [e^{-x/\sigma^2}]^{N-k+1} [1-e^{-x/\sigma^2}]^{k-1}\tag{33} \]

\(T=\alpha \hat{\sigma}^2\)可以得到

\[f_{\hat{T}}(\hat{T})=\frac{1}{\alpha}f_{x_k}(\hat{T}/\alpha)=\frac{k}{\alpha}C_N^k [e^{-\hat{T}/(\alpha\sigma^2)}]^{N-k+1} [1-e^{-\hat{T}/(\alpha\sigma^2)}]^{k-1}\tag{34} \]

所以

\[{\rm E}[P_{F_A}]=\int_0^\infty{e^{-\hat{T}/\sigma^2}f_{\hat{T}}(\hat{T})d\hat{T}}=\frac{k}{\alpha}C_N^k \int_0^\infty{[e^{-\hat{T}/(\alpha\sigma^2)}]^{N-k+1+\alpha}[1-e^{-\hat{T}/(\alpha\sigma^2)}]^{k-1}d\hat{T}} \tag{35} \]

在上式中令\(\lambda=1/\sigma^2,x=\frac{\lambda}{\alpha}\hat{T}\),则上面的积分比较难计算,则上式可以重新表示为

\[\begin{equation} \begin{aligned} {\rm E}[P_{F_A}]&=\frac{k}{\lambda}C_N^k \int_0^\infty{e^{-(N-k+1+\alpha)x}[1-e^{-x}]^{k-1}dx}=\frac{k}{\lambda}C_N^k B(N-k+1+\alpha,k)\\ &=\frac{k}{\lambda}C_N^k \frac{\Gamma(n-k+1+\alpha)\Gamma(k)}{\Gamma(N+1+\alpha)} \end{aligned} \end{equation}\tag{36} \]

上面的积分过程可以参考文献[2]。在上面公式中\(B\)表示\(\beta\)函数,\(\Gamma\)表示伽玛函数,当\(\alpha\)为整数时,(36)可进一步化简为

\[{\rm E}[P_{F_A}]=\sigma^2 \frac{N!(\alpha+N-k)!}{(N-k)!(\alpha+N)!}=\frac{\sigma^2 N!}{(N-k)!}\frac{(\alpha+N-k)!}{(\alpha+N)!}=\frac{\sigma^2 N!}{(N-k)!}\frac{1}{\prod_{i=0}^{k-1}(\alpha+N-i)}\tag{37} \]

在实际应用中,通常是设定\(N,k,P_{FA}\),然后首先估算得到\(\sigma^2\),最后根据(37)反推\(\alpha\)。从(37)可知,\(\alpha\)可以转化为求\(k\)次多项式方程的根。这边推导感觉有点问题,后面进一步学习后,再将其纠正!!!

(3)GO-CFAR推导

这种CFAR算法的具体做法是将参考窗平均划分为前参考窗和后参考窗,然后选择两个参考窗中噪声平均值较大的那个窗的结果作为噪声的估计值。具体地,设整个参考窗的长度为\(N\),则前、后参考窗的长度均为\(N/2\)。设第\(i\)个参考单元为\(x_i\),且\(x_i\)的概率密度函数为\(f_{x_i}=\frac{1}{\sigma^2}e^{-x_i/\sigma^2}\)。设前参考窗的充分统计量为\(Z_1=\frac{2}{N}\sum_{i=1}^{N/2}X_i\),后参考窗的充分统计量为\(Z_2=\frac{2}{N}\sum_{i=N/2+1}^N{X_i}\),则GO-CFAR的充分统计量为\(Z={\rm max}(Z_1,Z_2)\)。不难得到\(Z_1,Z_2,Z\)的概率密度函数可以分别表示为(用特征函数求解)

\[\begin{equation} \begin{aligned} f_{Z_1}(z_1)&=(\frac{N}{4 \sigma^2})^{N/2} \frac{z_1^{N/2-1}}{(N/2-1)!}e^{-N z_1/4 \sigma^2}\\ f_{Z_2}(z_2)&=(\frac{N}{4 \sigma^2})^{N/2} \frac{z_2^{N/2-1}}{(N/2-1)!}e^{-N z_2/4 \sigma^2}\\ f_Z(z)&=f_{Z_1}(z_1)F_{Z_2}(z_2)+F_{Z_1}(z_1)f_{Z_2}(z_2) \end{aligned} \end{equation} \]

(4)SO-CFAR推导

与GO-CFAR类似,这种算法也需要将参考窗划分为前后参考窗,不过SO-CFAR选择的是两个参考窗中估计噪声较小值作为噪声估计值。即\(Z={\rm min}(Z_1,Z_2)\)。此时\(Z\)的概率密度函数表达式如下

\[f_Z(z)=f_{Z_1}(z_1)+f_{Z_2}(z_2)-[f_{Z_1}(z_1)F_{Z_2}(z_2)+f_{Z_2}(z_2)F_{Z_1}(z_1)] \]

需要注意的是,对于GO-CFAR和SO-CFAR在推导过程中,充分统计量的计算相对复杂,本文不进行详细推导,相关内容将在后面的文档进行进一步说明。


参考文献

[1] Fundamentals of Radar Signal Processing. Chapter 6 and Chapter 7.

[2] Tables of Integrals, Series, and Products.

[3] 如何通俗地理解矩母函数


附录

(1)随机变量函数的概率密度函数

设已知随机变量\(X\)的pdf为\(f_X(x)\),cdf为\(F_X(x)\),若另外一个随机变量\(Y\)可以表示为\(X\)的函数,即\(y=h(x)\),则\(F_Y(y)=P(Y \leq y)=P(h(X)\leq y)=P(X \leq g(y))=F_X(g(y))\),所以\(f_Y(y)=F_Y(y)'=f_X(g(y))h'(y)\)。当然上面的推导只是为了展示基本思路,并没有考虑函数\(h(x)\)\(g(y)\)的详细情况,以下通过几个简单的例子来详细展示上述整个流程。

莱斯分布常用来表示正弦(余弦)信号加窄带高斯随机信号的包络分布。其概率密度函数可以表示为

\[f_X(x)= \begin{cases} \frac{x}{\sigma^2}e^{-\frac{x^2+m^2}{2\sigma^2}}I_0(\frac{x m}{\sigma^2}),&\text{$x \geq 0$}\\ 0,&\text{x < 0} \end{cases}\tag{1} \]

在(1)中\(x\)表示原始信号的包络,即原始采样信号的幅值,\(m\)表示主信号的峰值,\(I_0(x)\)表示修正贝塞尔函数。则若另一个随机变量\((Y=aX,a > 0)\),则\(Y\)的概率密度函数可以通过以下步骤就行求解:

\[F_Y(y)=P(Y \leq y)=p(aX \leq y)=p(X \leq y/a)=F_X(y/a) \\ \rightarrow f_Y(y)=F_Y'(y)=\frac{1}{a}f_x(y/a)\tag{2} \]

所以

\[f_Y(y)= \begin{cases} \frac{y}{a^2\sigma^2}e^{-\frac{y^2/a^2+m^2}{2\sigma^2}}I_0(\frac{y m/a} {\sigma^2}),&\text{$y \geq 0$}\\ 0,&\text{y < 0} \end{cases}\tag{3} \]

令举一例,若有一随机变量\((Z=X^2,X \geq 0)\),则\(Z\)的概率密度函数可以通过以下步骤求解:

\[F_Z(z)=P(Z \leq y)=p(X^2 \leq z)=p(X \leq \sqrt{z})=F_X(\sqrt{z}) \\ \rightarrow f_Z(z)=F_Z'(z)=\frac{1}{2\sqrt{z}}f_x(\sqrt{z})\tag{4} \]

所以

\[f_Z(z)= \begin{cases} \frac{1}{2\sigma^2}e^{-\frac{z+m^2}{2\sigma^2}}I_0(\frac{m \sqrt{z}} {\sigma^2}),&\text{$z \geq 0$}\\ 0,&\text{z < 0} \end{cases} \tag{5} \]

(2)随机变量的特征函数

随机变量的特征函数与其概率密度函数之间是两两对应的关系,设随便变量\(X\)的概率密度函数是\(f(x)\),其特征函数\(\phi(t)\)\(f(x)\)之间的关系可以表示为

\[\phi(t)={\rm E}(e^{itx})=\int_{-\infty}^{\infty}{e^{itx}f(x)dx} \\ f(x)=\frac{1}{2\pi}\int_{-\infty}^{\infty}{e^{-itx}\phi(t)dt}\tag{1} \]

以下通过试举指数函数的特征函数求解过程:

指数分的概率密度函数可以表示为\(f(x)=\lambda e^{-\lambda x},x \geq 0\),所以其特征函数为

\[\phi(t)=\lambda \int_0^{\infty}{[cos(tx)e^{-\lambda x}+i sin(tx)e^{-\lambda x}]dx}=\frac{\lambda^2}{\lambda^2+t^2}+i\frac{\lambda t}{\lambda^2+t^2}=(1-\frac{it}{\lambda})^{-1}\tag{2} \]

根据特征函数的定义,显然有如下结论:

\(X_1,X_2,...,X_n\)\(n\)个相互独立的随机变量,令\(Y=\sum_{i=1}^n{X_i}\)\(Z=\sum_{i=1}^n{a_i X_i+b_i}\),则\(Y\)\(Z\)的特征函数分别为

\[\phi_Y(t)=\prod_{i=1}^n{\phi_{X_i}(t)} \\ \phi_Z(t)=\prod_{i=1}^n{e^{i t b_i}\phi_{X_i}(a_i t)} \tag{3} \]

(3)随机变量的矩母函数

在随机变量的处理过程中我们经常需要计算它们的各阶矩。随机变量的矩具有一定的物理含义,能够反映随机变量在某些方面的特征,比如一阶矩\(E[x]\)表示随机变量的期望,二阶矩\(E[x^2]\)表示随机变量的方差,三阶矩\(E[x^3]\)表示随机变量分布的不确定性,....,各阶矩的定义式如下所示:

\[E[x]=\int{x f(x)dx}\\ E[x^2]=\int{x^2 f(x)dx}\\ \vdots\\ E[x^n]=\int{x^n f(x)dx}\tag{1} \]

在上面的表达式中\(f(x)\)为随机变量\(x\)的概率密度函数,可以看出求随机变量的矩涉及到函数积分,通常情况下积分比较难以计算。为此提出矩母函数(Moment Generating Function, MGF)的概念,用于辅助随机变量各阶矩的计算,具体地,MGF的定义如下:

\[{\rm MGF}(t)=E[e^{tx}]=\int{e^{tx}f(x)dx}\tag{2} \]

则有

\[E[x]=\frac{d}{dt}{\rm MGF}(t)|_{t=0}\\ E[x^2]=\frac{d^2}{dt^2}{\rm MGF}(t)|_{t=0}\\ \vdots\\ E[x^n]=\frac{d^n}{dt^n}{\rm MGF}(t)|_{t=0}\tag{3} \]

即要求随机变量的\(n\)阶矩,只需要对其MGF求\(n\)次导,并令\(t=0\)即可,显然MGF将随机变量的\(n\)阶矩由积分运算变换为微分运算,在很多情况下会极大地简化计算过程。

以下证明(3)式的正确性。根据泰勒展开公式有

\[e^{tx}=1+tx+\frac{(tx)^2}{2!}+...+\frac{(tx)^n}{n!},n \rightarrow \infty \tag{4} \]

所以

\[E[e^{tx}]=E[1]+tE[x]+\frac{t^2}{2!}E[x^2]+...+\frac{t^n}{n!}E[x^n],n \rightarrow \infty \tag{5} \]

分别对\(t\)\(n\)次导并令\(t=0\)可以得到:

\[\begin{equation} \begin{aligned} \frac{d}{dt}E[e^{tx}]&=E[x]+tE[x^2]+...+\frac{t^{n-1}}{(n-1)!}E[x^n]\rightarrow \frac{d}{dt}E[e^{tx}]|_{t=0}=E[x]\\ \frac{d^2}{dt^2}E[e^{tx}]&=E[x^2]+...+\frac{t^{n-2}}{(n-2)!}E[x^n]\rightarrow \frac{d}{dt}E[e^{tx}]|_{t=0}=E[x^2]\\ &\vdots\\ \frac{d^n}{dt^n}E[e^{tx}]&=E[x^n]\rightarrow \frac{d}{dt}E[e^{tx}]|_{t=0}=E[x^n] \end{aligned} \end{equation}\tag{6} \]

以下以指数分布函数的\(n\)阶矩求解过程进行示例说明:

指数分布的概率密度函数可以表示为\(f(x)=\lambda e^{-\lambda x},x>0\),所以

\[{\rm MGF}(t)=E[e^{tx}]=\lambda\int_0^{\infty}{e^{tx}e^{-\lambda x}dx}=\frac{\lambda}{\lambda-t}\tag{7} \]

显然(7)所示的简单表达式很容易进行求导,因此可以很方便地计算出\(n\)阶矩。

posted @ 2021-07-09 17:49  平和少年  阅读(18)  评论(0编辑  收藏  举报