现代法谱估计（2）Levinson-Durbin递推算法MATLAB及Python实现

原理

Levinson-Durbin递推算法是解Yule Walker方程的快速有效的算法，Yule Walker方程是p+1元线性方程组，它的一般解法是矩阵求逆或高斯消去法。通常都尽量避免使用矩阵求逆运算，因为它的运算量较大。如果利用高斯消去法直接求解线性方程组，其运算级约在p3数量级，而运用Levinson Durbin算法，可以将运算量减少到p2。
Levinson-Durbin算法是从一阶开始，由p-1阶模型的递推求解p阶模型的参数。用 ap(i)表示递推过程中在阶次为p时AR模型的第i个系数，i=1,2,…,p。G^2表示预测误差功率。将此功率用ρ表示，令 =ρ。再用ρp表示p阶预测时的最小预测误差功率，显然

该式很容易理解，因为零阶预测就相当于信号直通过去，所以信号功率有多大，预测误差功率就是多大。





从上式可以看出，用Levinson-Durbin算法求解AR参数的关键，是自相关函数序列的估计。为保证自相关函数的估计较准确，往往需要较多的样本数据，所以，在短记录的情况下，用Levinson-Durbin算法求解AR参数效果较差。反射系数相当于功率传输到终端接不匹配二端网络时所引起的功率反射程度。

程序和结果

（出于维护版权原因，此处只放截图）
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形式写的，公式是用Mathtype打出来的，故部分地方直接采用了截图形式）
下载链接：
现代法功率谱估计及MSE随SNR曲线