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

姓名： 学号： 日期：2015.12.14

1. 实验任务

信号为两个正弦信号加高斯白噪声，各正弦信号的信噪比均为10dB，长度为N，信号频率分别为f1和f2，初始相位120，取f1fs0.2，f1fs取不同的数值：0.3，0.25。fs为采样率。

（1）分别用 Levinson 递推法和 Burg 法进行功率谱估计，并分析改变数据长度、模型阶数对谱估计结果的影响。

（2）当正弦信号相位、频率、信噪比改变后，上述谱估计的结果有何变化？并作分析说明。

2. 原理分析

2.1 现代谱估计中的参数建模

根据参数模型来描述随机信号的方法，我们可以知道，如果能确定信号xn的信号模型，根据信号观测数据求出模型参数，系统函数用Hz表示，模型输入白噪声，其方差为

w2，信号的功率谱用下式求出：

PxxejwwHeiw22

按照这种求功率谱的思路，功率谱估计可分为三个步骤： （1）选择合适的信号模型；

（2）根据xn有限的观测数据，或者它的有限个自相关函数的估计值，估计模型的参数；

（3）计算墨香的输出功率谱。

其中以（1）、（2）两步最为关键。按照模型的不同，谱估计的方法有许多种，它们共同的特点是对信号观测区以外的数据不假设为0，而先根据信号观测数据估计模型参数，按照求模型输出功率的方法估计信号功率谱，回避了数据观测区以外的数据假设问题。

下面分析AR谱估计的两种方法：自相关法——列文森（Levenson）递推法和伯格（Burg）递推法。这两种方法均为已知信号观测数据，估计功率谱，两者共同特点是由信号观测数据

求模型系数时采用信号预测误差最小的原则。对于长记录数据，这些方法的估计质量是相似的，但对于短记录数据，不同方法之间存在差别。

2.2 自相关法——列文森（Levenson）递推法

自相关法的出发点是选择AR模型参数使预测误差功率最小，预测误差功率为

1Nnen21Nnxnaxnipii1p2

假设信号xn的数据区在0nN1范围，有P个预测系数，N个数据经过冲激响应为apii0,1,p的滤波器，输出预测误差en的长度为Np，因此应用下式计算：

1NNP1n0en21NNP1n0xnaxnipii1p2

上式中数据xn的两端需补充零点，相当于对无穷长的信号en的长度长于数据的长度，

加窗处理，得到长度为N的数据。上式对系数api的实部和虚部求微分使预测误差功率最小，得到

ˆxx0ˆxx1rrˆˆxx0rrxx1ˆxxp1rˆxxp2r式中自相关函数采用有偏自相关估计，即

ˆxxp1ap1ˆxx1rrˆaˆxxp2rr2p2xx （1）

aˆxx0ˆrrpppxx1N1m*xnxnmm0,1,2,,p^Nn0rxxm *^mp1,p2,,1rxxm对比上式，可知式（1）即为已推导出的Yule-Walker 方程，因此自相关法也是基于解Yule-Walker 方程的一种方法。但是直接解该方程，需要计算逆矩阵，不方便，因此，基于Yule-Walker 方程中自相关矩阵的性质，导出Levinson-Durbin递推法，这是一种高效的解方程的方法。Levinson-Durbin算法首先由一阶AR模型开始：

一阶AR模型p1的Yule-Walker方程为

rxx0rxx1112 r1r0xxxxa110由该方程解出

a11rxx1 rxx02121a1,1rxx0

然后令p2,3,4,，以此类推，可以得到一般递推公式如下：

2kp称为反射系数，kp1。1222，随着阶数增加，预测误差功率将减少p或不变。

由k=1开始递推，递推到k=p，依次得到各阶模型参数，

a112,12,a21,a22,2,ap1,ap2,2app,p

AR模型的各个系数及模型输入白噪声方差求出后，信号功率谱用下式计算

2pxxejwwHejw2p2w1/1apiejwi

i12这种方法计算简单，但需要预先估计出信号自相关函数，实际中只能按照信号的有限个观测数据估计自相关函数。当观测数据长度较短时，估计误差较大，会出现谱峰频率偏移和谱线（在信号谱峰附近产生虚假谱线）；如数据很长，估计自相关函数较准确，但计算量大，应适当选择数据长度。

2.3 伯格（Burg）递推法

Levinson-Durbin递推法需要由观测数据估计自相关函数，这是它的缺点。而伯格递推法则由信号观测数据直接计算AR模型参数。

伯格递推法利用Levinson-Durbin递推公式，导出前向预测误差与后向预测误差，并按照使它们最小的原则求出kp，从而实现不用估计自相关函数，直接用观测数据得出结果。

Burg递推法思想：借助格型预测误差滤波器，求前向、后向预测误差平均功率，选择kp使其最小，求出kp。之后，再利用Levinson-Durbin递推法求模型参数和输入噪声方差。

