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

数字信号处理实验报告姓名：学号：日期：2015.12.141.实验任务信号为两个正弦信号加高斯白噪声，各正弦信号的信噪比均为10dB，长度为N，信号频率分别为1f和2f，初始相位021，取2.01sff，sff1取不同的数值：0.3，0.25。sf为采样率。（1）分别用Levinson递推法和Burg法进行功率谱估计，并分析改变数据长度、模型阶数对谱估计结果的影响。（2）当正弦信号相位、频率、信噪比改变后，上述谱估计的结果有何变化？并作分析说明。2.原理分析2.1现代谱估计中的参数建模根据参数模型来描述随机信号的方法，我们可以知道，如果能确定信号xn的信号模型，根据信号观测数据求出模型参数，系统函数用zH表示，模型输入白噪声，其方差为2w，信号的功率谱用下式求出：22iwwjwxxeHeP按照这种求功率谱的思路，功率谱估计可分为三个步骤：（1）选择合适的信号模型；（2）根据nx有限的观测数据，或者它的有限个自相关函数的估计值，估计模型的参数；（3）计算墨香的输出功率谱。其中以（1）、（2）两步最为关键。按照模型的不同，谱估计的方法有许多种，它们共同的特点是对信号观测区以外的数据不假设为0，而先根据信号观测数据估计模型参数，按照求模型输出功率的方法估计信号功率谱，回避了数据观测区以外的数据假设问题。下面分析AR谱估计的两种方法：自相关法——列文森（Levenson）递推法和伯格（Burg）递推法。这两种方法均为已知信号观测数据，估计功率谱，两者共同特点是由信号观测数据求模型系数时采用信号预测误差最小的原则。对于长记录数据，这些方法的估计质量是相似的，但对于短记录数据，不同方法之间存在差别。2.2自相关法——列文森（Levenson）递推法自相关法的出发点是选择AR模型参数使预测误差功率最小，预测误差功率为21211npipininxanxNneN假设信号xn的数据区在01nN范围，有P个预测系数，N个数据经过冲激响应为0,1,piaip的滤波器，输出预测误差en的长度为Np，因此应用下式计算：210121011PNnpipiPNninxanxNneNen的长度长于数据的长度，上式中数据xn的两端需补充零点，相当于对无穷长的信号加窗处理，得到长度为N的数据。上式对系数pia的实部和虚部求微分使预测误差功率最小，得到12ˆˆˆˆ0111ˆˆˆˆ1022ˆˆˆˆ120pxxxxxxxxpxxxxxxxxppxxxxxxxxarrrprarrrprarprprrp（1）式中自相关函数采用有偏自相关估计，即1,,2,1,,2,1,01*^10*^ppmmrpmmnxnxNmrxxmNnxx对比上式，可知式（1）即为已推导出的Yule-Walker方程，因此自相关法也是基于解Yule-Walker方程的一种方法。但是直接解该方程，需要计算逆矩阵，不方便，因此，基于Yule-Walker方程中自相关矩阵的性质，导出Levinson-Durbin递推法，这是一种高效的解方程的方法。Levinson-Durbin算法首先由一阶AR模型开始：一阶AR模型1p的Yule-Walker方程为2111011100xxxxxxxxrrrra由该方程解出1110xxxxrar2211,110xxar然后令2,3,4,p，以此类推，可以得到一般递推公式如下：pk称为反射系数，1pk。22212p，随着阶数增加，预测误差功率将减少或不变。由k=1开始递推，递推到k=p，依次得到各阶模型参数，2221112122212,,,,,,,,pppppaaaaaaAR模型的各个系数及模型输入白噪声方差求出后，信号功率谱用下式计算222211/1pjwjwjwixxwwpiipeHeae这种方法计算简单，但需要预先估计出信号自相关函数，实际中只能按照信号的有限个观测数据估计自相关函数。当观测数据长度较短时，估计误差较大，会出现谱峰频率偏移和谱线分裂（在信号谱峰附近产生虚假谱线）；如数据很长，估计自相关函数较准确，但计算量大，应适当选择数据长度。2.3伯格（Burg）递推法Levinson-Durbin递推法需要由观测数据估计自相关函数，这是它的缺点。而伯格递推法则由信号观测数据直接计算AR模型参数。伯格递推法利用Levinson-Durbin递推公式，导出前向预测误差与后向预测误差，并按照使它们最小的原则求出pk，从而实现不用估计自相关函数，直接用观测数据得出结果。Burg递推法思想：借助格型预测误差滤波器，求前向、后向预测误差平均功率，选择pk使其最小，求出pk。之后，再利用Levinson-Durbin递推法求模型参数和输入噪声方差。设信号xn的观测数据区间：01nN,前向、后向预测误差功率分别用pf和pb表示，预测误差平均功率用p表示，公式分别为211NfpfpnpenNp2101NbpbpnenNp12ppfpb前向、后向观测误差公式分别为1pfppkknxnxnkaepkpkbpkpnxapnxne1*01pa上式中，信号项的自变量最大的是n,最小的是n-p，为了保证计算范围不超出给定的数据范围，在pf和pb计算公式中，选择求和范围为：1pnN。为求预测误差平均功率p最小时的反射系数pk，令0ppk,将前、后向预测误差的递推公式代入得121211*11112NpnbpfpNpnbpfppnenenenekBurg递推法求AR模型参数的递推公式总结如下：(1)12010,00NxxxxnrxnrN⌒⌒(2)00fbnxnenxne0,1,2,10,1,2,1nNnN…，…，(3)121211*11112NpnbpfpNpnbpfppnenenenek(4)121pppk(5),1,1,ppipippiaaka1,2,1ip…，(6),pppka(7)1,,2,111,,2,111*111NppnnekneneNppnneknenefppbpbpbppfpfp3.编程思想（1）编写程序产生题目要求的信号和噪声（2）然后分别用两种方法的递推流程进行谱估计（3）改变题目中要求的变量参数，分析结果的变化4.代码Levensonclc;clearall;fs=100;%采样频率Ts=1/fs;N=2^7;%数据长度p1=20;%阶数f1=0.2*fs;f2=0.25*fs;%设置信号频率pha1=0;pha2=0;%初始相位SNR=2;%设置信噪比%产生信号w=randn(1,N);Am=sqrt(2*10^(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('两个正弦信号波形');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%列文森递推法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-1a(p,k)=a(p-1,k)+a(p,p)*a(p-1,p-k);endcov(p)=(1-(abs(a(k,k)))^2)*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)).^2);%算出功率谱F=W*fs/(pi*2);%角频率坐标换算成频率坐标figure(2)plot(F,abs(PSD));xlabel('频率(Hz)');ylabel('功率谱');title('自相关法-列文森Levenson递推法的功率谱估计');grid;figure(3)p=1:p1;plot(p,cov(p),'b');xlabel('模型阶数');ylabel('预测误差功率大小');title('预测误差功率变化趋势');grid;Burgclc;clearall;fs=100;%采样频率Ts=1/fs;N=2^8;%数据长度p1=20;%阶数f1=0.2*fs;f2=0.25*fs;%设置信号频率pha1=0;pha2=0;SNR=2;%设置信噪比%产生信号w=randn(1,N);%噪声为高斯白噪声，功率为1.Am=sqrt(2*10^(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=denominator+(ef(p-1,n))^2+(eb(p-1,n-1))^2;endk(p)=numerator./(denominator+0.0001);cov(p)=(1-abs(k(p))^2).*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=1:0.01:pi;fori=1:pZ=Z+(a(p,i)*exp(-j*W*i));end;pxx=cov(p1+1)./(abs(1+Z).^2);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