基于Burg算法的最大熵谱估计

基于Burg 算法的最大熵谱估计

一、 实验目的

使用Matlab 平台实现基于Burg 算法的最大熵谱估计

二、 Burg 算法原理

现代谱估计是针对经典谱估计方差性能较差、分辨率较低的缺点提出并逐渐发展起来的，其分为参数模型谱估计和非参数模型谱估计。而参数模型谱估计主要有AR 模型、MA 模型、ARMA 模型等，其中AR 模型应用最多。

ARMA 模型功率谱的数学表达式为：

2

12121/1)(∑∑=-=-++=p i i j i q i i j i j e a e b e P ωωωσ

其中，P(e j ω

)为功率谱密度；s 2是激励白噪声的方差；a i 和b i 为模型参数。 若ARMA 模型中b i 全为0，就变成了AR 模型，又称线性自回归模型，其是一个全极点模型： 2

121/)(∑=-+=p i i j i j e a e P ωωσ

研究表明，ARMA 模型和MA 模型均可用无限阶的AR 模型来表示。且AR 模型的参数估计计算相对简单。同时，实际的物理系统通常是全极点系统。

要利用AR 模型进行功率谱估计，必须由Yule - Walker 方程求得AR 模型的参数。而目前求解Yule - Walker 方程主要有三种方法: Levinson-Durbin 递推算法、Burg 算法和协方差方法。其中Burg 算法计算结果较为准确，且对于短的时间序列仍能得到较正确的估计，因此应用广泛。

研究最大熵谱估计时，Levinson 递推一直受制于反射系数K m 的求出。而Burg 算法秉着使前、后向预测误差平均功率最小的基本思想，不直接估计AR 模型的参数，而是先估计反射系数K m ，再利用Levinson 关系式求得AR 模型的参数，继而得到功率谱估计。 Burg 定义m 阶前、后向预测误差为：

∑=-=m

i m m i n x i a n f 0)()()( (1)

∑=*--=m

i m i n x i m a n g m 0)()()( (2) 由式（1）和（2）又可得到前、后预测误差的阶数递推公式：

)1()()(11-+=--n g K n f n f m m m m (3)

)1()()(11-+=--*n g n f K n g m m m m (4)

定义m 阶前、后向预测误差平均功率为：

∑=+=N

m

n m m m n g n f P ])()([21

2

2

(5) 将阶数递推公式（3）和（4）代入（5），并令0=∂∂m

m

K P ，可得

∑∑+=--+=*---+--=N m n m m N m n m m m n g

n f n g n f K 1

2

1211

11])1()([21)

1()(

(6)

三、 Burg 算法递推步骤

Burg 算法的具体实现步骤：

步骤1 计算预测误差功率的初始值和前、后向预测误差的初始值，并令m = 1。

2

1

0)(1

∑==N n n x N P

)()()(00n x n g n f ==

步骤2 求反射系数

∑∑+=--+=*---+--

=N m n m m N m n m m m n g n f n g n f K 1

2

1211

11]

)1()([21)1()(

步骤3 计算前向预测滤波器系数

),()()(11i m a K i a i a m m m m -+=*

-- 1,...,1-=m i

m m K m a =)(

步骤4 计算预测误差功率

12

)1(--=m m m P K P

步骤5计算滤波器输出 )1()()(11-+=--n g K n f n f m m m m

)1()()(11-+=--*n g n f K n g m m m m

步骤6 令m ← m+1，并重复步骤2至步骤5，直到预测误差功率P m 不再明显减小。 最后，再利用Levinson 递推关系式估计AR 参数，继而得到功率谱估计。

四、 程序实现

%%%%%%%%%%%%基于Burg 算法的最大熵谱估计的Matlab 实现%%%%%%%%%%%% %%%设置两正弦小信号的归一化频率分别为0.175和0.20，信噪比SNR=30dB 、N=32%%% clear,clc; %清空内存及变量

N=32; %设置离散傅里叶变换点数，即最大阶数N 为32 SNR=30; %信噪比SNR 取为30dB

fs=1; %采样频率取为1Hz

t=1:N; %采样时间点从1变化到N

t=t/fs; %得到归一化频率采样点

y=sin(2*pi*0.175*t)+sin(2*pi*0.20*t); %信号归一化频率分别取为0.175和0.20 x=awgn(y,SNR); %在信号y 中加入高斯白噪声,信噪比为SNR 设定的数值 M=1; %设置起始计算的阶数M 为1

P(M)=0; %预测误差功率初值设为0

Rx(M)=0; %自相关函数初值设为0

for n=1:N %样本数从1变化到N

P(M)=P(M)+(abs(x(n)))^2; %计算预测误差功率和的初始值

ef(1,n)=x(n); %计算前向预测误差初值，令其等于此时的信号序列 eb(1,n)=x(n); %计算后向预测误差初值，令其等于此时的信号序列 end

P(M)=P(M)/N; %计算出预测误差功率的初始值

Rx(M)=P(M); %设定自相关函数初始值

M=2; %设置起始计算的阶数M 为2

A=0; %微分所得反射系数Km 的分子，初始值设为0 D=0; %微分所得反射系数Km 的分母，初始值设为0 for n=M:N %AR 阶数由M 变化到N