2019/12/26
DOA估计算法
较为详细的MUSIC算法原理及MATLAB实现_章子雎的博客-CSDN博客
DOA(Direction Of Arrival)波达方向定位技术主要有ARMA谱分析、最大似然法、熵谱分析法和特征分解法,特征分解法主要有MUSIC算法、
ESPRIT算法WSF算法等。
MUSIC (Multiple Signal Classification)算法,即多信号分类算法,由Schmidt等人于1979年提出。MUSIC算法是一种基于子空间分解的算法,
它利用信号子空间和噪声子空间的正交性,构建空间谱函数,通过谱峰搜索,估计信号的参数。对于声源定位来说,需要估计信号的DOA。
MUSIC算法对DOA的估计有很高的分辨率,且对麦克风阵列的形状没有特殊要求,因此应用十分广泛。
运用矩阵的定义,可得到更为简洁的表达式:
=
+
式中
https://blog.csdn.net/zhangziju/article/details/100730081
1/8
X
A
S
N
2019/12/26
= [
= [
= [
= [
(
(
1
1
(
1
),
),
),
1
),
(
(
(
(
(
2
(
2
(
2
), ...
), ...
), ...
2
), ...
(
)]
)]
,
,
)]
)]
,
。
(
较为详细的MUSIC算法原理及MATLAB实现_章子雎的博客-CSDN博客
为阵元的输出, 为方向响应向量, 是入射信号, 表示阵列噪声。
其中
=
2
有
= [
(
),
1
(
), ...
2
(
)] =
⎣⎢⎢⎢⎡
1
−
1⋮
1
−
2⋮
−
(
−1)
1
−
(
−1)
2
⋯⋯⋱⋯
⎦⎥⎥⎥⎤
1
−
D⋮
−1)
−
(
(
)
进行N点采样,要处理的问题就变成了通过输出信号
对
自然的将阵列信号看作是噪声干扰的若干空间谐波的叠加,从而将波达方向估计问题与谱估计联系起来。
对阵列输出X做相关处理,得到其协方差矩阵
的采样
{
(
) =
1, 2, ...,
(
)
}
估计信号源的波达方向角
,
1
...
2
,由此可以很
=
[
]H
其中 表示矩阵的共轭转置。
根据已假设信号与噪声互不相关、噪声为零均值白噪声,因此可得到:
=
其中
=
[(
=
2
+
[
+
)(
]H
称为信号相关矩阵
)
] =H
是噪声相关阵
[
]
+
[
] =H
+
2
是噪声功率
https://blog.csdn.net/zhangziju/article/details/100730081
2/8
X
x
t
x
t
x
t
M
T
S
S
t
S
t
S
t
D
T
A
a
θ
a
θ
a
θ
D
T
N
n
t
n
t
n
t
M
T
X
A
S
N
φ
k
s
i
n
θ
λ
π
d
k
A
a
θ
a
θ
a
θ
D
T
e
j
φ
e
j
M
φ
e
j
φ
e
j
M
φ
e
j
φ
e
j
M
φ
D
x
t
m
x
t
m
x
i
m
M
θ
θ
θ
D
R
x
E
X
X
H
R
x
E
A
S
N
A
S
N
A
E
S
S
A
H
H
E
N
N
A
R
A
S
H
R
N
R
s
E
S
S
R
N
σ
I
σ
2019/12/26
是
×
阶的单位矩阵
较为详细的MUSIC算法原理及MATLAB实现_章子雎的博客-CSDN博客
在实际应用中通常无法直接得到 ,能使用的只有样本的协方差矩阵:
^
, 是 的最大似然估计。
1 ∑
=
^
(
)
(
)
=1N
→ ∞
当采样数
先考虑理想情况,即无噪声的情况:
=
,对均匀线阵,矩阵A由
,他们是一致的,但实际情况将由于样本数有限而造成误差。根据矩阵特征分解的理论,可对阵列协方差矩阵进行特征分解,首
⎣⎢⎢⎢⎡
1
−
1⋮
⋯⋯⋱⋯
⎦⎥⎥⎥⎤
1
−
D⋮
−1)
−
(
1
−
2⋮
−
(
−1)
1
−
(
−1)
2
= [
(
),
1
(
), ...
2
(
)] =
所定义的范德蒙德矩阵,只要满足
=
,
=
,则他的各列相互独立。
若 为非奇异矩阵
(
) =
,各信号源两两不相干,且
>
,则
(
) =
,
由于
=
[
]H
,有:
=
即 为Hermite矩阵,它的特性是都是实数,又由于 为正定的,因此
… …
为半正定的,它有D个正特征值和
−
个零特征值。
再考虑有噪声存在的情况
=
+
2
由于
>2
0
, 为满秩阵,所以 有M个正实特征值
,
1
...
2
分别对应于M个特征向量
分别等于矩阵
的各特征值与 之和,其余
。又由于 为Hermite矩阵,所以各特征向量是正交的,即:
2
个特征值为 ,即 为 的最小特征值,它是
−
2
2
,
1
...
2
https://blog.csdn.net/zhangziju/article/details/100730081
=
0,
−
与信号有关的特征值只有D个,
=
维的,对应的特征向量
,
=
3/8
I
M
M
R
x
R
x
X
i
X
i
N
i
H
R
x
R
x
N
R
x
A
R
A
s
H
A
a
θ
a
θ
a
θ
D
T
e
j
φ
e
j
M
φ
e
j
φ
e
j
M
φ
e
j
φ
e
j
M
φ
D
θ
i
θ
i
j
j
R
s
R
a
n
k
R
s
D
M
D
r
a
n
d
A
R
A
s
H
D
R
x
E
X
X
R
s
H
R
x
R
s
R
s
A
R
A
H
s
M
D
R
x
A
R
A
s
H
σ
I
σ
R
x
R
x
λ
λ
λ
M
v
v
v
M
R
x
v
v
i
H
j
i
j
A
R
A
s
H
σ
M
D
σ
σ
R
M
D
v
i
i
2019/12/26
1, 2, ...,
中,也有D个是与信号有关的,另外
较为详细的MUSIC算法原理及MATLAB实现_章子雎的博客-CSDN博客
个是与噪声有关的,可利用特征分解的性质求出信号源的波达方向 。
−
MUSIC算法的原理及实现
通过对协方差矩阵的特征值分解,可得到如下结论:
将矩阵 的特征值进行从小到大的排序,即
≥1
≥2
... ≥
>
0
,其中D个较大的特征值对应于信号,
−
个较小的特征值对应于噪声。
矩阵 的属于这些特征值的特征向量也分别对应于各个信号和噪声,因此可把 的特征值(特征向量)划分为信号特征(特征向量)与噪声特征
(特征向量)。
设 为 的第 个特征值, 是与 个相对应的特征向量,有:
=
再设
=
2
是 的最小特征值
=
2
=
+ 1,
+ 2...
,
将
=
+
2
代入可得
=2
(
+
2
)
,
将其右边展开与左边比较得:
因
是
维的满秩矩阵,
(
)H
−1
存在;
=
0
∗
而 同样存在,则上式两边同乘以
−1
−1
(
)
−1
,
有:
−1
(
于是有
)
−1
=
0
https://blog.csdn.net/zhangziju/article/details/100730081
4/8
M
M
D
θ
k
R
x
λ
λ
λ
M
M
D
R
x
R
x
λ
i
R
x
i
v
i
λ
i
R
v
x
i
λ
v
i
i
λ
i
σ
R
x
R
v
x
i
σ
v
i
i
D
D
M
R
x
A
R
A
s
H
σ
I
σ
v
i
A
R
A
s
H
σ
I
v
i
A
R
A
v
s
H
i
A
A
H
D
D
A
A
R
s
R
A
A
A
s
H
H
R
A
A
A
A
R
A
v
s
H
H
s
H
i
2019/12/26
=
0,
=
+ 1,
+ 2, ...,
较为详细的MUSIC算法原理及MATLAB实现_章子雎的博客-CSDN博客
上式表明:噪声特征值所对应的特征向量(称为噪声特征向量) ,与矩阵 的列向量正交,而 的各列是与信号源的方向相对应的,这就是利用
噪声特征向量求解信号源方向的出发点。
用各噪声特征向量为例,构造一个噪声矩阵 :
=
[
,
+1
定义空间谱
]
, ...
)
:
+2
(
(
) =
1
(
)H
(
) =
∥
1
(
该式中分母是信号向量和噪声矩阵的内积,当
尖峰值,由该式,使 变化,通过寻找波峰来估计到达角。
2
)∥
(
)
和 的各列正交时,该分母为零,但由于噪声的存在,它实际上为一最小值,因此
(
)
有一
MUSIC算法实现的步骤
1.根据N个接收信号矢量得到下面协方差矩阵的估计值:
=
1 ∑
=1N
(
)
(
)
对上面得到的协方差矩阵进行特征分解
=
+
2
2.按特征值的大小排序将与信号个数 相等的特征值和对应的特征向量看做信号部分空间,将剩下的
间,得到噪声矩阵 :
−
个特征值和特征向量看做噪声部分空
=
=
[
0,
,
+1
=
+ 1,
+ 2, ...,
, ...
+2
]
https://blog.csdn.net/zhangziju/article/details/100730081
5/8
A
v
H
i
i
D
D
M
v
i
A
A
E
n
E
n
v
v
v
D
D
M
P
θ
m
u
P
θ
m
u
E
E
a
θ
a
θ
n
n
H
E
a
θ
n
H
a
θ
E
n
P
θ
m
u
θ
R
x
X
i
X
i
N
i
H
R
x
A
R
A
s
H
σ
I
D
M
D
E
n
A
v
H
i
i
D
D
M
E
n
v
v
v
D
D
M
2019/12/26
3.使 变化,按照式
(
) =
1
(
)
(
)
较为详细的MUSIC算法原理及MATLAB实现_章子雎的博客-CSDN博客
来计算谱函数,通过寻求峰值来得到波达方向的估计值。
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
clear; close all;
%%%%%%%% MUSIC for Uniform Linear Array%%%%%%%%
derad = pi/180; %角度->弧度
N = 8; % 阵元个数
M = 3; % 信源数目
theta = [-30 0 60]; % 待估计角度
snr = 10; % 信噪比
K = 512; % 快拍数
dd = 0.5; % 阵元间距
d=0:dd:(N-1)*dd;
A=exp(-1i*2*pi*d.'*sin(theta*derad)); %方向矢量
%%%%构建信号模型%%%%%
S=randn(M,K); %信源信号,入射信号
X=A*S; %构造接收信号
X1=awgn(X,snr,'measured'); %将白色高斯噪声添加到信号中
% 计算协方差矩阵
Rxx=X1*X1'/K;
% 特征值分解
[EV,D]=eig(Rxx); %特征值分解
EVA=diag(D)'; %将特征值矩阵对角线提取并转为一行
[EVA,I]=sort(EVA); %将特征值排序 从小到大
EV=fliplr(EV(:,I)); % 对应特征矢量排序
https://blog.csdn.net/zhangziju/article/details/100730081
6/8
θ
P
θ
m
u
a
θ
E
E
a
θ
H
n
n
H
2019/12/26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
较为详细的MUSIC算法原理及MATLAB实现_章子雎的博客-CSDN博客
% 遍历每个角度,计算空间谱
for iang = 1:361
angle(iang)=(iang-181)/2;
phim=derad*angle(iang);
a=exp(-1i*2*pi*d*sin(phim)).';
En=EV(:,M+1:N); % 取矩阵的第M+1到N列组成噪声子空间
Pmusic(iang)=1/(a'*En*En'*a);
end
Pmusic=abs(Pmusic);
Pmmax=max(Pmusic)
Pmusic=10*log10(Pmusic/Pmmax); % 归一化处理
h=plot(angle,Pmusic);
set(h,'Linewidth',2);
xlabel('入射角/(degree)');
ylabel('空间谱/(dB)');
set(gca, 'XTick',[-90:30:90]);
grid on;
实现结果
https://blog.csdn.net/zhangziju/article/details/100730081
7/8
2019/12/26
较为详细的MUSIC算法原理及MATLAB实现_章子雎的博客-CSDN博客
https://blog.csdn.net/zhangziju/article/details/100730081
8/8