信号处理---AR模型

a
k 1
N
0k
R xx ( l k ) R xx ( l ),
l 1, 2, , N (1)
a 01 a 02 R xx (0) R xx (1) R xx (2) R xx ( N ) a0 N 即 R xx ( l )
式中ak为自回归系数，称为AR系数；bk为滑动平均系数，称 为MA系数。
1
x(n)
p
b w(n l ) a
l l0
k
q
p
k
x(n k )
( a 0 1)
ak X ( z ) z
k 0

blW ( z ) z
l0
q
k 1
l
系统传输函数为：
H (z)
X (z) W (z)

bl z a
k 0 l 0 p k
q
l

k
B( z) A( z )
(极 点 )： A ( z ) （ 零 点 ） B( z )
a
k 0 q l 0
p
k
z
k
bz
l
1

z
2
如果除b0外其它的MA系数都等于0，即
AR（p）模型
p
可见，AR模型输出信号的自相关函数具有递推性质，即：
R xx ( m ) a k R xx ( m k ) 0 , m 0 k 1 p 2 R (0) a k R xx ( k ) , m 0 xx k 1
p
Yule-Walker方程（Y-W方程）
20
H p ( x 0， x1， ， x N ) ln p ( x 0， x1， ， x N ) dx

2 ln{(2 )

0
1
N
det[ R ( N )} p ( x 0， x1， ， x N ) dx ) dx
2 2
p(x ，x ， ，x
1
2R( N ) p(x ，x ， ，x
结论：
R xx (1) R xx (2) R xx (0) R xx (1) R xx ( N 1) R xx ( N 2) R xx (1) 0
R xx ( N 1) R xx ( N )
可以证明这种按最大熵外推自相关函数的结果与AR模型是等价的。
所以，上式实质为Yuler－walker方程，因此其同AR模
e
(z)
的串接
17
最大熵谱估计法（MESE）
Maximum Entropy Spectral Estimation） 熵的基本概念：熵代表一种不定度，最大熵代表最随 机，意味着对应的PSD最平坦。其定义为：
H p i ln p i
i
I i ln
1 pi
熵是信息量的一种量度，信息量与事件发生概率之间 有类似于反比例的关系，复合事件的信息量等于各独 立事件信息量之和。 对于随机过程，应该用联合概率密度函数来定义熵：
加权系数为-ak，那么 N x ( n ) a pk x ( n k )
k 1
N阶线性预测器：可看作用序列{x(n-N),x(n-N-1) , …… ,x(n-1)}激励一个冲击响应为-ak的线性时不变系 统的输出值。
{x(n-N),x(n-N-1) , …… ,x(n-1)}
一、谱估计的AR模型方法
数字系统的数学模型:有理分式传递函数的模 型的差分方程
w(n)
H ( z ) h(n)
x(n)
a 0 x ( t ) a1 x ( t 1) ... a p x ( t p ) b0 w ( t ) b1 w ( t 1) ... bq w ( t q )
RA 0 = R xx
H xx T
Pmin R xx (0) R A 0 = R xx (0) R xx A 0
R xx (0) R xx (1) R ... R xx ( N 1) R xx ( 1) R xx (0) ... R xx ( N 2) ... ... ... ... R xx (1 N ) R xx (2 N ) ... R xx (0)
H AR ( z )
X (z) W (z)

2
1 A( z )
1
1
1
S xx ( z )

a
k 1
p
k
z
k
A( z ) A( z
)
全极点模型
S xx ( )

A(e
2 j 2

p

2 2 j k
3
)
1 ak e
k 1
Yule-Walker方程（Y-W方程）推导 直接由模型差分方程推导，把模型的差
MATLAB里有专门实现L-D算法的函数可估计AR模 型参数： [a E]=aryule(x,wk.baidu.com),a为模型参数，E为噪声方差。
8
AR模型谱估计 功率谱密度
AR模型参数和 a1,a2，…，ap 激励源方差
将所估计的模型参数代入即可计算功率谱估计值：
ˆ S xx ( )

p
2
p
2 j k
P m in l 0 (2)
13
a
k 1
N
0k
R xx ( l k ) Pm in ,
将(1)和(2)两个关系式合并为一个式子：
Pm in， l 0 R xx ( l ) a 0 k R xx ( l k ) l 1, 2, , N k 1 0，
值的Rxx(N+1)：
R xx (0) det R ( N 1) R xx (1) R xx ( N 1) R xx (1) R xx (0) R xx ( N )
0
22

R xx ( N 1) R xx ( N ) R xx (0) 0
det R ( N 1) R xx ( N 1)
0 1
X
H
X
N
N
) dx 1及 x p ( x ) dx
N
H
1 2
ln{(2 ) det[ R ( N )]}
1 2
均值为0的高斯平稳随机过程的熵的表达式，它是R(N)的函数。
21
为了使得H取得最大值，应当使det[R(N)]取最大值。
根据外推或预测方法，求出使det[R(N+1)] 取最大
可以看出：N阶线性预测器的系数ak与AR模型中的
AR系数相等，预测误差概率最小值Pmin与AR模型 中的输入噪声方差 2 相等。所以，线性预测谱估
计与AR谱估计是等效的。
15
线性预测误差为
e(n) x(n)
a
k 1
p
pk
x(n k )
将上式进行 z 变换，得
p k E ( z ) 1 a pk z k 1
N
令
1, ak a0 k ,
N
k 0 k 1, 2, , N
有
Pm in， l 0 （ a k R xx (l k )= 0 ， l 1, 2, , N 3 ） k 1
14
将(3)写成矩阵展开形式为：
R xx (0) R (1) xx R xx ( N ) R xx ( 1) R xx (0) R xx ( N 1) ... ... ... 1 R xx [ ( N 1)] a 1 R xx (0) aN R xx ( N ) Pmin 0 0
分方程代入x(n)的自相关函数
x (n) ak x (n k ) w(n )
k 1 p
R xx ( m ) E x ( n ) x ( n m )
4
a k R xx ( m k ) 2 ....... m 0 k 1 R xx ( m ) p a k R xx ( m k ) .......... .... m 0 k 1
型等价 。
23
R xx (1) a1 R xx (0) a N R xx ( N 1) 0 R xx (2) a1 R xx (1) a N R xx ( N 2) 0 R xx ( N ) a1 R xx ( N 1) a N R xx (0) 0
1 a p ,k e
k 1
9
确定AR模型的阶的方法
• 不断增加阶数，观察预测误差功率，下降
到一定的程度； • 不断增加阶数，观察各阶模型预测误差序 列的周期图，最接近平坦时； • 发现后面的很多AR模型的参数接近于0；
10
线性预测谱估计
假设{x(n)}是一个N阶AR过程，现在时刻{x(n)}的 值 x ( n ) 可以由过去N个时刻的取样值的加权来预测，
12
将两个关系式写成矩阵展开式分别为：
R xx (0) R xx (1) R xx ( N 1) 即 R xx ( 1) R xx (0) R xx ( N 2) ... ... ... R xx (1 N ) a 01 R xx (2 N ) a 02 R xx (0) a 0 N R xx (1) R xx (2) R xx ( N )
H p ( x ) ln p ( x ) dx E ln( p ( x ))
18
对于零均值的高斯平稳随机过程则有：
p( x)
1

2
e
x / 2
2
2
p ( x1 , x 2 , , x N ) (2 )
-N /2
(det R ( N ))
1
2
1 exp X 2
5
写成单一正规矩阵的形式为：
R (0) R (1) ... R( p)
R (1) R (0) ... R ( p 1)
... ... ... ...
R ( p ) 1 2 a R ( p 1) 1 0 ... ... ... R (0) a p 0
H
R( N )
1
X
19
x1 x2 X xN
R xx (0) R xx (1) R(N ) R xx ( N )

R xx ( 1) R xx (0) R xx ( N 1) R xx ( N ) R xx ( ( N 1)) R xx (0)
以上利用了自相关函数的偶对称性。Y-W方程表明： 只要已知输出平稳随机信号的自相关函数，就能求出 AR模型中的参数{ak}，且需要的观测数据较少。
6
Yule-Walker方程的求解算法
• 算法，即Y-W方程的高效算法，运算量数 量级为p的二次方。
递推公式为：
a kk
a ki a k 1,i a kk a k 1, k i
X (z)
于是
H e ( z ) E (z) X (z) 1 a pk z
k 1 p k
预测误差滤波器
16
H e (z)
1 H (z)
1 a pk z
k 1
p
k
即预测误差滤波器是 x(n)的形成系统的逆滤波器，于是，有
e(n) (n)
AR 模型的 H(z)与其逆滤波器 H

2 k
2 R xx ( k ) a k 1, l R xx ( k l ) / k 1 l 1
(1 a kk )
2 2 k 1
k 1
,
2 0
xx ( 0 )
其中akk称为反射系数
7
给定初始值和AR模型的阶数p，可按照L-D算法流 程进行估计，流程终止规则为 m p 或 p
x (n)
-ak
11
预测误差为： e ( n ) x ( n ) x ( n ) 预测误差功率为： P E [| x ( n ) x ( n ) | 2 ]
确定系数ak的一个原则是使预测误差功率最小。根
据这一原则推导出的预测器系数-ak与x(n)的自相关
序列Rxx(m)之间的关系为：