功率谱估计Levinson递推法和Burg法

1、数字信号处理实验报告姓名：学号：日期：2015.12.141 .实验任务信号为两个正弦信号加高斯白噪声，各正弦信号的信噪比均为10dB,长度为N，信号频率分别为先和f2,初始相位%=52=。，取fjfs=0.2,fjfs取不同的数值：0.3,0.25。fs为采样率。(1)分别用Levinson递推法和Burg法进行功率谱估计，并分析改变数据长度、模型阶数对谱估计结果的影响。(2)当正弦信号相位、频率、信噪比改变后，上述谱估计的结果有何变化？并作分析说明。2 .原理分析2.1 现代谱估计中的参数建模根据参数模型来描述随机信号的方法，我们可以知道，如果能确定信号x(n)的信号模型，根据信号观测数据

2、求出模型参数，系统函数用H(z底示，模型输入白噪声，其方差为2仃w，信号的功率谱用下式求出：Pxxfejw)=Ow2He按照这种求功率谱的思路，功率谱估计可分为三个步骤：(1)选择合适的信号模型；(2)根据x(n府限的观测数据，或者它的有限个自相关函数的估计值，估计模型的参数；(3)计算墨香的输出功率谱。其中以(1)、(2)两步最为关键。按照模型的不同，谱估计的方法有许多种，它们共同的特点是对信号观测区以外的数据不假设为0,而先根据信号观测数据估计模型参数，按照求模型输出功率的方法估计信号功率谱，回避了数据观测区以外的数据假设问题。下面分析AR谱估计的两种方法：自相关法列文森(Levenson

3、)递推法和伯格(Burg)递推法。这两种方法均为已知信号观测数据，估计功率谱，两者共同特点是由信号观测数据求模型系数时采用信号预测误差最小的原则。对于长记录数据，这些方法的估计质量是相似的，但对于短记录数据，不同方法之间存在差别。2.2 自相关法列文森(Levenson递推法自相关法的出发点是选择AR模型参数使预测误差功率最小，预测误差功率为一2一18,1gee(n|=一£NNnRp2xnJapiXn-i=1假设彳t号x(n)的数据区在0MnMN1范围，有P个预测系数，N个数据经过冲激响应NZn=0pXn-二apiXn-ii1为api(i=0,1，p)的滤波器，输出预测误差e(n)的

4、长度为N+p,因此应用下式计算:NF'%ennz0e(n)的长度长于数据的长度，上式中数据x(n)的两端需补充零点，相当于对无穷长的信号加窗处理，得到长度为N的数据。上式对系数api的实部和虚部求微分使预测误差功率最小,p得到一朦(0)?x(-1)I朦(0)mm,0(p-1)l?x(p2)i?x(-p+1frap1&p+2)ap2iW)1app_一？(1)1"(2)rXx(p)(1)式中自相关函数采用有偏自相关估计，即i1N-4-ma-Zx*(nX(n+m)m=0,1,2,，pixx(m)=Nn田a*ixx(-m)m=_p+1,_p+2J',_1对比上式，可知

5、式(1)即为已推导出的Yule-Walker方程，因此自相关法也是基于解Yule-Walker方程的一种方法。但是直接解该方程，需要计算逆矩阵，不方便，因此，基于Yule-Walker方程中自相关矩阵的性质，导出Levinson-Durbin递推法，这是一种高效的解方程的方法。Levinson-Durbin算法首先由一阶AR模型开始：一阶AR模型(p=1)的Yule-Walker方程为rxx01rxx1rxx由该方程解出rxx1a)i=-rxx02/2-1-1_a1,1Irxx0.1然后令p=2,3,4,以此类推，可以得到一般递推公式如下：r“二一<%*=总卜口41左=L2.，夕7cr-

6、=(l-)cr:F、F*p-1口;二:(。)=修了5)222kp称为反射系数，kp<1o。1皂。2之之bp,随着阶数增加，预测误差功率将减少或不变。由k=1开始递推，递推到k=p,依次得到各阶模型参数，',a11,01J,1a21,a22,二2J,1ap1,ap2,app,二pJAR模型的各个系数及模型输入白噪声方差求出后，信号功率谱用下式计算2(pfpxx(ejw)=。：H(ejwJ=。"/1+£ape，1yJ这种方法计算简单，但需要预先估计出信号自相关函数，实际中只能按照信号的有限个观测数据估计自相关函数。当观测数据长度较短时，估计误差较大，会出现谱峰频率

7、偏移和谱线分裂(在信号谱峰附近产生虚假谱线)；如数据很长，估计自相关函数较准确，但计算量大，应适当选择数据长度。2.3伯格(Burg)递推法Levinson-Durbin递推法需要由观测数据估计自相关函数，这是它的缺点。而伯格递推法则由信号观测数据直接计算AR模型参数。伯格递推法利用Levinson-Durbin递推公式，导出前向预测误差与后向预测误差，并按照使它们最小的原则求出k,从而实现不用估计自相关函数，直接用观测数据得出结果。pBurg递推法思想：借助格型预测误差滤波器，求前向、后向预测误差平均功率，选才ikpp使其最小，求出kp。之后，再利用Levinson-Durbin递推法求模型

8、参数和输入噪声方差。设信号x(n)的观测数据区间：0wnwN1,前向、后向预测误差功率分别用Ppf和Ppb表示，预测误差平均功率用Pp表示，公式分别为np前向、后向观测误差公式分别为epnp十ppfpbfJ.epn)=xn工apkxn-kk1pebpn=xn-pVa*pkxn-pkk1ap0=1上式中，信号项的自变量最大的是n,最小的是n-p,为了保证计算范围不超出给定的数据范围，在Ppf和Ppb计算公式中，选择求和范围为：p<n<N-1odP为求预测误差平均功率Pp最小时的反射系数kp,令一p=0,将前、后向预测误差的递二kp推公式代入得kp=N4zn=pN1-2、efnepn-

