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模型参数的各个初始化值
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列的向量
p=zeros(1,P+1);%由于matlab中没有ρ，故用p来代替表示，ρ的范围是（0，P）共有P+1个
① ；
② ；
③ 。
3．AR模型的正则方程建立了参数 和 的自相关函数的关系，公式如下：
时， 时。
（二）两种AR模型阶次的算法
1．Yule-Walker算法(自相关法)
（1）算法主要思想
Yule-Walker算法通过解Yule-Walker方程获得AR模型参数。从低阶开始递推，直到阶次p，给出了在每一个阶次时的所有参数。公式如下：
参数法功率谱估计
一、信号的产生
（一）信号组成
在本实验中，需要事先产生待估计的信号，为了使实验结果较为明显，我产生了由两个不同频率的正弦信号（频率差相对较大）和加性高斯白噪声组成的信号。
（二）程序
N=1024;n=0:N-1;
xn=2*cos(2*pi*0.2*n)+ cos(2*pi*0.213*n)+randn(1,1024);
end
a(m,m)=(-2)*km1./km2;
for n=m:N
ef(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);
end
end
p(1)=Rx(1);
p(2)=Rx(1).*(1-abs(a(1,1).^2));
从波形图中可以十分清晰的分辨出两个不同频率的正弦波
2．Burg法
（1）算法主要思想
Burg法不是直接估计AR模型的参数，而是先估计反射系数。
使用线性预测的方法来计算不同阶数下的反射系数，其同时使用前向和后向线性预测，使前向和后向预测误差的平均功率相对各阶反射系数 最小，将反射系数代入Levinson-Durbin公式即可求解。
① ；
② ；
③ 。
（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
%下面计算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:N
km1=km1+ef(m-1,n).*eb(m-1,n-1);
km2=km2+(ef(m-1,n)).^2+(eb(m-1,n-1)).^2;
a=a(2:P+1,2:P+1);
for m=2:P %由于m=1时的各个值在上面已经给出，故从m=2开始求
for i=1:m-1
a(m,i)=a(m-1,i)+a(m,m)*a(m-1,m-i);
end
p(m)=p(m-1)*(1-abs(a(m,m)).^2);
end
z=[1,a(P,:)];
G=sqrt(p(P));
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));
（三）、两种算法的比较
由上面两图我们可以看出Burg算法得出的功率谱更平坦些。
（四）与古典法功率谱估计的比较
512点
在上面这两张图中我们就不能够很好的分辨两个靠的比较近的谱峰，在同样条件下参数法功率谱估计却能够实现谱峰的分辨，充分体现了参数法功率谱估计的优越性。
下面减小一下点数，再进行分析
256点古典法
[H, w]=freqz(G,z,512);%调用计算数字滤波器频响的函数
figure(2)
plot(w/(2*pi),10*log10(abs(H).^2));
ylabel('10log(PSD)')
title('(b)Burg法估计功率谱密度')
（4）结果分析
下面是程序运行后的结果
从上图中我们同样可以看到信号中有两个频率分量，一个在0.2处，一个在0.213左右处，与产生的信号相一致。
这样就产生了加有白噪声的两个正弦信号
其波形如下
二、参数模型法功率谱估计
（一）算法原理简介
1．参数模型法是现代谱估计的主要内容，思路如下：
①假定所研究的过程 是由一个白噪声序列 激励一个因果稳定的可逆线性系统 的输出；
②由已知的 ，或其自相关函数 估计 的参数；
③由 的参数来估计 的功率谱。
2．自回归模型，简称AR模型，它是一个全极点的模型。“自回归”的含义是：该模型现在的输出是现在的输入和过去p个输出的加权和。此模型可以表现为以下三式：
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);
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
[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）结果分析
（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列的向量
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
256点参数法功率谱估计
由上面两图的对比我们可以看到，参数法功率谱估计却还可以让人清晰地分辨两个挨得较近的谱峰。优点十分突出。
一些体会
在持续了半月的做作业期间内，我自己可以很清楚的看到自己的进步，从什么都不太懂到现在通过自己的努力可以编出相应的程序，虽然我知道自己编的程序还不是很完美，有很多地方都可以改得更简练，但这半月每天晚上查资料，调程序到1、2点的日子还是很让人怀念。大作业告一段落了，但我学习的脚步不会停下来，大作业恰恰给了我一个学习进步的机会，让我知道今后的路该怎么走。我想在以后的日子里，我还会继续修改这些程序，使它们更加完善。也让自己在编程中没有彻底弄明白的东西消化掉。
人生的路，且学且行。