AR模型功率谱估计及Matlab实现

南昌大学实验报告

学生姓名：学号：专业班级：

实验类型：□验证□综合□设计□创新实验日期：实验成绩：

一、实验名称

基于AR模型的功率谱估计及Matlab实现

二、实验目的

1.了解现代谱估计方法，深入研究AR模型法的功率谱估计

2.利用Matlab对AR模型法进行仿真

三、实验原理

1.现代谱估计

现代功率谱估计以信号模型为基础，如下图所示为x(n)的信号模型，输入白噪声ω(n)均值为0，方差为σω2，x(n)的功率谱可由下式计算：

P xx(e jω)=σω2|H(e jω)|2

如果通过观测数据估计出信号模型的参数，信号功率谱就可以按上式计算出来，这样估计功率谱的问题就变成由观测数据估计信号模型参数的问题。

2.功率谱估计的步骤：

（1）选择合适的信号模型；

（2）根据x(n)有限的观测数据，或者有限个自相关函数估计值，估计模型的参数；

（3）计算模型的输出功率谱。

3.模型选择

选择模型主要考虑是模型能够表示谱峰、谱谷和滚降的能力。对于尖峰的谱，选用具有极点的模型，如AR、ARMA模型；对于具有平坦的谱峰和深谷的信号，可以选用MA模型；既有极点又有零点的谱应选用ARMA模型，应该在选择模型合适的基础上，尽量减少模型的参数。

4.AR模型功率谱估计

在实际中，AR 模型的参数估计比较简单，对其有充分的研究，AR模型功率谱估计又称为自回归模型，它是一个全极点的模型，要利用AR模型进行功率谱估可以通过列文森(Levenson)递推算法由Yule-Walker 方程求AR模型的参数。

4.MATLAB中AR模型的谱估计的函数说明：

1．Pyulear函数：

功能：利用Yule--Walker方法进行功率谱估计.

格式：Pxx=Pyulear(x,ORDER,NFFT)

[Pxx,W]=Pyulear(x,ORDER,NFFT)

[Pxx,W]=Pyulear(x,ORDER,NFFT,Fs)

Pyulear(x,ORDER,NFFT,Fs,RANGE,MAGUNITS)

说明：Pxx =Pyulear(x,ORDER,NFFT)中，采用Yule--Walker方法估计序列x的功率谱，参数ORDER用来指定AR模型的阶数，NFFT为FFT算法的长度，默认值为256，若NFFT为偶数，则Pxx为（NFFT/2 + 1）维的列矢量，若NFFT为奇数，则Pxx为（NFFT + 1）/2维的列矢量；当x为复数时，Pxx长度为NFFT。[Pxx,W]=Pyulear(x,ORDER,NFFT)中，返回一个频率向量W.

[Pxx,W]=Pyulear(x,ORDER,NFFT,Fs)中，可以在F向量得到功率谱估计的频率点，Fs指定采样频率。

四、实验任务及要求

实验任务是利用Matlab采用Levenson算法完成对AR模型参数的求解，分析在不同模型阶数的下，功率谱估计的差异。

五、实验内容

输入一高斯白噪声干扰下的信号x(n)=cos(100πn)+sin(200πn)+w(n),其中w(n)为均值为0，方差为1的高斯白噪声。用自相关法求得AR模型参数的估计值并算得功率谱估计。并比较在不同阶次下功率谱估计的差别。

六、实验程序及波形图

%程序

fs=2000;

n=0:1/fs:0.4;

w=randn(size(n));

f1=100;f2=200;

x=cos(2*pi*f1*n)+sin(2*pi*f2*n)+w;%输入信号

figure(1)

plot(n,x);axis([0 0.4 -6 6]);grid;

xlabel('时间');ylabel('幅度');

title('输入信号');

N=3072;%取长度为3072

%使用自相关法得到功率谱估计；

[xpsd,f]=pyulear(x,10,N);%阶次P为10

pmax=max(xpsd);

xpsd=xpsd/pmax;

xpsd=10*log10(xpsd+0.000001);

figure(2)

plot(f,fftshift(xpsd));%将FFT的直流分量移到频谱中心

axis([0 3 -15 5]);

xlabel('频率');ylabel('功率谱密度');grid;

title('AR谱估计曲线（p=10）');

[xpsd,f]=pyulear(x,20,N);%阶次P为20

pmax=max(xpsd);

xpsd=xpsd/pmax;

xpsd=10*log10(xpsd+0.000001); figure(3)

plot(f,fftshift(xpsd));axis([0 3 -30 10]); xlabel('频率');ylabel('功率谱密度');grid; title('AR谱估计曲线（p=20）'); [xpsd,f]=pyulear(x,30,N);%阶次P为30 pmax=max(xpsd);

xpsd=xpsd/pmax;

xpsd=10*log10(xpsd+0.000001); figure(4)

plot(f,fftshift(xpsd));axis([0 3 -30 10]); xlabel('频率');ylabel('功率谱密度');grid; title('AR谱估计曲线（p=30）');