假设传输两个信号:
s
0
[
n
]
s_0[n]
s0[n]和
s
1
[
n
]
s_1[n]
s1[n],传输过程中叠加了高斯白噪声
w
[
n
]
w[n]
w[n],
w
[
n
]
∼
N
(
0
,
σ
2
)
w[n]\sim N(0,\sigma^2)
w[n]∼N(0,σ2)。在接收端,我们需要判断接收到的
x
[
n
]
x[n]
x[n]到底是哪个信号,就有了如下两个假设:
H
0
:
x
[
n
]
=
s
0
[
n
]
+
w
[
n
]
H
1
:
x
[
n
]
=
s
1
[
n
]
+
w
[
n
]
H_0:x[n]=s_0[n]+w[n]\\ H_1:x[n]=s_1[n]+w[n]
H0:x[n]=s0[n]+w[n]H1:x[n]=s1[n]+w[n]
假设发出两种信号的概率相等,当
Pr
(
x
⃗
∣
H
1
)
Pr
(
x
⃗
∣
H
0
)
>
γ
=
Pr
(
H
0
)
Pr
(
H
1
)
=
1
\frac{\Pr(\vec{x}|H_1)}{\Pr(\vec{x}|H_0)}>\gamma=\frac{\Pr(H_0)}{\Pr(H_1)}=1
Pr(x∣H0)Pr(x∣H1)>γ=Pr(H1)Pr(H0)=1
时,拒绝
H
0
H_0
H0,认为我们接收到了
s
1
s_1
s1信号。
为了进一步探讨该问题,将
Pr
(
x
⃗
∣
H
i
)
\Pr(\vec{x}|H_i)
Pr(x∣Hi)展开,因为接收到的信号
x
[
n
]
x[n]
x[n]每一位都满足
x
[
n
]
∼
N
(
s
i
[
n
]
,
σ
2
)
x[n]\sim N(s_i[n],\sigma^2)
x[n]∼N(si[n],σ2),故整个信号
x
⃗
\vec{x}
x有
Pr
(
x
⃗
∣
H
i
)
=
1
(
2
π
σ
2
)
N
2
exp
[
−
1
2
σ
2
∑
n
=
0
N
−
1
(
x
[
n
]
−
s
i
[
n
]
)
2
]
\Pr(\vec{x}|H_i)=\frac{1}{(2\pi\sigma^2)^{\frac{N}{2}}}\exp[-\frac{1}{2\sigma^2}\sum_{n=0}^{N-1}(x[n]-s_i[n])^2]
Pr(x∣Hi)=(2πσ2)2N1exp[−2σ21n=0∑N−1(x[n]−si[n])2]
注意到对于两个假设
H
0
H_0
H0和
H
1
H_1
H1,
Pr
(
x
⃗
∣
H
i
)
\Pr(\vec{x}|H_i)
Pr(x∣Hi)的区别在于
∑
n
=
0
N
−
1
(
x
[
n
]
−
s
i
[
n
]
)
2
\sum_{n=0}^{N-1}(x[n]-s_i[n])^2
∑n=0N−1(x[n]−si[n])2,故设
D
i
2
=
∑
n
=
0
N
−
1
(
x
[
n
]
−
s
i
[
n
]
)
2
=
∣
∣
x
⃗
−
s
i
⃗
∣
∣
2
D_i^2=\sum_{n=0}^{N-1}(x[n]-s_i[n])^2=||\vec{x}-\vec{s_i}||^2
Di2=n=0∑N−1(x[n]−si[n])2=∣∣x−si∣∣2
不难看出,我们做出判断的依据就在于接收信号构成的向量
x
⃗
\vec{x}
x与两已知信号向量
s
i
⃗
\vec{s_i}
si之差的模长的平方。
再展开
D
i
2
D_i^2
Di2,得
D
i
2
=
∑
n
=
0
N
−
1
x
2
[
n
]
−
2
∑
n
=
0
N
−
1
x
[
n
]
s
i
[
n
]
+
∑
n
=
0
N
−
1
s
i
2
[
n
]
D_i^2=\sum_{n=0}^{N-1}x^2[n]-2\sum_{n=0}^{N-1}x[n]s_i[n]+\sum_{n=0}^{N-1}s_i^2[n]
Di2=n=0∑N−1x2[n]−2n=0∑N−1x[n]si[n]+n=0∑N−1si2[n]
发现
D
1
2
D_1^2
D12和
D
2
2
D_2^2
D22之间的差别在于
−
2
∑
n
=
0
N
−
1
x
[
n
]
s
i
[
n
]
+
∑
n
=
0
N
−
1
s
i
2
[
n
]
-2\sum_{n=0}^{N-1}x[n]s_i[n]+\sum_{n=0}^{N-1}s_i^2[n]
−2∑n=0N−1x[n]si[n]+∑n=0N−1si2[n],故定义
T
i
(
x
)
=
∑
n
=
0
N
−
1
x
[
n
]
s
i
[
n
]
−
1
2
∑
n
=
0
N
−
1
s
i
2
[
n
]
T_i(x)=\sum_{n=0}^{N-1}x[n]s_i[n]-\frac{1}{2}\sum_{n=0}^{N-1}s_i^2[n]
Ti(x)=n=0∑N−1x[n]si[n]−21n=0∑N−1si2[n]
其中,
∑
n
=
0
N
−
1
s
i
2
[
n
]
\sum_{n=0}^{N-1}s_i^2[n]
∑n=0N−1si2[n]为信号
s
i
s_i
si的能量,故
T
i
T_i
Ti又可写为
T
i
(
x
)
=
∑
n
=
0
N
−
1
x
[
n
]
s
i
[
n
]
−
1
2
ε
i
T_i(x)=\sum_{n=0}^{N-1}x[n]s_i[n]-\frac{1}{2}\varepsilon_i
Ti(x)=n=0∑N−1x[n]si[n]−21εi
可得出误检概率
P
e
P_e
Pe的表达式为
P
e
=
Pr
(
H
1
∣
H
0
)
Pr
(
H
0
)
+
Pr
(
H
0
∣
H
1
)
Pr
(
H
1
)
=
1
2
[
Pr
(
H
1
∣
H
0
)
+
Pr
(
H
0
∣
H
1
)
]
=
1
2
[
Pr
(
T
1
(
x
⃗
)
−
T
0
(
x
⃗
)
>
0
∣
H
0
)
+
Pr
(
T
0
(
x
⃗
)
−
T
1
(
x
⃗
)
>
0
∣
H
1
)
]
\begin{aligned} P_e&=\Pr(H_1|H_0)\Pr(H_0)+\Pr(H_0|H_1)\Pr(H_1)\\ &=\frac{1}{2}[\Pr(H_1|H_0)+\Pr(H_0|H_1)]\\ &=\frac{1}{2}[\Pr(T_1(\vec{x})-T_0(\vec{x})>0|H_0)+\Pr(T_0(\vec{x})-T_1(\vec{x})>0|H_1)] \end{aligned}
Pe=Pr(H1∣H0)Pr(H0)+Pr(H0∣H1)Pr(H1)=21[Pr(H1∣H0)+Pr(H0∣H1)]=21[Pr(T1(x)−T0(x)>0∣H0)+Pr(T0(x)−T1(x)>0∣H1)]
故定义检验
T
(
x
⃗
)
=
T
1
(
x
⃗
)
−
T
0
(
x
⃗
)
T(\vec{x})=T_1(\vec{x})-T_0(\vec{x})
T(x)=T1(x)−T0(x),即
T
(
x
⃗
)
=
∑
n
=
0
N
−
1
x
[
n
]
(
s
1
[
n
]
−
s
0
[
n
]
)
−
1
2
(
ε
1
−
ε
0
)
T(\vec{x})=\sum_{n=0}^{N-1}x[n](s_1[n]-s_0[n])-\frac{1}{2}(\varepsilon_1-\varepsilon_0)
T(x)=n=0∑N−1x[n](s1[n]−s0[n])−21(ε1−ε0)
通过计算可得
T
T
T的期望与方差:
E
(
T
∣
H
0
)
=
∑
n
=
0
N
−
1
s
0
[
n
]
(
s
1
[
n
]
−
s
0
[
n
]
)
−
1
2
(
ε
1
−
ε
0
)
=
−
1
2
∣
∣
s
1
⃗
−
s
0
⃗
∣
∣
2
E
(
T
∣
H
1
)
=
1
2
∣
∣
s
1
⃗
−
s
0
⃗
∣
∣
2
V
a
r
(
T
∣
H
0
)
=
V
a
r
(
∑
n
=
0
N
−
1
x
[
n
]
(
s
1
[
n
]
−
s
0
[
n
]
)
∣
H
0
)
=
σ
2
∣
∣
s
1
⃗
−
s
0
⃗
∣
∣
2
V
a
r
(
T
∣
H
1
)
=
σ
2
∣
∣
s
1
⃗
−
s
0
⃗
∣
∣
2
\begin{aligned} &E(T|H_0)=\sum_{n=0}^{N-1}s_0[n](s_1[n]-s_0[n])-\frac{1}{2}(\varepsilon_1-\varepsilon_0)=-\frac{1}{2}||\vec{s_1}-\vec{s_0}||^2\\ &E(T|H_1)=\frac{1}{2}||\vec{s_1}-\vec{s_0}||^2\\ &Var(T|H_0)=Var(\sum_{n=0}^{N-1}x[n](s_1[n]-s_0[n])|H_0)=\sigma^2||\vec{s_1}-\vec{s_0}||^2\\ &Var(T|H_1)=\sigma^2||\vec{s_1}-\vec{s_0}||^2 \end{aligned}
E(T∣H0)=n=0∑N−1s0[n](s1[n]−s0[n])−21(ε1−ε0)=−21∣∣s1−s0∣∣2E(T∣H1)=21∣∣s1−s0∣∣2Var(T∣H0)=Var(n=0∑N−1x[n](s1[n]−s0[n])∣H0)=σ2∣∣s1−s0∣∣2Var(T∣H1)=σ2∣∣s1−s0∣∣2
所以,
(
T
∣
H
0
)
(T|H_0)
(T∣H0)和
(
T
∣
H
1
)
(T|H_1)
(T∣H1)均满足正态分布,且
(
T
∣
H
0
)
∼
N
(
−
1
2
∣
∣
s
1
⃗
−
s
0
⃗
∣
∣
2
,
σ
2
∣
∣
s
1
⃗
−
s
0
⃗
∣
∣
2
)
(
T
∣
H
1
)
∼
N
(
1
2
∣
∣
s
1
⃗
−
s
0
⃗
∣
∣
2
,
σ
2
∣
∣
s
1
⃗
−
s
0
⃗
∣
∣
2
)
(T|H_0)\sim N(-\frac{1}{2}||\vec{s_1}-\vec{s_0}||^2,\sigma^2||\vec{s_1}-\vec{s_0}||^2)\\ (T|H_1)\sim N(\frac{1}{2}||\vec{s_1}-\vec{s_0}||^2,\sigma^2||\vec{s_1}-\vec{s_0}||^2)
(T∣H0)∼N(−21∣∣s1−s0∣∣2,σ2∣∣s1−s0∣∣2)(T∣H1)∼N(21∣∣s1−s0∣∣2,σ2∣∣s1−s0∣∣2)
故误检概率可写为
P
e
=
Pr
(
T
>
0
∣
H
0
)
Pr
(
H
0
)
+
Pr
(
T
<
0
∣
H
1
)
Pr
(
H
1
)
=
1
2
[
Pr
(
T
>
0
∣
H
0
)
+
Pr
(
T
<
0
∣
H
1
)
]
=
Φ
(
−
1
2
∣
∣
s
1
⃗
−
s
0
⃗
∣
∣
2
σ
2
∣
∣
s
1
⃗
−
s
0
⃗
∣
∣
2
)
=
Φ
(
−
1
2
∣
∣
s
1
⃗
−
s
0
⃗
∣
∣
2
σ
2
)
\begin{aligned} P_e&=\Pr(T>0|H_0)\Pr(H_0)+\Pr(T<0|H_1)\Pr(H_1)\\ &=\frac{1}{2}[\Pr(T>0|H_0)+\Pr(T<0|H_1)]\\ &=\Phi(-\frac{\frac{1}{2}||\vec{s_1}-\vec{s_0}||^2}{\sqrt{\sigma^2||\vec{s_1}-\vec{s_0}||^2}})\\ &=\Phi(-\frac{1}{2}\sqrt{\frac{||\vec{s_1}-\vec{s_0}||^2}{\sigma^2}}) \end{aligned}
Pe=Pr(T>0∣H0)Pr(H0)+Pr(T<0∣H1)Pr(H1)=21[Pr(T>0∣H0)+Pr(T<0∣H1)]=Φ(−σ2∣∣s1−s0∣∣221∣∣s1−s0∣∣2)=Φ(−21σ2∣∣s1−s0∣∣2)
可以看出,误检概率与两信号向量之差的模长成正比,与噪音的标准差成反比。
我们可以用matlab验证该推导结果,每次检验均生成两个随机信号 s 1 , s 2 s_1,s_2 s1,s2,分别叠加噪音计算 T T T值,判断是否误检。
在这里只生成了长度为2的信号,一是因为长度太长每次检验的速度会增加,程序的总运行时间暴增;其次若信号长度太长,很难随机生成
∣
∣
s
1
⃗
−
s
0
⃗
∣
∣
2
||\vec{s_1}-\vec{s_0}||^2
∣∣s1−s0∣∣2较小的信号,使统计不准。
test.m
function res=test(sigma)
N=2;
s0=2*rand(1,N)-1;
s1=2*rand(1,N)-1;
w0=s0+sigma*randn(1,N);
w1=s1+sigma*randn(1,N);
T0=0;
T1=0;
E0=0;
E1=0;
for i=1:1:N
E0=E0+s0(i)^2;
E1=E1+s1(i)^2;
T0=T0+w0(i)*(s1(i)-s0(i));
T1=T1+w1(i)*(s1(i)-s0(i));
end
T0=T0-(E1-E0)/2;
T1=T1-(E1-E0)/2;
p=(s1-s0)*(s1-s0)';
%p=2*s1*s0'/(s1*s1'+s0*s0');
res=[p,(T0>0)+(T1<0)];
之后,统计两向量之差的模长与误检次数的关系,与推导结果比对,可得
在
∣
∣
s
1
⃗
−
s
0
⃗
∣
∣
2
||\vec{s_1}-\vec{s_0}||^2
∣∣s1−s0∣∣2偏大的部分,统计很不准确,这是因为对于随机生成的长度为2的信号,
∣
∣
s
1
⃗
−
s
0
⃗
∣
∣
2
||\vec{s_1}-\vec{s_0}||^2
∣∣s1−s0∣∣2那么大的样本很少,统计结果就会非常不准确。可以考虑使用蒙特卡洛法消除部分偏差。
统计及生成gif代码:
statistic.m
clc;clear;
cnt=1;
for sigma=0.2:0.01:0.6
str='σ=';
N=10000000;
d=1000;
y=zeros(1,d);
tot=zeros(1,d);
data=zeros(N,2);
hb=0;
lb=400;
for i=1:1:N
a=test(sigma);
data(i,1)=a(1);
data(i,2)=a(2);
hb=max(hb,a(1));
lb=min(lb,a(1));
end
hb=hb+0.01;
sep=(hb-lb)/d;
x=lb+sep/2:sep:hb-sep/2;
for i=1:1:N
idx=floor((data(i,1)-lb)/sep)+1;
y(idx)=y(idx)+data(i,2);
tot(idx)=tot(idx)+2;
end
y=y./tot;
figure
plot(x,y,"linewidth",2);
hold on;
pd=makedist("Normal","mu",0,"sigma",1);
y1=cdf(pd,trans(x,sigma));
plot(x,y1,"linewidth",2);
xlabel("$||s_1-s_0||^2$","interpreter","latex","fontsize",16);
ylabel("prE","fontsize",16);
str1=[str num2str(sigma)];
title(str1,"fontsize",18);
print(1,'-dbmp',sprintf('img/%d',cnt));
cnt=cnt+1;
close;
end
%xlim([0.1 0.7]);
%ylim([-0.0025 0.03]);
%依次读取生成的所有图片
for j=1:41
%获取当前图片
A=imread(sprintf('img/%d.bmp',j));
[I,map]=rgb2ind(A,256);
%生成gif,并保存
if(j==1)
imwrite(I,map,'movefig.gif','DelayTime',0.1,'LoopCount',Inf)
else
imwrite(I,map,'movefig.gif','WriteMode','append','DelayTime',0.1)
end
end
trans.m
function re=trans(x,sigma)
re=-sqrt(x/sigma^2)/2;