现代谱估计的原理及MATLAB仿真（二）（AR模型法、MVDR法、MUSIC法）
现代谱估计的原理及MATLAB仿真AR参数模型法（参数模型功率谱估计）、MVDR法（最小方差无失真响应法）、MUSIC法（多重信号分类法）
文章目录
- 前言
- 一、AR参数模型
- 1 原理
- 2 MATLAB仿真
- 二、MVDR法
- 1 原理
- 2 MATLAB仿真
- 三、MUSIC法
- 1 原理
- 2 MATLAB仿真
- 四、3种方法的对比
- 五、MATLAB源代码
前言
\;\;\;\;\; 现代功率谱估计方法包括AR参数模型法（参数模型功率谱估计）、MVDR法（最小方差无失真响应法）、MUSIC法（多重信号分类法）。本文在总结各种方法的原理后在MATLAB平台上完成了仿真，完成了对信号频率的估计，仿真不同大小的阶数对信号频率估计的影响以及这三种方法之间的对比。
提示：以下是本篇文章正文内容，转载请附上链接！
一、AR参数模型
1 原理
\;\;\;\;\;
AR的全称是auto-regressive（自回归）参数模型。该模型功率谱估计的基本思想是，认为随机过程
u
(
n
)
u(n)
u(n)就是白噪声通过LTI离散时间系统得到的响应，利用观测样本值估计出模型的参数（即得到频率响应
H
(
ω
)
H(\omega)
H(ω)）,也就得到了随机过程
u
(
n
)
u(n)
u(n)的功率谱
S
(
ω
)
=
∣
H
(
ω
)
∣
2
σ
2
S(\omega)=\mid H(\omega)\mid^2\sigma^2
S(ω)=∣H(ω)∣2σ2。
\;\;\;\;\;
设
v
(
n
)
v(n)
v(n)是零均值、方差为
σ
2
\sigma^2
σ2的白噪声。则
p
p
p 阶 AR(
p
p
p) 模型的输人
v
(
n
v(n
v(n)和输出
u
(
n
)
u(n)
u(n)满足如下差分方程：
u
(
n
)
=
−
∑
k
=
1
p
a
k
u
(
n
−
k
)
+
v
(
n
)
u(n)=-\sum_{k=1}^pa_ku(n-k)+v(n)
u(n)=−k=1∑paku(n−k)+v(n)那么AR模型的正则方程式可表示为
r
u
(
0
)
+
a
1
r
u
(
−
1
)
+
⋅
⋅
⋅
+
a
p
r
u
(
−
p
)
=
σ
2
r_u(0)+a_1r_u(-1)+\cdotp\cdotp\cdotp+a_pr_u(-p)=\sigma^2
ru(0)+a1ru(−1)+⋅⋅⋅+apru(−p)=σ2 和
[
r
u
(
0
)
r
u
(
−
1
)
⋯
r
u
(
−
p
+
1
)
r
u
(
1
)
r
u
(
0
)
⋯
r
u
(
−
p
+
2
)
⋮
⋮
⋱
⋮
r
u
(
p
−
1
)
r
u
(
p
−
2
)
⋯
r
u
(
0
)
]
[
a
1
a
2
⋮
a
p
]
=
[
−
r
u
(
1
)
−
r
u
(
2
)
⋮
−
r
u
(
p
)
]
\begin{bmatrix}r_u(0)&r_u(-1)&\cdots&r_u(-p+1)\\r_u(1)&r_u(0)&\cdots&r_u(-p+2)\\\vdots&\vdots&\ddots&\vdots\\r_u(p-1)&r_u(p-2)&\cdots&r_u(0)\end{bmatrix}\begin{bmatrix}a_1\\a_2\\\vdots\\a_p\end{bmatrix}=\begin{bmatrix}-r_u(1)\\-r_u(2)\\\vdots\\-r_u(p)\end{bmatrix}
ru(0)ru(1)⋮ru(p−1)ru(−1)ru(0)⋮ru(p−2)⋯⋯⋱⋯ru(−p+1)ru(−p+2)⋮ru(0)
a1a2⋮ap
=
−ru(1)−ru(2)⋮−ru(p)
上式可简单地表示为
R
p
θ
p
=
−
r
p
{\boldsymbol{R}}_p\boldsymbol{\theta}_p=-\boldsymbol{r}_p
Rpθp=−rp 其中，
r
p
=
[
r
u
(
1
)
r
u
(
2
)
⋯
r
u
(
p
)
]
T
\boldsymbol r_{p}=\begin{bmatrix}r_{u}(1)&r_{u}(2)&\cdots&r_{u}(p)\end{bmatrix}^{\mathrm{T}}
rp=[ru(1)ru(2)⋯ru(p)]T。因矩阵
R
ρ
{\boldsymbol{R}}_{\rho}
Rρ 是非奇异的，对上式两边左乘
R
p
−
1
{\boldsymbol{R}}_p^{-1}
Rp−1，有
θ
p
=
−
R
p
−
1
r
p
\boldsymbol{\theta}_p=-{\boldsymbol{R}}_p^{-1}\boldsymbol r_{p}
θp=−Rp−1rp 将
θ
p
\boldsymbol{\theta}_p
θp 代人含有
σ
2
\sigma^2
σ2的那个式子中即可得到
σ
2
\sigma^2
σ2。
\;\;\;\;\;
随机过程
u
(
n
)
u(n)
u(n)的功率谱可由下式给出：
S
A
R
(
ω
)
=
σ
2
∣
1
+
∑
k
=
1
p
a
k
e
−
j
ω
k
∣
2
S_{\mathrm{AR}}(\omega)=\frac{\sigma^2}{\left|1+\sum_{k=1}^pa_k\mathrm{e}^{-\mathrm{j}\omega k}\right|^2}
SAR(ω)=∣1+∑k=1pake−jωk∣2σ2
2 MATLAB仿真
\;\;\;\;\;
仿真采用的信号为复正弦信号。参数设置和仿真结果如下：
\;\;\;\;\;
从以上结果可看出，阶数为8时，20MHz和22MHz的两个信号不能被区分，但当阶数为16时，20MHz和22MHz的两个信号可以被区分，所以估计信号频率时阶数的选择十分重要。
二、MVDR法
1 原理
\;\;\;\;\;
最小方差无失真响应(MVDR , minimum variance distortionless response)算法，是另一类信号频率估计方法。考虑有
M
M
M 个权系数的横向滤波器，如下图所示。
\;\;\;\;\;
滤波器的输入为随机过程
x
(
n
)
x(n)
x(n) ，输出为
y
(
n
)
=
∑
i
=
0
M
−
1
w
i
∗
x
(
n
−
i
)
y(n)=\sum_{i=0}^{M-1}w_i^*x(n-i)
y(n)=i=0∑M−1wi∗x(n−i)其中
,
w
i
,w_i
,wi 表示横向滤波器的权系数。定义输人信号向量和权向量分别为
x
(
n
)
=
[
x
(
n
)
x
(
n
−
1
)
⋯
x
(
n
−
M
+
1
)
]
T
w
=
[
w
0
w
1
⋯
w
M
−
1
]
T
\begin{aligned}&\boldsymbol{x}(n)=\begin{bmatrix}x(n)&x(n-1)&\cdots&x(n-M+1)\end{bmatrix}^\mathrm{T}\\&\boldsymbol{w}=\begin{bmatrix}w_0&w_1&\cdots&w_{M-1}\end{bmatrix}^\mathrm{T}\end{aligned}
x(n)=[x(n)x(n−1)⋯x(n−M+1)]Tw=[w0w1⋯wM−1]T 则输出可表示为
y
(
n
)
=
w
H
x
(
n
)
=
x
T
(
n
)
w
∗
y(n)=\boldsymbol w^\mathrm{H}\boldsymbol x(n)=\boldsymbol x^\mathrm{T}(n)\boldsymbol w^*
y(n)=wHx(n)=xT(n)w∗ 信号
y
(
n
)
y(n)
y(n)的平均功率可以表示为
P
=
E
{
∣
y
(
n
)
∣
2
}
=
E
{
w
H
x
(
n
)
x
H
(
n
)
w
}
=
w
H
R
w
P= \mathbb{E} \{ \mid y( n) \mid ^2\} = \mathbb{E} \{ \boldsymbol w^\mathrm{H} \boldsymbol x( n) \boldsymbol x^\mathrm{H} ( n) \boldsymbol w\} = \boldsymbol w^\mathrm{H} \boldsymbol {Rw}
P=E{∣y(n)∣2}=E{wHx(n)xH(n)w}=wHRw其中，矩阵
R
∈
C
M
×
M
\boldsymbol R\in\mathbb{C}^{M\times M}
R∈CM×M为向量
x
(
n
)
\boldsymbol x(n)
x(n)的
M
M
M 维自相关矩阵，即
R
=
E
{
x
(
n
)
x
H
(
n
)
}
=
[
r
(
0
)
r
(
1
)
⋯
r
(
M
−
1
)
r
(
−
1
)
r
(
0
)
⋯
r
(
M
−
2
)
⋮
⋮
⋱
⋮
r
(
1
−
M
)
r
(
2
−
M
)
⋯
r
(
0
)
]
\begin{gathered}\boldsymbol R=E\{\boldsymbol x(n)\boldsymbol x^H(n)\}=\begin{bmatrix}r(0)&r(1)&\cdots&r(M-1)\\r(-1)&r(0)&\cdots&r(M-2)\\\vdots&\vdots&\ddots&\vdots\\r(1-M)&r(2-M)&\cdots&r(0)\end{bmatrix}\end{gathered}
R=E{x(n)xH(n)}=
r(0)r(−1)⋮r(1−M)r(1)r(0)⋮r(2−M)⋯⋯⋱⋯r(M−1)r(M−2)⋮r(0)
\;\;\;\;\;
当权向量取
w
^
M
V
D
R
=
R
^
−
1
a
(
ω
)
a
H
(
ω
)
R
^
−
1
a
(
ω
)
\hat{w}_{\mathrm{MVDR}}=\frac{\hat{\boldsymbol R}^{-1}\boldsymbol{a}(\omega)}{\boldsymbol{a}^{\mathrm{H}}(\omega)\boldsymbol{\hat{R}}^{-1}\boldsymbol{a}(\omega)}
w^MVDR=aH(ω)R^−1a(ω)R^−1a(ω) 则 MVDR 谱估计为
P
^
M
V
D
R
(
ω
)
=
1
a
H
(
ω
)
R
^
−
1
a
(
ω
)
\hat{P}_{\mathrm{MVDR}}(\omega)=\frac1{\boldsymbol a^{\mathrm{H}}(\omega)\hat{\boldsymbol R}^{-1}\boldsymbol a(\omega)}
P^MVDR(ω)=aH(ω)R^−1a(ω)1 其中
a
(
ω
)
=
[
1
e
−
j
ω
⋯
e
−
j
ω
(
M
−
1
)
]
T
\boldsymbol{a}(\omega)=\begin{bmatrix} 1 & \mathrm{e}^{-\mathrm{j}\omega} & \cdots & \mathrm{e}^{-\mathrm{j}\omega(M-1)} \end{bmatrix}^\mathrm{T}
a(ω)=[1e−jω⋯e−jω(M−1)]T 有信号的频率处会呈现出一个峰值。
2 MATLAB仿真
\;\;\;\;\;
仿真采用的信号为复正弦信号。参数设置和仿真结果如下：
\;\;\;\;\;
从以上结果可看出，阶数为16时，20MHz和22MHz的两个信号不能被区分，但当阶数为32时，20MHz和22MHz的两个信号可以被区分。
三、MUSIC法
1 原理
\;\;\;\;\;
信号频率估计的多重信号分类(MUSIC,multiple signal classification)算法，该算法于1979年由R.O.MUSIC算法利用了信号子空间和噪声子空间的正交性，构造空间谱函数，通过谱峰搜索，估计信号频率。
\;\;\;\;\;
根据
N
N
N 个观测样本值
x
(
0
)
,
x
(
1
)
,
.
.
.
,
x
(
N
−
1
)
x(0),x(1),...,x(N-1)
x(0),x(1),...,x(N−1)，估计自相关矩
R
^
∈
C
M
×
M
\hat{\boldsymbol{R}}\in\mathbb{C}^{M\times M}
R^∈CM×M。对
R
^
\hat{\boldsymbol R}
R^ 进行特征分解，得到
M
−
K
M-K
M−K 个最小特征值对应的特征向量，即得到噪声子空间的向量，构造矩阵
G
\boldsymbol G
G。
G
=
[
u
K
+
1
u
K
+
2
⋯
u
M
]
\boldsymbol G=\begin{bmatrix}\boldsymbol u_{K+1} & \boldsymbol u_{K+2} & \cdots & \boldsymbol u_M\end{bmatrix}
G=[uK+1uK+2⋯uM] 在
[
−
π
,
π
]
[-\pi,\pi]
[−π,π]内改变
ω
\omega
ω，计算下式的值， 峰值位置就是信号频率的估计值。
P
^
M
U
S
I
C
(
ω
)
=
1
a
H
(
ω
)
G
G
H
a
(
ω
)
=
1
∑
i
=
K
+
1
M
∣
a
H
(
ω
)
u
i
∣
2
\hat{P}_{\mathrm{MUSIC}}(\omega)=\frac{1}{\boldsymbol{a}^{\mathrm{H}}(\omega)\boldsymbol{G}\boldsymbol{G}^{\mathrm{H}}\boldsymbol{a}(\omega)}=\frac{1}{\sum_{i=K+1}^{M}\mid\boldsymbol{a}^{\mathrm{H}}(\omega)\boldsymbol{u}_{i}\mid^{2}}
P^MUSIC(ω)=aH(ω)GGHa(ω)1=∑i=K+1M∣aH(ω)ui∣21
2 MATLAB仿真
\;\;\;\;\;
仿真采用的信号为复正弦信号。参数设置和仿真结果如下：
\;\;\;\;\;
从以上结果可看出，阶数为8时，20MHz和22MHz的两个信号不能被区分，但当阶数为16时，20MHz和22MHz的两个信号可以被区分。
四、3种方法的对比
\;\;\;\;\;
参数设置如下，扫描点数和FFT点数相同，仿真结果如下。
\;\;\;\;\;
从以上结果可看出，MVDR信号频率估计的分辨率最低，其次是AR参数模型，MUSIC法的谱分辨率最高。
五、MATLAB源代码
现代谱估计AR参数模型法和MVDR法和MUSIC法超详细的MATLAB代码