设信号xn的观测数据区间：0nN1,前向、后向预测误差功率分别用pf和

pb表示，预测误差平均功率用p表示，公式分别为

f1en NpnppN12pfpb1ben Npn0pN12p

2pfpb1前向、后向观测误差公式分别为

fepnxnapkxnk

k1ppenxnpa*pkxnpk

bpk1 ap01

上式中，信号项的自变量最大的是n,最小的是n-p，为了保证计算范围不超出给定的数据范围，在pf和pb计算公式中，选择求和范围为：pnN1 。

为求预测误差平均功率p最小时的反射系数kp，令推公式代入得

*2epf1nebp1n1N1pkp0,将前、后向预测误差的递

kpnpfp1ennpN12ebp1n12

Burg递推法求AR模型参数的递推公式总结如下： (1) rxx0⌒⌒21N1xn,0rxx0 Nn0(2)

e0fnxneb0nxn

n0,1,2,…，N1

n0,1,2,…，N1*2epf1nebp1n1N1(3) kpnpfp1ennpN12ebp1n1p1

2

(4) p1kp2(5) ap,iap1,ikpap1,pi i1,2,…，p1 (6) kpap,p

(7)

epfnepf1nkpebp1n1np1,p2,,N1enebpbp1n1k*fpp1ennp1,p2,,N1

3. 编程思想

（1）编写程序产生题目要求的信号和噪声 （2）然后分别用两种方法的递推流程进行谱估计 （3）改变题目中要求的变量参数，分析结果的变化

4. 代码

Levenson

clc; clear all;

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('两个正弦信号波形');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 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)); end

PSD=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;

Burg

clc; clear all;

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); 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).*eb(p-1,n-1); eb(p,n)=eb(p-1,n-1)+k(p).*ef(p-1,n); end end

