参数法功率谱估计

参数法功率谱估计

一、信号的产生

（一）信号组成

在本实验中，需要事先产生待估计的信号，为了使实验结果较为明显，我产生了由两个不同频率的正弦信号（频率差相对较大）和加性高斯白噪声组成的信号。

（二）程序

N=1024;n=0:N-1;

xn=2*cos(2*pi*0.2*n)+ cos(2*pi*0.213*n)+randn(1,1024);

这样就产生了加有白噪声的两个正弦信号

其波形如下

0100200300400500600

-8-6

-4

-2

2

4

6

8

10

(a) 两个正弦信号与白噪声叠加的时域波形

二、参数模型法功率谱估计

（一）算法原理简介

1．参数模型法是现代谱估计的主要内容，思路如下：

① 假定所研究的过程)(n x 是由一个白噪声序列)(n 激励一个因果稳定的可逆线性系统)(z H 的输出；

② 由已知的)(n x ，或其自相关函数)(m r x 估计)(z H 的参数；

③ 由)(z H 的参数来估计)(n x 的功率谱。

2．自回归模型，简称AR 模型，它是一个全极点的模型。“自回归”的含义是：该模型现在的输出是现在的输入和过去p 个输出的加权和。此模型可以表现

为以下三式：

① ∑=+--=p k k n u k n x a n x 1

)()()(；

② ∑=-+==p k k

k z a z A z H 111)(1)(；

③ 212

1)(∑=-+=p k jwk

k jw x e a e P σ。

3．AR 模型的正则方程建立了参数k a 和)(n x 的自相关函数的关系，公式如下：

=)(m r x ∑=--p k x k k m r a 1)( 1≥m 时，=)(m r x 21)(σ+-∑=k r a p

k x k 0=m 时。

（二）两种AR 模型阶次的算法

1．Yule-Walker 算法(自相关法)

（1）算法主要思想

Yule-Walker 算法通过解Yule-Walker 方程获得AR 模型参数。从低阶开始递推，直到阶次p ，给出了在每一个阶次时的所有参数。公式如下：

① 11

11/])()()([--=-∑+--=m m k x x m m m r k m r k a k ρ；

② )()()(11k m a k k a k a m m m m -+=--；

③ ]1[21m

m m k -=-ρρ。

（2）运算简要框图

Yule-Walker 法谱估计运算简要框图

(3) 程序示例

clear all;close all;

N=512;n=0:N-1;

xn=2*cos(2*pi*0.2*n)+ cos(2*pi*0.213*n)+2*randn(1,512);

Rx=zeros(1,N+1);

%从课本上的公式来看，Rx （m ）中的m 属于（0，m ），即共有m+1个，故在这里设Rx 是一个一行，N+1列的向量

figure(1)

plot(n,xn);

title('（a ）两正弦信号加白噪声波形')

%下面用书中所讲自相关函数估计中的渐进无偏估计来估计自相关函数

for m=1:N+1;%由于在matlab中，下角标不能是0，m属于（0，m），在此只能从1到N+1

sum=0;

for n=1:(N+1-m);%同样道理，把书中公式里m换成m-1,N换成N+1，求和下限变为1

sum=sum+xn(n).*xn(n+m-1);

end

Rx(m)= sum/N;%切记，这里的Rx(1)才是自相关函数在0点的取值。Rx（m）只是一个存储数据的代号，为了跟书中公式一致，才叫Rx end

%下面估计各参数

P=50;

a=zeros(P,P);%a中有两个变量m,i，所以设a是P行P列的向量

km=zeros(1,P);%因为阶次是P，故反射系数有P个

p=zeros(1,P+1);%由于matlab中没有ρ，故用p来代替表示，ρ的范围是（0，P）共有P+1个

%下面计算AR模型参数的各个初始化值

p(1)=Rx(1);

a(1,1)=-Rx(2)/Rx(1);km(1)=a(1,1);

p(2)=Rx(1).*(1-abs(a(1,1).^2));

for m=2:P %由于m=1时的各个值在上面已经给出，故从m=2开始求

sum1=0;

for i=1:m-1

sum1=sum1+a(m-1,i).*Rx(m-i+1);

end

a(m,m)=-(Rx(m+1)+sum1)/p(m);%km(m)=a(m,m);求出km

for i=1:m-1

a(m,i)=a(m-1,i)+a(m,m)*a(m-1,m-i);

end

p(m+1)=p(m)*(1-abs(a(m,m)).^2);

end

z=[1,a(P,:)];

G=sqrt(p(P));

[H w]=freqz(G,z,512);%调用计算数字滤波器频响的函数

figure(2)

plot(w/(2*pi),10*log10(abs(H).^2));

title('自相关法')

ylabel('10log(PSD)')

title('(b) yule-walker法估计功率谱密度')

（4）结果分析