9、1n-pBurg递推法求AR莫型参数的递推公式总结如下:crxx1N2小(0户工x(n),P(0尸rxx(0)n=0(2)fn)=x(n)n=0,1,2,，N-16b(n)=x(n)n=0,1,2,，N-1N_1-2%e；ne；*n-1Pp=1_kp2限(5) ap,i=a；,ipapAp-Li=1,2,p-1(6) kp=ap,pepn=epnkpe"n-1n=p1,p2,N-1epn)=e"n一1kjepnn=p1,p2,N一13.编程思想(1)编写程序产生题目要求的信号和噪声(2)然后分别用两种方法的递推流程进行谱估计(3)改变题目中要求的变量参数，分析结果的变化4.

10、代码Levensonclc;clearall;fs=100;%采样频率Ts=1/fs;N=2A7;%数据长度p1=20;%数f1=0.2*fs;f2=0.25*fs;%设置信号频率pha1=0;pha2=0;%J始相位SNR=2;%设置信噪比%产生信号w=randn(1,N);Am=sqrt(2*10A(SNR/10);k=1:N;x1=Am*sin(2*pi*f1*k*Ts+pha1);x2=Am*sin(2*pi*f2*k*Ts+pha2);xx=x1+x2;x=w+x1+x2;figure(1)subplot(2,1,1);plot(xx);ylim(-20,20);title('

11、;两个正弦信号波形);gridon;subplot(2,1,2);plot(x);ylim(-20,20);title('加噪信号波形);gridon;%计算自相关函数R=;form=1:Ns=0;forn=1:N-(m-1)s=s+x(n)*x(n+m-1);endR(m)=s/N;end发U文森递推法a(1,1)=-R(2)/R(1);cov(1)=R(1)+a(1,1)*R(2);forp=2:p1sum=0;fork=1:p-1sum=sum+a(p-1,k)*R(p-k+1);enda(p,p)=-(R(p+1)+sum)/cov(p-1);%#算反射系数fork=1:p-1

12、a(p,k)=a(p-1,k)+a(p,p)*a(p-1,p-k);endcov(p)=(1-(abs(a(k,k)A2)*cov(p-1);end%计算功率谱，功率谱以2*pi为周期，又信号为实信号，只需输出0到P1即可；W=0.01:0.01:pi;Z=0;fork=1:p1Z=Z+(a(p1,k).*exp(-j*k*W);endPSD=1./(abs(1+Z)A2);%?出功率谱F=W*fs/(pi*2);%角频率坐标换算成频率坐标figure(2)plot(F,abs(PSD);xlabel('频率(Hz)');ylabel('功率谱);title('

13、自相关法-列文森Levenson递推法的功率谱估计);grid;figure(3)p=1:p1;plot(p,cov(p),'b');xlabel('模型阶数');ylabel('预测误差功率大小);title('预测误差功率变化趋势');grid;Burgclc;clearall;fs=100;%采样频率Ts=1/fs;N=2A8;煨据长度p1=20;%数f1=0.2*fs;f2=0.25*fs;%设置信号频率pha1=0;pha2=0;SNR=2;婕置信噪比%产生信号w=randn(1,N);%噪声为高斯白噪声，功率为1.Am=sqr

14、t(2*10A(SNR/10);k=1:N;x1=Am*sin(2*pi*f1*k*Ts+pha1);x2=Am*sin(2*pi*f2*k*Ts+pha2);x=w+x1+x2;%递推ef=zeros(p1,N);eb=zeros(p1,N);ef(1,：)=x;eb(1,:)=x;cov(1)=x*x'/N;k(1)=0;a=zeros(p1+1,p1+1);forp=2:p1+1numerator=0;denominator=0;forn=p:Nnumerator=numerator+(-2)*ef(p-1,n)*eb(p-1,n-1);denominator=denominat

15、or+(ef(p-1,n)A2+(eb(p-1,n-1)A2;endk(p)=numerator./(denominator+0.0001);cov(p)=(1-abs(k(p)A2).*cov(p-1);a(p,p)=k(p);forn=p:Nef(p,n)=ef(p-1,n)+k(p).*eb(p-1,n-1);eb(p,n)=eb(p-1,n-1)+k(p).*ef(p-1,n);endendk=k(1,2:p1+1);a=a(2:p1+1,2:p1+1);forp=2:p1fori=1:p-1a(p,i)=a(p-1,i)+k(p)*a(p-1,p-i);endend%功率谱Z=0;W

16、=1:0.01:pi;fori=1:pZ=Z+(a(p,i)*exp(-j*W*i);end;pxx=cov(p1+1)./(abs(1+Z).A2);F=W*fs/(pi*2);%角频率坐标换算成频率坐标figure(1)plot(F,pxx);xlabel('Frequency(Hz),);ylabel('PowerSpectralDensity');title('伯格(Burg)递推法的功率谱估计);grid;%figure(2)%p=1:p1;%plot(p,cov(p),'b');%xlabel('模型阶次');%yla

17、bel('预测误差功率大小');%title('预测误差功率的变化趋势');%grid;5.实验结果及分析5.1产生信号两个正弦信号波形加噪信号波形5.2Levensond的推法1.数据长度的影响（阶数设置为20阶）N=32N=64160140120港100¥功8060402000111-j1,5101520253035404550频率（Hz）N=128600自相关法-列文森Levenson递推法的功率谱估计500400-rife组率300功2001000011J115101520253035404550频率（Hz）N=1024300025002000

18、率1500功100050000自相关法-列文森Levenson递推法的功率谱估计,jL5101520253035404550频率（Hz）N=8192600050004000埴率3000功2000100025频率（Hz）N=163846000自相关法-列文森Levenson递推法的功率谱估计5000400030002000100010051015202530频率（Hz）354045502,模型阶数的影响（数据长度设置为N=1024）预测误差功率的变化（从1阶到50阶）预测误差功率变化趋势20太15I105005101520253035404550模型阶数不同阶数（pl表示）功率谱的图像如下:P1

19、=570自相关法-列文森Levenson递推法的功率谱估计6050JI、W4011t30201000J1510152025频率30(Hz)35404550P1=10300250自相关法-列文森Levenson递推法的功率谱估计200缶150功100I00510152025频率3035404550(Hz)P1=201000用举功500005101520253035404550频率（Hz）P1=30自相关法-列文森Levenson递推法的功率谱估计频率（Hz）P1=50频率（Hz）3.相位的影响（设置N=128p1=20）0自相关法-列文森Levenson递推法的功率谱估计600500400遭率3

20、00面2001000:J1105101520253035404550频率(Hz)Pi/6自相关法-列文森Levenson递推法的功率谱估计450400350300加250至功20015010050005101520253035404550频率(Hz)Pi/31/450400350300250200150100500自相关法-列文森Levenson递推法的功率谱估计05101520253035404550频率(Hz)Pi/23503002502001501005010015101520253035404550400自相关法-列文森Levenson递推法的功率谱估计频率(Hz)Pi频率（Hz）4.

21、频率的影响(N=128p1=20初始相位为0)Fs=100,f1=0.2*fsf2=0.25*fs频率(Hz)Fs=100,f1=0.2*fsf2=0.3*fs400350300250200150100501105101520253035404550自相关法-列文森Levenson递推法的功率谱估计频率(Hz)Fs=1000,f1=0.2*fsf2=0.25*fs4002001J50100150200250300350400450500频率（Hz）50010000Fs=1000,f1=0.2*fsf2=0.3*fs自相关法-列文森Levenson递推法的功率谱估计频率（Hz）5.信噪比的影响S

22、NR=20dB12001000逆800¥功八60040020000111JI5101520253035404550频率(Hz)SNR=10dB600500自相关法-列文森Levenson递推法的功率谱估计400埴率300面200100J1/11051015202530频率(Hz)35404550SNR=6dB300自相关法-列文森Levenson递推法的功率谱估计250200霆150场100501111I/10510152025频率(Hz)3035404550SNR=2dB80604020111115101520253035404550频率(Hz)5.3Burg递推法20010001

23、5I/20253035404550Frequency(Hz)N=64500伯格(Burg)递推法的功率谱估计450400y350D300De250W2001501100500151J20253035404550Frequency(Hz)N=2567000151111600500400300200100202530354045伯格（Burg）递推法的功率谱估计50Frequency(Hz)N=102412000100008000600040002000015伯格（Burg）递推法的功率谱估计2025403035Frequency(Hz)4550N=819245004000350030002500

24、200015001000500015202540伯格（Burg）递推法的功率谱估计453035Frequency(Hz)50N=1638460005000400030002000100001530Frequency(Hz)2 .模型阶数白影响（N=128）模型阶次不同阶数的功率谱015/i1ii1/f20253035404550Frequency(Hz)10QO521P1=10计估谱率功的法推递UB格伯101010'0'0'0'0ooQQQoO876543201512025404550P1=20伯格(Burg)递推法的功率谱估计400350300250200150100503035Frequency(Hz)15010050Frequency(Hz)3 .相位的影响(N=256p1=20)0Pi/6Frequency(Hz)5O10000008642yrncouaJpHJBO-Frequency(Hz)5o1oooooooooooooo7654321jjy/mDnuapJgasL/3piFrequency(Hz)5o