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

i 0

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

实验目的

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

Burg 算法原理

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

型等，其中AR 模型应用最多。

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

其中，P(e j 「为功率谱密度；s 2是激励白噪声的方差；a i 和b i 为模型参数。

若ARMA 模型中b i 全为0,就变成了 AR 模型，又称线性自回归模型， 其是一个全极点 模型：

P(e j )

研究表明，ARMA 模型和MA 模型均

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

要利用AR 模型进行功率谱估计，必须由

Yule - Walker 方程求得AR 模型的参数。而目

前求解Yule - Walker 方程主要有三种方法：Levinson-Durbin 递推算法、Burg 算法和协方差方 法。其中Burg 算法计算结果较为准确，且对于短的时间序列仍能得到较正确的估计，因此 应用广泛。

研究最大熵谱估计时，Levin so n 递推一直受制于反射系数 K m 的求出。而Burg 算法秉着 使前、后向预测误差平均功率最小的基本思想，不直接估计

AR 模型的参数，而是先估计反

射系数K m ,再利用Levin son 关系式求得AR 模型的参数，继而得到功率谱估计。

Burg 定义m 阶前、后向预测误差为：

AR 模型、MA

P(e j )

2

1

b i e

i 1

a i e j

a i e j

f m ( n)

m

a m (i)x( n i)

(1)

m

g m( n) a m(m i)x( n i)

i 0

由式(1 )和(2)又可得到前、后预测误差的阶数递推公式:

f m(n) f m i (n) K m

g m i(n 1)

g m( n) K m f mi( n)g m 1(n 1)⑷定义m阶前、后向预测误差平均功率为：

1N

2 2

P m - [|f m(门)g m(门)]

2n m

将阶数递推公式(3)和(4)代入(5)，并令一也0，可得

K m

N

f m 1( n)

g m 1( n 1)

n m 1

K m 1 N

-[f m 1( n)2 g m 1(n 1)2]

2 n m 1

三、Burg算法递推步骤

Burg算法的具体实现步骤:

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

2

x(n)

f°( n) g°( n) x(n)

步骤2求反射系数

N

f m 1( n)

g m 1( n 1)

n m 1 ______________________________

K m 1 N

-[f m 1(n)2 g m 1(n 1)2]

2 n m 1

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

P。

a m(i) a m』)K m a m 1(m i), i 1,...,m 1

a m(m) K m

步骤4计算预测误差功率

步骤5计算滤波器输出

f m ( n) f m 1( n) K m

g mi (n 1) g m ( n) K m f m 1( n) g m 1(n 1)

步骤6令m J m+1,并重复步骤2至步骤5,直到预测误差功率 P m 不再明显减小。 最后，再利用Levin so n 递推关系式估计 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=awg n(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

P m (1 K m 2

)P m 1