自回归（AR）模型的功率谱估计

随机信号的参数模型

假定随机信号 $x(n)$ 是由白噪声 $w(n)$ 激励某一确定系统的响应。如下图所示：

随机信号 $x(n)$、白噪声 $w(n)$和系统的冲击响应 $h(n)$ 之间的关系为：$$x(n)=h(n)\ast w(n)=\sum _{k=-\infty }^{+\infty }h(k)w(n-k)$$其中，$\ast$ 为卷积操作。如果确定白噪声 $w(n)$ 的参数（假设为$a_k$），也就相当于确定系统中的参数 $h(k)$。这样，就能将研究随机信号转化成研究产生随机信号的系统。

自回归滑动平均模型

在自回归滑动平均（Auto Regression Moving Average，ARMA）模型中，随机噪声 $x(n)$ 由白噪声 $w(n)$、信号的若干过去值 $x(n-k)$ 和若干激励值 $w(n-k)$线性组合而成。用线性差分方程表示为$$x(n)+\sum _{k=1}^p{a}_{k}x(n-k)=w(n)+\sum _{k=1}^p{b}_{k}w(n-k)$$对两边进行z变换，有：$$\sum _{k=0}^p{a}_{k}X(z)z^{-k}=\sum _{k=0}^p{b}_{k}W(z)z^{-k}$$其中，$a_k$ 为模型的参数，p为模型的阶数。系统函数为：$$H(z)=\frac{X(z)}{W(z)}=\frac{\sum _{k=0}^p{b}_k{z}^{-k}}{1+\sum _{k=1}^p{a}_k{z}^{-k}}$$

滑动平均模型

在滑动平均（Moving Average，MR）模型中，随机信号 $x(n)$ 由白噪声 $w(n)$ 和若干激励值 $w(n-k)$ 线性组合而成。用线性差分方程表示为：$$x(n)=w(n)+\sum _{k=1}^p{b}_{k}w(n-k)$$对两边进行z变换，有：$$X(z)=\sum _{k=0}^p{b}_{k}W(z)z^{-k}$$其中，$a_k$ 为模型的参数，p为模型的阶数。系统函数为：$$H(z)=\frac{X(z)}{W(z)}=X(z)=\sum _{k=0}^p{b}_k{z}^{-k}$$

自回归模型

在自回归（Auto Regression，AR）模型中，随机信号 $x(n)$ 由白噪声 $w(n)$ 和信号本身的若干次过去值 $x(n-k)$ 线性组合而成。用线性差分方程表示为：$$x(n)+\sum _{k=1}^p{a}_{k}x(n-k)=w(n)$$对两边进行z变换，有：$$\sum _{k=0}^p{a}_{k}X(z)z^{-k}=W(z)$$其中，$a_k$ 为模型的参数，p为模型的阶数。系统函数为：$$H(z)=\frac{X(z)}{W(z)}=\frac{1}{W(z)}=\frac{1}{1+\sum _{k=1}^p{a}_k{z}^{-k}}$$

自回归模型的参数估计

信号处理中的自相关函数，反映同一个信号在不同时刻取值间的相关程度。信号的自相关函数和信号的功率谱密度互为傅里叶变换，通过信号的自相关函数可以求得信号的功率谱密度。

随机噪声 $x(n)$ 的自相关函数 $R_{xx}(m)$ 为：$$R_{xx}(m)=E[x(n)x(n+m)]$$将随机信号 $x(n)$ 的表达式代入，得： R x x ( m ) = E { x ( n ) [ w ( n + m ) − ∑ k = 1 p a k x ( n + m − k ) ] } = E [ x ( n ) w ( n + m ) ] − ∑ k = 1 p a k E [ x ( n ) x ( n + m − k ) ] = R x w ( m ) − ∑ k = 1 p a k R x x ( m − k )