k=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); end end %功率谱

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('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 产生信号

两个正弦信号波形20100-10-20020406080100120140加噪信号波形20100-10-20020406080100120140

5.2 Levenson递推法

1. 数据长度的影响（阶数设置为20阶）

N=32

自相关法-列文森Levenson递推法的功率谱估计45403530功率谱25201510500510152025频率(Hz)3035404550

N=

自相关法-列文森Levenson递推法的功率谱估计180160140120功率谱100806040200051015202530频率(Hz)35404550

N=128

自相关法-列文森Levenson递推法的功率谱估计600500400功率谱30020010000510152025频率(Hz)3035404550

N=1024

自相关法-列文森Levenson递推法的功率谱估计300025002000功率谱1500100050000510152025频率(Hz)3035404550

N=8192

自相关法-列文森Levenson递推法的功率谱估计600050004000功率谱30002000100000510152025频率(Hz)3035404550

N=16384

自相关法-列文森Levenson递推法的功率谱估计600050004000功率谱30002000100000510152025频率(Hz)3035404550

2. 模型阶数的影响（数据长度设置为N=1024）

预测误差功率的变化（从1阶到50阶）

预测误差功率变化趋势2520预测误差功率大小1510500510152025模型阶数3035404550

不同阶数（p1表示）功率谱的图像如下：

P1=5

自相关法-列文森Levenson递推法的功率谱估计706050功率谱4030201000510152025频率(Hz)3035404550

P1=10

自相关法-列文森Levenson递推法的功率谱估计300250200功率谱1501005000510152025频率(Hz)3035404550

P1=20

自相关法-列文森Levenson递推法的功率谱估计15001000功率谱50000510152025频率(Hz)3035404550

P1=30

自相关法-列文森Levenson递推法的功率谱估计4500400035003000功率谱25002000150010005000051015202530频率(Hz)35404550

P1=50

自相关法-列文森Levenson递推法的功率谱估计350030002500功率谱20001500100050000510152025频率(Hz)3035404550

3. 相位的影响（设置N=128 p1=20）

0

自相关法-列文森Levenson递推法的功率谱估计600500400功率谱30020010000510152025频率(Hz)3035404550

Pi/6

自相关法-列文森Levenson递推法的功率谱估计450400350300功率谱2502001501005000510152025频率(Hz)3035404550

Pi/3

自相关法-列文森Levenson递推法的功率谱估计400350300250功率谱2001501005000510152025频率(Hz)3035404550

Pi/2

自相关法-列文森Levenson递推法的功率谱估计450400350300功率谱2502001501005000510152025频率(Hz)3035404550

Pi

自相关法-列文森Levenson递推法的功率谱估计450400350300功率谱2502001501005000510152025频率(Hz)3035404550

4. 频率的影响（N=128 p1=20 初始相位为0）

Fs=100, f1=0.2*fs f2=0.25*fs

自相关法-列文森Levenson递推法的功率谱估计700600500功率谱40030020010000510152025频率(Hz)3035404550

Fs=100, f1=0.2*fs f2=0.3*fs

自相关法-列文森Levenson递推法的功率谱估计400350300250功率谱2001501005000510152025频率(Hz)3035404550

Fs=1000,f1=0.2*fs f2=0.25*fs

自相关法-列文森Levenson递推法的功率谱估计600500400功率谱3002001000050100150200250300频率(Hz)350400450500

Fs=1000,f1=0.2*fs f2=0.3*fs

自相关法-列文森Levenson递推法的功率谱估计350300250功率谱200150100500050100150200250300频率(Hz)350400450500

5. 信噪比的影响

SNR =20dB

自相关法-列文森Levenson递推法的功率谱估计140012001000功率谱80060040020000510152025频率(Hz)3035404550

SNR =10dB

自相关法-列文森Levenson递推法的功率谱估计600500400功率谱30020010000510152025频率(Hz)3035404550

SNR =6dB

自相关法-列文森Levenson递推法的功率谱估计300250200功率谱1501005000510152025频率(Hz)3035404550

SNR=2dB

自相关法-列文森Levenson递推法的功率谱估计12010080功率谱6040200051015202530频率(Hz)35404550

5.3 Burg递推法

1. 数据长度的影响

N=32

伯格（Burg）递推法的功率谱估计1000900800Power Spectral Density70060050040030020010001520253035Frequency(Hz)404550

N=

伯格（Burg）递推法的功率谱估计500450400Power Spectral Density3503002502001501005001520253035Frequency(Hz)404550

N=256

伯格（Burg）递推法的功率谱估计700600500Power Spectral Density40030020010001520253035Frequency(Hz)404550

N=1024

伯格（Burg）递推法的功率谱估计1200010000Power Spectral Density800060004000200001520253035Frequency(Hz)404550

N=8192

伯格（Burg）递推法的功率谱估计450040003500Power Spectral Density3000250020001500100050001520253035Frequency(Hz)404550

N=16384

伯格（Burg）递推法的功率谱估计60005000Power Spectral Density400030002000100001520253035Frequency(Hz)404550

2.

模型阶数的影响（N=128）

预测误差功率的变化趋势

预测误差功率的变化趋势2520预测误差功率大小1510500510152025模型阶次3035404550

不同阶数的功率谱 P1=5

伯格（Burg）递推法的功率谱估计250200Power Spectral Density1501005001520253035Frequency(Hz)404550

P1=10

伯格（Burg）递推法的功率谱估计800700600Power Spectral Density50040030020010001520253035Frequency(Hz)404550

P1=20

伯格（Burg）递推法的功率谱估计400350300Power Spectral Density2502001501005001520253035Frequency(Hz)404550

P1=50

伯格（Burg）递推法的功率谱估计250200Power Spectral Density1501005001520253035Frequency(Hz)404550

3. 相位的影响（N=256 p1=20）

0

伯格（Burg）递推法的功率谱估计300250Power Spectral Density2001501005001520253035Frequency(Hz)404550

Pi/6

伯格（Burg）递推法的功率谱估计12001000Power Spectral Density80060040020001520253035Frequency(Hz)404550

Pi/3

伯格（Burg）递推法的功率谱估计700600500Power Spectral Density40030020010001520253035Frequency(Hz)404550

Pi/2

伯格（Burg）递推法的功率谱估计600500Power Spectral Density40030020010001520253035Frequency(Hz)404550

Pi

伯格（Burg）递推法的功率谱估计700600500Power Spectral Density40030020010001520253035Frequency(Hz)404550

4. 频率的影响

Fs=100, f1=0.2*fs f2=0.25*fs

伯格（Burg）递推法的功率谱估计300250Power Spectral Density2001501005001520253035Frequency(Hz)404550

Fs=100, f1=0.2*fs f2=0.3*fs

伯格（Burg）递推法的功率谱估计350300250Power Spectral Density2001501005001520253035Frequency(Hz)404550

Fs=1000, f1=0.2*fs f2=0.25*fs

伯格（Burg）递推法的功率谱估计400350300Power Spectral Density250200150100500150200250300350Frequency(Hz)400450500

Fs=1000, f1=0.2*fs f2=0.3*fs

伯格（Burg）递推法的功率谱估计800700600Power Spectral Density5004003002001000150200250300350Frequency(Hz)400450500

分析：

5. 信噪比的影响

SNR=20dB

伯格（Burg）递推法的功率谱估计900800700Power Spectral Density60050040030020010001520253035Frequency(Hz)404550

SNR=10dB

伯格（Burg）递推法的功率谱估计600500Power Spectral Density40030020010001520253035Frequency(Hz)404550

SNR=6dB

伯格（Burg）递推法的功率谱估计450400350Power Spectral Density3002502001501005001520253035Frequency(Hz)404550

SNR=2dB

伯格（Burg）递推法的功率谱估计350300250Power Spectral Density2001501005001520253035Frequency(Hz)404550

分析：

信噪比越大，功率谱估计效果越好。信噪比较小的时候，可以看到，谱估计图像中，频率在20Hz和25Hz的地方，频谱不够尖锐，