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

基于Burg算法的最大熵谱估计一、实验目的使用Matlab平台实现基于Burg算法的最大熵谱估计二、Burg算法原理现代谱估计是针对经典谱估计方差性能较差、分辨率较低的缺点提出并逐渐发展起来的，其分为参数模型谱估计和非参数模型谱估计。而参数模型谱估计主要有AR模型、MA模型、ARMA模型等，其中AR模型应用最多。ARMA模型功率谱的数学表达式为：212121/1)(piijiqiijijeaebeP其中，P(ejω)为功率谱密度；s2是激励白噪声的方差；ai和bi为模型参数。若ARMA模型中bi全为0，就变成了AR模型，又称线性自回归模型，其是一个全极点模型：2121/)(piijijeaeP研究表明，ARMA模型和MA模型均可用无限阶的AR模型来表示。且AR模型的参数估计计算相对简单。同时，实际的物理系统通常是全极点系统。要利用AR模型进行功率谱估计，必须由Yule-Walker方程求得AR模型的参数。而目前求解Yule-Walker方程主要有三种方法:Levinson-Durbin递推算法、Burg算法和协方差方法。其中Burg算法计算结果较为准确，且对于短的时间序列仍能得到较正确的估计，因此应用广泛。研究最大熵谱估计时，Levinson递推一直受制于反射系数Km的求出。而Burg算法秉着使前、后向预测误差平均功率最小的基本思想，不直接估计AR模型的参数，而是先估计反射系数Km，再利用Levinson关系式求得AR模型的参数，继而得到功率谱估计。Burg定义m阶前、后向预测误差为：mimminxianf0)()()((1)miminximangm0)()()((2)由式（1）和（2）又可得到前、后预测误差的阶数递推公式：)1()()(11ngKnfnfmmmm(3))1()()(11ngnfKngmmmm(4)定义m阶前、后向预测误差平均功率为：NmnmmmngnfP])()([2122(5)将阶数递推公式（3）和（4）代入（5），并令0mmKP，可得NmnmmNmnmmmngnfngnfK12121111])1()([21)1()((6)三、Burg算法递推步骤Burg算法的具体实现步骤：步骤1计算预测误差功率的初始值和前、后向预测误差的初始值，并令m=1。210)(1NnnxNP)()()(00nxngnf步骤2求反射系数NmnmmNmnmmmngnfngnfK12121111])1()([21)1()(步骤3计算前向预测滤波器系数),()()(11imaKiaiammmm1,...,1mimmKma)(步骤4计算预测误差功率12)1(mmmPKP步骤5计算滤波器输出)1()()(11ngKnfnfmmmm)1()()(11ngnfKngmmmm步骤6令m←m+1，并重复步骤2至步骤5，直到预测误差功率Pm不再明显减小。最后，再利用Levinson递推关系式估计AR参数，继而得到功率谱估计。四、程序实现%%%%%%%%%%%%基于Burg算法的最大熵谱估计的Matlab实现%%%%%%%%%%%%%%%设置两正弦小信号的归一化频率分别为0.175和0.20，信噪比SNR=30dB、N=32%%%clear,clc;%清空内存及变量N=32;%设置离散傅里叶变换点数，即最大阶数N为32SNR=30;%信噪比SNR取为30dBfs=1;%采样频率取为1Hzt=1:N;%采样时间点从1变化到Nt=t/fs;%得到归一化频率采样点y=sin(2*pi*0.175*t)+sin(2*pi*0.20*t);%信号归一化频率分别取为0.175和0.20x=awgn(y,SNR);%在信号y中加入高斯白噪声,信噪比为SNR设定的数值M=1;%设置起始计算的阶数M为1P(M)=0;%预测误差功率初值设为0Rx(M)=0;%自相关函数初值设为0forn=1:N%样本数从1变化到NP(M)=P(M)+(abs(x(n)))^2;%计算预测误差功率和的初始值ef(1,n)=x(n);%计算前向预测误差初值，令其等于此时的信号序列eb(1,n)=x(n);%计算后向预测误差初值，令其等于此时的信号序列endP(M)=P(M)/N;%计算出预测误差功率的初始值Rx(M)=P(M);%设定自相关函数初始值M=2;%设置起始计算的阶数M为2A=0;%微分所得反射系数Km的分子，初始值设为0D=0;%微分所得反射系数Km的分母，初始值设为0forn=M:N%AR阶数由M变化到NA=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阶前、后向预测误差平均功率)endKm=-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);forn=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的分母，初始值设为0forn=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;endKm=-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);form=1:M-2%AR阶数m由1变化到（M-2）a(M-1,m)=a(M-2,m)+Km*a(M-2,M-1-m);%递推计算各阶前向预测滤波器的系数endwhileFPE(M-1)TH%比较此刻阶数的误差平均功率与之前设置的平均功率的大小TH=FPE(M-1);%小于之前数值时，覆盖得到新的最小平均功率，并进行递推运算forn=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))^2);%预测误差功率FPE(M-1)=P(M)*(N+M)/(N-M);form=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;sum1=0;%采样周期T赋值；功率谱初值设为0f=0.01:0.01:0.5;%选取数据采样点，归一化频率0.01～0.5，间隔为0.01form=1:M-1;%AR阶数为1～(M-1)sum1=sum1+a(M-1,m)*exp(-j*2*pi*m*f*T);%傅里叶变换，得到AR参数的估计ends1=(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算法的基本思想是使前、后向预测误差的平均功率最小，且没有使用自相关估计法，因而，结果与真实值更接近。