基于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

21211

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

212

11

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

A=A+ef(M-1,n)*eb(M-1,n-1); %计算分子的和

D=D+(abs(ef(M-1,n)))^2+(abs(eb(M-1,n-1)))^2; %计算分母的和（即M阶前、

后向预测误差平均功率) end

Km=-2*A/D; %计算反射系数Km（此时起始阶数为2）

a(M-1,M-1)=-2*A/D; %计算前向预测滤波器系数

P(M)=P(M-1)*(1-(abs(Km))^2); %计算预测误差功率

FPE(M-1)=P(M)*(N+M)/(N-M); %设置最大预测误差平均功率

TH=FPE(M-1);

for n=M:N %AR阶数由M变化到N

ef(M,n)=ef(M-1,n)+Km*eb(M-1,n-1); %计算滤波器输出的前向预测误差

eb(M,n)=eb(M-1,n-1)+Km*ef(M-1,n); %计算滤波器输出的后向预测误差end

M=M+1; %阶数叠加，以便递推计算下一阶数据

A=0; %反射系数Km的分子，初始值设为0

D=0; %反射系数Km的分母，初始值设为0

for n=M:N %同前，进行递推运算

A=A+ef(M-1,n)*eb(M-1,n-1);

D=D+(abs(ef(M-1,n)))^2+(abs(eb(M-1,n-1)))^2;

end

Km=-2*A/D;

a(M-1,M-1)=-2*A/D;

P(M)=P(M-1)*(1-(abs(Km))^2);

FPE(M-1)=P(M)*(N+M)/(N-M);

for m=1:M-2 %AR阶数m由1变化到（M-2）a(M-1,m)=a(M-2,m)+Km*a(M-2,M-1-m); %递推计算各阶前向预测滤波器的系数end

while FPE(M-1)

ef(M,n)=ef(M-1,n)+Km*eb(M-1,n-1); %递推计算滤波器输出的各阶前向预测误差

eb(M,n)=eb(M-1,n-1)+Km*ef(M-1,n); %递推计算滤波器输出的各阶后向预测误差end

Km=-2*A/D; %反射系数

a(M-1,M-1)=-2*A/D; %前向预测滤波器系数

P(M)=P(M-1)*(1-(abs(Km))^2); %预测误差功率

FPE(M-1)=P(M)*(N+M)/(N-M);

for m=1:M-2 %AR阶数m由1变化到（M-2）a(M-1,m)=a(M-2,m)+Km*a(M-2,M-1-m); %递推得到各阶前向预测滤波器的系数end

end

T=1/fs;sum1=0; %采样周期T赋值；功率谱初值设为0

f=0.01:0.01:0.5; %选取数据采样点，归一化频率0.01～0.5，间隔为0.01 for m=1:M-1; %AR阶数为1～(M-1)

sum1=sum1+a(M-1,m)*exp(-j*2*pi*m*f*T); %傅里叶变换，得到AR参数的估计end

s1=(abs(1+sum1)).^2; %由Fejer-Riesz定理，得到最大熵谱估计，即ARMA功率谱s=P(M)*T./s1; %求得各阶功率谱的矩阵

plot(f,10*log10(s),'k'); %画出功率谱随频率变化的曲线图

xlabel('f/fs'); %X轴坐标名称

ylabel('功率谱/dB'); %Y轴坐标名称

五、仿真结果及分析

结果分析：如上，是阶数分别为16和32时所得的功率谱曲线。如图可知，Burg算法得到的谱线分辨率很高，曲线平滑而且干扰很少，可清晰分辨出功率谱的峰值及其对应的频率值（处于初始设定的频率值0.175～0.20之间）；并且随着阶数的增加，预测误差功率随之减小，但变化不大，对于频率分辨率影响也很小。相较于Levinson算法，Burg算法的基本思想是使前、后向预测误差的平均功率最小，且没有使用自相关估计法，因而，结果与真实值更接近。

时间序列分析——最经典的

【时间简“识”】 说明：本文摘自于经管之家(原人大经济论坛) 作者：胖胖小龟宝。原版请到经管之家(原人大经济论坛) 查看。 1.带你看看时间序列的简史 现在前面的话—— 时间序列作为一门统计学，经济学相结合的学科，在我们论坛，特别是五区计量经济学中是热门讨论话题。本月楼主推出新的系列专题——时间简“识”，旨在对时间序列方面进行知识扫盲（扫盲，仅仅扫盲而已……），同时也想借此吸引一些专业人士能够协助讨论和帮助大家解疑答惑。 在统计学的必修课里，时间序列估计是遭吐槽的重点科目了，其理论性强，虽然应用领域十分广泛，但往往在实际操作中会遇到很多“令人发指”的问题。所以本帖就从基础开始，为大家絮叨絮叨那些关于“时间”的故事！ Long long ago，有多long估计大概7000年前吧，古埃及人把尼罗河涨落的情况逐天记录下来，这一记录也就被我们称作所谓的时间序列。记录这个河流涨落有什么意义当时的人们并不是随手一记，而是对这个时间序列进行了长期的观察。结果，他们发现尼罗河的涨落非常有规律。掌握了尼罗河泛滥的规律，这帮助了古埃及对农耕和居所有了规划，使农业迅速发展，从而创建了埃及灿烂的史前文明。

