原理

AR模型的系统函数可以表示为：

如果在白噪声 激励下模型的输出为x(n)，则模型输入、输出关系的时域表达式为：

此式为AR模型的差分方程。将白噪声 激励AR模型产生的输出x(n)叫做AR过程。
根据相关卷积定理，若y(n)=x(n)*h(n)，则有

即卷积的相关等于相关的卷积。如果对上式两边求傅里叶变换，根据维纳辛钦定理和相关定理，有

即输出自功率谱等于输入自功率谱与系统能量谱的乘积。
根据谱分解定理，任何平稳随机信号x(n)都可以看成是由高斯白噪声激励一个因果稳定的可逆系统H(z)产生的输出，如下图所示

因此上述公式可以变换为

上式说明，平稳随机信号x(n)的功率谱可以用模型H(z)的参数来表示。
对于p阶AR模型的系统函数

可以看出其由p+1个待定参数：a1-ap和G。
系统输出x（n）可表示为：

上式可以解释为：用n时刻之前的p个值的线性组合来预测n时刻的值x(n)，预测误差为Gw(n) 。在均方误差最小准则下，组合系数的选择应该使预测误差的均方误差最小。令均方误差为e(n)，有

则有

由于

可以计算出

在最小均方误差（MMSE）准则下，要使预测值最佳地逼近x(n)，参数的选择应使

也即

可得到

即

由上式可得p个方程，写成矩阵式为

上式用到自相关函数的偶对称性质。由这p个方程，可以求出p个参数ai。有了这些参数就可以根据自相关函数和参数ai求解系统增益G。联立可得

上式即为AR模型的正则方程，也叫Yule-Walker方程。
利用以上求出的p+1个参数，我们可以可以估计出该随机信号的功率谱。

程序及结果

（出于维护版权原因，此次只放截图）
MATLAB
程序：


结果：

Python
程序：






结果：

分析

由上图可见，我给程序输入的N为256，取的阶数p为128，信号中f1=0.1，f2=0.13，在图中我们可以看到对应的位置出现了峰值，由于给定幅度大小不同，故峰值大小有很大差别。与经典法谱估计不同的是，图中只出现了一半，而不是像经典法中的出现四个两两对称的峰，这是由于程序中使用了freqz函数（Python中为scipy.signal.freqz）求其频率响应，MATLAB中我们可以打这个函数的帮助手册：

h :Complex, n-element frequency response vector.
h:复n元频率响应矢量。
w: Frequency vector of length n, in radians/sample. w consists of n points equally spaced around the upper half of the unit circle (from 0 to π radians/sample).
w：长度n的频率矢量，单位为弧度/样本。w由n个点组成，这些点平均分布在单位圆的上半部分（从0到π弧度/样本）。

可见w的范围是0到π弧度，所以用freqz函数求出的频率响应只有0到π部分。由于h是复频率响应矢量，因此要求其功率谱要对频率响应求模平方。然后我利用了寻峰函数来得出峰值对应的横坐标（Python中是自己写了一个寻峰估计频率的方法），程序中将所有的峰值及其索引分别放入到两个数组中，且根据峰值大小降序排列，确定最大的两个峰值和其索引位置，然后由索引位置除以2N（freqz函数得出的是单位圆上半部分,频率范围对应为0-0.5，而对应的序列长度为N，因此根据比例关系可知，索引值/N=频率值/0.5也即，频率值=索引值/2N）。

（因原博客是以word形式写的，CSDN不支持Mathtype公式，故部分地方直接采用了截图形式）