本文主要内容摘自下面文献: [1] Fadel F. Digham (2007). “On the Energy Detection of Unknown Signals Over Fading Channels.” IEEE TRANSACTIONS ON COMMUNICATION 55(1): 4. [2] H.Urkowitz (1967). “Energy Detection of Unknown Deterministic Signals.” Proceedings of the IEEE 55(4): 9.【参见博文《合作频谱感知:未知确定信号的能量检测》】 [3] Zhi Quan, S. C., and Ali H. Sayed (2008). “Optimal Linear Cooperation for Spectrum Sensing in Cognitive Radio Networks.” IEEE JOURNAL OF SELECTED TOPICS IN SIGNAL PROCESSING, 2(1): 13. [4] 博文《卡方分布》
1、系统模型
Urkowitz (1967)一文中,分别讨论了低通和带通两种情况。这里主要集中讨论带通情况。设发送的低通信号为
S
L
P
(
t
)
=
s
c
(
t
)
+
j
s
s
(
t
)
S_{\rm LP}(t)=s_c(t)+js_s(t)
S L P ( t ) = s c ( t ) + j s s ( t ) ,信道为
h
=
α
e
j
θ
h=\alpha e^{j\theta}
h = α e j θ ,则带通接收信号为
r
(
t
)
=
{
R
{
[
h
s
L
P
(
t
)
+
n
L
P
(
t
)
]
e
j
2
π
f
c
t
}
:
H
1
R
{
n
L
P
(
t
)
e
j
2
π
f
c
t
}
:
H
0
r(t)={\Large\{ }\begin{aligned} {\mathcal R}\{[hs_{\rm LP}(t)+n_{\rm LP}(t)]e^{j2\pi f_ct}\}:H_1\\ {\mathcal R}\{n_{\rm LP}(t)e^{j2\pi f_ct}\}:H_0 \end{aligned}
r ( t ) = { R { [ h s L P ( t ) + n L P ( t ) ] e j 2 π f c t } : H 1 R { n L P ( t ) e j 2 π f c t } : H 0 此时,定义信噪比为
γ
=
α
2
E
s
N
0
\gamma =\frac{\alpha^2 E_s}{N_0}
γ = N 0 α 2 E s 。 根据Urkowitz (1967),我们可以得到判决变量为
r
=
2
N
0
∫
0
T
r
2
(
t
)
d
t
.
r=\frac{2}{N_0}\int_{0}^{T}r^2(t)dt.
r = N 0 2 ∫ 0 T r 2 ( t ) d t . 下面我们分别针对
H
0
H_0
H 0 和
H
1
H_1
H 1 进行推导。
H
0
H_0
H 0 带通噪声可以表示为
n
(
t
)
=
n
c
(
t
)
cos
(
2
π
f
c
t
)
−
n
s
(
t
)
sin
(
2
π
f
c
t
)
,
n(t)=n_c(t)\cos(2\pi f_ct)-n_s(t)\sin(2\pi f_ct),
n ( t ) = n c ( t ) cos ( 2 π f c t ) − n s ( t ) sin ( 2 π f c t ) , 因此,判决变量为
r
=
2
N
0
∫
0
T
n
2
(
t
)
d
t
=
1
N
0
∫
0
T
[
n
c
2
(
t
)
+
n
s
2
(
t
)
]
d
t
\begin{aligned} r&=\frac{2}{N_0}\int_{0}^{T}n^2(t)dt\\ &=\frac{1}{N_0}\int_{0}^{T}\left[ n_c^2(t)+n_s^2(t)\right]dt\\ \end{aligned}
r = N 0 2 ∫ 0 T n 2 ( t ) d t = N 0 1 ∫ 0 T [ n c 2 ( t ) + n s 2 ( t ) ] d t 进一步,我们把
n
c
(
t
)
n_c(t)
n c ( t ) 和
n
s
(
t
)
n_s(t)
n s ( t ) 分别用有限个抽样样值表示,注意这时的抽样速率为
W
W
W 而非
2
W
2W
2 W ,若
N
=
2
W
T
N=2WT
N = 2 W T ,则有
n
c
(
t
)
=
∑
i
=
1
N
/
2
n
c
i
s
i
n
c
(
W
t
−
i
)
n
s
(
t
)
=
∑
i
=
1
N
/
2
n
s
i
s
i
n
c
(
W
t
−
i
)
,
\begin{aligned} n_c(t)&=\sum_{i=1}^{N/2}n_{ci}{\rm sinc}\left(Wt-i \right)\\ n_s(t)&=\sum_{i=1}^{N/2}n_{si}{\rm sinc}\left(Wt-i \right), \end{aligned}
n c ( t ) n s ( t ) = i = 1 ∑ N / 2 n c i s i n c ( W t − i ) = i = 1 ∑ N / 2 n s i s i n c ( W t − i ) , 利用sinc函数的性质,我们可以得到
r
=
1
N
0
W
∑
i
=
1
N
/
2
[
n
c
i
2
+
n
s
i
2
]
=
∑
i
=
1
N
/
2
[
(
n
c
i
N
0
W
)
2
+
(
n
s
i
N
0
W
)
2
]
r=\frac{1}{N_0W}\sum_{i=1}^{N/2}\left[ n_{ci}^2+n_{si}^2\right]=\sum_{i=1}^{N/2}\left[ \left(\frac{n_{ci}}{\sqrt{N_0W}}\right)^2+\left(\frac{n_{si}}{\sqrt{N_0W}}\right)^2\right]
r = N 0 W 1 i = 1 ∑ N / 2 [ n c i 2 + n s i 2 ] = i = 1 ∑ N / 2 [ ( N 0 W
n c i ) 2 + ( N 0 W
n s i ) 2 ] 满足中心
χ
2
\chi^2
χ 2 分布,自由度为
2
T
W
2TW
2 T W 。
H
1
H_1
H 1 低通等效信号为
s
L
P
(
t
)
=
s
c
(
t
)
+
j
s
s
(
t
)
,
s_{\rm LP}(t)=s_c(t)+js_s(t),
s L P ( t ) = s c ( t ) + j s s ( t ) , 等效低通信道为
h
L
P
=
α
c
+
j
α
s
h_{\rm LP}=\alpha_c+j\alpha_s
h L P = α c + j α s ,则经过衰落信道后的低通等效接收信号(不含噪声)为
s
L
P
′
(
t
)
=
α
c
s
c
(
t
)
+
j
α
c
s
s
(
t
)
+
j
α
s
s
c
(
t
)
−
α
s
s
s
(
t
)
,
s'_{\rm LP}(t)=\alpha_cs_c(t)+j\alpha_c s_s(t)+j\alpha_s s_c(t)-\alpha_ss_s(t),
s L P ′ ( t ) = α c s c ( t ) + j α c s s ( t ) + j α s s c ( t ) − α s s s ( t ) , 则带通信号可以表示为
s
′
(
t
)
=
R
[
s
L
P
′
(
t
)
e
j
2
π
f
c
t
]
=
[
α
c
s
c
(
t
)
−
α
s
s
s
(
t
)
]
cos
(
2
π
f
c
t
)
−
[
α
c
s
s
(
t
)
+
α
s
s
c
(
t
)
]
sin
(
2
π
f
c
t
)
\begin{aligned} s'(t)&={\mathcal R}\left[ s'_{\rm LP}(t)e^{j2\pi f_ct} \right]\\ &=\left[ \alpha_cs_c(t)-\alpha_ss_s(t)\right]\cos(2\pi f_ct)-\left[\alpha_c s_s(t)+\alpha_s s_c(t)\right]\sin(2\pi f_ct)\\ \end{aligned}
s ′ ( t ) = R [ s L P ′ ( t ) e j 2 π f c t ] = [ α c s c ( t ) − α s s s ( t ) ] cos ( 2 π f c t ) − [ α c s s ( t ) + α s s c ( t ) ] sin ( 2 π f c t ) 由此,可以得到接收信号为
r
(
t
)
=
[
α
c
s
c
(
t
)
−
α
s
s
s
(
t
)
+
n
c
(
t
)
]
cos
(
2
π
f
c
t
)
−
[
α
c
s
s
(
t
)
+
α
s
s
c
(
t
)
+
n
s
(
t
)
]
sin
(
2
π
f
c
t
)
=
r
c
(
t
)
cos
(
2
π
f
c
t
)
−
r
s
(
t
)
sin
(
2
π
f
c
t
)
\begin{aligned} r(t)&=\left[ \alpha_cs_c(t)-\alpha_ss_s(t)+n_c(t)\right]\cos(2\pi f_ct)-\left[\alpha_c s_s(t)+\alpha_s s_c(t)+n_s(t)\right]\sin(2\pi f_ct)\\ &=r_c(t)\cos(2\pi f_ct)-r_s(t)\sin(2\pi f_ct) \end{aligned}
r ( t ) = [ α c s c ( t ) − α s s s ( t ) + n c ( t ) ] cos ( 2 π f c t ) − [ α c s s ( t ) + α s s c ( t ) + n s ( t ) ] sin ( 2 π f c t ) = r c ( t ) cos ( 2 π f c t ) − r s ( t ) sin ( 2 π f c t ) 同理,可以用抽样样值表示
r
c
(
t
)
=
∑
i
=
1
N
/
2
[
α
c
s
c
i
−
α
s
s
s
i
+
n
c
i
]
s
i
n
c
(
W
t
−
i
)
r
s
(
t
)
=
∑
i
=
1
N
/
2
[
α
c
s
s
i
+
α
s
s
c
i
+
n
s
i
]
s
i
n
c
(
W
t
−
i
)
,
\begin{aligned} r_c(t)&=\sum_{i=1}^{N/2}[\alpha_c s_{ci}-\alpha_s s_{si}+n_{ci}]{\rm sinc}\left(Wt-i \right)\\ r_s(t)&=\sum_{i=1}^{N/2}[\alpha_c s_{si}+\alpha_s s_{ci}+n_{si}]{\rm sinc}\left(Wt-i \right), \end{aligned}
r c ( t ) r s ( t ) = i = 1 ∑ N / 2 [ α c s c i − α s s s i + n c i ] s i n c ( W t − i ) = i = 1 ∑ N / 2 [ α c s s i + α s s c i + n s i ] s i n c ( W t − i ) , 由此,我们可以得到
(2)
r
=
2
N
0
∫
0
T
r
2
(
t
)
d
t
=
1
N
0
∫
0
T
[
r
c
2
(
t
)
+
r
s
2
(
t
)
]
d
t
,
\tag{2} \begin{aligned} r=&\frac{2}{N_0}\int_{0}^{T}r^2(t)dt\\ &=\frac{1}{N_0}\int_{0}^{T}\left[ r_c^{2}(t)+r_s^{2}(t)\right]dt\\ \end{aligned},
r = N 0 2 ∫ 0 T r 2 ( t ) d t = N 0 1 ∫ 0 T [ r c 2 ( t ) + r s 2 ( t ) ] d t , ( 2 ) 即
H
1
:
r
=
1
N
0
W
[
∑
i
=
1
N
/
2
(
α
c
s
c
i
−
α
s
s
s
i
+
n
c
i
)
2
+
∑
i
=
1
N
/
2
(
α
c
s
s
i
+
α
s
s
c
i
+
n
s
i
)
2
]
=
∑
i
=
1
N
/
2
(
α
c
s
c
i
−
α
s
s
s
i
+
n
c
i
N
0
W
)
2
+
∑
i
=
1
N
/
2
(
α
c
s
s
i
+
α
s
s
c
i
+
n
s
i
N
0
W
)
2
\begin{aligned} H_1:\qquad r&=\frac{1}{N_0W}\left[ \sum_{i=1}^{N/2}(\alpha_c s_{ci}-\alpha_s s_{si}+n_{ci})^2+\sum_{i=1}^{N/2}(\alpha_c s_{si}+\alpha_s s_{ci}+n_{si})^2\right]\\ &=\sum_{i=1}^{N/2}\left(\frac{\alpha_c s_{ci}-\alpha_s s_{si}+n_{ci}}{\sqrt{N_0W}}\right)^2+\sum_{i=1}^{N/2}\left(\frac{\alpha_c s_{si}+\alpha_s s_{ci}+n_{si}}{\sqrt{N_0W}}\right)^2 \end{aligned}
H 1 : r = N 0 W 1 ⎣ ⎡ i = 1 ∑ N / 2 ( α c s c i − α s s s i + n c i ) 2 + i = 1 ∑ N / 2 ( α c s s i + α s s c i + n s i ) 2 ⎦ ⎤ = i = 1 ∑ N / 2 ( N 0 W
α c s c i − α s s s i + n c i ) 2 + i = 1 ∑ N / 2 ( N 0 W
α c s s i + α s s c i + n s i ) 2 满足非中心
χ
2
\chi^2
χ 2 分布,方差
σ
2
=
1
\sigma^2=1
σ 2 = 1 ,非中心参数为
μ
=
2
γ
\mu=2\gamma
μ = 2 γ ,
γ
=
α
E
s
N
0
\gamma=\frac{\alpha E_s}{N_0}
γ = N 0 α E s 。
我们来推导下非中心参数:
μ
=
∑
i
=
1
N
/
2
(
α
c
2
s
c
i
2
N
0
W
+
α
s
2
s
s
i
2
N
0
W
)
+
∑
i
=
1
N
/
2
(
α
c
2
s
s
i
2
N
0
W
+
α
s
2
s
c
i
2
N
0
W
)
=
α
c
2
N
0
∑
i
=
1
N
/
2
s
c
i
2
W
+
α
s
2
N
0
∑
i
=
1
N
/
2
s
s
i
2
W
+
α
c
2
N
0
∑
i
=
1
N
/
2
s
s
i
2
W
+
α
s
2
N
0
∑
i
=
1
N
/
2
s
c
i
2
W
=
2
α
2
E
s
N
0
\begin{aligned} \mu&=\sum_{i=1}^{N/2}\left(\frac{\alpha_c^2 s^2_{ci}}{{N_0W}}+\frac{\alpha_s^2 s^2_{si}}{{N_0W}}\right)+\sum_{i=1}^{N/2}\left(\frac{\alpha_c^2 s^2_{si}}{{N_0W}}+\frac{\alpha_s^2 s^2_{ci}}{{N_0W}}\right)\\ &=\frac{\alpha_c^2}{N_0} \sum_{i=1}^{N/2}\frac{s^2_{ci}}{{W}}+\frac{\alpha_s^2}{N_0} \sum_{i=1}^{N/2}\frac{s^2_{si}}{{W}}+\frac{\alpha_c^2}{N_0} \sum_{i=1}^{N/2}\frac{s^2_{si}}{{W}}+\frac{\alpha_s^2}{N_0} \sum_{i=1}^{N/2}\frac{s^2_{ci}}{{W}}\\ &=\frac{2\alpha^2 E_s}{N_0} \end{aligned}
μ = i = 1 ∑ N / 2 ( N 0 W α c 2 s c i 2 + N 0 W α s 2 s s i 2 ) + i = 1 ∑ N / 2 ( N 0 W α c 2 s s i 2 + N 0 W α s 2 s c i 2 ) = N 0 α c 2 i = 1 ∑ N / 2 W s c i 2 + N 0 α s 2 i = 1 ∑ N / 2 W s s i 2 + N 0 α c 2 i = 1 ∑ N / 2 W s s i 2 + N 0 α s 2 i = 1 ∑ N / 2 W s c i 2 = N 0 2 α 2 E s
因此,
r
r
r 的概率密度函数可以表示为
f
R
(
r
)
=
{
r
N
2
−
1
e
−
r
2
σ
2
σ
N
2
N
2
Γ
(
N
2
)
,
H
0
1
2
σ
2
(
r
μ
)
N
−
2
4
e
−
μ
+
r
2
σ
2
I
N
2
−
1
(
μ
r
σ
2
)
,
H
1
,
f_R(r)=\left\{\begin{aligned} \frac{r^{\frac{N}{2}-1}e^{-\frac{r}{2\sigma^2}}}{\sigma^N2^{\frac{N}{2}}\Gamma(\frac{N}{2})},&\quad H_0\\ \frac{1}{2\sigma^2}(\frac{r}{\mu})^{\frac{N-2}{4}}e^{-\frac{\mu+r}{2\sigma^2}}I_{\frac{N}{2}-1}(\frac{\sqrt{\mu r}}{\sigma^2}),&\quad H_1, \end{aligned}\right.
f R ( r ) = ⎩ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎧ σ N 2 2 N Γ ( 2 N ) r 2 N − 1 e − 2 σ 2 r , 2 σ 2 1 ( μ r ) 4 N − 2 e − 2 σ 2 μ + r I 2 N − 1 ( σ 2 μ r
) , H 0 H 1 , 这里的
I
v
(
⋅
)
I_v(\cdot)
I v ( ⋅ ) 为
v
v
v 阶第一类修正贝塞尔函数。
2、AWGN信道上的检测概率与虚警概率
为了计算检测与虚警概率,有很多近似方法。下面我们主要讨论高斯分布来近似。我们知道,随着样本数
N
N
N 的增大,根据中心极限定理,可以认为
y
y
y 的分布接近高斯分布。 我们知道检测变量
r
r
r 满足如下分布:
r
∼
{
χ
N
2
,
H
0
χ
N
2
(
μ
)
,
H
1
r\sim \left\{\begin{aligned} \chi_N^2,&\quad H_0\\ \chi_N^2(\mu),&\quad H_1 \end{aligned} \right.
r ∼ { χ N 2 , χ N 2 ( μ ) , H 0 H 1 这里
μ
=
2
E
s
α
2
N
0
=
2
γ
\mu=\frac{2E_s\alpha^2}{N_0}=2\gamma
μ = N 0 2 E s α 2 = 2 γ ,显然
γ
\gamma
γ 为信噪比。
用高斯分布近似卡方分布: 若
m
=
∑
i
=
1
N
r
i
2
m=\sum_{i=1}^{N}r_i^2
m = ∑ i = 1 N r i 2 ,其中
r
i
∼
N
(
0
,
σ
2
)
r_i\sim {\mathcal N}(0,\sigma^2)
r i ∼ N ( 0 , σ 2 ) ,显然
m
σ
2
\frac{m}{\sigma^2}
σ 2 m 为自由度等于
N
N
N 的卡方分布。当
N
N
N 足够大时,我们可以把卡方分布近似成高斯分布。由于
(
r
i
σ
2
)
2
(\frac{r_i}{\sigma^2})^2
( σ 2 r i ) 2 的均值1,方差为2,因此
r
i
r_i
r i 的均值为
σ
2
\sigma^2
σ 2 ,方差为
2
σ
4
2\sigma^4
2 σ 4 ,因此有
m
∼
N
(
N
σ
2
,
2
N
σ
4
)
.
m\sim {\mathcal N}(N\sigma^2,2N\sigma^4).
m ∼ N ( N σ 2 , 2 N σ 4 ) . 若
m
=
∑
i
=
1
N
r
i
2
m=\sum_{i=1}^{N}r_i^2
m = ∑ i = 1 N r i 2 ,其中
r
i
∼
N
(
s
,
σ
2
)
r_i\sim {\mathcal N}(s,\sigma^2)
r i ∼ N ( s , σ 2 ) ,显然
m
σ
2
\frac{m}{\sigma^2}
σ 2 m 为自由度等于
N
N
N 的非中心卡方分布,其均值为
N
s
+
N
Ns+N
N s + N ,方差为
2
(
N
+
2
N
s
2(N+2Ns
2 ( N + 2 N s ),这里
N
s
=
μ
Ns=\mu
N s = μ 为非中心参数。同样当
N
N
N 足够大时,我们将其近似为高斯分布。此时
r
i
r_i
r i 的均值为
N
(
1
+
s
)
σ
2
=
(
N
+
μ
)
σ
2
N(1+s)\sigma^2=(N+\mu)\sigma^2
N ( 1 + s ) σ 2 = ( N + μ ) σ 2 ,方差为
2
N
(
1
+
2
s
)
σ
4
=
2
(
N
+
2
μ
)
σ
4
2N(1+2s)\sigma^4=2(N+2\mu)\sigma^4
2 N ( 1 + 2 s ) σ 4 = 2 ( N + 2 μ ) σ 4 ,因此有
m
∼
N
(
(
N
+
μ
)
σ
2
,
2
(
N
+
2
μ
)
σ
4
)
)
.
m\sim {\mathcal N}\left((N+\mu)\sigma^2,2(N+2\mu)\sigma^4)\right).
m ∼ N ( ( N + μ ) σ 2 , 2 ( N + 2 μ ) σ 4 ) ) .
当
N
N
N 足够大时,我们可以把卡方分布近似成均值为
m
r
=
{
N
σ
2
;
H
0
(
N
+
μ
)
σ
2
;
H
1
m_r=\left\{ \begin{aligned} N\sigma^2;&\quad H_0\\ (N+\mu)\sigma^2;&\quad H_1 \end{aligned}\right.
m r = { N σ 2 ; ( N + μ ) σ 2 ; H 0 H 1 方差为
σ
r
2
=
{
2
N
σ
4
;
H
0
2
(
N
+
2
μ
)
σ
4
;
H
1
\sigma^2_r=\left\{ \begin{aligned} 2N\sigma^4;&\quad H_0\\ 2(N+2\mu)\sigma^4;&\quad H_1 \end{aligned}\right.
σ r 2 = { 2 N σ 4 ; 2 ( N + 2 μ ) σ 4 ; H 0 H 1 的高斯分布,即
r
∼
N
(
m
r
,
σ
r
2
)
r\sim {\mathcal N}(m_r,\sigma^2_r)
r ∼ N ( m r , σ r 2 ) 。 若判决门限设为
η
\eta
η ,则检测概率为
P
d
=
P
r
(
r
>
η
∣
H
1
)
=
Q
[
η
−
E
(
r
∣
H
1
)
V
a
r
(
r
∣
H
1
)
]
,
P_d={\rm Pr}(r>\eta|H_1)=Q\left[\frac{\eta-{\rm E}(r|H_1)}{\sqrt{{\rm Var}(r|H_1)}} \right],
P d = P r ( r > η ∣ H 1 ) = Q [ V a r ( r ∣ H 1 )
η − E ( r ∣ H 1 ) ] , 虚警概率为
P
f
=
P
r
(
r
>
η
∣
H
0
)
=
Q
[
η
−
E
(
r
∣
H
0
)
V
a
r
(
r
∣
H
0
)
]
.
P_f={\rm Pr}(r>\eta|H_0)=Q\left[\frac{\eta-{\rm E}(r|H_0)}{\sqrt{{\rm Var}(r|H_0)}} \right].
P f = P r ( r > η ∣ H 0 ) = Q [ V a r ( r ∣ H 0 )
η − E ( r ∣ H 0 ) ] . 显然,不同门限值大小,会影响
P
d
P_d
P d 与
P
f
P_f
P f 。此外,我们还有漏检概率为
P
m
=
1
−
P
d
.
P_m=1-P_d.
P m = 1 − P d . 图1给出了漏检概率
P
m
P_m
P m 随虚检概率(
P
f
P_f
P f )变化的曲线,这里没有考虑衰落,即
α
=
1
\alpha=1
α = 1 。显然,
P
m
P_m
P m 随着
P
f
P_f
P f 的增加而减少。我们还可以发现,随着信噪比
γ
\gamma
γ 的增大,曲线下面积在减少,因此检测性能变好。 图1 漏检概率(
P
m
P_m
P m )与虚检概率(
P
f
P_f
P f )关系曲线