burg法实现功率谱估计

用burg 法实现功率谱估计
参数模型法是现代谱估计中的主要内容，AR 模型参数的求解有三种方法：自相关法、Burg 递推算法和改进协方差法。

Burg 算法不是直接估计AR 模型的参数，而是先估计反射系数Km,再利用Levinson 关系式求得AR 模型的参数。

Burg 算法采用的数据加窗方法是协方差法，不含有对已知数据段之外的数据做人为的假设。

1.其原理如下：
Burg 算法是使前向预测误差和后向预测误差均方误差之和最小来求取Km 的，它不对已知数据段之外的数据做认为假设。

计算m 阶预测误差的递推表示公式如下：
x(n)
(n)(n)(n)1)-(n (n)1)
-(n (n)(n)0f 0f 1-m m 1-b m 1-m f 1-m m e e e e ==+=+=e k e e k e b b m b m f
求取反射系数的公式如下：
}
1)]-(n [(n)]{[1)]
-(n (n)[2-2b 1-m 2f 1-m b 1-m f 1-m m e e e e +=E E k 对于平稳随机过程，可以用时间平均代替集合平均，因此上式可写成：
[][][][]{}p ,2,1,1)-(n (n)1)-(n (n)2-1-2
1-2
1-1-m
n 1-1-，⋯=+=∑∑==m N m n b m f m N b m f m m e e e e k 这样便可求得AR 模型的反射系数。

将m 阶AR 模型的反射系数和m-1阶AR 模型的系数代入到Levinson 关系式中，可以求得AR 模型其他的p-1个参数。

Levinson 关系式如下：
1-m 1,2,i i),-(m (i)(i)1-m 1-m m ，⋯=+=a k a a m
m 阶AR 模型的第m+1个参数G ，
ρm 2
G =,
其中ρm 是预测误差功率，可由递推公式
)-(12m 1-m m k ρρ= 求得。

易知为进行该式的递推，必须知道0阶AR 模型误差功率
ρ0 ρ0=[]
(0)(n)E x 2R x = 可知该式由给定序列易于求得。

完成上述过程，即最终求得了表征该随机信号的AR 模型的p+1个参数 。

然后根据
)(e j x ωS =2
j 2)H(e ωωσ 即可求得该随机信号的功率谱密度。

一.数据为：
12()1*2*()j n j n x n A e A e v n ωω=++
实现功率谱估计的matlan 源代码：
clear
Fs=500;
n=0:1/Fs:0.5;
w1=200*pi;A1=5;
w2=300*pi;A2=12;
xn=A1*(cos(w1*n)+j*sin(w1*n))+A2*(j*sin(w2*n)+cos(w2*n))+randn(size(n));
%xn=A1*exp(jw1n)+A2*exp(jw2n);
subplot(211);
plot(n,xn);
xlabel('n');
ylabel('xn');
title('xn=A1*exp(jw1n)+A2*exp(jw2n)+e(n)');
ymax_xn=max(xn)+0.2;
ymin_xn=min(xn)-0.2;
axis([0 0.5 ymin_xn ymax_xn]);
p=floor(length(xn)/3)+2;
nfft=1024;
[xpsd,f]=pburg(xn,p,nfft,Fs);
pmax=max(xpsd);
xpsd=xpsd/pmax;
xpsd=10*log10(xpsd);
subplot(212);
plot(f,xpsd);
title('Power Spectral estimate with burg');
ylabel('Power Spectral estimate(dB)');
xlabel('f(Hz)');
grid on;
ymin_psd=min(xpsd)-1;
ymax_psd=max(xpsd)+1;
axis([0 Fs/2 ymin_psd ymax_psd]);
实验结果：
（1）下图依次是阶数为N/3，
N/2
二.数据为：
=-+-+-
x n x n x n x n p ()(1)(2)......() Matlab源代码：
clear
clf'
clc
N=100;
Fs=500;
%产生x1信号即AR信号
vn=rands(1,N);
xn=zeros(1,N);
xn(1)=vn(1);xn(2)=vn(2);
a1=0.78;
a2=-0.96;
for n=3:N
xn(n)=vn(n)+a1*xn(n-1)+a2*xn(n-2) end
subplot(211);
plot(xn);
title('xn(n)=noise(n)+a1*xn(n-1)+a2*xn(n-2)'); p=floor(length(xn)/5);
nfft=1024;
[xpsd,f]=pburg(xn,p,nfft,Fs);
pmax=max(xpsd);
xpsd=xpsd/pmax;
xpsd=10*log10(xpsd);
subplot(212);
plot(f,xpsd);
title('Power Spectral estimate with burg'); ylabel('Power Spectral estimate(dB)'); xlabel('f(Hz)');
grid on;
ymin_psd=min(xpsd)-1;
ymax_psd=max(xpsd)+1;
axis([0 Fs/2 ymin_psd ymax_psd]);。