功率谱估计的参数模型方法---AR模型谱估计：自相关法&协方差法

第一篇自己的博客，这篇文章将参考文献[1]中的AR模型方法与参考文献[2]中的方法联系起来，算是个笔记吧。看到这个代码有点懵acm.m，这样简单的问题都思考了好一会儿了，LaTex也不怎么会用。。。我的知识真是太匮乏了

信号的线性模型（零极点模型/Pole-Zero Modeling）

常系数线性差分方程：
$Y[n]+a_{1}Y[n-1]+\cdots +a_{p}Y[n-p]=b_{0}X[n]+b_{1}X[n-1]+\cdots +b_{q}X[n-q]$
其中，$X[n]$是输入随机序列，$Y[n]$是输出序列。

（1）AR模型（自回归/全零点模型）

如果则如果$q=0$（右端仅有一项）
$Y[n]+a_{1}Y[n-1]+\cdots +a_{p}Y[n-p]=b_{0}X[n]$
记为$AR(p)$，称$Y[n]$为$p$阶自回归模型

（2）MA模型（移动平均/全极点模型）

如果则如果$p=0$（左端仅有一项）
$Y[n]=b_{0}X[n]+b_{1}X[n-1]+\cdots +b_{q}X[n-q]$
记为$MA(q)$，称$Y[n]$为$q$阶移动平均模型

（3）ARMA模型（自回归/全零点模型）

如果则如果$p,q0$记为$ARMA(p，q)$，称$Y[n]$为$(p,q)$阶自回归移动平均模型

常系数线性差分方程对应的系统函数$H(z)$为
$H(z)=\frac{b_0+b_1{z}^{-1}+\cdots +b_q{z}^{-q}}{1+a_1{z}^{-1}+\cdots +a_p{z}^{-p}})$
将随机信号建模为零均值单位方差白噪声激励$H(z)$产生：假设输入随机序列$X[n]$是零均值单位方差白噪声，功率谱为常数$1$，自相关函数为单位脉冲函数$\delta (k)$,则输出功率谱为
$P(e^{j\omega })=\frac{|b_0+b_1{e}^{-j\omega }+\cdots +b_q{e}^{-jq\omega }{|}^2}{|1+a_1{e}^{-j\omega }+\cdots +a_p{e}^{-jp\omega }{|}^2}$
一旦模型选定，下一步是由给定数据估计参数模型$a_i,b_j$,对于AR模型需要估计$b_0,a_1,\cdots ,a_p$。基于最小均方误差准则的四种参数估计方法：自相关法（Yule-Walk 法）、协方差法、修正协方差法、Burg法，后两种方法最小化前向加反向预测误差的平方和来求解参数。
$arg\underset{a_k}{\max }{\xi }_{MSE}$
其中 ξ M S E = E { ∣ e ( n ) ∣ 2 } \xi_{MSE}=E\{|{e(n)|}^2\} ξMSE​=E{∣e(n)∣2},$e(n)=x(n)-\hat{x}(n)$

自相关法（Yule-Walk 法）

为求解$a_k$,取${\xi }_{MSE}$对$a_k$求偏导 等于$0$,化简得
${\sum }_{l=1}^p{a}_l{r}_x(k-l)=-r_x(k)k=1,2,\cdots ,p$
写成矩阵形式为：
${\mathbf{R}}_{\mathbf{x}}\mathbf{a}=-{\mathbf{r}}_{\mathbf{x}}$ \quad 参考文献[1]中的式(3.7)
此式被称为线性预测的正则方程（Normal equations），其中${\mathbf{R}}_{\mathbf{x}}$是数据$x(n)$的自相关矩阵，是$N\ast N$的Hermite、Toeplitz矩阵：
${\mathbf{R}}_{\mathbf{x}}=\left[\begin{array}{cccc}{r}_x(0)& r_x^{\ast }(1)& \cdots & r_x^{\ast }(p-1)\\ r_x(1)& r_x(0)& \cdots & r_x^{\ast }(p-2)\\ \cdots & \cdots & \cdots & \cdots \\ r_x(p-1)& r_x(p-2)& \cdots & r_x(0)\end{array}\right]$
其中：
$r_x(k)=\frac{1}{N}{\sum }_{n=k}^{N}x(n)x^{\ast }(n-k)$$N$为数据长度，$whenn<0orn>N,x(n)=0$
$a={[a_1,a_2,\cdots ,a_p]}^T$
式中${\mathbf{R}}_{\mathbf{x}}$又可以写成
${\mathbf{R}}_{\mathbf{x}}={\mathbf{X}}^H\mathbf{X}$
$\mathbf{X}=\left[\begin{array}{ccccc}x(0)& 0& 0& \cdots & 0\\ x(1)& x(1)& 0& \cdots & 0\\ x(2)& x(1)& x(0)& \cdots & 0\\ \cdots & \cdots & \cdots & \cdots & \cdots \\ x(N)& x(N-1)& x(N-2)& \cdots & x(N-p+1)\\ 0& x(N)& x(N-1)& \cdots & x(N-p+2)\\ \cdots & \cdots & \cdots & \cdots & \cdots \\ 0& 0& 0& \cdots & x(N)\end{array}\right]$
${\mathbf{r}}_{\mathbf{x}}$可以写成${\mathbf{R}}_{\mathbf{x}}={\mathbf{X}}^H{\mathbf{X}}_2$，${\mathbf{X}}_2$是$\mathbf{X}$的第二列。
因此正则方程可写为：
${\mathbf{X}}^H\mathbf{X}\mathbf{a}=-{\mathbf{X}}^H{\mathbf{X}}_2$
即：
$\mathbf{X}\mathbf{a}=-{\mathbf{X}}_2$ \quad 参考文献[2]中的式(4.21)(4.81)(4.123)
利用自相关法估计参数的代码acm.m在这里

协方差法

Reference

1.杨绿溪. 现代数字信号处理[M]. 2007.
2.Hayes M H. Statistical Digital Signal Processing and Modeling[J]. 1996. Chapter 4, 式(4.40) \quad 书中MATLAB代码在这里
3.陈明. 信息与通信工程中的随机过程(第2版)[M]. 2005.