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

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

2、噪声，其方差为，信号的功率谱用下式求出：按照这种求功率谱的思路，功率谱估计可分为三个步骤：（1） 选择合适的信号模型；（2） 根据有限的观测数据，或者它的有限个自相关函数的估计值，估计模型的参数；（3） 计算墨香的输出功率谱。其中以（1）、（2）两步最为关键。按照模型的不同，谱估计的方法有许多种，它们共同的特点是对信号观测区以外的数据不假设为0，而先根据信号观测数据估计模型参数，按照求模型输出功率的方法估计信号功率谱，回避了数据观测区以外的数据假设问题。下面分析AR谱估计的两种方法：自相关法列文森（Levenson）递推法和伯格（Burg）递推法。这两种方法均为已知信号观测数据，估计功率谱，两

3、者共同特点是由信号观测数据求模型系数时采用信号预测误差最小的原则。对于长记录数据，这些方法的估计质量是相似的，但对于短记录数据，不同方法之间存在差别。2.2 自相关法列文森（Levenson）递推法自相关法的出发点是选择AR模型参数使预测误差功率最小，预测误差功率为假设信号的数据区在范围，有个预测系数，个数据经过冲激响应为的滤波器，输出预测误差的长度为，因此应用下式计算：的长度长于数据的长度，上式中数据的两端需补充零点，相当于对无穷长的信号加窗处理，得到长度为N的数据。上式对系数的实部和虚部求微分使预测误差功率最小，得到 （1）式中自相关函数采用有偏自相关估计，即对比上式，可知式（1）即为已推

4、导出的Yule-Walker 方程，因此自相关法也是基于解Yule-Walker 方程的一种方法。但是直接解该方程，需要计算逆矩阵，不方便，因此，基于Yule-Walker 方程中自相关矩阵的性质，导出Levinson-Durbin递推法，这是一种高效的解方程的方法。Levinson-Durbin算法首先由一阶AR模型开始：一阶AR模型的Yule-Walker方程为由该方程解出然后令，以此类推，可以得到一般递推公式如下：称为反射系数，。，随着阶数增加，预测误差功率将减少或不变。由k=1开始递推，递推到k=p，依次得到各阶模型参数，AR模型的各个系数及模型输入白噪声方差求出后，信号功率谱用下式计

5、算这种方法计算简单，但需要预先估计出信号自相关函数，实际中只能按照信号的有限个观测数据估计自相关函数。当观测数据长度较短时，估计误差较大，会出现谱峰频率偏移和谱线分裂（在信号谱峰附近产生虚假谱线）；如数据很长，估计自相关函数较准确，但计算量大，应适当选择数据长度。2.3 伯格（Burg）递推法Levinson-Durbin递推法需要由观测数据估计自相关函数，这是它的缺点。而伯格递推法则由信号观测数据直接计算AR模型参数。伯格递推法利用Levinson-Durbin递推公式，导出前向预测误差与后向预测误差，并按照使它们最小的原则求出，从而实现不用估计自相关函数，直接用观测数据得出结果。Burg递

6、推法思想：借助格型预测误差滤波器，求前向、后向预测误差平均功率，选择使其最小，求出。之后，再利用Levinson-Durbin递推法求模型参数和输入噪声方差。设信号的观测数据区间：,前向、后向预测误差功率分别用和表示，预测误差平均功率用表示，公式分别为前向、后向观测误差公式分别为 上式中，信号项的自变量最大的是n,最小的是n-p，为了保证计算范围不超出给定的数据范围，在和计算公式中，选择求和范围为： 。为求预测误差平均功率最小时的反射系数，令,将前、后向预测误差的递推公式代入得Burg递推法求AR模型参数的递推公式总结如下：(1) (2) (3) (4) (5) (6) (7) 3. 编程思想

7、（1） 编写程序产生题目要求的信号和噪声（2） 然后分别用两种方法的递推流程进行谱估计（3） 改变题目中要求的变量参数，分析结果的变化4. 代码Levensonclc;clear all;fs=100; %采样频率Ts=1/fs;N=27; %数据长度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)

8、;xx=x1+x2;x=w+x1+x2;figure(1)subplot(2,1,1);plot(xx);ylim(-20,20);title('两个正弦信号波形');grid on;subplot(2,1,2);plot(x);ylim(-20,20);title('加噪信号波形');grid on;%计算自相关函数R=;for m=1:N s=0; for n=1:N-(m-1) s=s+x(n)*x(n+m-1); end R(m)=s/N;end %列文森递推法a(1,1)=-R(2)/R(1); cov(1)=R(1)+a(1,1)*R(2);for

9、p=2:p1 sum=0; for k=1:p-1 sum=sum+a(p-1,k)*R(p-k+1); end a(p,p)=-(R(p+1)+sum)/cov(p-1);%计算反射系数 for k=1:p-1 a(p,k)=a(p-1,k)+a(p,p)*a(p-1,p-k); end cov(p)=(1-(abs(a(k,k)2)*cov(p-1);end %计算功率谱，功率谱以2*pi为周期，又信号为实信号，只需输出0到P1即可；W=0.01:0.01:pi;Z=0;for k=1:p1 Z=Z+(a(p1,k).*exp(-j*k*W);endPSD=1./(abs(1+Z).2);

10、 %算出功率谱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;clear all;f

11、s=100; %采样频率Ts=1/fs;N=28; %数据长度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

12、; k(1)=0; a=zeros(p1+1,p1+1);for p=2:p1+1 numerator=0; denominator=0; for n=p:N numerator= numerator+ (-2)*ef(p-1,n)*eb(p-1,n-1); denominator=denominator+(ef(p-1,n)2+(eb(p-1,n-1)2; end k(p)=numerator./(denominator+0.0001); cov(p)=(1-abs(k(p)2).*cov(p-1); a(p,p)=k(p); for n=p:N ef(p,n)=ef(p-1,n)+k(p)

13、.*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);for p=2:p1 for i=1:p-1 a(p,i)=a(p-1,i)+k(p)*a(p-1,p-i); endend%功率谱Z=0;W=1:0.01:pi;for i=1:p Z=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('

14、;Frequency(Hz)');ylabel('Power Spectral Density');title('伯格（Burg）递推法的功率谱估计');grid;% figure(2)% p=1:p1;% plot(p,cov(p),'b');% xlabel('模型阶次');% ylabel('预测误差功率大小');% title('预测误差功率的变化趋势');% grid;5. 实验结果及分析5.1 产生信号5.2 Levenson递推法1. 数据长度的影响（阶数设置为20阶）N=32

15、N=64N=128N=1024N=8192N=163842. 模型阶数的影响（数据长度设置为N=1024）预测误差功率的变化（从1阶到50阶）不同阶数（p1表示）功率谱的图像如下：P1=5P1=10P1=20P1=30P1=503. 相位的影响（设置N=128 p1=20）0Pi/6Pi/3Pi/2Pi4. 频率的影响（N=128 p1=20 初始相位为0）Fs=100, f1=0.2*fs f2=0.25*fsFs=100, f1=0.2*fs f2=0.3*fsFs=1000,f1=0.2*fs f2=0.25*fsFs=1000,f1=0.2*fs f2=0.3*fs5. 信噪比的影响SNR =20dB SNR =10dBSNR =6dBSNR=2dB5.3 Burg递推法1. 数据长度的影响N=32N=64N=256N=1024N=8192N=16384 2. 模型阶数的影响（N=128）预测误差功率的变化趋势不同阶数的功率谱P1=5P1=10P1=20P1=503. 相位的影响（N=256 p