基于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 模型、MAP(e j )21b i ei 1a i e ja i e jf m ( n)ma m (i)x( n i)(1)mg m( n) a m(m i)x( n i)i 0由式(1 )和(2)又可得到前、后预测误差的阶数递推公式:f m(n) f m i (n) K mg m i(n 1)g m( n) K m f mi( n)g m 1(n 1)⑷定义m阶前、后向预测误差平均功率为：1N2 2P m - [|f m(门)g m(门)]2n m将阶数递推公式(3)和(4)代入(5)，并令一也0，可得K mNf m 1( n)g m 1( n 1)n m 1K m 1 N-[f m 1( n)2 g m 1(n 1)2]2 n m 1三、Burg算法递推步骤Burg算法的具体实现步骤:步骤1计算预测误差功率的初始值和前、后向预测误差的初始值，并令2x(n)f°( n) g°( n) x(n)步骤2求反射系数Nf 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 1a m(m) K m步骤4计算预测误差功率步骤5计算滤波器输出f m ( n) f m 1( n) K mg 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 为32SNR=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变化到NP(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 为2A=0; %微分所得反射系数 Km 的分子，初始值设为 0 D=0; %微分所得反射系数 Km 的分母，初始值设为 0 for n=M:N%AR 阶数由M 变化到NP m (1 K m 2)P m 1A=A+ef(M-1,n)*eb(M-1,n-1);%计算分子的和D=D+(abs(ef(M-1,n))F2+(abs(eb(M-1,n-1))F2; %计算分母的和(即M 阶前、后向预测误差平均功率) end Km=-2*A/D;%计算反射系数Km (此时起始阶数为2)a(M-1,M-1)=-2*A/D;%计算前向预测滤波器系数P(M)=P(M-1)*(1-(abs(Km))A2);% 计算预测误差功率% 设置最大预测误差平均功率FPE(M-1)=P(M)*(N+M)/(N-M);TH=FPE(M-1);for n=M:N%AR 阶数由M 变化到Nef(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); %计算滤波器输出的后向预测误差endM=M+1;% 阶数叠加，以便递推计算下一阶数据A=0;% 反射系数Km 的分子，初始值设为0D=0;%反射系数Km 的分母，初始值设为0for n=M:N%同前，进行递推运算A=A+ef(M-1,n)*eb(M-1,n-1);D=D+(abs(ef(M-1, n))F2+(abs(eb(M-1,门-1)))人2; endKm=-2*A/D;a(M-1,M-1)=-2*A/D;P(M)=P(M-1)*(1-(abs(Km))A2);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); %递推计算各阶前向预测滤波器的系数endwhile FPE(M-1)<TH %比较此刻阶数的误差平均功率与之前设置的平均功率的大小TH=FPE(M-1); %小于之前数值时，覆盖得到新的最小平均功率，并进行递推运算for n=M:Nef(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); %递推计算滤波器输出的各阶后向预测误差endKm=-2*A/D; % 反射系数a(M-1,M-1)=-2*A/D;P(M)=P(M-1)*(1-(abs(Km))A2);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); %递推得到各阶前向预测滤波器的系数endendT=1/fs;sum仁0; %采样周期T赋值；功率谱初值设为0f=0.01:0.01:0.5; %选取数据采样点，归一化频率0.01-0.5，间隔为0.01 for m=1:M-1; %AR 阶数为1 〜(M-1)sum仁sum1+a(M-1,m)*exp(-j*2*pi*m*f*T); %傅里叶变换，得到AR 参数的估计end%由Fejer-Riesz定理，得到最大熵谱估计，即ARMA功率谱xlabel('f/fs'); %X 轴坐标名称五、仿真结果及分析%前向预测滤波器系数%预测误差功率s1= (abs(1+sum1)).A2;s=P(M)*T./s1; %求得各阶功率谱的矩阵plot(f,10*log10(s),'k'); %画出功率谱随频率变化的曲线图ylabel('功率谱/dB'); %Y轴坐标名称卑悭I4E5结果分析：如上，是阶数分别为16和32时所得的功率谱曲线。

如图可知，Burg算法得到的谱线分辨率很高，曲线平滑而且干扰很少，可清晰分辨出功率谱的峰值及其对应的频率值（处于初始设定的频率值0.175〜0.20之间）;并且随着阶数的增加，预测误差功率随之减小，但变化不大，对于频率分辨率影响也很小。

相较于Levinson算法，Burg算法的基本思想是使前、后向预测误差的平均功率最小，且没有使用自相关估计法，因而，结果与真实值更接近。