参数法功率谱估计一、信号的产生（一）信号组成在本实验中，需要事先产生待估计的信号，为了使实验结果较为明显，我产生了由两个不同频率的正弦信号（频率差相对较大）和加性高斯白噪声组成的信号。
（二）程序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-2246810(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 kk z a z A z H 111)(1)(；③ 2121)(∑=-+=p k jwkk 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 pk x k 0=m 时。
（二）两种AR 模型阶次的算法1．Yule-Walker 算法(自相关法)（1）算法主要思想Yule-Walker 算法通过解Yule-Walker 方程获得AR 模型参数。
从低阶开始递推，直到阶次p ，给出了在每一个阶次时的所有参数。
公式如下：① 1111/])()()([--=-∑+--=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[21mm 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+1sum=0;for n=1:(N+1-m);%同样道理，把书中公式里m换成m-1,N换成N+1，求和下限变为1sum=sum+xn(n).*xn(n+m-1);endRx(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-1sum1=sum1+a(m-1,i).*Rx(m-i+1);enda(m,m)=-(Rx(m+1)+sum1)/p(m);%km(m)=a(m,m);求出kmfor i=1:m-1a(m,i)=a(m-1,i)+a(m,m)*a(m-1,m-i);endp(m+1)=p(m)*(1-abs(a(m,m)).^2);endz=[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）结果分析00.050.10.150.20.250.30.350.40.450.5051015202530(b) yule-walker 法估计功率谱密度10l o g (P S D )从波形图中可以十分清晰的分辨出两个不同频率的正弦波2．Burg 法（1）算法主要思想Burg法不是直接估计AR模型的参数，而是先估计反射系数。
使用线性预测的方法来计算不同阶数下的反射系数，其同时使用前向和后向线性预测，使前向和后向预测误差的平均功率相对各阶反射系数m k最小，将反射系数代入Levinson-Durbin公式即可求解。
（2）运算简要框图X(n)（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)+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+1sum=0;for n=1:(N+1-m);%同样道理，把书中公式里m换成m-1,N换成N+1，求和下限变为1sum=sum+xn(n).*xn(n+m-1);endRx(m)= sum/N;%切记，这里的Rx(1)才是自相关函数在0点的取值。
Rx（m）只是一个存储数据的代号，为了跟书中公式一致，才叫RxendP=50;a=zeros(P,P);%a中有两个变量m,i，所以设a是P行P列的向量p=zeros(1,P+1);%由于matlab中没有ρ，故用p来代替表示，ρ的范围是（0，P）共有P+1个%下面计算AR模型参数a(1,1)=-Rx(2)/Rx(1);ef=zeros(P,N);eb=zeros(P,N);ef(1,:)=xn;eb(1,:)=xn;for m=2:P+1;km1=0;km2=0;for n=m:Nkm1=km1+ef(m-1,n).*eb(m-1,n-1);km2=km2+(ef(m-1,n)).^2+(eb(m-1,n-1)).^2;enda(m,m)=(-2)*km1./km2;for n=m:Nef(m,n)=ef(m-1,n)+a(m,m).*eb(m-1,n-1);eb(m,n)=eb(m-1,n-1)+a(m,m).*ef(m-1,n);endendp(1)=Rx(1);p(2)=Rx(1).*(1-abs(a(1,1).^2));a=a(2:P+1,2:P+1);for m=2:P %由于m=1时的各个值在上面已经给出，故从m=2开始求for i=1:m-1a(m,i)=a(m-1,i)+a(m,m)*a(m-1,m-i);endp(m)=p(m-1)*(1-abs(a(m,m)).^2);endz=[1,a(P,:)];G=sqrt(p(P));[H, w]=freqz(G,z,512);%调用计算数字滤波器频响的函数figure(2)plot(w/(2*pi),10*log10(abs(H).^2));ylabel('10log(PSD)')title('(b)Burg 法估计功率谱密度 ')（4）结果分析下面是程序运行后的结果00.050.10.150.20.250.30.350.40.450.5-10-5510152025303510l o g (P S D )(b)Burg 法估计功率谱密度从上图中我们同样可以看到信号中有两个频率分量，一个在0.2处，一个在0.213左右处，与产生的信号相一致。
（三）、两种算法的比较00.050.10.150.20.250.30.350.40.450.5051015202530(b) yule-walker 法估计功率谱密度10l o g (P S D )00.050.10.150.20.250.30.350.40.450.5-10-5510152025303510l o g (P S D )(b)Burg 法估计功率谱密度由上面两图我们可以看出Burg 算法得出的功率谱更平坦些。
（四）与古典法功率谱估计的比较512点00.10.20.30.40.50.60.70.80.910100200300S k00.10.20.30.40.50.60.70.80.91-5005010010l o g (P S D )(c) 周期图法估计功率谱密度00.10.20.30.40.50.60.70.80.910100200300S k00.10.20.30.40.50.60.70.80.91-5005010010l o g (P S D )(d) Bartlett 法估计功率谱密度 L=2在上面这两张图中我们就不能够很好的分辨两个靠的比较近的谱峰，在同样条件下参数法功率谱估计却能够实现谱峰的分辨，充分体现了参数法功率谱估计的优越性。
下面减小一下点数，再进行分析256点古典法00.10.20.30.40.50.60.70.80.910100200300400S k 00.10.20.30.40.50.60.70.80.91-5005010010l o g (P S D )(c) 周期图法估计功率谱密度256点参数法功率谱估计00.050.10.150.20.250.30.350.40.450.5051015202510l o g (P S D )(b) yule-walker 法估计功率谱密度由上面两图的对比我们可以看到，当采样的点数减少的时候，古典法的估计就十分不准确了，而参数法功率谱估计却还可以让人清晰地分辨两个挨得较近的谱峰。