Rxx​(m)​=E{x(n)[w(n+m)−k=1∑p​ak​x(n+m−k)]}=E[x(n)w(n+m)]−k=1∑p​ak​E[x(n)x(n+m−k)]=Rxw​(m)−k=1∑p​ak​Rxx​(m−k)​由于 R x w ( m ) = { σ p 2 m=0 0 m>0 R_{xw}(m)=

{σp2m=00m>0" role="presentation">{σ2p0m=0m>0$$\{\begin{array}{ll}{\sigma }_p^2& \text{m=0}\\ 0& \text{m>0}\end{array}$$

Rxw​(m)={σp2​0​m=0m>0​则自相关函数 $R_{xx}(m)$ 可以表示为： R x x ( m ) = { σ p 2 − ∑ k = 1 p a k R x x ( m − k ) m=0 − ∑ k = 1 p a k R x x ( m − k ) m>0 R_{xx}(m)=

{σp2−∑k=1pakRxx(m−k)m=0−∑k=1pakRxx(m−k)m>0" role="presentation">⎧⎩⎨⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪σ2p−∑k=1pakRxx(m−k)−∑k=1pakRxx(m−k)m=0m>0$$\{\begin{array}{ll}{\sigma }_p^2-{\displaystyle \sum _{k=1}^p{a}_k{R}_{xx}(m-k)}& \text{m=0}\\ -{\displaystyle \sum _{k=1}^p{a}_k{R}_{xx}(m-k)}& \text{m>0}\end{array}$$

Rxx​(m)=⎩⎪⎪⎪⎪⎨⎪⎪⎪⎪⎧​σp2​−k=1∑p​ak​Rxx​(m−k)−k=1∑p​ak​Rxx​(m−k)​m=0m>0​ 其中，$m=0,1,...,p$。用矩阵（Yule-Walker方程）可以表示为： [ R x x ( 0 ) R x x ( − 1 ) ⋯ R x x ( − p ) R x x ( 1 ) R x x ( 0 ) ⋯ R x x ( 1 − p ) ⋮ ⋮ ⋱ ⋮ R x x ( p ) R x x ( p − 1 ) ⋯ R x x ( 0 ) ] [ 1 a 1 ⋮ a p ] = [ σ p 2 0 ⋮ 0 ]

[Rxx(0)Rxx(−1)⋯Rxx(−p)Rxx(1)Rxx(0)⋯Rxx(1−p)⋮⋮⋱⋮Rxx(p)Rxx(p−1)⋯Rxx(0)]" role="presentation">⎡⎣⎢⎢⎢⎢⎢Rxx(0)Rxx(1)⋮Rxx(p)Rxx(−1)Rxx(0)⋮Rxx(p−1)⋯⋯⋱⋯Rxx(−p)Rxx(1−p)⋮Rxx(0)⎤⎦⎥⎥⎥⎥⎥$$\left[\begin{array}{cccc}{R}_{xx}(0)& R_{xx}(-1)& \cdots & R_{xx}(-p)\\ R_{xx}(1)& R_{xx}(0)& \cdots & R_{xx}(1-p)\\ ⋮& ⋮& \ddots & ⋮\\ R_{xx}(p)& R_{xx}(p-1)& \cdots & R_{xx}(0)\end{array}\right]$$

[1a1⋮ap]" role="presentation">⎡⎣⎢⎢⎢⎢1a1⋮ap⎤⎦⎥⎥⎥⎥$$\left[\begin{array}{c}1\\ a_1\\ ⋮\\ a_p\end{array}\right]$$

=

[σp20⋮0]" role="presentation">⎡⎣⎢⎢⎢⎢σ2p0⋮0⎤⎦⎥⎥⎥⎥$$\left[\begin{array}{c}{\sigma }_p^2\\ 0\\ ⋮\\ 0\end{array}\right]$$

⎣⎢⎢⎢⎡​Rxx​(0)Rxx​(1)⋮Rxx​(p)​Rxx​(−1)Rxx​(0)⋮Rxx​(p−1)​⋯⋯⋱⋯​Rxx​(−p)Rxx​(1−p)⋮Rxx​(0)​⎦⎥⎥⎥⎤​⎣⎢⎢⎢⎡​1a1​⋮ap​​⎦⎥⎥⎥⎤​=⎣⎢⎢⎢⎡​σp2​0⋮0​⎦⎥⎥⎥⎤​目前，求解Yule-Walker方程的方法主要有Yule-Walker法、Levinson-Durbin法、Burg法、协方差法和修正协方差法这五种，这五种方法有各自的优缺点。下面分别介绍这五种方法：

Yule-Walker法（预测误差最小）

对于预测模型，假设第n个数据 $x(n)$ 是由前p个数据 $x(n-p),x(n-p+1),...,x(n-1)$ 预测得到的。即：$$\hat{x}(n)=-\sum _{k=1}^p{a}_{k}x(n-k)$$那么预测误差 $e(n)$ 为：$$e(n)=x(n)+\hat{x}(n)=x(n)+\sum _{k=1}^p{a}_{k}x(n-k)$$而对于自回归（AR）模型而言，有：$$w(n)=x(n)+\sum _{k=1}^p{a}_{k}x(n-k)$$如果自回归模型和预测模型的系数 $a_k$ 相同，那么预测误差 $e(n)$ 就等于白噪声 $w(n)$。因此，Yule-Walker法的基本思想就是计算预测误差功率 ${\rho }_e$ 最小时模型的预测系数 $a_k$：$${\rho }_e=\frac{1}{N}E[e^2(n)]=\frac{1}{N}E[x(n)+\sum _{k=1}^p{a}_{k}x(n-k){]}^2$$为了使预测误差功率最小，有：$$\frac{\partial {\rho }_e}{\partial a_k}=\frac{\partial E[e^2(n)]}{\partial a_k}=0$$而 E [ e 2 ( n ) ] = E { [ x ( n ) + ∑ k = 1 p a k x ( n − k ) ] 2 } = E [ x 2 ( n ) ] + 2 ∑ k = 1 p a k E [ x ( n ) x ( n − k ) ] + ∑ i = 1 p a i ∑ k = 1 p a k E [ x ( n − i ) x ( n − k ) ] = R x x ( 0 ) + 2 ∑ k = 1 p a k R x x ( k ) + ∑ i = 1 p a i ∑ k = 1 p a k R x x ( i − k )

E[e2(n)]​=E{[x(n)+k=1∑p​ak​x(n−k)]2}=E[x2(n)]+2k=1∑p​ak​E[x(n)x(n−k)]+i=1∑p​ai​k=1∑p​ak​E[x(n−i)x(n−k)]=Rxx​(0)+2k=1∑p​ak​Rxx​(k)+i=1∑p​ai​k=1∑p​ak​Rxx​(i−k)​则有：$$\frac{\partial E[e^2(n)]}{\partial a_k}=2R_{xx}(k)+2\sum _{i=1}^p{a}_i{R}_{xx}(i-k)=0$$即：$$\sum _{i=1}^p{a}_i{R}_{xx}(i-k)=-R_{xx}(k)$$其中，$k=1,2,...,p$。写成矩阵形式： [ R x x ( 0 ) R x x ( 1 ) ⋯ R x x ( p − 1 ) R x x ( − 1 ) R x x ( 0 ) ⋯ R x x ( p − 2 ) ⋮ ⋮ ⋱ ⋮ R x x ( 1 − p ) R x x ( 2 − p ) ⋯ R x x ( 0 ) ] [ a 1 a 2 ⋮ a p ] = [ − R x x ( 1 ) − R x x ( 2 ) ⋮ − R x x ( p ) ]

[Rxx(0)Rxx(1)⋯Rxx(p−1)Rxx(−1)Rxx(0)⋯Rxx(p−2)⋮⋮⋱⋮Rxx(1−p)Rxx(2−p)⋯Rxx(0)]" role="presentation">⎡⎣⎢⎢⎢⎢⎢Rxx(0)Rxx(−1)⋮Rxx(1−p)Rxx(1)Rxx(0)⋮Rxx(2−p)⋯⋯⋱⋯Rxx(p−1)Rxx(p−2)⋮Rxx(0)⎤⎦⎥⎥⎥⎥⎥$$\left[\begin{array}{cccc}{R}_{xx}(0)& R_{xx}(1)& \cdots & R_{xx}(p-1)\\ R_{xx}(-1)& R_{xx}(0)& \cdots & R_{xx}(p-2)\\ ⋮& ⋮& \ddots & ⋮\\ R_{xx}(1-p)& R_{xx}(2-p)& \cdots & R_{xx}(0)\end{array}\right]$$

[a1a2⋮ap]" role="presentation">⎡⎣⎢⎢⎢⎢a1a2⋮ap⎤⎦⎥⎥⎥⎥$$\left[\begin{array}{c}{a}_1\\ a_2\\ ⋮\\ a_p\end{array}\right]$$

=

[−Rxx(1)−Rxx(2)⋮−Rxx(p)]" role="presentation">⎡⎣⎢⎢⎢⎢⎢−Rxx(1)−Rxx(2)⋮−Rxx(p)⎤⎦⎥⎥⎥⎥⎥$$\left[\begin{array}{c}-R_{xx}(1)\\ -R_{xx}(2)\\ ⋮\\ -R_{xx}(p)\end{array}\right]$$

⎣⎢⎢⎢⎡​Rxx​(0)Rxx​(−1)⋮Rxx​(1−p)​Rxx​(1)Rxx​(0)⋮Rxx​(2−p)​⋯⋯⋱⋯​Rxx​(p−1)Rxx​(p−2)⋮Rxx​(0)​⎦⎥⎥⎥⎤​⎣⎢⎢⎢⎡​a1​a2​⋮ap​​⎦⎥⎥⎥⎤​=⎣⎢⎢⎢⎡​−Rxx​(1)−Rxx​(2)⋮−Rxx​(p)​⎦⎥⎥⎥⎤​
则预测系数 $a_k$ 可以表示成： [ a 1 a 2 ⋮ a p ] = [ R x x ( 0 ) R x x ( 1 ) ⋯ R x x ( p − 1 ) R x x ( − 1 ) R x x ( 0 ) ⋯ R x x ( p − 2 ) ⋮ ⋮ ⋱ ⋮ R x x ( 1 − p ) R x x ( 2 − p ) ⋯ R x x ( 0 ) ] − 1 [ − R x x ( 1 ) − R x x ( 2 ) ⋮ − R x x ( p ) ]

[a1a2⋮ap]" role="presentation">⎡⎣⎢⎢⎢⎢a1a2⋮ap⎤⎦⎥⎥⎥⎥$$\left[\begin{array}{c}{a}_1\\ a_2\\ ⋮\\ a_p\end{array}\right]$$

=

[Rxx(0)Rxx(1)⋯Rxx(p−1)Rxx(−1)Rxx(0)⋯Rxx(p−2)⋮⋮⋱⋮Rxx(1−p)Rxx(2−p)⋯Rxx(0)]" role="presentation">⎡⎣⎢⎢⎢⎢⎢Rxx(0)Rxx(−1)⋮Rxx(1−p)Rxx(1)Rxx(0)⋮Rxx(2−p)⋯⋯⋱⋯Rxx(p−1)Rxx(p−2)⋮Rxx(0)⎤⎦⎥⎥⎥⎥⎥$$\left[\begin{array}{cccc}{R}_{xx}(0)& R_{xx}(1)& \cdots & R_{xx}(p-1)\\ R_{xx}(-1)& R_{xx}(0)& \cdots & R_{xx}(p-2)\\ ⋮& ⋮& \ddots & ⋮\\ R_{xx}(1-p)& R_{xx}(2-p)& \cdots & R_{xx}(0)\end{array}\right]$$

^{-1}

[−Rxx(1)−Rxx(2)⋮−Rxx(p)]" role="presentation">⎡⎣⎢⎢⎢⎢⎢−Rxx(1)−Rxx(2)⋮−Rxx(p)⎤⎦⎥⎥⎥⎥⎥$$\left[\begin{array}{c}-R_{xx}(1)\\ -R_{xx}(2)\\ ⋮\\ -R_{xx}(p)\end{array}\right]$$

⎣⎢⎢⎢⎡​a1​a2​⋮ap​​⎦⎥⎥⎥⎤​=⎣⎢⎢⎢⎡​Rxx​(0)Rxx​(−1)⋮Rxx​(1−p)​Rxx​(1)Rxx​(0)⋮Rxx​(2−p)​⋯⋯⋱⋯​Rxx​(p−1)Rxx​(p−2)⋮Rxx​(0)​⎦⎥⎥⎥⎤​−1⎣⎢⎢⎢⎡​−Rxx​(1)−Rxx​(2)⋮−Rxx​(p)​⎦⎥⎥⎥⎤​
然后根据模型的系数 $a_k$求 ${\sigma }_p^2$：$${\sigma }_p^2=R_{xx}(0)+\sum _{k=1}^p{a}_k{R}_{xx}(-k)$$

Levinson-Durbin法（递归）

假设 $p$ 阶自回归模型的系数为 $a_{pk}$，其中 $p$ 表示此时自回归模型的阶数，也就是当前阶数情况下总共有几个系数，$k$ 表示当前系数的序列索引。对于Yule-Walker方程，

当 $p=1$ 时，有： [ R x x ( 0 ) R x x ( 1 ) R x x ( 1 ) R x x ( 0 ) ] [ 1 a 11 ] = [ σ 1 2 0 ]

= [Rxx​(0)Rxx​(1)​Rxx​(1)Rxx​(0)​][1a11​​]=[σ12​0​]解得：$$a_{11}=-\frac{R_{xx}(1)}{R_{xx}(0)}$$$${\sigma }_1^2=(1-|a_{11}{|}^2)R_{xx}(0)$$当 $p=2$ 时，有： [ R x x ( 0 ) R x x ( 1 ) R x x ( 2 ) R x x ( 1 ) R x x ( 0 ) R x x ( 1 ) R x x ( 2 ) R x x ( 1 ) R x x ( 0 ) ] [ 1 a 21 a 22 ] = [ σ 2 2 0 0 ]

[Rxx(0)Rxx(1)Rxx(2)Rxx(1)Rxx(0)Rxx(1)Rxx(2)Rxx(1)Rxx(0)]" role="presentation">⎡⎣⎢Rxx(0)Rxx(1)Rxx(2)Rxx(1)Rxx(0)Rxx(1)Rxx(2)Rxx(1)Rxx(0)⎤⎦⎥$$\left[\begin{array}{ccc}{R}_{xx}(0)& R_{xx}(1)& R_{xx}(2)\\ R_{xx}(1)& R_{xx}(0)& R_{xx}(1)\\ R_{xx}(2)& R_{xx}(1)& R_{xx}(0)\end{array}\right]$$

[1a21a22]" role="presentation">⎡⎣⎢1a21a22⎤⎦⎥$$\left[\begin{array}{c}1\\ a_{21}\\ a_{22}\end{array}\right]$$

=

[σ2200]" role="presentation">⎡⎣⎢σ2200⎤⎦⎥$$\left[\begin{array}{c}{\sigma }_2^2\\ 0\\ 0\end{array}\right]$$

⎣⎡​Rxx​(0)Rxx​(1)Rxx​(2)​Rxx​(1)Rxx​(0)Rxx​(1)​Rxx​(2)Rxx​(1)Rxx​(0)​⎦⎤​⎣⎡​1a21​a22​​⎦⎤​=⎣⎡​σ22​00​⎦⎤​解得：$$a_{22}=-\frac{R_{xx}(2)-a_{11}{R}_{xx}(1)}{{\sigma }_1^2}$$$$a_{21}=a_{11}+a_{22}{a}_{11}$$$${\sigma }_2^2=(1-|a_{22}{|}^2{\sigma }_1^2)$$推广到 $k$ 阶时，有：$$a_{kk}=-\frac{R_{xx}(k)-\sum _{l=1}^{k-1}{a}_{k-1,l}{R}_{xx}(k-l)}{{\sigma }_{k-1}^2}$$$$a_{ki}=a_{k-1,i}+a_{kk}{a}_{k-1,k-i},i=1,2,...,k-1$$$${\sigma }_k^2=(1-|a_{kk}{|}^2){\sigma }_{k-1}^2,{\sigma }_0^2=R_{xx}(0)$$Levinson-Durbin法的基本思想就是通过递归到所需要的阶数时，求得所估计的模型系数 $a_{pi},i=1,2,...,k$ 和 ${\sigma }_p^2$。

Burg法（前后向平均预测误差最小+递归）

对于预测模型，假设随机信号 $x(n)$ 是由前 $p$ 个数据 $x(n-p),x(n-p+1),...,x(n-1)$ 预测得到的：$$\hat{x}(n)=-\sum _{k=1}^p{a}_{pk}x(n-k)$$则预测误差 $e_p(n)$ 为：$$e_p(n)=x(n)-\hat{x}(n)=x(n)+\sum _{k=1}^p{a}_{pk}x(n-k)$$而随机信号 $x(n-p)$ 是由后 $p$ 个数据 $x(n-p+1),x(n-p+2),...,x(n)$ 预测得到的：$$\hat{x}(n-p)=-\sum _{k=1}^p{a}_{pk}x(n-p+k)$$则预测误差 $b_p(n)$ 为:$$b_p(n)=x(n-p)-\hat{x}(n-p)=x(n-p)+\sum _{k=1}^p{a}_{pk}x(n-p+k)$$称 $e_p(n)$ 和 $b_p(n)$ 分别为前向预测误差和后向预测误差。Burg法的基本思想就是递归地求前后向预测误差平均功率 ${\rho }_p$ 最小时的反射系数 $K_p=a_{pp}$，然后再根据反射系数 $K_p$ 求模型的参数 $a_{pk}$ 和 ${\sigma }_p^2$。

假设随机信号 $x(n)$ 的观测区间为 $0\le n\le N-1$，则前向预测功率 ${\rho }_{pe}$、后向预测功率 ${\rho }_{pb}$ 和平均预测功率 ${\rho }_p$ 分别为：$${\rho }_{pe}=\frac{1}{N-p}\sum _{n=p}^{N-1}|e_p(n){|}^2=\frac{1}{N-p}\sum _{n=p}^{N-1}|x(n)+\sum _{k=1}^p{a}_{pk}x(n-k){|}^2$$$${\rho }_{pb}=\frac{1}{N-p}\sum _{n=p}^{N-1}|b_p(n){|}^2=\frac{1}{N-p}\sum _{n=p}^{N-1}|x(n-p)+\sum _{k=1}^p{a}_{pk}x(n-p+k){|}^2$$ $${\rho }_p=\frac{1}{2}({\rho }_{pe}+{\rho }_{pb})$$当反射系数最小时，有：$$\frac{\partial {\rho }_p}{\partial K_p}=0$$解得：$$K_p=a_{pp}=\frac{-2\sum _{n=p}^{N-1}[e_{p-1}(n)b_{p-1}(n-1)]}{\sum _{n=p}^{N-1}[e_{p-1}^2(n)+b_{p-1}^2(n-1)]}$$由于 $e_0(n)=b_0(n)=x(n)$，且前向预测误差和后向预测误差分别存在递推公式：$$e_p(n)=e_{p-1}(n)+K_p{b}_{p-1}(n-1)$$$$b_p(n)=b_{p-1}(n-1)+K_p{e}_{p-1}(n)$$故可以根据反射系数 $K_p$ 递推地求出模型的参数$a_{pi}$ 和 ${\sigma }_p^2$，有：$$a_{pi}=a_{p-1,i}+K_p{a}_{p-1,p-i},i=1,2,...,p-1$$$${\sigma }_p^2=(1-K_p^2){\sigma }_{p-1}^2$$

协方差法（预测误差最小）

协方差法和Yule-Walker法都是通过最小化预测功率求模型的参数。但不同于Yule-Walker法，协方差法的预测功率为：$${\rho }_e=\frac{1}{N-p}\sum _{n=p}^{N-1}|e^2(n)|=\frac{1}{N-p}\sum _{n=p}^{N-1}|x(n)+\sum _{k=1}^p{a}_{pk}x(n-k){|}^2$$当预测功率最小时，有：$$\frac{\partial {\rho }_e}{\partial a_{pk}}=\frac{1}{N-p}\sum _{n=p}^{N-1}[x(n)+\sum _{k=1}^p{a}_{pk}x(n-k)]x(n-l)=0$$其中，$l=1,2,...,p$。用 ${\hat{c}}_{xx}(k,l)$ 来表示上式，得：$$\sum _{k=1}^p{a}_{pk}{\hat{c}}_{xx}(l,k)=-{\hat{c}}_{xx}(l,0)$$其中： c ^ x x ( k , l ) = { 1 N − p ∑ n = p N − 1 [ x ( n − k ) x ( n − l ) ] c ^ x x ( l , k ) \hat{c}_{xx}(k,l)=

c^xx​(k,l)=⎩⎪⎨⎪⎧​N−p1​n=p∑N−1​[x(n−k)x(n−l)]c^xx​(l,k)​写成矩阵形式为： [ c ^ x x ( 1 , 1 ) c ^ x x ( 1 , 2 ) ⋯ c ^ x x ( 1 , p ) c ^ x x ( 2 , 1 ) c ^ x x ( 2 , 2 ) ⋯ c ^ x x ( 2 , p ) ⋮ ⋮ ⋱ ⋮ c ^ x x ( p , 1 ) c ^ x x ( p , 2 ) ⋯ c ^ x x ( p , p ) ] [ a p 1 a p 2 ⋮ a p p ] = [ − c ^ x x ( 1 , 0 ) − c ^ x x ( 2 , 0 ) ⋮ − c ^ x x ( p , 0 ) ]

[c^xx(1,1)c^xx(1,2)⋯c^xx(1,p)c^xx(2,1)c^xx(2,2)⋯c^xx(2,p)⋮⋮⋱⋮c^xx(p,1)c^xx(p,2)⋯c^xx(p,p)]" role="presentation">⎡⎣⎢⎢⎢⎢⎢c^xx(1,1)c^xx(2,1)⋮c^xx(p,1)c^xx(1,2)c^xx(2,2)⋮c^xx(p,2)⋯⋯⋱⋯c^xx(1,p)c^xx(2,p)⋮c^xx(p,p)⎤⎦⎥⎥⎥⎥⎥$$\left[\begin{array}{cccc}{\hat{c}}_{xx}(1,1)& {\hat{c}}_{xx}(1,2)& \cdots & {\hat{c}}_{xx}(1,p)\\ {\hat{c}}_{xx}(2,1)& {\hat{c}}_{xx}(2,2)& \cdots & {\hat{c}}_{xx}(2,p)\\ ⋮& ⋮& \ddots & ⋮\\ {\hat{c}}_{xx}(p,1)& {\hat{c}}_{xx}(p,2)& \cdots & {\hat{c}}_{xx}(p,p)\end{array}\right]$$

[ap1ap2⋮app]" role="presentation">⎡⎣⎢⎢⎢⎢⎢ap1ap2⋮app⎤⎦⎥⎥⎥⎥⎥$$\left[\begin{array}{c}{a}_{p1}\\ a_{p2}\\ ⋮\\ a_{pp}\end{array}\right]$$

=

[−c^xx(1,0)−c^xx(2,0)⋮−c^xx(p,0)]" role="presentation">⎡⎣⎢⎢⎢⎢⎢−c^xx(1,0)−c^xx(2,0)⋮−c^xx(p,0)⎤⎦⎥⎥⎥⎥⎥$$\left[\begin{array}{c}-{\hat{c}}_{xx}(1,0)\\ -{\hat{c}}_{xx}(2,0)\\ ⋮\\ -{\hat{c}}_{xx}(p,0)\end{array}\right]$$

⎣⎢⎢⎢⎡​c^xx​(1,1)c^xx​(2,1)⋮c^xx​(p,1)​c^xx​(1,2)c^xx​(2,2)⋮c^xx​(p,2)​⋯⋯⋱⋯​c^xx​(1,p)c^xx​(2,p)⋮c^xx​(p,p)​⎦⎥⎥⎥⎤​⎣⎢⎢⎢⎡​ap1​ap2​⋮app​​⎦⎥⎥⎥⎤​=⎣⎢⎢⎢⎡​−c^xx​(1,0)−c^xx​(2,0)⋮−c^xx​(p,0)​⎦⎥⎥⎥⎤​则可以求出模型的参数 $a_{pk}$： [ a p 1 a p 2 ⋮ a p p ] = [ c ^ x x ( 1 , 1 ) c ^ x x ( 1 , 2 ) ⋯ c ^ x x ( 1 , p ) c ^ x x ( 2 , 1 ) c ^ x x ( 2 , 2 ) ⋯ c ^ x x ( 2 , p ) ⋮ ⋮ ⋱ ⋮ c ^ x x ( p , 1 ) c ^ x x ( p , 2 ) ⋯ c ^ x x ( p , p ) ] − 1 [ − c ^ x x ( 1 , 0 ) − c ^ x x ( 2 , 0 ) ⋮ − c ^ x x ( p , 0 ) ]

[ap1ap2⋮app]" role="presentation">⎡⎣⎢⎢⎢⎢⎢ap1ap2⋮app⎤⎦⎥⎥⎥⎥⎥$$\left[\begin{array}{c}{a}_{p1}\\ a_{p2}\\ ⋮\\ a_{pp}\end{array}\right]$$

=

[c^xx(1,1)c^xx(1,2)⋯c^xx(1,p)c^xx(2,1)c^xx(2,2)⋯c^xx(2,p)⋮⋮⋱⋮c^xx(p,1)c^xx(p,2)⋯c^xx(p,p)]" role="presentation">⎡⎣⎢⎢⎢⎢⎢c^xx(1,1)c^xx(2,1)⋮c^xx(p,1)c^xx(1,2)c^xx(2,2)⋮c^xx(p,2)⋯⋯⋱⋯c^xx(1,p)c^xx(2,p)⋮c^xx(p,p)⎤⎦⎥⎥⎥⎥⎥$$\left[\begin{array}{cccc}{\hat{c}}_{xx}(1,1)& {\hat{c}}_{xx}(1,2)& \cdots & {\hat{c}}_{xx}(1,p)\\ {\hat{c}}_{xx}(2,1)& {\hat{c}}_{xx}(2,2)& \cdots & {\hat{c}}_{xx}(2,p)\\ ⋮& ⋮& \ddots & ⋮\\ {\hat{c}}_{xx}(p,1)& {\hat{c}}_{xx}(p,2)& \cdots & {\hat{c}}_{xx}(p,p)\end{array}\right]$$

^{-1}

[−c^xx(1,0)−c^xx(2,0)⋮−c^xx(p,0)]" role="presentation">⎡⎣⎢⎢⎢⎢⎢−c^xx(1,0)−c^xx(2,0)⋮−c^xx(p,0)⎤⎦⎥⎥⎥⎥⎥$$\left[\begin{array}{c}-{\hat{c}}_{xx}(1,0)\\ -{\hat{c}}_{xx}(2,0)\\ ⋮\\ -{\hat{c}}_{xx}(p,0)\end{array}\right]$$

⎣⎢⎢⎢⎡​ap1​ap2​⋮app​​⎦⎥⎥⎥⎤​=⎣⎢⎢⎢⎡​c^xx​(1,1)c^xx​(2,1)⋮c^xx​(p,1)​c^xx​(1,2)c^xx​(2,2)⋮c^xx​(p,2)​⋯⋯⋱⋯​c^xx​(1,p)c^xx​(2,p)⋮c^xx​(p,p)​⎦⎥⎥⎥⎤​−1⎣⎢⎢⎢⎡​−c^xx​(1,0)−c^xx​(2,0)⋮−c^xx​(p,0)​⎦⎥⎥⎥⎤​

${\sigma }_p^2$ 的估计值： σ p 2 = ρ m i n = 1 N − p ∑ n = p N − 1 [ x ( n ) + ∑ k = 1 p a p k x ( n − k ) ] [ x ( n ) + ∑ l = 1 p a p k x ( n − l ) ] = 1 N − p ∑ n = p N − 1 [ x ( n ) + ∑ k = 1 p a p k x ( n − k ) ] x ( n ) = c ^ x x ( 0 , 0 ) + ∑ k = 1 p a p k c ^ x x ( 0 , k )

σp2​=ρmin​​=N−p1​n=p∑N−1​[x(n)+k=1∑p​apk​x(n−k)][x(n)+l=1∑p​apk​x(n−l)]=N−p1​n=p∑N−1​[x(n)+k=1∑p​apk​x(n−k)]x(n)=c^xx​(0,0)+k=1∑p​apk​c^xx​(0,k)​可以将参数 $a_{pk}$ 和 ${\sigma }_p^2$ 的求解合并在一个矩阵中，有： [ c ^ x x ( 0 , 0 ) c ^ x x ( 0 , 1 ) ⋯ c ^ x x ( 0 , p ) c ^ x x ( 1 , 0 ) c ^ x x ( 1 , 1 ) ⋯ c ^ x x ( 1 , p ) ⋮ ⋮ ⋱ ⋮ c ^ x x ( p , 0 ) c ^ x x ( p , 1 ) ⋯ c ^ x x ( p , p ) ] [ 1 a p 1 ⋮ a p p ] = [ σ p 2 0 ⋮ 0 ]

[c^xx(0,0)c^xx(0,1)⋯c^xx(0,p)c^xx(1,0)c^xx(1,1)⋯c^xx(1,p)⋮⋮⋱⋮c^xx(p,0)c^xx(p,1)⋯c^xx(p,p)]" role="presentation">⎡⎣⎢⎢⎢⎢⎢c^xx(0,0)c^xx(1,0)⋮c^xx(p,0)c^xx(0,1)c^xx(1,1)⋮c^xx(p,1)⋯⋯⋱⋯c^xx(0,p)c^xx(1,p)⋮c^xx(p,p)⎤⎦⎥⎥⎥⎥⎥$$\left[\begin{array}{cccc}{\hat{c}}_{xx}(0,0)& {\hat{c}}_{xx}(0,1)& \cdots & {\hat{c}}_{xx}(0,p)\\ {\hat{c}}_{xx}(1,0)& {\hat{c}}_{xx}(1,1)& \cdots & {\hat{c}}_{xx}(1,p)\\ ⋮& ⋮& \ddots & ⋮\\ {\hat{c}}_{xx}(p,0)& {\hat{c}}_{xx}(p,1)& \cdots & {\hat{c}}_{xx}(p,p)\end{array}\right]$$

[1ap1⋮app]" role="presentation">⎡⎣⎢⎢⎢⎢⎢1ap1⋮app⎤⎦⎥⎥⎥⎥⎥$$\left[\begin{array}{c}1\\ a_{p1}\\ ⋮\\ a_{pp}\end{array}\right]$$

=

[σp20⋮0]" role="presentation">⎡⎣⎢⎢⎢⎢σ2p0⋮0⎤⎦⎥⎥⎥⎥$$\left[\begin{array}{c}{\sigma }_p^2\\ 0\\ ⋮\\ 0\end{array}\right]$$

⎣⎢⎢⎢⎡​c^xx​(0,0)c^xx​(1,0)⋮c^xx​(p,0)​c^xx​(0,1)c^xx​(1,1)⋮c^xx​(p,1)​⋯⋯⋱⋯​c^xx​(0,p)c^xx​(1,p)⋮c^xx​(p,p)​⎦⎥⎥⎥⎤​⎣⎢⎢⎢⎡​1ap1​⋮app​​⎦⎥⎥⎥⎤​=⎣⎢⎢⎢⎡​σp2​0⋮0​⎦⎥⎥⎥⎤​

修正协方差法（前后向平均预测误差最小）

修正的协方差法与Burg法类似，使用前后向预测平均误差最小的方法来估计自回归模型的参数。但不同于Burg法，修正的协方差法的后向预测误差 $b_p(n)$ 与Burg法中的后向预测误差不一样。

假设随机信号 $x(n)$ 是由前 $p$ 个数据 $x(n-p),x(n-p+1),...,x(n-1)$ 预测得到的，则有：$$\hat{x}(n)=-\sum _{k=1}^p{a}_{pk}x(n-k)$$在修正的协方差法中，前向预测误差 $e_p(n)$ 与Burg法中的前向预测误差一致：$$e_p(n)=x(n)-\hat{x}(n)=x(n)+\sum _{k=1}^p{a}_{pk}x(n-k)$$对于后向预测误差 $b_p(n)$，假设随机信号 $x(n)$ 是由后 $p$ 个数据 $x(n+1),x(n+2),...,x(n+p)$ 预测得到的：$$\hat{x}(n)=-\sum _{k=1}^p{a}_{pk}x(n+k)$$那么，后向预测误差 $b_p(n)$ 可以表示为：$$b_p(n)=x(n)-\hat{x}(n)=x(n)+\sum _{k=1}^p{a}_{pk}x(n+k)$$假设随机信号 $x(n)$ 的观测区间为 $0\le n\le N-1$，则前向预测功率 ${\rho }_{pe}$、后向预测功率 ${\rho }_{pb}$ 和平均预测功率 ${\rho }_p$ 可以分别表示为：$${\rho }_{pe}=\frac{1}{N-p}\sum _{n=p}^{N-1}|e_p(n){|}^2=\frac{1}{N-p}\sum _{n=p}^{N-1}|x(n)+\sum _{k=1}^p{a}_{pk}x(n-k){|}^2$$$${\rho }_{pb}=\frac{1}{N-p}\sum _{n=0}^{N-1-p}|b_p(n){|}^2=\frac{1}{N-p}\sum _{n=0}^{N-1-p}|x(n)+\sum _{k=1}^p{a}_{pk}x(n+k){|}^2$$$${\rho }_p=\frac{1}{2}({\rho }_{pe}+{\rho }_{pb})$$要使平均误差功率 ${\rho }_p$ 最小，则： ∂ ρ p ∂ a p k = 1 N − p { ∑ n = p N − 1 [ x ( n ) + ∑ k = 1 p a p k x ( n − k ) ] x ( n − l ) + ∑ n = 0 N − 1 − p [ x ( n ) + ∑ k = 1 p a p k x ( n + k ) ] x ( n + l ) } = 0 \frac{\partial{\rho_p}}{\partial{a_{pk}}}=\frac{1}{N-p}\{\sum^{N-1}_{n=p}[x(n)+\sum^p_{k=1}a_{pk}x(n-k)]x(n-l)+\sum^{N-1-p}_{n=0}[x(n)+\sum^p_{k=1}a_{pk}x(n+k)]x(n+l)\}=0 ∂apk​∂ρp​​=N−p1​{n=p∑N−1​[x(n)+k=1∑p​apk​x(n−k)]x(n−l)+n=0∑N−1−p​[x(n)+k=1∑p​apk​x(n+k)]x(n+l)}=0其中，$l=1,2,...,p$。经过简化后，有：$$\sum _{k=1}^p{a}_{pk}[\sum _{n=p}^{N-1}x(n-k)x(n-l)+\sum _{n=0}^{N-1-p}x(n+k)x(n+l)]=-[\sum _{n=p}^{N-1}x(n)x(n-l)+\sum _{n=0}^{N-1-p}x(n)x(n+l)]$$则上式可以表示为：$$\sum _{k=1}^p{a}_{pk}{\hat{c}}_{xx}(l,k)=-{\hat{c}}_{xx}(l,0)$$其中： c ^ x x ( l , k ) = { 1 2 ( N − p ) ∑ n = p N − 1 x ( n − k ) x ( n − l ) + ∑ n = 0 N − 1 − p x ( n + k ) x ( n + l ) c ^ x x ( k , l ) \hat{c}_{xx}(l,k)=

c^xx​(l,k)=⎩⎪⎪⎨⎪⎪⎧​2(N−p)1​n=p∑N−1​x(n−k)x(n−l)+n=0∑N−1−p​x(n+k)x(n+l)c^xx​(k,l)​写成矩阵形式： [ c ^ x x ( 1 , 1 ) c ^ x x ( 1 , 2 ) ⋯ c ^ x x ( 1 , p ) c ^ x x ( 2 , 1 ) c ^ x x ( 2 , 2 ) ⋯ c ^ x x ( 2 , p ) ⋮ ⋮ ⋱ ⋮ c ^ x x ( p , 1 ) c ^ x x ( p , 2 ) ⋯ c ^ x x ( p , p ) ] [ a p 1 a p 2 ⋮ a p p ] = − [ c ^ x x ( 1 , 0 ) c ^ x x ( 2 , 0 ) ⋮ c ^ x x ( p , 0 ) ]

[c^xx(1,1)c^xx(1,2)⋯c^xx(1,p)c^xx(2,1)c^xx(2,2)⋯c^xx(2,p)⋮⋮⋱⋮c^xx(p,1)c^xx(p,2)⋯c^xx(p,p)]" role="presentation">⎡⎣⎢⎢⎢⎢⎢c^xx(1,1)c^xx(2,1)⋮c^xx(p,1)c^xx(1,2)c^xx(2,2)⋮c^xx(p,2)⋯⋯⋱⋯c^xx(1,p)c^xx(2,p)⋮c^xx(p,p)⎤⎦⎥⎥⎥⎥⎥$$\left[\begin{array}{cccc}{\hat{c}}_{xx}(1,1)& {\hat{c}}_{xx}(1,2)& \cdots & {\hat{c}}_{xx}(1,p)\\ {\hat{c}}_{xx}(2,1)& {\hat{c}}_{xx}(2,2)& \cdots & {\hat{c}}_{xx}(2,p)\\ ⋮& ⋮& \ddots & ⋮\\ {\hat{c}}_{xx}(p,1)& {\hat{c}}_{xx}(p,2)& \cdots & {\hat{c}}_{xx}(p,p)\end{array}\right]$$

[ap1ap2⋮app]" role="presentation">⎡⎣⎢⎢⎢⎢⎢ap1ap2⋮app⎤⎦⎥⎥⎥⎥⎥$$\left[\begin{array}{c}{a}_{p1}\\ a_{p2}\\ ⋮\\ a_{pp}\end{array}\right]$$

=-

[c^xx(1,0)c^xx(2,0)⋮c^xx(p,0)]" role="presentation">⎡⎣⎢⎢⎢⎢⎢c^xx(1,0)c^xx(2,0)⋮c^xx(p,0)⎤⎦⎥⎥⎥⎥⎥$$\left[\begin{array}{c}{\hat{c}}_{xx}(1,0)\\ {\hat{c}}_{xx}(2,0)\\ ⋮\\ {\hat{c}}_{xx}(p,0)\end{array}\right]$$

⎣⎢⎢⎢⎡​c^xx​(1,1)c^xx​(2,1)⋮c^xx​(p,1)​c^xx​(1,2)c^xx​(2,2)⋮c^xx​(p,2)​⋯⋯⋱⋯​c^xx​(1,p)c^xx​(2,p)⋮c^xx​(p,p)​⎦⎥⎥⎥⎤​⎣⎢⎢⎢⎡​ap1​ap2​⋮app​​⎦⎥⎥⎥⎤​=−⎣⎢⎢⎢⎡​c^xx​(1,0)c^xx​(2,0)⋮c^xx​(p,0)​⎦⎥⎥⎥⎤​则模型的参数 $a_{pk}$ 估计为： [ a p 1 a p 2 ⋮ a p p ] = [ c ^ x x ( 1 , 1 ) c ^ x x ( 1 , 2 ) ⋯ c ^ x x ( 1 , p ) c ^ x x ( 2 , 1 ) c ^ x x ( 2 , 2 ) ⋯ c ^ x x ( 2 , p ) ⋮ ⋮ ⋱ ⋮ c ^ x x ( p , 1 ) c ^ x x ( p , 2 ) ⋯ c ^ x x ( p , p ) ] − 1 [ − c ^ x x ( 1 , 0 ) − c ^ x x ( 2 , 0 ) ⋮ − c ^ x x ( p , 0 ) ]

[ap1ap2⋮app]" role="presentation">⎡⎣⎢⎢⎢⎢⎢ap1ap2⋮app⎤⎦⎥⎥⎥⎥⎥$$\left[\begin{array}{c}{a}_{p1}\\ a_{p2}\\ ⋮\\ a_{pp}\end{array}\right]$$

=

[c^xx(1,1)c^xx(1,2)⋯c^xx(1,p)c^xx(2,1)c^xx(2,2)⋯c^xx(2,p)⋮⋮⋱⋮c^xx(p,1)c^xx(p,2)⋯c^xx(p,p)]" role="presentation">⎡⎣⎢⎢⎢⎢⎢c^xx(1,1)c^xx(2,1)⋮c^xx(p,1)c^xx(1,2)c^xx(2,2)⋮c^xx(p,2)⋯⋯⋱⋯c^xx(1,p)c^xx(2,p)⋮c^xx(p,p)⎤⎦⎥⎥⎥⎥⎥$$\left[\begin{array}{cccc}{\hat{c}}_{xx}(1,1)& {\hat{c}}_{xx}(1,2)& \cdots & {\hat{c}}_{xx}(1,p)\\ {\hat{c}}_{xx}(2,1)& {\hat{c}}_{xx}(2,2)& \cdots & {\hat{c}}_{xx}(2,p)\\ ⋮& ⋮& \ddots & ⋮\\ {\hat{c}}_{xx}(p,1)& {\hat{c}}_{xx}(p,2)& \cdots & {\hat{c}}_{xx}(p,p)\end{array}\right]$$

^{-1}

[−c^xx(1,0)−c^xx(2,0)⋮−c^xx(p,0)]" role="presentation">⎡⎣⎢⎢⎢⎢⎢−c^xx(1,0)−c^xx(2,0)⋮−c^xx(p,0)⎤⎦⎥⎥⎥⎥⎥$$\left[\begin{array}{c}-{\hat{c}}_{xx}(1,0)\\ -{\hat{c}}_{xx}(2,0)\\ ⋮\\ -{\hat{c}}_{xx}(p,0)\end{array}\right]$$

⎣⎢⎢⎢⎡​ap1​ap2​⋮app​​⎦⎥⎥⎥⎤​=⎣⎢⎢⎢⎡​c^xx​(1,1)c^xx​(2,1)⋮c^xx​(p,1)​c^xx​(1,2)c^xx​(2,2)⋮c^xx​(p,2)​⋯⋯⋱⋯​c^xx​(1,p)c^xx​(2,p)⋮c^xx​(p,p)​⎦⎥⎥⎥⎤​−1⎣⎢⎢⎢⎡​−c^xx​(1,0)−c^xx​(2,0)⋮−c^xx​(p,0)​⎦⎥⎥⎥⎤​白噪声的方差 ${\sigma }_p^2$ 估计为： σ p 2 = ρ m i n = 1 2 ( N − p ) { ∑ n = p N − 1 [ x ( n ) + ∑ k = 1 p a p k x ( n − k ) ] [ x ( n ) + ∑ l = 1 p a p k x ( n − l ) ] + ∑ n = 0 N − 1 − p [ x ( n ) + ∑ k = 1 p a p k x ( n + k ) ] [ x ( n ) + ∑ l = 1 p a p k x ( n + l ) ] } = 1 2 ( N − p ) { ∑ n = p N − 1 [ x ( n ) + ∑ k = 1 p a p k x ( n − k ) ] x ( n ) + ∑ n = 0 N − 1 − p [ x ( n ) + ∑ k = 1 p a p k x ( n + k ) ] x ( n ) } = 1 2 ( N − p ) [ ∑ n = p N − 1 x 2 ( n ) + ∑ n = 0 N − 1 − p x 2 ( n ) ] + ∑ k = 1 p a p k { 1 2 ( N − p ) [ ∑ n = p N − 1 x ( n ) x ( n − k ) + ∑ n = 0 N − 1 − p x ( n ) x ( n + k ) ] } = c ^ x x ( 0 , 0 ) + ∑ k = 1 p a p k c ^ x x ( 0 , k )

σp2=ρmin=12(N−p){∑n=pN−1[x(n)+∑k=1papkx(n−k)][x(n)+∑l=1papkx(n−l)]+∑n=0N−1−p[x(n)+∑k=1papkx(n+k)][x(n)+∑l=1papkx(n+l)]}=12(N−p){∑n=pN−1[x(n)+∑k=1papkx(n−k)]x(n)+∑n=0N−1−p[x(n)+∑k=1papkx(n+k)]x(n)}=12(N−p)[∑n=pN−1x2(n)+∑n=0N−1−px2(n)]+∑k=1papk{12(N−p)[∑n=pN−1x(n)x(n−k)+∑n=0N−1−px(n)x(n+k)]}=c^xx(0,0)+∑k=1papkc^xx(0,k)" role="presentation">σ2p=ρmin=12(N−p){∑n=pN−1[x(n)+∑k=1papkx(n−k)][x(n)+∑l=1papkx(n−l)]+∑n=0N−1−p[x(n)+∑k=1papkx(n+k)][x(n)+∑l=1papkx(n+l)]}=12(N−p){∑n=pN−1[x(n)+∑k=1papkx(n−k)]x(n)+∑n=0N−1−p[x(n)+∑k=1papkx(n+k)]x(n)}=12(N−p)[∑n=pN−1x2(n)+∑n=0N−1−px2(n)]+∑k=1papk{12(N−p)[∑n=pN−1x(n)x(n−k)+∑n=0N−1−px(n)x(n+k)]}=c^xx(0,0)+∑k=1papkc^xx(0,k)$$\begin{array}{rl}{\sigma }_p^2& ={\rho }_{min}\\ & =\frac{1}{2(N-p)}\{\sum _{n=p}^{N-1}[x(n)+\sum _{k=1}^p{a}_{pk}x(n-k)][x(n)+\sum _{l=1}^p{a}_{pk}x(n-l)]+\sum _{n=0}^{N-1-p}[x(n)+\sum _{k=1}^p{a}_{pk}x(n+k)][x(n)+\sum _{l=1}^p{a}_{pk}x(n+l)]\}\\ & =\frac{1}{2(N-p)}\{\sum _{n=p}^{N-1}[x(n)+\sum _{k=1}^p{a}_{pk}x(n-k)]x(n)+\sum _{n=0}^{N-1-p}[x(n)+\sum _{k=1}^p{a}_{pk}x(n+k)]x(n)\}\\ & =\frac{1}{2(N-p)}[\sum _{n=p}^{N-1}{x}^2(n)+\sum _{n=0}^{N-1-p}{x}^2(n)]+\sum _{k=1}^p{a}_{pk}\{\frac{1}{2(N-p)}[\sum _{n=p}^{N-1}x(n)x(n-k)+\sum _{n=0}^{N-1-p}x(n)x(n+k)]\}\\ & ={\hat{c}}_{xx}(0,0)+\sum _{k=1}^p{a}_{pk}{\hat{c}}_{xx}(0,k)\end{array}$$

σp2​​=ρmin​=2(N−p)1​{n=p∑N−1​[x(n)+k=1∑p​apk​x(n−k)][x(n)+l=1∑p​apk​x(n−l)]+n=0∑N−1−p​[x(n)+k=1∑p​apk​x(n+k)][x(n)+l=1∑p​apk​x(n+l)]}=2(N−p)1​{n=p∑N−1​[x(n)+k=1∑p​apk​x(n−k)]x(n)+n=0∑N−1−p​[x(n)+k=1∑p​apk​x(n+k)]x(n)}=2(N−p)1​[n=p∑N−1​x2(n)+n=0∑N−1−p​x2(n)]+k=1∑p​apk​{2(N−p)1​[n=p∑N−1​x(n)x(n−k)+n=0∑N−1−p​x(n)x(n+k)]}=c^xx​(0,0)+k=1∑p​apk​c^xx​(0,k)​合并成一个矩阵来求模型参数 $a_{pk}$ 和 ${\sigma }_p^2$，则有： [ c ^ x x ( 0 , 0 ) c ^ x x ( 0 , 1 ) ⋯ c ^ x x ( 0 , p ) c ^ x x ( 1 , 0 ) c ^ x x ( 1 , 1 ) ⋯ c ^ x x ( 1 , p ) ⋮ ⋮ ⋱ ⋮ c ^ x x ( p , 0 ) c ^ x x ( p , 1 ) ⋯ c ^ x x ( p , p ) ] [ 1 a p 1 ⋮ a p p ] = − [ σ p 2 0 ⋮ 0 ]

[c^xx(0,0)c^xx(0,1)⋯c^xx(0,p)c^xx(1,0)c^xx(1,1)⋯c^xx(1,p)⋮⋮⋱⋮c^xx(p,0)c^xx(p,1)⋯c^xx(p,p)]" role="presentation">⎡⎣⎢⎢⎢⎢⎢c^xx(0,0)c^xx(1,0)⋮c^xx(p,0)c^xx(0,1)c^xx(1,1)⋮c^xx(p,1)⋯⋯⋱⋯c^xx(0,p)c^xx(1,p)⋮c^xx(p,p)⎤⎦⎥⎥⎥⎥⎥$$\left[\begin{array}{cccc}{\hat{c}}_{xx}(0,0)& {\hat{c}}_{xx}(0,1)& \cdots & {\hat{c}}_{xx}(0,p)\\ {\hat{c}}_{xx}(1,0)& {\hat{c}}_{xx}(1,1)& \cdots & {\hat{c}}_{xx}(1,p)\\ ⋮& ⋮& \ddots & ⋮\\ {\hat{c}}_{xx}(p,0)& {\hat{c}}_{xx}(p,1)& \cdots & {\hat{c}}_{xx}(p,p)\end{array}\right]$$

[1ap1⋮app]" role="presentation">⎡⎣⎢⎢⎢⎢⎢1ap1⋮app⎤⎦⎥⎥⎥⎥⎥$$\left[\begin{array}{c}1\\ a_{p1}\\ ⋮\\ a_{pp}\end{array}\right]$$

=-

[σp20⋮0]" role="presentation">⎡⎣⎢⎢⎢⎢σ2p0⋮0⎤⎦⎥⎥⎥⎥$$\left[\begin{array}{c}{\sigma }_p^2\\ 0\\ ⋮\\ 0\end{array}\right]$$

⎣⎢⎢⎢⎡​c^xx​(0,0)c^xx​(1,0)⋮c^xx​(p,0)​c^xx​(0,1)c^xx​(1,1)⋮c^xx​(p,1)​⋯⋯⋱⋯​c^xx​(0,p)c^xx​(1,p)⋮c^xx​(p,p)​⎦⎥⎥⎥⎤​⎣⎢⎢⎢⎡​1ap1​⋮app​​⎦⎥⎥⎥⎤​=−⎣⎢⎢⎢⎡​σp2​0⋮0​⎦⎥⎥⎥⎤​

自回归模型的功率谱估计

当求出 $p$ 阶模型的参数 $a_k$ 和白噪声方差 ${\sigma }_p^2$，就可以根据这些已知条件来计算模型的功率谱 $P_{xx}(e^{jw})$： P x x ( e j w ) = σ p 2 ∣ H ( e j w ) ∣ 2 = σ p 2 ∣ 1 1 + ∑ k = 1 p a k e − j w k ∣ 2 P_{xx}(e^{jw})=\sigma^2_p|H(e^{jw})|^2=\sigma^2_p

^2 Pxx​(ejw)=σp2​∣H(ejw)∣2=σp2​∣∣∣∣∣∣​1+k=1∑p​ak​e−jwk1​​∣∣∣∣∣∣​2