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)0f0f1-mm1-bm1-mf1-mmeeeeekeekebbmbmf求取反射系数的公式如下：}1)]-(n[(n)]{[1)]-(n(n)[2-2b1-m2f1-mb1-mf1-mmeeeeEEk对于平稳随机过程，可以用时间平均代替集合平均，因此上式可写成：p,2,1,1)-(n(n)1)-(n(n)2-1-21-21-1-mn1-1-，mNmnbmfmNbmfmmeeeek这样便可求得AR模型的反射系数。将m阶AR模型的反射系数和m-1阶AR模型的系数代入到Levinson关系式中，可以求得AR模型其他的p-1个参数。Levinson关系式如下：1-m1,2,ii),-(m(i)(i)1-m1-mm，akaamm阶AR模型的第m+1个参数G，m2G,其中m是预测误差功率，可由递推公式)-(12m1-mmk求得。易知为进行该式的递推，必须知道0阶AR模型误差功率00=(0)(n)Ex2Rx可知该式由给定序列易于求得。完成上述过程，即最终求得了表征该随机信号的AR模型的p+1个参数。然后根据)(ejxS=2j2)H(e即可求得该随机信号的功率谱密度。一.数据为：12()1*2*()jnjnxnAeAevn实现功率谱估计的matlan源代码：clearFs=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([00.5ymin_xnymax_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('PowerSpectralestimatewithburg');ylabel('PowerSpectralestimate(dB)');xlabel('f(Hz)');gridon;ymin_psd=min(xpsd)-1;ymax_psd=max(xpsd)+1;axis([0Fs/2ymin_psdymax_psd]);实验结果：（1）下图依次是阶数为N/3，N/2二.数据为：()(1)(2)......()xnxnxnxnpMatlab源代码：clearclf'clcN=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;forn=3:Nxn(n)=vn(n)+a1*xn(n-1)+a2*xn(n-2)endsubplot(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('PowerSpectralestimatewithburg');ylabel('PowerSpectralestimate(dB)');xlabel('f(Hz)');gridon;ymin_psd=min(xpsd)-1;ymax_psd=max(xpsd)+1;axis([0Fs/2ymin_psdymax_psd]);