好~~从上面那个故事我们看到了 1、时间序列的定义——按照时间的顺序把随机事件变化发展的过程记录下来就构成了一个时间序列。 2、时间序列分析的定义——对时间序列进行观察、研究，找寻它变化发展的规律，预测它将来的走势就是时间序列分析。 既然有了序列，那怎么拿来分析呢 时间序列分析方法分为描述性时序分析和统计时序分析。 1、描述性时序分析——通过直观的数据比较或绘图观测，寻找序列中蕴含的发展规律，这种分析方法就称为描述性时序分析 描述性时序分析方法具有操作简单、直观有效的特点，它通常是人们进行统计时序分析的第一步。 2、统计时序分析 （1）频域分析方法 原理：假设任何一种无趋势的时间序列都可以分解成若干不同频率的周期波动 发展过程： 1）早期的频域分析方法借助富里埃分析从频率的角度揭示时间序列的规律 2）后来借助了傅里叶变换，用正弦、余弦项之和来逼近某个函数 3）20世纪60年代，引入最大熵谱估计理论，进入现代谱分析阶段 特点：非常有用的动态数据分析方法，但是由于分析方法复杂，结果抽象，有一定的使用局限性 （2）时域分析方法

基于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)

功率谱分析(Matlab)

功率谱分析 由题目内容，设采样频率fs=1000HZ，数据长度为256，模型阶数为14，f1=200,f2=300、250。 （1）用最大熵法进行谱估计 运行程序后，观察图像f1和f2相差较小时，功率谱变化更剧烈；模型的阶数越高，图像中能够获得的信息就越多，但同时计算量也就越大；增加数据长度可以获得更多的信息，提高了谱分析的分辨率，这是因为AR模型的谱估计隐含着对数据和自相关函数的外推，其长度可能会超过给定长度，分辨率不受信源信号的限制。 （2）分别用Levinson递推法和Burg法进行功率谱分析 ①Levinson递推法 运行程序后，观察图像，f1和f2相差较小时，功率谱变化更剧烈；模型的阶数越高，图像中能够获得的信息就越多，但同时计算量也就越大；增加数据长度可以获得更多的信息，提高了谱分析的分辨率，但本题中信号为正弦信号加白噪声，故图像观察不明显。 ②Burg法 运行程序后，观察图像，f1和f2相差较小时，功率谱变化更剧烈；模型的阶数越高，图像中能够获得的信息就越多，但同时计算量也就越大；增加数据长度可以获得更多的信息，提高了谱分析的分辨率。 (3)改变信号的相位、频率、信噪比，上述谱分析结果有何变化 如果正弦信号的频率过大，超过fs/2，会产生频率混叠现象，输入f1=600HZ，会在400HZ处产生一个波峰；降低信噪比会导致谱分辨率下降；信号起始相位的变动可导致谱线的偏移和分裂（我的图像观察不到）。 最大熵法估计 N=1024; Nfft=256; Fs=1000; n=0:N-1; t=n/Fs; x1=sin(2*pi*200*t); x2=sin(2*pi*300*t); %0.3 xn=x1+awgn(x1,10)+x2+awgn(x2,10); [Pxx1,f]=pmem(xn,14,Nfft,Fs); subplot(4,1,1) plot(f,10*log10(Pxx1)); xlabel('Frequency (Hz)');ylabel('Power Spectrum (dB)');

基于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)

最大熵原理及其应用

论文名称：最大熵原理及其应用班级：13级通信工程班 专业：通信工程 学号： 学生姓名：指导老师： 时间：2015年11月8日 摘要 熵是源于物理学的基本概念，后来Shannon在信息论中引入了信息熵的概念，它在统计

物理中的成功使人们对熵的理论和应用有了广泛和高度的重视。最大熵原理是一种在实际问题中已得到广泛应用的信息论方法。本文从信息熵的概念出发，对最大熵原理做了简要介绍，并论述了最大熵原理的合理性，最后提及它在一些领域的应用，通过在具体例子当中应用最大熵原理，展示该原理的适用场合，以期对最大熵原理及其应用有更深刻的理解。 关键词：熵；信息熵；最大熵原理；不适定性问题 引言 科学技术的发展使人类跨入了高度发展的信息化时代。在政治、军事、经济等各个领域，信息的重要性不言而喻，有关信息理论的研究正越来越受到重视，信息论方法也逐渐被广泛应用于各个领域。 信息论一般指的是香农信息论，主要研究在信息可以度量的前提下如何有效地、可靠地、安全地传递信息，涉及消息的信息量、消息的传输以及编码问题。1948年C.E.Shannon 为解决通信工程中不确定信息的编码和传输问题创立信息论，提出信息的统计定义和信息熵、互信息概念，解决了信息的不确定性度量问题，并在此基础上对信息论的一系列理论和方法进行了严格的推导和证明，使以信息论为基础的通信工程获得了巨大的发展。信息论从它诞生的那时起就吸引了众多领域学者的注意，他们竞相应用信息论的概念和方法去理解和解决本领域中的问题。近年来，以不确定性信息为研究对象的信息论理论和方法在众多领域得到了广泛应用，并取得了许多重要的研究成果。迄今为止，较为成熟的研究成果有:A.N.Kolmogorov在1956年提出的关于信息量度定义的三种方法——概率法，组合法，计算法；A.N.Kolmogorov在1968年阐明并为J.Chaitin在1987年系统发展了的关于算法信息的理论。这些成果大大丰富了信息理论的概念、方法和应用范围。 在信息论中，最大熵的含义是最大的不确定性，它解决的一大类问题是在先验知识不充分的条件下进行决策或推断等。熵方法在谱估计、图象滤波、图象重建、天文信号处理、专家系统等中都有广泛的应用。最大熵原理在实际问题中的应用近年来一直在不断地发展。 1.信息熵的概念 信息熵是将熵概念成功地扩展到信息科学领域。熵是描述客观事物无序性的参数，它最早是由R.Clausius于1865年引入热力学中的一个物理概念，通常称之为热力学熵。后来L.Boltzmann赋予熵统计意义上的解释，称之为统计热力学熵。1929年，匈牙利科学家