递推与线性预测 Recursive Linear Prediction
1. 线性预测引入
线性预测是要解决这样的问题。假设我们有n个随机变量，事实上，如果我们把随机变量的角标看做是时间，这就是一个随机过程，否则的话，就是一组随机变量
Z 1 , Z 2 , . . . Z n ， Z n + 1 Z_1,Z_2,...Z_n，Z_{n+1} Z1,Z2,...Zn，Zn+1
现在我们希望能够利用前面的n个变量，通过线性的方法，去估计第n+1个变量，这就是线性预测问题
Z
=
(
Z
1
,
Z
2
,
.
.
.
Z
n
)
T
=
>
Z
n
+
1
Z=(Z_1,Z_2,...Z_n)^T => Z_{n+1}
Z=(Z1,Z2,...Zn)T=>Zn+1
事实上，研究随机过程，很大程度上就是为了进行线性预测。
我们定义一个参数α
α
=
(
α
1
,
.
.
.
,
α
n
)
T
\alpha = (\alpha_1,...,\alpha_n)^T
α=(α1,...,αn)T
我们通过对Z进行线性组合来估计Zn+1
α T Z = ∑ k = 1 n α k Z n + 1 − k → Z n + 1 \alpha^TZ = \sum_{k=1}^n \alpha_k Z_{n+1-k} \rightarrow Z_{n+1} αTZ=k=1∑nαkZn+1−k→Zn+1
我们的优化条件这样的
m i n E ( Z n + 1 − ∑ k = 1 n α k Z n + 1 − k ) 2 minE(Z_{n+1}-\sum_{k=1}^n \alpha_k Z_{n+1-k})^2 minE(Zn+1−k=1∑nαkZn+1−k)2
注意这里只是将Zk的角标进行了倒转，不影响问题的本质
我们之前学过了正交性原理，这里的最优估计，就等价于所有的原材料与残差之间都是正交的
m i n E ( Z n + 1 − ∑ k = 1 n α k Z n + 1 − k ) 2 = > E ( Z n + 1 − k ∗ ( Z n + 1 − ∑ k = 1 n α k Z n + 1 − k ) ) = 0 k = 1 , 2 , . . . , n minE(Z_{n+1}-\sum_{k=1}^n \alpha_k Z_{n+1-k})^2 \\ =>E(Z_{n+1-k}*(Z_{n+1}-\sum_{k=1}^n \alpha_k Z_{n+1-k})) = 0 \\ k = 1,2,...,n minE(Zn+1−k=1∑nαkZn+1−k)2=>E(Zn+1−k∗(Zn+1−k=1∑nαkZn+1−k))=0k=1,2,...,n
因此，我们就有
= > E ( Z n + 1 − k ∗ Z n + 1 ) − ∑ i = 1 n α i E ( Z n + 1 − i ∗ Z n + 1 − k ) = 0 =>E(Z_{n+1-k}*Z_{n+1})-\sum_{i=1}^n \alpha_i E(Z_{n+1-i}*Z_{n+1-k}) = 0 =>E(Zn+1−k∗Zn+1)−i=1∑nαiE(Zn+1−i∗Zn+1−k)=0
这里我们假设Z是个平稳过程，因此Z的自相关就仅仅与下标的差(时间差)有关系
= > R Z ( − k ) − ∑ i = 1 n α i R Z ( − i + k ) = 0 => R_Z(-k)- \sum_{i=1}^n \alpha_i R_Z(-i+k) =0 =>RZ(−k)−i=1∑nαiRZ(−i+k)=0
因为相关函数是对称的
R Z ( − k ) = R Z ( k ) R_Z(-k) = R_Z(k) RZ(−k)=RZ(k)
所以有
R
Z
(
k
)
−
∑
i
=
1
n
α
i
R
Z
(
k
−
i
)
=
0
(
k
=
1
,
.
.
.
,
n
)
R_Z(k)- \sum_{i=1}^n \alpha_i R_Z(k-i) =0 \quad (k=1,...,n)
RZ(k)−i=1∑nαiRZ(k−i)=0(k=1,...,n)
这个方程就Wiener-Hopf方程
我们来表示一下这个方程
( R ( 0 ) R ( 1 ) . . . R ( n − 1 ) R ( 1 ) R ( 0 ) . . . R ( n − 2 ) . . . . . . . . . . . . R ( n − 1 ) R ( n − 2 ) . . . R ( 0 ) ) ∗ ( α 1 . . . . . . α n ) = ( R ( 1 ) . . . . . . R ( n ) ) ( i ) \begin{pmatrix} R(0) & R(1) &...&R(n-1) \\ R(1) & R(0) &...&R(n-2) \\ ... & ... &... & ... \\ R(n-1)& R(n-2)& ... & R(0) \end{pmatrix} * \begin{pmatrix} \alpha_1 \\ ... \\ ... \\ \alpha_n \end{pmatrix} = \begin{pmatrix} R(1) \\ ... \\ ... \\ R(n) \end{pmatrix} \quad\quad\quad(i) ⎝⎜⎜⎛R(0)R(1)...R(n−1)R(1)R(0)...R(n−2)............R(n−1)R(n−2)...R(0)⎠⎟⎟⎞∗⎝⎜⎜⎛α1......αn⎠⎟⎟⎞=⎝⎜⎜⎛R(1)......R(n)⎠⎟⎟⎞(i)
2. 线性预测与LPC编码
其实线性预测问题与我们是息息相关的。比如在语音信号的编码上。
可能我们会感到很奇怪，因为语音信号并不是个平稳的信号，如何使用线性预测的方法进行计算呢?
事实上，虽然长时间的语音信号不具备平稳性，但是短时间还是可以的。人们经过研究，发现如果在时间轴上把语音信号分成不同的段，每段数据200-300ms，在每一段时间内，语音信号是近似平稳的。
那么这一段200-300ms的数据大约有多长呢?
人声音频率大概是200-3400HZ，奈奎斯特之后是400-6800，再加入冗余，语音采样率大概是8k/s，因此一段数据大概有1600个左右的样本数据。有了线性预测之后，我们在进行语音信号传递的时候，就不需要发送这1600个样本数据了，我们可以发送20-30个系数，然后再传一段20-30个数字长度的引子即可，就可以使用前面的数据对后面需要的数据进行递推了。
这种语音编码方式有极高的压缩比，这种编码技术叫做线性预测编码，也就是LPC。事实上,我们每次打电话都在用线性预测编码。这种技术现在已经内置于每一支手机上了。这是成熟技术。
3. 递归与线性预测求解
3.1 问题分析
我们回过来看看方程式(i)，事实上，对这个方程式的求解具有里程碑式的意义。
我们思考一下线性方程式的求解有什么方法，其中最简单的方法就是高斯-若尔当消元法。但是高斯-若尔当消元法时间复杂度有O(n3)
如果我们要把算法落地到计算机上，这个复杂度是非常高的。想想如果我们要使用LPC编码，就必须要求解这个方程，但是一般手机处理性能并没有那么高，算法的时间复杂度问题，成了求解线性预测问题的拦路虎
我们如果降低线性方程组运算的复杂度。这是一个领域，叫做数值线性代数。很多人在研究这个问题。
人们研究的关键都是矩阵的结构，人们希望在高斯消去的基础上往前走。某种意义上，如果我们不对矩阵结构做特殊约定，高斯消元法就是效率最高的算法了。
但是一旦矩阵有特殊性质了，比如对称、正交等，线性方程组就会有更加高效的算法。
3.2 Toeplitz矩阵
3.2.1 结构
( R ( 0 ) R ( 1 ) . . . R ( n − 1 ) R ( 1 ) R ( 0 ) . . . R ( n − 2 ) . . . . . . . . . . . . R ( n − 1 ) R ( n − 2 ) . . . R ( 0 ) ) \begin{pmatrix} R(0) & R(1) &...&R(n-1) \\ R(1) & R(0) &...&R(n-2) \\ ... & ... &... & ... \\ R(n-1)& R(n-2)& ... & R(0) \end{pmatrix} ⎝⎜⎜⎛R(0)R(1)...R(n−1)R(1)R(0)...R(n−2)............R(n−1)R(n−2)...R(0)⎠⎟⎟⎞
仔细观察我们的矩阵，我们发现它是对称的，如果矩阵有对称性，可以把数量级降到平方O(n2)，但是还是不够。对于求解计算还是太复杂，我们希望要把它降到O(n)，这是因为这个矩阵除了对称性，还有其他的结构。
因为这个矩阵，看起了是nxn，实际上只有n个数字，并且行与行之间是循环挪动的结构,这样的矩阵叫Toeplitz矩阵。
正是因为我们能够对这个矩阵进行高效处理，才能让语音线性编码技术走入千家万户。
3.2.2 性质
Toeplitz矩阵具有这样的性质：如果我们把Toeplitz矩阵的线性方程组的变量逆转一个方向，我们得到方程的解也会逆转一个方向
T ∈ R n x n ， T is Toeplistz If T ∗ x = y = > T ∗ x ~ = y ~ T \in R^{nxn} ，\text{T is Toeplistz} \\ \text{If } T*x = y \\ =>T*\widetilde{x} = \widetilde{y} T∈Rnxn，T is ToeplistzIf T∗x=y=>T∗x =y
其中
x = ( x 1 , x 2 , . . . , x n ) T y = ( y 1 , y 2 , . . . , y n ) T x ~ = ( x n , . . . , x 1 ) T y ~ = ( y n , . . . , y 1 ) T x = (x_1,x_2,...,x_n)^T\\ y = (y_1,y_2,...,y_n)^T \\ \widetilde{x} = (x_n,...,x_1)^T \\ \widetilde{y} = (y_n,...,y_1)^T x=(x1,x2,...,xn)Ty=(y1,y2,...,yn)Tx =(xn,...,x1)Ty =(yn,...,y1)T
3.2.3 证明
我们来证明一下
假设
y i = ∑ j = 1 n T i j ∗ x j , ( i = 1 , . . . , n ) ( i i ) y_i = \sum_{j=1}^n T_{ij}*x_j ,(i=1,...,n) \quad\quad(ii) yi=j=1∑nTij∗xj,(i=1,...,n)(ii)
y n + 1 − i = ∑ j = 1 n T i j ∗ x n + 1 − j , ( i = 1 , . . . , n ) ( i i i ) y_{n+1-i} = \sum_{j=1}^n T_{ij}*x_{n+1-j} ,(i=1,...,n) \quad\quad(iii) yn+1−i=j=1∑nTij∗xn+1−j,(i=1,...,n)(iii)
我们对(ii)式子进行变量代换，进行逆向表示
y n + 1 − i = ∑ j = 1 n T ( n + 1 − i , n + 1 − j ) ∗ x n + 1 − j , ( i = 1 , . . . , n ) ( i v ) y_{n+1-i} = \sum_{j=1}^n T_{(n+1-i,n+1-j)}*x_{n+1-j} ,(i=1,...,n) \quad\quad(iv) yn+1−i=j=1∑nT(n+1−i,n+1−j)∗xn+1−j,(i=1,...,n)(iv)
对比(iv)与(iii)，也就是我们需要证明
T ( n + 1 − i , n + 1 − j ) = T i j ( v ) T_{(n+1-i,n+1-j)} = T_{ij} \quad\quad(v) T(n+1−i,n+1−j)=Tij(v)
Toeplitz矩阵具有Tij仅仅依赖于|i-j|的特性
T i j = T ∣ i − j ∣ T_{ij} = T_{|i-j|} Tij=T∣i−j∣
因此显然(v)式成立
3.3 方程求解
3.3.1 递推思路
现在我们希望能利用O(n)的时间复杂度去求解线性方程组，我们可以利用递推的思想。
假设我们的方程是n阶的。如果我们能拿到n-1阶方程的解，然后我们可以利用简单的计算得到n阶的方程，然后我们就可以不断的递推，直到找到1阶方程的解。
α ( n ) < = α ( n − 1 ) < = α ( n − 2 ) < = . . . < = α ( 2 ) α ( 1 ) \alpha^{(n)} <= \alpha^{(n-1)} <= \alpha^{(n-2)} <= ... <= \alpha^{(2)} \alpha^{(1)} α(n)<=α(n−1)<=α(n−2)<=...<=α(2)α(1)
如果这种做法成立，我们就需要解1阶方程就行
( R ( 0 ) R ( 1 ) . . . R ( n − 1 ) R ( 1 ) R ( 0 ) . . . R ( n − 2 ) . . . . . . . . . . . . R ( n − 1 ) R ( n − 2 ) . . . R ( 0 ) ) ∗ ( α 1 . . . . . . α n ) = ( R ( 1 ) . . . . . . R ( n ) ) \begin{pmatrix} R(0) & R(1) &...&R(n-1) \\ R(1) & R(0) &...&R(n-2) \\ ... & ... &... & ... \\ R(n-1)& R(n-2)& ... & R(0) \end{pmatrix} * \begin{pmatrix} \alpha_1 \\ ... \\ ... \\ \alpha_n \end{pmatrix} = \begin{pmatrix} R(1) \\ ... \\ ... \\ R(n) \end{pmatrix} ⎝⎜⎜⎛R(0)R(1)...R(n−1)R(1)R(0)...R(n−2)............R(n−1)R(n−2)...R(0)⎠⎟⎟⎞∗⎝⎜⎜⎛α1......αn⎠⎟⎟⎞=⎝⎜⎜⎛R(1)......R(n)⎠⎟⎟⎞
当n=1时
n = 1 = > R ( 0 ) α ( 1 ) = R ( 1 ) α ( 1 ) = R ( 1 ) R ( 0 ) n=1 => R(0) \alpha^{(1)} = R(1) \\ \alpha^{(1)} = \frac{ R(1)}{ R(0)} n=1=>R(0)α(1)=R(1)α(1)=R(0)R(1)
这里进行符号解释
α k ( n ) \alpha^{(n)}_{k} αk(n)
表示的是n阶线性方程组的αk的解
3.3.2 寻找递推关系
3.3.2.1 augmented normal equation
下面我们要找到n阶和n+1阶方程之间的关系
我们可以对方程进行如下的变形
( R ( 0 ) R ( 1 ) . . . R ( n − 1 ) R ( 1 ) R ( 0 ) . . . R ( n − 2 ) . . . . . . . . . . . . R ( n − 1 ) R ( n − 2 ) . . . R ( 0 ) ) ∗ ( α 1 ( n ) . . . . . . α n ( n ) ) − ( R ( 1 ) . . . . . . R ( n ) ) = 0 \begin{pmatrix} R(0) & R(1) &...&R(n-1) \\ R(1) & R(0) &...&R(n-2) \\ ... & ... &... & ... \\ R(n-1)& R(n-2)& ... & R(0) \end{pmatrix} * \begin{pmatrix} \alpha^{(n)}_1 \\ ... \\ ... \\ \alpha^{(n)}_n \end{pmatrix} - \begin{pmatrix} R(1) \\ ... \\ ... \\ R(n) \end{pmatrix} = 0 ⎝⎜⎜⎛R(0)R(1)...R(n−1)R(1)R(0)...R(n−2)............R(n−1)R(n−2)...R(0)⎠⎟⎟⎞∗⎝⎜⎜⎛α1(n)......αn(n)⎠⎟⎟⎞−⎝⎜⎜⎛R(1)......R(n)⎠⎟⎟⎞=0
把后面的向量提取进前面的矩阵中，变成了n*(n+1)的矩阵
= > ( R ( 1 ) R ( 0 ) . . . R ( n − 1 ) . . . . . . . . . . . . R ( n ) R ( n − 1 ) . . . R ( 0 ) ) ∗ ( − 1 α 1 ( n ) . . . α n ( n ) ) = 0 => \begin{pmatrix} R(1) & R(0) &...&R(n-1) \\ ... & ... &...&... \\ R(n) & R(n-1) &... & R(0) \\ \end{pmatrix} * \begin{pmatrix} -1 \\ \alpha^{(n)}_1 \\ ... \\ \alpha^{(n)}_n \end{pmatrix} = 0 =>⎝⎛R(1)...R(n)R(0)...R(n−1).........R(n−1)...R(0)⎠⎞∗⎝⎜⎜⎛−1α1(n)...αn(n)⎠⎟⎟⎞=0
这一步叫做augmented normal equation (修正过的标准方程）
3.3.2.2 n阶矩阵补全(n+1)阶矩阵
得到修正过的标准方程还是有点问题，因为他是n*(n+1)的矩阵，我们希望把它补全为(n+1)的矩阵,我们发现这个方程与(n+1)阶的Toeplitz方程第一行就差一个{R(0),…,R(n)}，我们就补上这个行向量
( R ( 0 ) R ( 1 ) . . . R ( n ) R ( 1 ) R ( 0 ) . . . R ( n − 1 ) . . . . . . . . . . . . R ( n ) R ( n − 1 ) . . . R ( 0 ) ) ∗ ( − 1 α 1 ( n ) . . . α n ( n ) ) \begin{pmatrix} R(0) & R(1) &...&R(n) \\ R(1) & R(0) &...&R(n-1) \\ ... & ... &...&... \\ R(n) & R(n-1) &... & R(0) \\ \end{pmatrix} * \begin{pmatrix} -1 \\ \alpha^{(n)}_1 \\ ... \\ \alpha^{(n)}_n \end{pmatrix} ⎝⎜⎜⎛R(0)R(1)...R(n)R(1)R(0)...R(n−1)............R(n)R(n−1)...R(0)⎠⎟⎟⎞∗⎝⎜⎜⎛−1α1(n)...αn(n)⎠⎟⎟⎞
补全之后，我们需要重新审视一下右边的向量，这个乘积就已经不是一个全0的了,发现得到的右边第一行就是预测误差
− R ( 0 ) + ∑ k = 1 n α k ( n ) ∗ R ( k ) -R(0)+\sum_{k=1}^{n}\alpha^{(n)}_k*R(k) −R(0)+k=1∑nαk(n)∗R(k)
我们来解释一下，先把预测误差写出来
E ( Z n + 1 − ∑ k = 1 n α k ( n ) ∗ Z n + 1 − k ) 2 = E [ ( Z n + 1 − ∑ k = 1 n α k ( n ) ∗ Z n + 1 − k ) ∗ ( Z n + 1 − ∑ k = 1 n α k ( n ) ∗ Z n + 1 − k ) ] = E [ Z n + 1 ∗ ( Z n + 1 − ∑ k = 1 n α k ( n ) ∗ Z n + 1 − k ) − ( ∑ k = 1 n α k ( n ) ∗ Z n + 1 − k ) ∗ ( Z n + 1 − ∑ k = 1 n α k ( n ) ∗ Z n + 1 − k ) ] E(Z_{n+1} - \sum_{k=1}^{n}\alpha^{(n)}_k*Z_{n+1-k})^2 \\ = E[(Z_{n+1} - \sum_{k=1}^{n}\alpha^{(n)}_k*Z_{n+1-k})*(Z_{n+1} - \sum_{k=1}^{n}\alpha^{(n)}_k*Z_{n+1-k})] \\ =E[Z_{n+1}*(Z_{n+1} - \sum_{k=1}^{n}\alpha^{(n)}_k*Z_{n+1-k})-(\sum_{k=1}^{n}\alpha^{(n)}_k*Z_{n+1-k})*(Z_{n+1} - \sum_{k=1}^{n}\alpha^{(n)}_k*Z_{n+1-k})] E(Zn+1−k=1∑nαk(n)∗Zn+1−k)2=E[(Zn+1−k=1∑nαk(n)∗Zn+1−k)∗(Zn+1−k=1∑nαk(n)∗Zn+1−k)]=E[Zn+1∗(Zn+1−k=1∑nαk(n)∗Zn+1−k)−(k=1∑nαk(n)∗Zn+1−k)∗(Zn+1−k=1∑nαk(n)∗Zn+1−k)]
根据正交性原理，所有的原材料与残差的内积为0，所以，后一个减号里面的东西就是0，也就是
E ( Z n + 1 − ∑ k = 1 n α k ( n ) ∗ Z n + 1 − k ) 2 = E [ Z n + 1 ∗ ( Z n + 1 − ∑ k = 1 n α k ( n ) ∗ Z n + 1 − k ) ] = R ( 0 ) − ∑ k = 1 n α k ( n ) ∗ R ( k ) E(Z_{n+1} - \sum_{k=1}^{n}\alpha^{(n)}_k*Z_{n+1-k})^2 \\ = E[Z_{n+1}*(Z_{n+1} - \sum_{k=1}^{n}\alpha^{(n)}_k*Z_{n+1-k})] \\ = R(0)-\sum_{k=1}^{n}\alpha^{(n)}_k*R(k) E(Zn+1−k=1∑nαk(n)∗Zn+1−k)2=E[Zn+1∗(Zn+1−k=1∑nαk(n)∗Zn+1−k)]=R(0)−k=1∑nαk(n)∗R(k)
所以，我们用n阶矩阵补全的(n+1)*(n+1)阶矩阵如下
(
R
(
0
)
R
(
1
)
.
.
.
R
(
n
)
R
(
1
)
R
(
0
)
.
.
.
R
(
n
−
1
)
.
.
.
.
.
.
.
.
.
.
.
.
R
(
n
)
R
(
n
−
1
)
.
.
.
R
(
0
)
)
∗
(
−
1
α
1
(
n
)
.
.
.
α
n
(
n
)
)
=
(
−
e
n
0
.
.
.
0
)
\begin{pmatrix} R(0) & R(1) &...&R(n) \\ R(1) & R(0) &...&R(n-1) \\ ... & ... &...&... \\ R(n) & R(n-1) &... & R(0) \\ \end{pmatrix} * \begin{pmatrix} -1 \\ \alpha^{(n)}_1 \\ ... \\ \alpha^{(n)}_n \end{pmatrix} =\begin{pmatrix} -e_n \\ 0 \\ ... \\ 0 \end{pmatrix}
⎝⎜⎜⎛R(0)R(1)...R(n)R(1)R(0)...R(n−1)............R(n)R(n−1)...R(0)⎠⎟⎟⎞∗⎝⎜⎜⎛−1α1(n)...αn(n)⎠⎟⎟⎞=⎝⎜⎜⎛−en0...0⎠⎟⎟⎞
3.3.2.3 (n+1)阶矩阵补全
n阶矩阵补全
( R ( 0 ) . . . R ( n ) . . . . . . . . . R ( n ) . . . R ( 0 ) ) ∗ ( − 1 α 1 ( n ) . . . α n ( n ) ) = ( − e n 0 . . . 0 ) ( 1 ) \begin{pmatrix} R(0) &...&R(n) \\ ... &...&... \\ R(n) &... & R(0) \\ \end{pmatrix} * \begin{pmatrix} -1 \\ \alpha^{(n)}_1 \\ ... \\ \alpha^{(n)}_n \end{pmatrix} =\begin{pmatrix} -e_n \\ 0 \\ ... \\ 0 \end{pmatrix} \quad\quad(1) ⎝⎛R(0)...R(n).........R(n)...R(0)⎠⎞∗⎝⎜⎜⎛−1α1(n)...αn(n)⎠⎟⎟⎞=⎝⎜⎜⎛−en0...0⎠⎟⎟⎞(1)
我们用同样方法对(n+1)阶矩阵补全
( R ( 0 ) . . . R ( n ) R ( n + 1 ) . . . . . . . . . . . . R ( n ) . . . R ( 0 ) R ( 1 ) R ( n + 1 ) . . . R ( 1 ) R ( 0 ) ) ∗ ( − 1 α 1 ( n + 1 ) . . . α n + 1 ( n + 1 ) ) = ( − e n + 1 0 . . . 0 ) ( 2 ) \begin{pmatrix} R(0) &... &R(n) &R(n+1) \\ ... &... &...& ... \\ R(n) & ... &R(0)& R(1) \\ R(n+1) & ... &R(1) & R(0) \\ \end{pmatrix} * \begin{pmatrix} -1 \\ \alpha^{(n+1)}_1 \\ ... \\ \alpha^{(n+1)}_{n+1} \end{pmatrix} =\begin{pmatrix} -e_{n+1} \\ 0 \\ ... \\ 0 \end{pmatrix} \quad\quad(2) ⎝⎜⎜⎛R(0)...R(n)R(n+1)............R(n)...R(0)R(1)R(n+1)...R(1)R(0)⎠⎟⎟⎞∗⎝⎜⎜⎛−1α1(n+1)...αn+1(n+1)⎠⎟⎟⎞=⎝⎜⎜⎛−en+10...0⎠⎟⎟⎞(2)
3.3.2.4 递推
我们发现(n)维补全矩阵就是(n+1)维补全矩阵的左上角，他们之间就差一行一列。
我们基于(n)维矩阵补全的矩阵构造一个与(n+1)维矩矩阵补全的矩阵结构类似的矩阵
(
R
(
0
)
.
.
.
R
(
n
)
R
(
n
+
1
)
.
.
.
.
.
.
.
.
.
.
.
.
R
(
n
)
.
.
.
R
(
0
)
R
(
1
)
R
(
n
+
1
)
.
.
.
R
(
1
)
R
(
0
)
)
∗
(
−
1
α
1
(
n
)
.
.
.
α
n
(
n
)
0
)
=
(
−
e
n
0
.
.
.
.
.
.
A
n
)
(
3
)
\begin{pmatrix} R(0) &... &R(n) &R(n+1) \\ ... &... &...& ... \\ R(n) & ... &R(0)& R(1) \\ R(n+1) & ... &R(1) & R(0) \\ \end{pmatrix} * \begin{pmatrix} -1 \\ \alpha^{(n)}_1 \\ ... \\ \alpha^{(n)}_n\\ 0 \end{pmatrix} =\begin{pmatrix} -e_{n} \\ 0 \\ ...\\ ... \\ A_n \end{pmatrix} \quad\quad(3)
⎝⎜⎜⎛R(0)...R(n)R(n+1)............R(n)...R(0)R(1)R(n+1)...R(1)R(0)⎠⎟⎟⎞∗⎝⎜⎜⎜⎜⎛−1α1(n)...αn(n)0⎠⎟⎟⎟⎟⎞=⎝⎜⎜⎜⎜⎛−en0......An⎠⎟⎟⎟⎟⎞(3)
这里注意，因为(n)维矩阵补全的矩阵是(n+1)*(n+1)的，而(n+1)维矩阵的补全矩阵是(n+2)*(n+2)的,所以α(n)的列向量会少一个维度，为了不让绿色框中的向量对结果产生影响，我们对最后一行补个0，这样的话我们需要在最后面补一个An
A n = − R ( n + 1 ) + ∑ k = 1 n α k ( n ) R ( n + 1 − k ) A_n = -R(n+1) +\sum_{k=1}^n \alpha^{(n)}_k R(n+1-k) An=−R(n+1)+k=1∑nαk(n)R(n+1−k)
我们使用Toeplitz矩阵的性质
T ∈ R n x n ， T is Toeplistz If T ∗ x = y = > T ∗ x ~ = y ~ T \in R^{nxn} ，\text{T is Toeplistz} \\ \text{If } T*x = y \\ =>T*\widetilde{x} = \widetilde{y} T∈Rnxn，T is ToeplistzIf T∗x=y=>T∗x =y
其中
x = ( x 1 , x 2 , . . . , x n ) T y = ( y 1 , y 2 , . . . , y n ) T x ~ = ( x n , . . . , x 1 ) T y ~ = ( y n , . . . , y 1 ) T x = (x_1,x_2,...,x_n)^T\\ y = (y_1,y_2,...,y_n)^T \\ \widetilde{x} = (x_n,...,x_1)^T \\ \widetilde{y} = (y_n,...,y_1)^T x=(x1,x2,...,xn)Ty=(y1,y2,...,yn)Tx =(xn,...,x1)Ty =(yn,...,y1)T
我们对(3)式使用这个性质
( R ( 0 ) . . . R ( n ) R ( n + 1 ) . . . . . . . . . . . . R ( n ) . . . R ( 0 ) R ( 1 ) R ( n + 1 ) . . . R ( 1 ) R ( 0 ) ) ∗ ( 0 α n ( n ) . . . α 1 ( n ) − 1 ) = ( A n 0 . . . . . . − e n ) ( 4 ) \begin{pmatrix} R(0) &... &R(n) &R(n+1) \\ ... &... &...& ... \\ R(n) & ... &R(0)& R(1) \\ R(n+1) & ... &R(1) & R(0) \\ \end{pmatrix} * \begin{pmatrix} 0 \\ \alpha^{(n)}_n \\ ... \\ \alpha^{(n)}_1\\ -1 \end{pmatrix} =\begin{pmatrix} A_n \\ 0 \\ ...\\ ... \\ -e_n \end{pmatrix} \quad\quad(4) ⎝⎜⎜⎛R(0)...R(n)R(n+1)............R(n)...R(0)R(1)R(n+1)...R(1)R(0)⎠⎟⎟⎞∗⎝⎜⎜⎜⎜⎛0αn(n)...α1(n)−1⎠⎟⎟⎟⎟⎞=⎝⎜⎜⎜⎜⎛An0......−en⎠⎟⎟⎟⎟⎞(4)
我们利用(3)(4)两个方程进行线性组合
( 3 ) + k ∗ ( 4 ) (3) + k*(4) (3)+k∗(4)
( R ( 0 ) . . . R ( n ) R ( n + 1 ) . . . . . . . . . . . . R ( n ) . . . R ( 0 ) R ( 1 ) R ( n + 1 ) . . . R ( 1 ) R ( 0 ) ) ∗ ( ( − 1 α 1 ( n ) . . . α n ( n ) 0 ) + k ( 0 α n ( n ) . . . α 1 ( n ) − 1 ) ) = ( − e n + k ∗ A n 0 . . . . . . − k ∗ e n + A n ) ( 5 ) \begin{pmatrix}R(0) &... &R(n) &R(n+1) \\ ... &... &...& ... \\ R(n) & ... &R(0)& R(1) \\ R(n+1) & ... &R(1) & R(0) \\ \end{pmatrix} *(\begin{pmatrix} -1 \\ \alpha^{(n)}_1 \\ ... \\ \alpha^{(n)}_n\\ 0 \end{pmatrix} + k\begin{pmatrix} 0 \\ \alpha^{(n)}_n \\ ... \\ \alpha^{(n)}_1\\ -1 \end{pmatrix} ) =\begin{pmatrix} -e_n+k*A_n \\ 0 \\ ...\\ ... \\ -k*e_n +A_n \end{pmatrix} \quad\quad(5) ⎝⎜⎜⎛R(0)...R(n)R(n+1)............R(n)...R(0)R(1)R(n+1)...R(1)R(0)⎠⎟⎟⎞∗(⎝⎜⎜⎜⎜⎛−1α1(n)...αn(n)0⎠⎟⎟⎟⎟⎞+k⎝⎜⎜⎜⎜⎛0αn(n)...α1(n)−1⎠⎟⎟⎟⎟⎞)=⎝⎜⎜⎜⎜⎛−en+k∗An0......−k∗en+An⎠⎟⎟⎟⎟⎞(5)
令
− k ∗ e n + A n = 0 -k*e_n +A_n = 0 −k∗en+An=0
即
k
=
A
n
e
n
k = \frac{A_n}{e_n}
k=enAn
− e n + k ∗ A n = − ( 1 − k 2 ) ∗ e n ( 6 ) -e_n + k*A_n = -(1-k^2)*e_n \quad\quad(6) −en+k∗An=−(1−k2)∗en(6)
(6)代入(5)中得
(
R
(
0
)
.
.
.
R
(
n
)
R
(
n
+
1
)
.
.
.
.
.
.
.
.
.
.
.
.
R
(
n
)
.
.
.
R
(
0
)
R
(
1
)
R
(
n
+
1
)
.
.
.
R
(
1
)
R
(
0
)
)
∗
(
−
1
α
1
(
n
)
+
k
∗
α
n
(
n
)
.
.
.
α
n
(
n
)
+
k
∗
α
1
(
n
)
−
k
)
=
(
−
(
1
−
k
2
)
∗
e
n
0
.
.
.
.
.
.
0
)
=
(
−
e
n
+
1
0
.
.
.
.
.
.
0
)
(
7
)
\begin{pmatrix} R(0) &... &R(n) &R(n+1) \\ ... &... &...& ... \\ R(n) & ... &R(0)& R(1) \\ R(n+1) & ... &R(1) & R(0) \\ \end{pmatrix} * \begin{pmatrix} -1 \\ \alpha^{(n)}_1 +k*\alpha^{(n)}_n \\ ... \\ \alpha^{(n)}_n +k*\alpha^{(n)}_1\\ -k \end{pmatrix} =\begin{pmatrix} -(1-k^2)*e_n \\ 0 \\ ...\\ ... \\ 0 \end{pmatrix} = \begin{pmatrix} -e_{n+1} \\ 0 \\ ...\\ ... \\ 0 \end{pmatrix} \quad\quad(7)
⎝⎜⎜⎛R(0)...R(n)R(n+1)............R(n)...R(0)R(1)R(n+1)...R(1)R(0)⎠⎟⎟⎞∗⎝⎜⎜⎜⎜⎛−1α1(n)+k∗αn(n)...αn(n)+k∗α1(n)−k⎠⎟⎟⎟⎟⎞=⎝⎜⎜⎜⎜⎛−(1−k2)∗en0......0⎠⎟⎟⎟⎟⎞=⎝⎜⎜⎜⎜⎛−en+10......0⎠⎟⎟⎟⎟⎞(7)
其中
k = A n e n A n = − R ( n + 1 ) + ∑ k = 1 n α k ( n ) R ( n + 1 − k ) e n = R ( 0 ) − ∑ k = 1 n α k ( n ) ∗ R ( k ) k = \frac{A_n}{e_n} \\ A_n = -R(n+1) +\sum_{k=1}^n \alpha^{(n)}_k R(n+1-k) \\ e_n =R(0)-\sum_{k=1}^{n}\alpha^{(n)}_k*R(k) k=enAnAn=−R(n+1)+k=1∑nαk(n)R(n+1−k)en=R(0)−k=1∑nαk(n)∗R(k)
与(n+1)阶补全矩阵比较
( R ( 0 ) . . . R ( n ) R ( n + 1 ) . . . . . . . . . . . . R ( n ) . . . R ( 0 ) R ( 1 ) R ( n + 1 ) . . . R ( 1 ) R ( 0 ) ) ∗ ( − 1 α 1 ( n + 1 ) . . . α n + 1 ( n + 1 ) ) = ( − e n + 1 0 . . . 0 ) ( 2 ) \begin{pmatrix} R(0) &... &R(n) &R(n+1) \\ ... &... &...& ... \\ R(n) & ... &R(0)& R(1) \\ R(n+1) & ... &R(1) & R(0) \\ \end{pmatrix} * \begin{pmatrix} -1 \\ \alpha^{(n+1)}_1 \\ ... \\ \alpha^{(n+1)}_{n+1} \end{pmatrix} =\begin{pmatrix} -e_{n+1} \\ 0 \\ ... \\ 0 \end{pmatrix} \quad\quad(2) ⎝⎜⎜⎛R(0)...R(n)R(n+1)............R(n)...R(0)R(1)R(n+1)...R(1)R(0)⎠⎟⎟⎞∗⎝⎜⎜⎛−1α1(n+1)...αn+1(n+1)⎠⎟⎟⎞=⎝⎜⎜⎛−en+10...0⎠⎟⎟⎞(2)
3.3.2.5 Levinson-Durbin Recursive
我们就得到的递推叫做Levinson-Durbin Recursive，是Levinson1947年的工作，durbin发表与1960年。这篇论文目前只被引用了800次
下面整理一下结论
- n=1
α ( 1 ) = R ( 1 ) R ( 0 ) e 1 = R ( 0 ) − α 1 ( 1 ) R ( 1 ) = R ( 0 ) − R ( 1 ) 2 R ( 0 ) k 1 = − R ( 2 ) + α 1 ( 1 ) R ( 1 ) R ( 0 ) − R ( 1 ) 2 R ( 0 ) \alpha^{(1)} = \frac{R(1)}{R(0)} \\ e_1 = R(0) - \alpha_1^{(1)}R(1) = R(0) - \frac{R(1)^2}{R(0)} \\ k_1 = \frac{-R(2)+\alpha_1^{(1)}R(1)}{ R(0) - \frac{R(1)^2}{R(0)}} α(1)=R(0)R(1)e1=R(0)−α1(1)R(1)=R(0)−R(0)R(1)2k1=R(0)−R(0)R(1)2−R(2)+α1(1)R(1)
- n=i
α ( i ) = ( α 1 ( i ) , . . . , α i ( i ) ) T k i = − R ( i + 1 ) + ∑ j = 1 i α j ( i ) R ( i + 1 − j ) R ( 0 ) − ∑ j = 1 i α j ( i ) ∗ R ( j ) \alpha^{(i)} = (\alpha_1^{(i)},...,\alpha_i^{(i)})^T \\ k_i = \frac{-R(i+1)+\sum_{j=1}^{i}\alpha_j^{(i)}R(i+1-j)}{R(0)-\sum_{j=1}^{i}\alpha^{(i)}_j*R(j)} α(i)=(α1(i),...,αi(i))Tki=R(0)−∑j=1iαj(i)∗R(j)−R(i+1)+∑j=1iαj(i)R(i+1−j)
- n=i+1
α j ( i + 1 ) = α j ( i ) + k i ∗ α i + 1 − j ( i ) j = 1 , 2 , . . . , i α i + 1 ( i + 1 ) = − k i \alpha_j^{(i+1)} = \alpha_j^{(i)} + k_i * \alpha_{i+1-j}^{(i)} \quad\quad j = 1,2,...,i \\ \alpha_{i+1}^{(i+1)} = -k_i αj(i+1)=αj(i)+ki∗αi+1−j(i)j=1,2,...,iαi+1(i+1)=−ki
4. Levinson迭代的探讨
4.1 前向预测与后向预测
我们首先来介绍一下什么是前向预测与后向预测
我们假设有这样的数据(Z1,…,Zn+1)
我们用(Z1,…,Zn)去估计Zn+1被称为前向预测。而如果我们用(Z2,…,Zn+1)去预测Z1那么就被称之为后向预测。
前向预测的目的是做数据预测，而后向预测的目的是做了做数据平滑。因此我们有更多的数据之后，我们可以用后面的数据来对时间更靠前的数据进行估计，估计的目的在于，通过获得的更多的知识用在之前的时间点上，从而消除前面噪声的影响，让我们的数据更加光滑。
F o r w a r d − B a c k w o r d P r e d i c t i o n ( Z 1 , . . . , Z n + 1 ) = > ( Z 1 , . . . , Z n ) = > Z n + 1 ( F o r w a r d ) = > P r e d i c t i o n ( Z 1 , . . . , Z n + 1 ) = > ( Z 2 , . . . , Z n + 1 ) = > Z 1 ( B a c k w a r d ) = > S m o o t h i n g Forward-Backword \quad Prediction \\ (Z_1,...,Z_{n+1}) => (Z_1,...,Z_n) => Z_{n+1} (Forward)=>Prediction \\ (Z_1,...,Z_{n+1}) => (Z_2,...,Z_{n+1}) => Z_{1} (Backward)=>Smoothing \\ Forward−BackwordPrediction(Z1,...,Zn+1)=>(Z1,...,Zn)=>Zn+1(Forward)=>Prediction(Z1,...,Zn+1)=>(Z2,...,Zn+1)=>Z1(Backward)=>Smoothing
4.2 Levinson迭代的本质
我们来继续品味一下Levinson迭代中的翻转操作到底是什么意思。
假设我们要做n阶前向预测，预测Zn+1。这个预测本质上就是投影，把Zn+1投影到前面n个数据中
F o r w a r d : f p ( n ) ( n + 1 ) = P r o j ( Z 1 , . . . , Z n ) Z n + 1 Forward: f_p^{(n)}(n+1) = Proj_{(Z_1,...,Z_n)}Z_{n+1} Forward:fp(n)(n+1)=Proj(Z1,...,Zn)Zn+1
我们希望能够进行递推,通过更加低阶的预测来预测Zn+1。
( Z 2 , . . . , Z n ) → Z n + 1 f p ( n − 1 ) ( n + 1 ) (Z_2,...,Z_n) \rightarrow Z_{n+1} \\ f_p^{(n-1)}(n+1) (Z2,...,Zn)→Zn+1fp(n−1)(n+1)
这里先思考一个问题，同样是n-1阶预测，使用(Z1,…,Zn-1)预测Zn和用(Z2,…,Zn)预测Zn+1的系数有没有区别呢?
( Z 1 , . . . , Z n − 1 ) → Z n = > f p ( n − 1 ) ( n ) ( Z 2 , . . . , Z n ) → Z n + 1 = > f p ( n − 1 ) ( n + 1 ) (Z_1,...,Z_{n-1}) \rightarrow Z_{n} => \quad f_p^{(n-1)}(n) \\ (Z_2,...,Z_n) \rightarrow Z_{n+1} => \quad f_p^{(n-1)}(n+1) \\ (Z1,...,Zn−1)→Zn=>fp(n−1)(n)(Z2,...,Zn)→Zn+1=>fp(n−1)(n+1)
事实上，这两组系数完全相同。因为Z是平稳的，可以不计起点，只考虑预测的阶数。
f p ( n − 1 ) ( n + 1 ) = f p ( n − 1 ) ( n ) \quad f_p^{(n-1)}(n+1) = \quad f_p^{(n-1)}(n) fp(n−1)(n+1)=fp(n−1)(n)
让我们回到n阶预测，我们假设已经知道了(Z2,…,Zn)对Zn+1的预测,也就是n-1阶的前向预测。现在要进行n阶预测
P r o j ( Z 1 , . . . , Z n ) Z n + 1 = P r o j { Z 1 } ∪ { Z 2 , . . . , Z n } Z n + 1 Proj_{(Z_1,...,Z_n)}Z_{n+1} = Proj_{\{Z_1\} \cup \{Z_2,...,Z_n\}}Z_{n+1} Proj(Z1,...,Zn)Zn+1=Proj{Z1}∪{Z2,...,Zn}Zn+1
我们思考一下下面的式子是否成立
P r o j { Z 1 } ∪ { Z 2 , . . . , Z n } Z n + 1 ? = P r o j { Z 1 } Z n + 1 + P r o j { Z 2 , . . . , Z n } Z n + 1 Proj_{\{Z_1\} \cup \{Z_2,...,Z_n\}}Z_{n+1}\quad ?= \quad Proj_{\{Z_1\}}Z_{n+1} + Proj_{\{Z_2,...,Z_n\}}Z_{n+1} Proj{Z1}∪{Z2,...,Zn}Zn+1?=Proj{Z1}Zn+1+Proj{Z2,...,Zn}Zn+1
这就向我们前面说明维纳滤波时候一样，因为Z1与{Z2,…,Zn}不一定是正交的，所以这个结论一般不成立。但是我们可以先做正交化，让Z1y与{Z2,…,Zn}张成的向量空间正交
{ Z 1 , . . . , Z n } = { Z 2 , . . . , Z n } ⊕ { Z 1 − P r o j { Z 2 , . . . , Z n } Z 1 } \{Z_1,...,Z_n\} = \{Z_2,...,Z_n\} \oplus \{ Z_1 - Proj_{\{ Z_2,...,Z_n\}} Z_{1}\} {Z1,...,Zn}={Z2,...,Zn}⊕{Z1−Proj{Z2,...,Zn}Z1}
我们可以看到，后面部分就是后向预测的预测误差。
因此，我们令
e b ( n − 1 ) = Z 1 − P r o j { Z 2 , . . . , Z n } Z 1 e^{(n-1)}_b = Z_1 - Proj_{\{ Z_2,...,Z_n\}} Z_{1} eb(n−1)=Z1−Proj{Z2,...,Zn}Z1
所以，我们知道了，一个高阶的前向预测，就是一个低阶的前向预测和一个低阶的后向预测之间的某种组合，这个就和看Levinson迭代的思想非常相近了
{ Z 1 , . . . , Z n } = { Z 2 , . . . , Z n } ⊕ { Z 1 − P r o j { Z 2 , . . . , Z n } Z 1 } = > { Z 2 , . . . , Z n } ⊕ { e b ( n − 1 ) } \{Z_1,...,Z_n\} = \{Z_2,...,Z_n\} \oplus \{ Z_1 - Proj_{\{ Z_2,...,Z_n\}} Z_{1}\} \\ =>\{Z_2,...,Z_n\} \oplus \{ e^{(n-1)}_b \} {Z1,...,Zn}={Z2,...,Zn}⊕{Z1−Proj{Z2,...,Zn}Z1}=>{Z2,...,Zn}⊕{eb(n−1)}
我们来计算一下
P r o j { Z 1 , . . . , Z n } Z n + 1 = P r o j { Z 2 , . . . , Z n } Z n + 1 + P r o j { e b ( n − 1 ) } Z n + 1 = f p ( n − 1 ) ( n + 1 ) + E ( Z n + 1 ∗ e b ( n − 1 ) ) [ E ( e b ( n − 1 ) ) ] 2 ∗ e b ( n − 1 ) Proj_{\{Z_1,...,Z_n\}}Z_{n+1} = Proj_{\{Z_2,...,Z_n\}}Z_{n+1} + Proj_{\{ e^{(n-1)}_b \}}Z_{n+1} \\ = f_p^{(n-1)}(n+1) + \frac{E(Z_{n+1}* e^{(n-1)}_b)}{[E(e_b^{(n-1)})]^2}* e^{(n-1)}_b Proj{Z1,...,Zn}Zn+1=Proj{Z2,...,Zn}Zn+1+Proj{eb(n−1)}Zn+1=fp(n−1)(n+1)+[E(eb(n−1))]2E(Zn+1∗eb(n−1))∗eb(n−1)
因为后向预测误差等于Z1减去后向预测
e b ( n − 1 ) = Z 1 − b p ( n − 1 ) e^{(n-1)}_b = Z_1 - b_p^{(n-1)} eb(n−1)=Z1−bp(n−1)
同时令k
k = E ( Z n + 1 ∗ e b ( n − 1 ) ) [ E ( e b ( n − 1 ) ) 2 ] [ 1 ] k = \frac{E(Z_{n+1}* e^{(n-1)}_b)}{[E(e_b^{(n-1)})^2]} \quad\quad[1] k=[E(eb(n−1))2]E(Zn+1∗eb(n−1))[1]
则原式可化为
f p ( n − 1 ) ( n + 1 ) + E ( Z n + 1 ∗ e b ( n − 1 ) ) [ E ( e ( n − 1 ) ) ] 2 ∗ e b ( n − 1 ) = f p ( n − 1 ) ( n + 1 ) + k ∗ ( Z 1 − b p ( n − 1 ) ) = f p ( n − 1 ) ( n + 1 ) − k ∗ b p ( n − 1 ) + k ∗ Z 1 f_p^{(n-1)}(n+1) + \frac{E(Z_{n+1}* e^{(n-1)}_b)}{[E(e^{(n-1)})]^2}* e^{(n-1)}_b = f_p^{(n-1)}(n+1) + k*( Z_1 - b_p^{(n-1)}) \\ = f_p^{(n-1)}(n+1) - k*b_p^{(n-1)} + k*Z_1 fp(n−1)(n+1)+[E(e(n−1))]2E(Zn+1∗eb(n−1))∗eb(n−1)=fp(n−1)(n+1)+k∗(Z1−bp(n−1))=fp(n−1)(n+1)−k∗bp(n−1)+k∗Z1
这跟Levinson迭代是一样的
( R ( 0 ) . . . R ( n ) R ( n + 1 ) . . . . . . . . . . . . R ( n ) . . . R ( 0 ) R ( 1 ) R ( n + 1 ) . . . R ( 1 ) R ( 0 ) ) ∗ ( ( − 1 α 1 ( n ) . . . α n ( n ) 0 ) + k ( 0 α n ( n ) . . . α 1 ( n ) − 1 ) ) = ( − e n + k ∗ A n 0 . . . . . . − k ∗ e n + A n ) \begin{pmatrix} R(0) &... &R(n) &R(n+1) \\ ... &... &...& ... \\ R(n) & ... &R(0)& R(1) \\ R(n+1) & ... &R(1) & R(0) \\ \end{pmatrix} * ( \begin{pmatrix} -1 \\ \alpha^{(n)}_1 \\ ... \\ \alpha^{(n)}_n\\ 0 \end{pmatrix} + k\begin{pmatrix} 0 \\ \alpha^{(n)}_n \\ ... \\ \alpha^{(n)}_1\\ -1 \end{pmatrix} ) =\begin{pmatrix} -e_n+k*A_n \\ 0 \\ ...\\ ... \\ -k*e_n +A_n \end{pmatrix} ⎝⎜⎜⎛R(0)...R(n)R(n+1)............R(n)...R(0)R(1)R(n+1)...R(1)R(0)⎠⎟⎟⎞∗(⎝⎜⎜⎜⎜⎛−1α1(n)...αn(n)0⎠⎟⎟⎟⎟⎞+k⎝⎜⎜⎜⎜⎛0αn(n)...α1(n)−1⎠⎟⎟⎟⎟⎞)=⎝⎜⎜⎜⎜⎛−en+k∗An0......−k∗en+An⎠⎟⎟⎟⎟⎞
(
α
1
(
n
+
1
)
,
.
.
.
,
α
n
(
n
+
1
)
,
α
n
+
1
(
n
+
1
)
)
=
(
α
1
(
n
)
+
k
∗
α
n
(
n
)
,
.
.
.
,
α
n
(
n
)
+
k
∗
α
1
(
n
)
,
−
k
)
(\alpha_1^{(n+1)},...,\alpha_n^{(n+1)},\alpha_{n+1}^{(n+1)}) = (\alpha_1^{(n)} + k* \alpha_n^{(n)},...,\alpha_n^{(n)} + k* \alpha_1^{(n)},-k)
(α1(n+1),...,αn(n+1),αn+1(n+1))=(α1(n)+k∗αn(n),...,αn(n)+k∗α1(n),−k)
这里的k跟我们用前向预测+后向预测中的k就差一个符号，这个不本质。
我们来计算一下[1]
k = E ( Z n + 1 ∗ e b ( n − 1 ) ) [ E ( e b ( n − 1 ) ) 2 ] k = \frac{E(Z_{n+1}* e^{(n-1)}_b)}{[E(e_b^{(n-1)})^2]} \\ k=[E(eb(n−1))2]E(Zn+1∗eb(n−1))
计算一下分子和分母
E ( Z n + 1 ∗ e b ( n − 1 ) ) = E ( Z n + 1 ∗ ( Z 1 − ∑ k = 1 n − 1 α k ( n − 1 ) ∗ Z ( k + 1 ) ) = R ( n ) − ∑ k = 1 n − 1 α k ( n − 1 ) ∗ R ( n − k ) E(Z_{n+1}* e^{(n-1)}_b) = E(Z_{n+1}*(Z_1 - \sum_{k=1}^{n-1}\alpha^{(n-1)}_k*Z_{(k+1})) \\ = R(n) - \sum_{k=1}^{n-1}\alpha^{(n-1)}_k*R(n-k) E(Zn+1∗eb(n−1))=E(Zn+1∗(Z1−k=1∑n−1αk(n−1)∗Z(k+1))=R(n)−k=1∑n−1αk(n−1)∗R(n−k)
[ E ( e b ( n − 1 ) ) ] 2 = E [ ( Z 1 − ∑ k = 1 n − 1 α k ( n − 1 ) ∗ Z ( k + 1 ) 2 ] = E [ Z 1 ∗ ( Z 1 − ∑ k = 1 n − 1 α k ( n − 1 ) ∗ Z ( k + 1 ) ] = R ( 0 ) − ∑ k = 1 n − 1 α k ( n − 1 ) ∗ R ( k ) [E(e_b^{(n-1)})]^2 = E[(Z_1 - \sum_{k=1}^{n-1}\alpha^{(n-1)}_k*Z_{(k+1})^2] \\ =E[Z_1*(Z_1 - \sum_{k=1}^{n-1}\alpha^{(n-1)}_k*Z_{(k+1})] = R(0) - \sum_{k=1}^{n-1}\alpha_k^{(n-1)}*R(k) [E(eb(n−1))]2=E[(Z1−k=1∑n−1αk(n−1)∗Z(k+1)2]=E[Z1∗(Z1−k=1∑n−1αk(n−1)∗Z(k+1)]=R(0)−k=1∑n−1αk(n−1)∗R(k)
最终结果
k = R ( n ) − ∑ k = 1 n − 1 α k ( n − 1 ) ∗ R ( n − k ) R ( 0 ) − ∑ k = 1 n − 1 α k ( n − 1 ) ∗ R ( k ) k = \frac{R(n) - \sum_{k=1}^{n-1}\alpha^{(n-1)}_k*R(n-k)}{ R(0) - \sum_{k=1}^{n-1}\alpha_k^{(n-1)}*R(k)} k=R(0)−∑k=1n−1αk(n−1)∗R(k)R(n)−∑k=1n−1αk(n−1)∗R(n−k)
我们可以看出，对前向预测和后向预测线性组合的k，实际上和Levinson迭代中翻转步骤构造的k是一样的。因此我们对这个迭代有了新认识。其实那个翻转就是前向后向预测交互作用的影响。
一开始，我们觉得k是为了销项而凑出来的，现在我们觉得，这个k是一个前向预测的投影系数。
4.3 计算高阶后向预测
现在我们再来想，如何利用低阶前向预测和低阶后向预测去计算高阶后向预测呢。
P r o j ( Z 2 , . . . , Z n + 1 ) Z 1 = P r o j { Z n + 1 } ∪ { Z 2 , . . . , Z n } Z 1 { Z 2 , . . . , Z n + 1 } = { Z 2 , . . . , Z n } ⊕ { Z n + 1 − P r o j { Z 2 , . . . , Z n } Z n + 1 } e f ( n − 1 ) ( n + 1 ) = Z n + 1 − P r o j { Z 2 , . . . , Z n } Z n + 1 Proj_{(Z_2,...,Z_{n+1})}Z_{1} = Proj_{\{Z_{n+1}\} \cup \{Z_2,...,Z_n\}}Z_{1} \\ \{Z_2,...,Z_{n+1}\} = \{Z_2,...,Z_n\} \oplus \{ Z_{n+1} - Proj_{\{ Z_2,...,Z_n\}} Z_{n+1}\} \\ e_f^{(n-1)}(n+1) = Z_{n+1} - Proj_{\{ Z_2,...,Z_n\}} Z_{n+1} Proj(Z2,...,Zn+1)Z1=Proj{Zn+1}∪{Z2,...,Zn}Z1{Z2,...,Zn+1}={Z2,...,Zn}⊕{Zn+1−Proj{Z2,...,Zn}Zn+1}ef(n−1)(n+1)=Zn+1−Proj{Z2,...,Zn}Zn+1
因此
P r o j ( Z 2 , . . . , Z n + 1 ) Z 1 = P r o j { e f ( n − 1 ) ( n + 1 ) } Z 1 + P r o j { Z 2 , . . . , Z n } Z 1 Proj_{(Z_2,...,Z_{n+1})}Z_{1} = Proj_{\{e_f^{(n-1)}(n+1)\}} Z_1 + Proj_{\{Z_2,...,Z_n\} }Z_1 Proj(Z2,...,Zn+1)Z1=Proj{ef(n−1)(n+1)}Z1+Proj{Z2,...,Zn}Z1
后面的计算就完全相同了
4.4 小结
从这里，我们看Levinson迭代，其实并不仅仅是代数上的小技巧，而是蕴含了很深的信号处理的本质。每一个高阶的预测，都是一个低阶的前向预测与一个低阶的后向预测的线性组合。
8801
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
