Y-W方程的L-D递推算法和伯格递推算法的比较

**AR模型参数与自相关函数的关系Y-W方程的矩阵表示pmmnmRnamRyppny,,1,0][][)(][21例:p=3时的Y-W方程000)3()2()1(1]0[]1[]2[]3[]1[]0[]1[]2[]2[]1[]0[]1[]3[]2[]1[]0[2333aaaRRRRRRRRRRRRRRRRyyyyyyyyyyyyyyyy若已知Ry[n]，由Y-W方程解出各参数a3(1),a3(2),a3(3),2。则可由AR模型参数获得功率谱Py()的估计值。Y-W方程的L-D递推算法一阶Y-W方程的解0)1(1]0[]1[]1[]0[211aRRRRyyyy解此方程得]0[]1[)1(1yyRRa)1(]1[]0[121aRRyy))1(1](0[21aRyY-W方程的L-D递推算法二阶Y-W方程的解00)2()1(1]0[]1[]2[]1[]0[]1[]2[]1[]0[2222aaRRRRRRRRRyyyyyyyyy]1[]0[]1[]2[]0[)2(2222yyyyyRRRRRa211]1[)1(]2[yyRaR]1[]0[]2[]1[]1[]0[)1(222yyyyyyRRRRRRa)1()2()1(121aaa)2(]2[)1(]1[]0[2222aRaRRyyy2122])2(1[aY-W方程的L-D递推算法p阶Y-W方程的递推解21111][)(][)(pyppnypnpRnapRpa)1,,2,1()()()()(11pnnpapananapppp2122])(1[ppppaY-W方程的L-D递推算法(1)计算自相关函数的估计值(2)求解一阶模型参数关函数的估计值]0[]1[)1(1yyRRa))1(1](0[2121aRyL-D算法估计功率谱的步骤(4)求出功率谱估计21j2ARe)(1)(ˆpnnpnaPY-W方程的L-D递推算法(3)由递推算法求解p阶模型参数L-D算法估计功率谱的步骤21111][)(][)(pyppnypnpRnapRpa)1,,2,1()()()()(11pnnpapananapppp2122])(1[ppppa利用L-D算法估计频谱的MATLAB函数x：进行功率谱估计的输入有限长序列；p:模型的阶数NFFT：DFT的点数；Fs：绘制功率谱曲线的抽样频率，默认值为1；Pxx：功率谱估计值；f：Pxx值所对应的频率点[Pxx,f]=pyulear(x,p,NFFT,Fs)例：利用L-D算法进行谱估计一序列含有白噪声和两个频率的余弦信号，][)π4.0cos()π3.0cos(][kkkkx利用L-D算法估计该序列的功率谱。N=128;p=40;NFFT=2048;Fs=2;n=0:N-1;randn('state',0);x=cos(0.3*pi*n)+cos(0.4*pi*n)+randn(size(n));[P,f]=periodogram(x,[],NFFT,2);[Py,fy]=pyulear(x,p,NFFT,2);subplot(211);plot(f,10*log(P));grid;title('Periodogram');axis([01-6060]);subplot(212);plot(fy,10*log(Py));grid;title('Yule');axis([01-6060]);谱估计结果——p=40,N=12800.10.20.30.40.50.60.70.80.91-60-40-200204060Periodogram00.10.20.30.40.50.60.70.80.91-60-40-200204060Yule伯格(Burg)递推算法L-D算法缺点：在计算相关函数估计时，对N个观测数据以外的数据作零的假设，故谱估计误差较大。伯格(Burg)递推算法基本思想：直接从观测的数据利用线性预测器的前向和后向预测的总均方误差之和为最小的准则来估计反射系数，进而通过L-D算法的递推公式求出AR模型的优化参数。伯格(Burg)递推算法前向预测误差的递推公式2阶前向预测误差]2[)2(]1[)1(][][22f2kyakyakyke1阶后向预测误差][)1(]1[][1b1kyakykeL-D算法的递推公式)1()2()1()1(1212aaaa]2[)2(]1[)1()2()1(][][2121f2kyakyaaakyke]1[)1(]2[)2(]1[)1(][121kyakyakyaky]1[][][b12f1f2keKkekeK2=a2(2)2阶预测器的反射系数伯格(Burg)递推算法预测误差的递推公式一般地]1[][][b1f1fkeKkekepppp同理可得后向预测误差的递推公式][]1[][f1b1bkeKkekeppppKp=ap(p)为p阶预测器的反射系数。伯格(Burg)递推算法预测误差滤波器的格形结构]1[][][b1f1fkeKkekepppp][]1[][f1b1bkeKkekepppp1z1z2K1z1K][ky][0kef][1kef][2kef][kefp][kebp][2keb][1keb][0keb1K2KpKpK伯格(Burg)递推算法反射系数Kp的确定前向和后向均方预测误差的总和为}][][{2b2f1kekeEppNpkp由0}][][2][][2{1pbpbppfpfpNpkppKkekeKkekeKE可得}]1[][{]1[][221211111kekekekeKbpfpNpkbpfpNpkp伯格(Burg)递推算法步骤(1)确定初始条件][][][b0f0kykeke][121020kyNNk(2)从p=1开始迭代计算:计算AR模型参数}]1[][{]1[][22b12f11b1f11kekekekeKppNpkppNpkp递推p阶均方误差2122)1(pppK伯格(Burg)递推算法(3)递推高一阶前、后向预测误差]1[][][11keKkekebppfpfp][]1[][11keKkekefppbpbp(4)若阶数小于p，则阶数加1，回到步骤(2)进行下一次迭代，直到达到预定阶数p。(5)估计功率谱21j2ARe)(1)(ˆpnnpnaP伯格(Burg)递推算法步骤伯格(Burg)递推算法[Pxx,f]=pburg(x,p,NFFT,Fs)Burg算法估计频谱的MATLAB函数x：进行功率谱估计的输入有限长序列；p:模型的阶数NFFT：DFT的点数；Fs：绘制功率谱曲线的抽样频率，默认值为1；Pxx：功率谱估计值；f：Pxx值所对应的频率点利用Burg法进行谱估计程序N=512;NFFT=1024;Fs=2;p=40;n=0:N-1;randn('state',0);x=cos(0.3*pi*n)+cos(0.32*pi*n)+randn(size(n));[P,f]=pyulear(x,p,NFFT,2);[Pw,f2]=pburg(x,p,NFFT,2);subplot(211);plot(f,10*log(P));grid;title('L-D');axis([01-3060]);subplot(212);plot(f2,10*log(Pw));grid;title('Burg');axis([01-3060]);AR模型阶数p=50的谱估计结果00.10.20.30.40.50.60.70.80.91-200204060L-D00.10.20.30.40.50.60.70.80.91-200204060BurgAR模型阶数p=80的谱估计结果00.10.20.30.40.50.60.70.80.91-200204060L-D00.10.20.30.40.50.60.70.80.91-200204060Burg