原理
前面的Yule Walker方程和Levinson Durbin算法都用到了信号的自相关序列,但是这样可能会存在自相关估计不准的问题(默认为序列长度为N,N的取值以外取不到的点都默认为0)。而Burg是一种不需要自相关函数,直接由观测数据求解反射系数的方法,是一种与预测误差格型滤波器密切相关的算法。
首先将预测的概念予以推广,不再局限于由过去估计现在或未来,也可以由现在估计过去。为便于区分,将前者叫做前向预测,后者叫做后向预测,将“由随机序列的一些已知值的线性组合去估计序列的未知值”称作线性预测。将p阶前向预测表示为:
式中上角标“f”表示“前向”(forward),ap(k)表示前向预测系数。将p阶前向预测误差(FPE)表示为:
相应的,p阶后向预测表示为:
式中上角标“b”表示“后向”(backward),bp(k)表示预测系数。将p阶后向预测误差(BPE)表示为:
以证明,在实系数情况下,前向预测系数与后向预测系数相等,即
所以上式可以写为:
借助Levinson-Durbin关系式,可以得到:
递推的初始值为
即零阶预测时预测误差等于信号值,因为零阶相当于信号直接通过。
一般将递推过程中不断变化的阶次用m表示。于是可将计算m阶预测误差的递推公式表示如下:
Burg算法是使前向预测均方误差和后向预测均方误差之和最小来求取km的,令
于是有
解得
对于平稳随机过程,可以用时间平均代替集合平均,因此上式可以写成
从严格意义上说,“时间平均”应除以求和项的总项数,但该因式对分子和分母是相同的,所以可以不考虑。
程序及结果
MATLAB
(出于维护版权原因,此处只放截图)
程序:
结果:
Python
程序:
结果:
分析
由上图可见,我给程序输入的N为256,取的阶数p为128,信号中f1=0.1,f2=0.13,在图中我们可以看到对应的位置出现了峰值。w的范围是0到π弧度,也即用freqz函数求出的频率响应只有0到π部分。由于h是复频率响应矢量,因此要求其功率谱要对频率响应求模平方。利用寻峰函数findpeaks来得出峰值对应的横坐标(Python中是自己写了一个寻峰估计频率的方法),程序中将所有的峰值及其索引分别放入到两个数组中,且根据峰值大小降序排列,确定最大的两个峰值和其索引位置,然后由索引位置除以2N (freqz函数得出的是单位圆上半部分,频率范围对应为0-0.5,而对应的序列长度为N,因此根据比例关系可知,索引值/N=频率值/0.5也即,频率值=索引值/2N)。
(因原博客用word编写,CSDN不支持Mathtype公式,故部分采用截图形式)
下载链接:
现代法功率谱估计及MSE随SNR曲线