利用协方差法估计AR模型参数

随机信号分析基础大作业

利用协方差法估计AR模型参数进而估计功率谱

严奎（学号：**********）

陈韬（学号：**********）

朱燕豪（学号：**********）

2011年01月15日

作业综述：

本作业中采用面向对象的程序设计方法，将用到的子程序封装在一个类中，防止其他函数的干扰，具有良好的信息内聚性。类中定义的有获得（0,1）分布随机数的函数uniform()，产生高斯分布随机数的函数gauss()，产生自回归滑动平均模型ARMA （p,q ）数据的函数arma()，用乔布斯基（Cholesky ）算法求解对称正定方程组的函数cholesky （），计算ARMA 模型的功率谱密度的函数psd （），用协方差方法估计AR 模型参数，进而实现功率谱估计的函数covar （）。采用的编程工具是V C++6.0以及VS2010，用MATLAB 对生成的数据进行画图。

一．题目要求

给定一段信号数据及采样率，利用现代谱估计理论编程估计信号的功率谱。

二．基本原理及方法

现代谱估计是通过观测数据估计参数模型再按照求参数模型输出功率的方法估计信号功率谱，主要是针对经典谱估计的分辨率低和方差性能不好等问题提出的，应用最广的是AR 参数模型。现代谱估计的参数模型有自回归滑动平均（ARMA ）模型、自回归（AR ）模型、滑动平均（MA ）模型，Wold 分解定理阐明了三者之间的关系：任何有限方差的ARMA 或MA 模型的平稳随机过程可以用无限阶的AR 模型表示，任何有限方差的ARMA 或MA 模型的平稳随机过程可以用无限阶的AR 模型表示。但是由于只有AR 模型参数估计是一组线性方程，而实际的物理系统往往是全极点系统，因而AR 应用最广。

我们用协方差法估计AR 模型参数，进而实现功率谱估计。若已知平稳随机序列x （n ）的AR 模型为

∑==-+p

i n w i n x i a n x 1

)

()()()(

其中a （i ）是AR 系数，w(n)是均值为零，均方差为σ的白噪声。

1. 计算协方差

∑-==---=1

,,1,0,),()(1),(N p

n xx p k j k n x j n x P N k j c

2. 用乔布斯基算法解对称正定方程组

N 阶对称正定方程组的矩阵形式为AX=B ，即

⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡-=⎥⎥

⎥⎥⎦⎤⎢⎢⎢

⎢⎣⎡⎥⎥⎥⎥⎥⎥⎦

⎤

⎢⎢⎢

⎢⎢⎢⎣⎡)0,()0,2()0,1()()2()1(),()3,()2,()1,(),3()3,3()2,3()1,3(),2()

3,2()2,2()1,2(),1()3,1()2,1()1,1(p c c c p a a a p p c p c p c p c p c c c c p c c c c p c c c c xx xx xx xx xx xx xx xx xx xx xx xx xx xx xx

xx xx xx xx

矩阵A 的乔布斯基分解

T LDL A =

这里D 是主对角元素都为正实数的对角阵，即D=diag(d1,d2,…,dn)，L 为主对角元素是1的

下三角矩阵。用乔布斯基算法解对称正定方程组的方法是，先用回代法求解方程组LY=B ，

得到Y 之后，再用回代法求解方程

Y

D L -1T =

3.计算激励白噪声的方差

∑=+=r

k xx xx k c k a c 1

2)

,0()()0,0(σ

4．用AR 模型参数的估计值，可以计算功率谱密度

2

1

2

)(1)(∑=-+

=

p i jwi

e

i a w S σ

三．算法设计与实现

1．程序流程图

采用协方差的方法进行功率谱估计。如下图所示

图1算法流程图1 2.主要模块的设计：

1.产生随机序列的函数uniform（），

采用线性同余法由种子seed产生随机数。

2.产生高斯白噪声的函数gauss（），

gauss(double mean,double sigma,long int * s)

{

int i;double x,y;

for(x=0,i=0;i<12;i++)

x+=uniform(0.0,1.0,s);

x=x-6.0;

y=mean+x*sigma;

return(y);

}

3.ARMA模型数据的生成函数为arma（）略

4.乔里斯基算法解对称正定方程组的函数cholesky（）略

5.由协方差函数covar（）求AR参数；

6.再根据AR参数求出功率谱的函数psd()略；

7.最终用MA TLAB的画图工具给出直观的功率谱图形，

四．结果分析

输入平稳随机序列x（n）的AR模型为

x

x

n

n

x

n

x

n

n

+

w

-

-

-

-

-

+

-

x=

)1

.0

)3

(n

924

(

)

(

)4

654

(

.2

.2

)

(

76

.3

(

)2

809

其中1，-2.76,3.809,-2.654,0.924为AR系数，

根据要求产生W(n)是均值为零，方差为1的白噪声。

根据均匀分布产生（0,1）分布的随机序列，再由均值和方差生成高斯白噪声如下图所示：