一、MUSIC常用方法介绍
在讲MUSIC方法之前 ,老师提到了阵列处理,是一种传感器在空域上的扩展,利用阵列接收器做DOA有两个关键性假设:
1、远场平面波假设
2、窄带信号(信号带宽
≪
\ll
≪载频,ps:在《随机信号分析教程》这本书中提到,窄带信号的带宽近似小于等于载频的10%,20%以下也行)
二、MUSIC方法推导
1、给出均匀线元阵列图像如下:
假设输入共有M个信号且互不相关,分别记为
S
1
(
t
)
、
S
2
(
t
)
.
.
.
S
M
(
t
)
S_1(t)、S_2(t)...S_M(t)
S1(t)、S2(t)...SM(t),阵列接收器(阵元)共有N个,分别记为
X
0
(
t
)
、
X
1
(
t
)
、
.
.
.
X
N
−
1
(
t
)
X_0(t)、X_1(t)、...X_{N-1}(t)
X0(t)、X1(t)、...XN−1(t)。我们给出M个信号到第一个阵元上的角度为:
θ
1
、
θ
2
、
.
.
.
θ
M
\theta_1、\theta_2、...\theta_M
θ1、θ2、...θM。
已知基本上都有了,现在我们需要在推导一下每个阵元得到的信号值
X
i
(
t
)
X_i(t)
Xi(t)与输入信号之间的关系。
假设输入信号为
S
(
t
)
S(t)
S(t),其在一个阵元的入射角为
θ
\theta
θ,那么其到达第二个信号的输入则为
S
(
t
−
T
)
S(t-T)
S(t−T),T为一个时延,参考下图:
可以得到
T
=
d
∗
s
i
n
θ
c
T=\frac{d*sin\theta}{c}
T=cd∗sinθ,其中c代表光速,那么对应
S
(
t
)
S(t)
S(t)在第k个阵元上的输入则为
S
(
t
−
k
∗
d
∗
s
i
n
θ
c
)
S(t-\frac{k*d*sin\theta}{c})
S(t−ck∗d∗sinθ),
利用实部加虚部表示可以得到以及窄带信号对时延不敏感可以得到如下,
S
(
t
−
k
∗
d
∗
s
i
n
θ
c
)
=
S
(
t
)
e
x
p
(
j
2
π
f
(
k
∗
d
∗
s
i
n
θ
c
)
)
S(t-\frac{k*d*sin\theta}{c})=S(t)exp(j2{\pi}f(\frac{k*d*sin\theta}{c}))
S(t−ck∗d∗sinθ)=S(t)exp(j2πf(ck∗d∗sinθ))
证明过程如下:
S
(
t
)
=
S
~
(
t
)
e
x
p
(
−
j
2
π
f
t
)
S(t)=\tilde{S}(t)exp(-j2{\pi}ft)
S(t)=S~(t)exp(−j2πft)
S
(
t
−
T
)
=
S
~
(
t
−
T
)
e
x
p
(
−
j
2
π
f
t
)
e
x
p
(
j
2
π
f
T
)
S(t-T)=\tilde{S}(t-T)exp(-j2{\pi}ft)exp(j2{\pi}fT)
S(t−T)=S~(t−T)exp(−j2πft)exp(j2πfT)
由于窄带信号对时延不敏感:
S
~
(
t
−
T
)
≈
S
~
(
t
)
\tilde{S}(t-T)\approx\tilde{S}(t)
S~(t−T)≈S~(t)故:
S
(
t
−
T
)
=
S
~
(
t
−
T
)
e
x
p
(
−
j
2
π
f
t
)
e
x
p
(
j
2
π
f
T
)
=
S
~
(
t
)
e
x
p
(
−
j
2
π
f
t
)
e
x
p
(
j
2
π
f
T
)
=
S
(
t
)
e
x
p
(
j
2
π
f
T
)
S(t-T)=\tilde{S}(t-T)exp(-j2{\pi}ft)exp(j2{\pi}fT)\\=\tilde{S}(t)exp(-j2{\pi}ft)exp(j2{\pi}fT)\\=S(t)exp(j2{\pi}fT)
S(t−T)=S~(t−T)exp(−j2πft)exp(j2πfT)=S~(t)exp(−j2πft)exp(j2πfT)=S(t)exp(j2πfT)
利用
c
=
λ
∗
f
c=\lambda*f
c=λ∗f可以进一步化简为
S
(
t
)
e
x
p
(
j
2
π
(
k
∗
d
∗
s
i
n
θ
λ
)
)
S(t)exp(j2{\pi}(\frac{k*d*sin\theta}{\lambda}))
S(t)exp(j2π(λk∗d∗sinθ))
综合上面的已知我们给出N个阵元接收到的信号公式如下:
X
0
(
t
)
=
S
1
(
t
)
+
S
2
(
t
)
.
.
.
+
S
M
(
t
)
X_0(t)=S_1(t)+S_2(t)...+S_M(t)
X0(t)=S1(t)+S2(t)...+SM(t)
X
1
(
t
)
=
S
1
(
t
)
e
x
p
(
j
2
π
d
λ
s
i
n
θ
1
)
+
S
2
(
t
)
e
x
p
(
j
2
π
d
λ
s
i
n
θ
2
)
.
.
.
+
S
M
(
t
)
e
x
p
(
j
2
π
d
λ
s
i
n
θ
M
)
X_1(t)=S_1(t)exp(j2{\pi}\frac{d}{\lambda}sin\theta_1)+S_2(t)exp(j2{\pi}\frac{d}{\lambda}sin\theta_2)...+S_M(t)exp(j2{\pi}\frac{d}{\lambda}sin\theta_M)
X1(t)=S1(t)exp(j2πλdsinθ1)+S2(t)exp(j2πλdsinθ2)...+SM(t)exp(j2πλdsinθM)
⋮
\vdots
⋮
X
N
−
1
(
t
)
=
S
1
(
t
)
e
x
p
(
j
(
N
−
1
)
2
π
d
λ
s
i
n
θ
1
)
+
S
2
(
t
)
e
x
p
(
(
N
−
1
)
j
2
π
d
λ
s
i
n
θ
2
)
.
.
.
+
S
M
(
t
)
e
x
p
(
j
(
N
−
1
)
2
π
d
λ
s
i
n
θ
M
)
X_{N-1}(t)=S_1(t)exp(j(N-1)2{\pi}\frac{d}{\lambda}sin\theta_1)+S_2(t)exp((N-1)j2{\pi}\frac{d}{\lambda}sin\theta_2)...+S_M(t)exp(j(N-1)2{\pi}\frac{d}{\lambda}sin\theta_M)
XN−1(t)=S1(t)exp(j(N−1)2πλdsinθ1)+S2(t)exp((N−1)j2πλdsinθ2)...+SM(t)exp(j(N−1)2πλdsinθM)
我们用
Φ
i
\Phi_i
Φi代表
2
π
d
λ
s
i
n
θ
i
2{\pi}\frac{d}{\lambda}sin\theta_i
2πλdsinθi,对上式进行化简,并利用向量形式表示,可以得到如下结果:
X=
(
X
1
(
t
)
X
2
(
t
)
.
.
.
X
N
−
1
(
t
)
)
T
(X_1(t){\ }X_2(t){\ }...{\ }X_{N-1}(t))^T
(X1(t) X2(t) ... XN−1(t))T
A
(
θ
)
=
(
1
1
.
.
.
1
e
x
p
(
j
Φ
1
)
e
x
p
(
j
Φ
2
)
.
.
.
e
x
p
(
j
Φ
M
)
⋮
⋮
.
.
.
⋮
e
x
p
(
j
(
N
−
1
)
Φ
1
)
e
x
p
(
j
(
N
−
1
)
Φ
2
)
.
.
.
e
x
p
(
j
(
N
−
1
)
Φ
M
)
)
\bold{A(\theta)}=\left(\begin{matrix}1&1&...&1\\exp(j\Phi_1)&exp(j\Phi_2)&...&exp(j\Phi_M)\\{\vdots}&{\vdots}&...&\vdots\\exp(j(N-1)\Phi_1)&exp(j(N-1)\Phi_2)&...&exp(j(N-1)\Phi_M)\end{matrix}\right)
A(θ)=
1exp(jΦ1)⋮exp(j(N−1)Φ1)1exp(jΦ2)⋮exp(j(N−1)Φ2)............1exp(jΦM)⋮exp(j(N−1)ΦM)
S
=
(
S
1
(
t
)
S
2
(
t
)
.
.
.
S
M
(
t
)
)
T
\bold{S}=(S_1(t){\ }S_2(t){\ }...S_M(t))^T
S=(S1(t) S2(t) ...SM(t))T
构造模型:
X
=
A
(
θ
)
∗
S
+
ε
\bold{X}=\bold{A(\theta)}*\bold{S}+\bold{\varepsilon}
X=A(θ)∗S+ε
目前来讲是“一缺三”,左侧数据是我们可以获取的,右侧阵列流行是我们要找的,输入信号未知,噪声未知,MUSIC方法将借用PCA的观点来求出最终的阵列流行
A
(
θ
)
\bold{A(\theta)}
A(θ)。
对模型两侧采取相关,并进行特征分解可以得到如下结果:
R
x
=
E
(
A
(
θ
)
S
+
ε
)
(
A
(
θ
)
S
+
ε
)
H
=
A
(
θ
)
R
s
A
(
θ
)
H
+
σ
2
I
=
U
Λ
U
H
R_x=E(\bold{A(\theta)}\bold{S}+\bold{\varepsilon})(\bold{A(\theta)}\bold{S}+\bold{\varepsilon})^H=\bold{A(\theta)}R_s\bold{A(\theta)}^H+\sigma^2I=U\Lambda{U^H}
Rx=E(A(θ)S+ε)(A(θ)S+ε)H=A(θ)RsA(θ)H+σ2I=UΛUH
(
U
s
U
ε
)
(
Λ
s
0
0
Λ
ε
)
(
U
s
H
U
ε
H
)
=
A
(
θ
)
R
s
A
(
θ
)
H
+
σ
2
I
(U_s{\ }U_{\varepsilon})\left(\begin{matrix}\Lambda_s&0\\0&\Lambda_{\varepsilon}\end{matrix}\right)\left(\begin{matrix}U_s^H\\U_{\varepsilon}^H\end{matrix}\right)=\bold{A(\theta)}R_s\bold{A(\theta)}^H+\sigma^2I
(Us Uε)(Λs00Λε)(UsHUεH)=A(θ)RsA(θ)H+σ2I
U
s
Λ
s
U
s
H
+
U
ε
Λ
ε
U
ε
H
=
U
s
Λ
s
U
s
H
+
σ
2
I
2
=
A
(
θ
)
R
s
A
(
θ
)
H
+
σ
2
I
1
U_s\Lambda_sU_s^H+U_{\varepsilon}\Lambda_{\varepsilon}U_{\varepsilon}^H=U_s\Lambda_sU_s^H+\sigma^2I_2=\bold{A(\theta)}R_s\bold{A(\theta)}^H+\sigma^2I_1
UsΛsUsH+UεΛεUεH=UsΛsUsH+σ2I2=A(θ)RsA(θ)H+σ2I1
注意上边的
I
1
I_1
I1和
I
2
I_2
I2的行列是不一样的,因此进一步化简可得到如下结果:
(
U
s
U
ε
)
(
Λ
s
−
σ
2
I
0
0
0
)
(
U
s
H
U
ε
H
)
=
A
(
θ
)
R
s
A
(
θ
)
H
(U_s{\ }U_{\varepsilon})\left(\begin{matrix}\Lambda_s-\sigma^2I&0\\0&0\end{matrix}\right)\left(\begin{matrix}U_s^H\\U_{\varepsilon}^H\end{matrix}\right)=\bold{A(\theta)}R_s\bold{A(\theta)}^H
(Us Uε)(Λs−σ2I000)(UsHUεH)=A(θ)RsA(θ)H
U
s
(
Λ
s
−
σ
2
I
)
U
s
H
=
A
(
θ
)
R
s
A
(
θ
)
H
U_s(\Lambda_s-\sigma^2I)U_s^H=\bold{A(\theta)}R_s\bold{A(\theta)}^H
Us(Λs−σ2I)UsH=A(θ)RsA(θ)H
由于我们给出的信号源是假定互不相关的,因此
R
s
R_s
Rs满秩,同样的左侧
Λ
s
−
σ
2
I
\Lambda_s-\sigma^2I
Λs−σ2I满秩,利用子空间方法,下面论证
S
p
a
c
e
(
U
s
)
=
S
p
a
c
e
(
A
(
θ
)
)
Space(U_s)=Space(A(\theta))
Space(Us)=Space(A(θ))
利用二次型,我们得到如下公式:
y
H
U
s
(
Λ
s
−
σ
2
I
)
U
s
H
y
=
y
H
A
(
θ
)
R
s
A
(
θ
)
H
y
y^HU_s(\Lambda_s-\sigma^2I)U_s^Hy=y^H\bold{A(\theta)}R_s\bold{A(\theta)}^Hy
yHUs(Λs−σ2I)UsHy=yHA(θ)RsA(θ)Hy
由于中间部分满秩,如果
U
s
H
y
=
0
U_s^Hy=0
UsHy=0那么
A
(
θ
)
H
y
=
0
\bold{A(\theta)}^Hy=0
A(θ)Hy=0,因此可以得到
k
e
r
(
A
H
(
θ
)
)
=
k
e
r
(
U
s
)
ker(\bold{A}^H(\theta))=ker(U_s)
ker(AH(θ))=ker(Us),故
s
p
a
n
(
A
H
(
θ
)
)
=
s
p
a
n
(
U
s
)
span(\bold{A}^H(\theta))=span(U_s)
span(AH(θ))=span(Us)
上面用到的是零空间+列空间是整个空间维数为n固定不变导出的。
现在我们知道了
S
p
a
c
e
(
U
s
)
=
S
p
a
c
e
(
A
(
θ
)
)
Space(U_s)=Space(A(\theta))
Space(Us)=Space(A(θ)),同时由于信号空间和噪声空间不相关即正交,那么我们容易得到,从信号空间中取得的向量均正交与噪声空间,又因为信号空间和阵列流行空间是同一个空间,因此我们可以得到阵列流行空间中的人一个向量都正交于噪声空间,即有
A
(
θ
)
\bold{A}(\theta)
A(θ)的每一组列向量都正交于
U
ε
U_{\varepsilon}
Uε,具体公式如下:
G
(
θ
)
=
0
∑
k
=
1
M
∣
∣
A
k
(
θ
)
H
∗
U
ε
∣
∣
=
0
G(\theta)=0\sum_{k=1}^M||A_k(\theta)^H*U_{\varepsilon}||=0
G(θ)=0∑k=1M∣∣Ak(θ)H∗Uε∣∣=0
MUSIC方法是利用
1
G
(
θ
)
\frac{1}{G(\theta)}
G(θ)1在
θ
∈
(
0
,
180
)
\theta\in(0,180)
θ∈(0,180)上进行扫描,形成伪谱,峰值高的地方对应的就是信号源的角度,这种方法的分辨率很高,但是仍然有一定的缺点。
首先由于我们并不知道信号源的数量,因此我们需要通过数据相关阵的特征值进行划分(含有信号能量的特征值大,只有噪声的特征值小),这种划分对信噪比要求很高,倘若信噪比较低,就会导致特征值没有明显的大小区别,没法进行信号源数量的确定,一般情况下,MUSIC方法要求5-6db至少,效果好的话最好能到8-10db,这就给实验环境形成了很多限制 。
以上就是这学期第一节课的课程内容,希望期末自己能再瞅瞅。