随机过程的线性预测
1. 随机过程线性预测问题
1.1 问题引入
前面介绍了一些线性预测的方法,比如levension迭代。但是,现在我们想要知道,对随机信号进行线性预测,究竟能够达到一个怎么样的程度。
我们先把线性预测问题重申一下,假设我们要预测一个宽平稳的随机过程。
{ Z ( k ) } \{ Z(k) \} {Z(k)}
我们现在有P个数据,全部都发生在k时刻之前,我们希望用这P个数据去预测k时刻的值
Z ( k ) ← { Z ( k − 1 ) , . . . , Z ( k − P ) } Z(k) \leftarrow \{ Z(k-1),...,Z(k-P) \} Z(k)←{Z(k−1),...,Z(k−P)}
预测的方法是对这P个数据进行线性组合。
∑ m = 1 P α m Z ( k − m ) \sum_{m=1}^P \alpha_m Z(k-m) m=1∑PαmZ(k−m)
我们假设线性的系数已经是最优的了
{ α 1 , . . , α P } = a r g m i n α E ∣ Z ( k ) − ∑ m = 1 P α m Z ( k − m ) ∣ 2 \{ \alpha_1,..,\alpha_P \} = argmin_{\alpha} E|Z(k) - \sum_{m=1}^P \alpha_m Z(k-m) |^2 {α1,..,αP}=argminαE∣Z(k)−m=1∑PαmZ(k−m)∣2
ϵ P = Z ( k ) − ∑ m = 1 P α m Z ( k − m ) \epsilon_P = Z(k) - \sum_{m=1}^P \alpha_m Z(k-m) ϵP=Z(k)−m=1∑PαmZ(k−m)
这个方程的求解之前介绍过,通过levension迭代可以解这个wiener hopf方程
Wiener Hopf R Z α = r \text{ Wiener Hopf} \\ R_Z \alpha = r Wiener HopfRZα=r
一般来说,随着P的增加,估计的误差会不断减小。但是我们想知道,随着P的无穷增加,误差会不会存在一个极限。我们之前探讨过使用Cramér–Rao来判断估计问题的极限,不过这里是线性估计问题,下界可能会有所差异
E ∣ ϵ P ∣ 2 = σ P 2 σ P 2 ≥ σ P + 1 2 ≥ . . . P → ∞ → σ ∞ 2 E|\epsilon_P|^2 = \sigma_P^2 \\ \sigma_P^2 \geq \sigma_{P+1}^2 \geq ... \underrightarrow{P \rightarrow \infty} \quad \sigma_{\infty}^2 E∣ϵP∣2=σP2σP2≥σP+12≥...P→∞σ∞2
我们现在希望能够在有限的条件下,找到这个估计的下界。
这个估计的下界是存在的,我们先说结论
kolmogorov-Szago Identity σ ∞ 2 = e x p ( 1 2 π ∫ − π π l o g ( S Z ( ω ) ) d ω ) \text{kolmogorov-Szago Identity} \sigma_{\infty}^2 = exp(\frac{1}{2\pi} \int _{-\pi}^{\pi} log( S_Z(\omega))d\omega) kolmogorov-Szago Identityσ∞2=exp(2π1∫−ππlog(SZ(ω))dω)
其中SZ(ω)是Z(k)的功率谱密度
1.2 公式解读
我们可以发现,估计的下界与功率谱有关系。因此,我们一旦知道了一个随机过程的功率谱密度,我们就能够知道预测可以做到怎么样的地步
为什么预测会与功率谱密度相关呢
假设功率谱密度比较集中在零点附近,log之后,对功率谱做积分,会得到负值非常大的数字,积分就会趋近于负无穷,取指数之后,就接近与0。说明我们能够很准确的预测。这是因为,如果功率谱比较集中,相关函数就会比较发散,比较胖。就说明不同时刻随机过程的取值,相关度比较高,说明数据之间具有密切的关系,预测就比较好做
这里加一些自己的注解。因为相关函数和功率谱是傅里叶变换对的关系,相关函数是时域,功率谱是频域。如果功率谱比较集中,极端的看就是频谱上只有有限的几根线,时域上就是个标准的正弦函数曲线的叠加。因为相关函数是两个时间差的函数,一旦时域是无限长的函数,就说明任意两个时刻的点都具有相关性。相关性一旦强了,预测就好做
反之,我们知道,如果功率谱很宽,积分会收敛,指数值就能够求出来,就说明误差就会有下界。因为功率谱一旦宽了,相关函数就会比较密集,不同时刻的相关度就会比较低。既然数据相关性不大,那么预测也难做。
1.3 白噪声的预测误差分析
我们以白噪声为例,白噪声的功率谱是个常数。我们求解一下误差下界的这个式子,最后求出来的结果就是功率谱密度
White Noise σ ∞ P = e x p ( 1 2 π ∫ − π π l o g ( S Z ( ω ) ) d ω ) = e x p ( 1 2 π ∫ − π π l o g ( C ) d ω ) = C = S Z ( ω ) \text{White Noise} \sigma_{\infty}^P = exp(\frac{1}{2\pi} \int _{-\pi}^{\pi} log( S_Z(\omega))d\omega) \\ = exp(\frac{1}{2\pi} \int _{-\pi}^{\pi} log( C)d\omega) \\ = C = S_Z(\omega) White Noiseσ∞P=exp(2π1∫−ππlog(SZ(ω))dω)=exp(2π1∫−ππlog(C)dω)=C=SZ(ω)
再来看一下相关函数,如果功率谱密度是常数,则相关函数的这个积分,只有t=0的时候有值
R Z ( τ ) = 1 2 π ∫ − ∞ + ∞ S Z ( ω ) e x p ( j ω τ ) d ω R Z ( τ ) = C δ ( τ ) R Z ( 0 ) = C R_Z(\tau) = \frac{1}{2\pi} \int_{-\infty}^{+\infty} S_Z(\omega) exp(j\omega \tau)d\omega \\ R_Z(\tau) = C \delta(\tau) \\ R_Z(0) = C RZ(τ)=2π1∫−∞+∞SZ(ω)exp(jωτ)dωRZ(τ)=Cδ(τ)RZ(0)=C
因此
σ ∞ P = R Z ( 0 ) = E ( Z ( t ) Z ( t ) ) = E ( ( Z ( t ) − 0 ) 2 ) \sigma_{\infty}^P = R_Z(0) = E(Z(t)Z(t)) = E((Z(t)-0)^2) σ∞P=RZ(0)=E(Z(t)Z(t))=E((Z(t)−0)2)
我们可以看出来,对白噪声进行估计,其估计误差的下界就是方差
2. 相关矩阵的正定性
我们下面要对随机过程是否能够进行完美预测进行分析。我们首先要从相关矩阵出发,看看根据对相关矩阵的解读能够解读出什么信息
2.1 问题描述与结论
2.1.1 相关矩阵的建立
我们首先要建立这样一个问题
假设我们有一个P阶的随机矢量
Z ( P ) = ( Z ( k − 1 ) , . . . , Z ( k − P ) ) T Z_{(P)} = (Z_{(k-1)},...,Z_{(k-P)})^T Z(P)=(Z(k−1),...,Z(k−P))T
我们假设这是一个复向量,我们要求这个复向量的相关矩阵。因为复数自身的共轭乘积一定是非负的,所以这个相关矩阵的对角元一定都不是负的,因此矩阵一定是非负定的。
这个P阶的相关矩阵表示为
R Z ( P ) = E ( Z ( P ) Z ( P ) H ) ≥ 0 R_Z^{(P)} = E(Z_{(P)}Z_{(P)}^H) \geq 0 RZ(P)=E(Z(P)Z(P)H)≥0
可以得到一个P阶的toplitz矩阵
R Z ( P ) = ( r 0 r 1 . . . r ( P − 1 ) r − 1 r 0 . . . r ( P − 2 ) . . . . . . . . . . . . r ( − P + 1 ) r ( − P + 2 ) . . . r 0 ) R_Z^{(P)} = \begin{pmatrix} r_0 & r_1 &...& r_{(P-1)} \\ r_{-1} & r_0 &...& r_{(P-2)} \\ ... & ... &...& ... \\ r_{(-P+1)}& r_{(-P+2)}&...&r_0 \end{pmatrix} RZ(P)=⎝⎜⎜⎛r0r−1...r(−P+1)r1r0...r(−P+2)............r(P−1)r(P−2)...r0⎠⎟⎟⎞
2.1.2 相关矩阵的性质
这个矩阵具有这样的性质:如果下标是相反数,对应的值是共轭的
r
−
k
=
r
k
∗
r_{-k} = r_k^*
r−k=rk∗
证明
r − k = Z 0 ∗ Z k ∗ r k = Z k ∗ Z 0 ∗ r − k = r k ∗ r_{-k} = Z_0 * Z_k^* \\ r_k = Z_k * Z_0^* \\ r_{-k} = r_k^* r−k=Z0∗Zk∗rk=Zk∗Z0∗r−k=rk∗
因此,这个矩阵的转置等于其共轭,其共轭转置等于其本身
R Z ( p ) ‾ = ( R Z ( p ) ) T R Z ( p ) = ( R Z ( p ) ) H \overline{R_Z^{(p)}} = (R_Z^{(p)})^T \\ R_Z^{(p)} = (R_Z^{(p)} )^H RZ(p)=(RZ(p))TRZ(p)=(RZ(p))H
2.1.3 相关矩阵正定性的结论
在这里我们给出一个结论
假设构成P阶toplitz矩阵的元素r§
r ( P ) = ( r 0 , . . . , r P − 1 ) r_{(P)} = (r_0,...,r_{P-1}) r(P)=(r0,...,rP−1)
结论:如果这个序列做离散傅里叶变换是非负的,那么这个序列构成的toplitz矩阵一定是非负定的。如果这个序列任意阶toplitz矩阵是非负定的,那么这个序列的离散傅里叶变换一定是非负的
Herglotz { ∀ P , R Z ( P ) ≥ 0 } ⇔ { S ( ω ) = ∑ k = − ∞ + ∞ r k e x p ( − j ω k ) ≥ 0 } \text{Herglotz} \\ \{\forall P, R_Z^{(P)} \geq 0 \} \Leftrightarrow \{ S(\omega) = \sum_{k=-\infty}^{+\infty} r_k exp(-j \omega k) \geq 0\} Herglotz{∀P,RZ(P)≥0}⇔{S(ω)=k=−∞∑+∞rkexp(−jωk)≥0}
后面的式子也是功率谱密度是数学定义
Mathematical S ( ω ) = ∑ k = − ∞ + ∞ r k e x p ( − j ω k ) \text{Mathematical} S(\omega) = \sum_{k=-\infty}^{+\infty} r_k exp(-j \omega k) MathematicalS(ω)=k=−∞∑+∞rkexp(−jωk)
其实相关函数的傅里叶变换是非负的是显然的,因为我们可以从物理定义看出来。只不过我们现在要从数学定义出发去证明
物理定义
Physical 连续时间 S Z ( ω ) = l i m T → ∞ 1 T E ∣ ∫ − T 2 T 2 Z ( t ) e x p ( − j ω t ) d t ∣ 2 ≥ 0 \text{Physical} \\ \text{连续时间} S_Z(\omega) = lim_{T \rightarrow \infty} \frac{1}{T} E |\int _{-\frac{T}{2}}^{\frac{T}{2}} Z(t) exp(-j \omega t)dt|^2 \geq 0 \\ Physical连续时间SZ(ω)=limT→∞T1E∣∫−2T2TZ(t)exp(−jωt)dt∣2≥0
- 先短时傅里叶变换,规避积分不收敛
- 求模取平方
- 取均值
- 时间归一化,并且T趋近于无穷大
离散时间 = l i m P → ∞ 1 2 P + 1 E ∣ ∑ k = − P P Z ( k ) e x p ( − j ω k ) ∣ 2 ≥ 0 \text{离散时间} = lim_{P \rightarrow \infty} \frac{1}{2P+1} E|\sum_{k=-P}^P Z(k) exp(-j\omega k)|^2 \geq 0 离散时间=limP→∞2P+11E∣k=−P∑PZ(k)exp(−jωk)∣2≥0
2.2 证明
我们假设知道了一个复序列的傅里叶变换是非负的,我们要证明这个复序列构成的toplitz是非负定的
{ ∀ P , R Z ( P ) ≥ 0 } ⇐ { S ( ω ) = ∑ k = − ∞ + ∞ r k e x p ( − j ω k ) ≥ 0 } \{\forall P, R_Z^{(P)} \geq 0 \} \Leftarrow \{ S(\omega) = \sum_{k=-\infty}^{+\infty} r_k exp(-j \omega k) \geq 0\} {∀P,RZ(P)≥0}⇐{S(ω)=k=−∞∑+∞rkexp(−jωk)≥0}
证明一个矩阵是非负定的,只要他的二次型不小于0即可。
假设我们有一个复序列Z
∀
Z
=
(
Z
1
,
.
.
,
Z
P
)
T
\forall Z = (Z_1,..,Z_P)^T
∀Z=(Z1,..,ZP)T
Z H R Z ( P ) Z = ∑ k = 0 P − 1 ∑ m = 0 P − 1 Z k Z m ∗ R Z ( P ) ( k , m ) = ∑ k = 0 P − 1 ∑ m = 0 P − 1 Z k Z m ∗ r k − m = ∑ k = 0 P − 1 ∑ m = 0 P − 1 Z k Z m ∗ ( 1 2 π ∫ − π + π S ( ω ) e x p ( j ω ( k − m ) ) ) d ω = 1 2 π ∫ − π + π ( ∑ k = 0 P − 1 ∑ m = 0 P − 1 Z k Z m ∗ ) S ( ω ) e x p ( j ω ( k − m ) ) ) d ω = 1 2 π ∫ − π + π ( ∑ k = 0 P − 1 Z k e x p ( j ω k ) ) ( ∑ m = 0 P − 1 Z m ∗ e x p ( − j ω m ) ) S ( ω ) d ω = 1 2 π ∫ − π + π ( ∑ k = 0 P − 1 Z k e x p ( j ω k ) ) ( ∑ m = 0 P − 1 Z m e x p ( j ω m ) ‾ ) S ( ω ) d ω = 1 2 π ∫ − π + π ∣ ∑ k = 0 P − 1 Z k e x p ( j ω k ) ∣ 2 S ( ω ) d ω Z^H R_Z^{(P)} Z = \sum_{k=0}^{P-1} \sum_{m=0}^{P-1} Z_k Z_m^* R_Z^{(P)}(k,m) \\ = \sum_{k=0}^{P-1} \sum_{m=0}^{P-1} Z_k Z_m^* r_{k-m} \\ = \sum_{k=0}^{P-1} \sum_{m=0}^{P-1} Z_k Z_m^* (\frac{1}{2 \pi}\int_{-\pi}^{+\pi} S(\omega) exp(j \omega (k-m)) )d \omega \\ = \frac{1}{2 \pi}\int_{-\pi}^{+\pi} (\sum_{k=0}^{P-1} \sum_{m=0}^{P-1} Z_k Z_m^* )S(\omega) exp(j \omega (k-m)) )d \omega \\ = \frac{1}{2 \pi}\int_{-\pi}^{+\pi} (\sum_{k=0}^{P-1} Z_k exp(j\omega k) )( \sum_{m=0}^{P-1} Z_m^* exp(-j\omega m)) S(\omega) d \omega \\ = \frac{1}{2 \pi}\int_{-\pi}^{+\pi} (\sum_{k=0}^{P-1} Z_k exp(j\omega k) )(\overline{ \sum_{m=0}^{P-1} Z_m exp(j\omega m)}) S(\omega) d \omega \\ = \frac{1}{2 \pi}\int_{-\pi}^{+\pi} |\sum_{k=0}^{P-1} Z_k exp(j\omega k) |^2 S(\omega) d \omega ZHRZ(P)Z=k=0∑P−1m=0∑P−1ZkZm∗RZ(P)(k,m)=k=0∑P−1m=0∑P−1ZkZm∗rk−m=k=0∑P−1m=0∑P−1ZkZm∗(2π1∫−π+πS(ω)exp(jω(k−m)))dω=2π1∫−π+π(k=0∑P−1m=0∑P−1ZkZm∗)S(ω)exp(jω(k−m)))dω=2π1∫−π+π(k=0∑P−1Zkexp(jωk))(m=0∑P−1Zm∗exp(−jωm))S(ω)dω=2π1∫−π+π(k=0∑P−1Zkexp(jωk))(m=0∑P−1Zmexp(jωm))S(ω)dω=2π1∫−π+π∣k=0∑P−1Zkexp(jωk)∣2S(ω)dω
因为模的平方一定是大于等于0的,并且我们已知条件中,复序列rk的傅里叶变换也是大于等于0的,功率谱密度恒大于等于0
∣
∑
k
=
0
P
−
1
Z
k
e
x
p
(
j
ω
k
)
∣
2
≥
0
S
(
ω
)
=
∑
k
=
−
∞
+
∞
r
k
e
x
p
(
−
j
ω
k
)
≥
0
|\sum_{k=0}^{P-1} Z_k exp(j\omega k) |^2 \geq 0 \\ S(\omega) = \sum_{k=-\infty}^{+\infty} r_k exp(-j \omega k) \geq 0
∣k=0∑P−1Zkexp(jωk)∣2≥0S(ω)=k=−∞∑+∞rkexp(−jωk)≥0
因此我们可以证明出,toplitz矩阵的二次型一定是非负的。
Z H R Z ( P ) Z = 1 2 π ∫ − π + π ∣ ∑ k = 0 P − 1 Z k e x p ( j ω k ) ∣ 2 S ( ω ) d ω ≥ 0 Z^H R_Z^{(P)} Z = \frac{1}{2 \pi}\int_{-\pi}^{+\pi} |\sum_{k=0}^{P-1} Z_k exp(j\omega k) |^2 S(\omega) d \omega \geq 0 ZHRZ(P)Z=2π1∫−π+π∣k=0∑P−1Zkexp(jωk)∣2S(ω)dω≥0
因此toplitz矩阵一定是非负定的。
R Z ( P ) ≥ 0 R_Z^{(P)} \geq 0 RZ(P)≥0
但是相关矩阵是非负定的有两种情况,一种是相关矩阵可能是奇异的,另一种是相关矩阵恒大于0。这两种情况能够得到不同的线性预测结果。下面我们会针对相关矩阵的奇异性,进行继续讨论
3. 完美预测
3.1 相关矩阵的奇异性
相关矩阵是非负定的,就说明,相关矩阵可能是奇异的,也就是相关矩阵的二次型可能是0
R Z ( P ) ≥ 0 R_Z^{(P)} \geq 0 RZ(P)≥0
3.2 相关矩阵完全奇异
因此,我们就取降一阶的toplitz矩阵来进行分析
如果我们发现降一阶的矩阵还不是正定的,就再降一阶,直到变成正定矩阵
低一阶的toplitz矩阵
R Z ( P ) ≥ 0 → R Z ( P − 1 ) ≥ 0 → R Z ( P − 2 ) ≥ 0 → . . . R_Z^{(P)} \geq 0 \rightarrow R_Z^{(P-1)} \geq 0 \rightarrow R_Z^{(P-2)} \geq 0 \rightarrow ... RZ(P)≥0→RZ(P−1)≥0→RZ(P−2)≥0→...
如果矩阵一直有奇异性,降低到头,就是r0
如果r0也是0,那么这个随机过程整个就是0
r ( 0 ) = 0 → Z ( k ) = 0 r(0) = 0 \rightarrow Z(k) = 0 r(0)=0→Z(k)=0
那么这个随机过程就已经是确定的了,完全没有随机性了
3.3 相关矩阵部分奇异
如果r0大于0,那么一定能够找到一个临界点,使得相关矩阵没有奇异性
R Z ( n ) ≥ 0 and R Z ( n − 1 ) > 0 R_Z^{(n)} \geq 0 \text{ and } R_Z^{(n-1)} >0 RZ(n)≥0 and RZ(n−1)>0
如果存在这样的n,那么功率谱一定能够写成这样的形式
S ( ω ) = 2 π ∑ i = 1 n β i δ ( ω − ω i ) ⇔ r k = ∑ i = 1 n β i e x p ( j ω i k ) S(\omega) = 2 \pi \sum_{i=1}^n \beta_i \delta(\omega - \omega_i) \Leftrightarrow r_k = \sum_{i=1}^n \beta_i exp(j \omega_i k) S(ω)=2πi=1∑nβiδ(ω−ωi)⇔rk=i=1∑nβiexp(jωik)
这是一个谐波信号。谐波信号一定是能够被完美预测的
我们来看一个典型的谐波信号
Z ( t ) = A c o s ( ω t + ϕ ) A . r . v . ϕ ∼ U ( 0 , 2 π ) R Z ( τ ) = E ( A 2 ) 2 c o s ( ω τ ) Z(t) = A cos(\omega t +\phi) \\ A.r.v. \quad \phi \sim U(0,2\pi) \\ R_Z(\tau) = \frac{E(A^2)}{2} cos(\omega \tau) Z(t)=Acos(ωt+ϕ)A.r.v.ϕ∼U(0,2π)RZ(τ)=2E(A2)cos(ωτ)
其中A是随机变量。φ是(0,2π)的均匀分布
这个谐波信号几乎没有随机性,只有A和φ两个变量。t和φ、A是解耦的。对余弦波形进行预测,只要在一个周期采样四个点就能够预测
现在信号是这样的
Z ( t ) = ∑ i = 1 n A i c o s ( ω i t + ϕ i ) ⇒ R Z ( τ ) = ∑ i = 1 n E ( A 2 ) 2 c o s ( ω i τ ) Z(t) = \sum_{i=1}^n A_i cos(\omega _i t +\phi _i) \Rightarrow R_Z(\tau) = \sum_{i=1}^n \frac{E(A^2)}{2} cos(\omega_i \tau) Z(t)=i=1∑nAicos(ωit+ϕi)⇒RZ(τ)=i=1∑n2E(A2)cos(ωiτ)
3.4 小结
我们得到了与随机过程的线性预测有关的第一个结论。如果相关矩阵有奇异性。无论在哪一阶发生了奇异,这个随机过程就是有这个阶数构成的谐波。谐波里面就有这么多的分量。
谐波信号得到的功率谱是线谱,仅仅有有限个数的线构成的谱。线谱是没有随机性的
4. 非完美预测
4.1 非完美预测的条件
刚才我们了解到了,如果我们的功率谱是线谱,随机信号是几乎没有随机性的,就可以做到完美预测。因此,如果我们想要分析线性预测问题的误差,该随机信号的预测问题就必须是个非完美预测问题。
如果我们希望随机信号无法做到完美预测,我们需要有这样的条件
- 功率谱必须是连续谱
- paley-wiener条件
paley-wiener 1 2 π ∫ − π + π l o g S Z ( ω ) d ω > − ∞ \text{paley-wiener} \frac{1}{2 \pi} \int _{-\pi}^{+\pi} log S_Z(\omega) d\omega > -\infty paley-wiener2π1∫−π+πlogSZ(ω)dω>−∞
我们刚才说过,如果这个积分等于-∞,误差极限会趋于0,因此是能够被完美预测的。如果要达到非完美预测,这个积分就要收敛
如果随机信号能够满足paley-wiener条件,它的功率谱一定能够表示成下面的样子
S Z ( ω ) = ∣ B ( ω ) ∣ 2 S_Z(\omega) = |B(\omega)|^2 SZ(ω)=∣B(ω)∣2
其中,B(Z)是单位圆内的解析函数
B ( Z ) is Analytic at |Z| <1 B ( Z ) = ∑ k = 0 ∞ b k Z k , b 0 > 0 B(Z) \text{is Analytic at |Z| <1} \\ B(Z) = \sum_{k=0}^{\infty} b_k Z^k, b_0 >0 B(Z)is Analytic at |Z| <1B(Z)=k=0∑∞bkZk,b0>0
所谓解析函数,就是在解析区间内任意点都能够进行泰勒级数展开,并且没有奇点。并且函数可以转化为幂级数展开,并且幂级数收敛,幂级数中不能有复指数。
这个解析函数还有一个要求,就是b0一定是大于0的
B(ω)就是在单位圆上的取值
B ( ω ) = B ( e x p ( j ω ) ) B(\omega) = B(exp(j \omega)) B(ω)=B(exp(jω))
总结一下。维纳做的工作说明,如果随机过程是非完美预测的,一定满足paley-wiener条件。其功率谱一定能够表示为B的形式。并且B要满足两个条件:在单位圆内解析,零点展开的系数大于0。
4.2 谱分解
我们注意到,功率谱是个二阶量,没有相位信号,要拆分成两个一阶量,必定涉及到相位恢复的问题。
而且这个解析函数B有很多个,我们必须从中进行挑选。这个拆分、挑选的过程就叫做谱分解。在前面维纳滤波中介绍过
Spectral Factorization \text{Spectral Factorization} Spectral Factorization
我们挑选的条件是相位极小。并且要求在单位圆中没有极点。在单位圆外面没有零点。这样的话B和1/B就都是稳定的
我们是需要把B求解出来的,但是我们需要做点准备,因此后面要解决一下其他的问题
4.3 线性代数分析法与预测误差
事实上,我们还是在研究随机过程线性预测的最优误差问题。我们首先通过线性代数的分析对误差进行表示。在这个过程中,我们会先通过讨论相关矩阵的LU分解,得到一些关系式,用于后面的误差分析。
4.3.1 相关矩阵的LU分解
4.3.1.1 Cholesky分解
我们首先要对相关矩阵做一些线性代数的分析。我们要对相关矩阵进行LU分解。LU分解就是高斯消去。同时在这里也叫做Cholesky
因为这里相关矩阵具有对称性,其共轭转置等于其本身,因此L和U矩阵也具有对称性
R Z ( P ) = ( R Z ( P ) ) H R_Z^{(P)} = (R_Z^{(P)})^H RZ(P)=(RZ(P))H
因此,可以得到L和U矩阵是共轭转置的关系
可得
R Z ( P ) = ( L ( p ) ) H ( L ( p ) ) ( a ) R_Z^{(P)} = (L^{(p)})^H (L^{(p)}) \quad\quad(a) RZ(P)=(L(p))H(L(p))(a)
4.3.1.2 递推关系
现在我们希望P+1阶的相关矩阵也能做Cholesky分解。但是,我们想要通过寻找和P阶相关矩阵之间的递推关系来实现。也就是找到P阶相关矩阵和P+1阶相关矩阵之间的关系
R Z ( P + 1 ) = ( L ( p + 1 ) ) H ( L ( p + 1 ) ) R_Z^{(P+1)} =(L^{(p+1)})^H (L^{(p+1)}) RZ(P+1)=(L(p+1))H(L(p+1))
L ( P ) → L ( P + 1 ) L^{(P)} \rightarrow L^{(P+1)} L(P)→L(P+1)
下三角矩阵可以写成这样的形式
L ( P + 1 ) = ( c 0 b A ) L^{(P+1)} = \begin{pmatrix} c & 0 \\ b & A \end{pmatrix} L(P+1)=(cb0A)
下三角矩阵用分块矩阵的方式进行表示。c表示最左上角的元素,是个常数。b是一个向量,A是一个矩阵。
R Z ( P + 1 ) = ( L ( p + 1 ) ) H ( L ( p + 1 ) ) = ( c ‾ b H 0 A H ) ∗ ( c 0 b A ) = ( ∣ c ∣ 2 + b H b b H A A H b A H A ) = ( r 0 r 1 . . . r P r − 1 . . . R Z ( P ) r − P ) = ( r 0 r 1 . . . r P r 1 ∗ . . . R Z ( P ) r P ∗ ) R_Z^{(P+1)}=(L^{(p+1)})^H (L^{(p+1)}) = \begin{pmatrix} \overline{c} & b^H \\ 0 & A^H \end{pmatrix} *\begin{pmatrix} c & 0 \\ b & A \end{pmatrix}\\ = \begin{pmatrix} |c|^2 +b^H b & b^H A \\ A^H b & A^H A \end{pmatrix} \\ = \begin{pmatrix} r_0 & r_1 &...& r_{P} \\ r_{-1} & \\ ... & & R_Z^{(P)} \\ r_{-P}& \end{pmatrix} = \begin{pmatrix} r_0 & r_1 &...& r_{P} \\ r_{1}^* & \\ ... & & R_Z^{(P)} \\ r_{P}^*& \end{pmatrix} RZ(P+1)=(L(p+1))H(L(p+1))=(c0bHAH)∗(cb0A)=(∣c∣2+bHbAHbbHAAHA)=⎝⎜⎜⎛r0r−1...r−Pr1...RZ(P)rP⎠⎟⎟⎞=⎝⎜⎜⎛r0r1∗...rP∗r1...RZ(P)rP⎠⎟⎟⎞
我们得到了三个等式
r 0 = ∣ c ∣ 2 + b H b ( 1 ) A H A = R Z ( P ) ( 2 ) A H b = ( r 1 ∗ . . . r P ∗ ) ( 3 ) r_0 = |c|^2 +b^H b \quad\quad (1) \\ A^HA = R_Z^{(P)} \quad\quad (2) \\ A^H b = \begin{pmatrix} r_{1}^* \\ ... \\ r_{P}^* \end{pmatrix} \quad\quad (3) r0=∣c∣2+bHb(1)AHA=RZ(P)(2)AHb=⎝⎛r1∗...rP∗⎠⎞(3)
这里需要说明的是LU分解并不是唯一的
H = L H L L ~ = d i a g ( ± 1 , . . . , ± 1 ) L H = L ~ H L ~ H = L^H L \\ \widetilde {L} = diag(±1,...,±1) L \\ H = \widetilde {L}^H \widetilde {L} H=LHLL =diag(±1,...,±1)LH=L HL
只要L乘以一个对角阵,对角阵元都是正负1,那么得到的新的L矩阵,依旧可以用来做LU分解
所以,这里我们要加一个限定条件,要求L矩阵所有的对角元都是非负的。那么结果就是唯一的了
对LU分解做了限定之后,我们就可以继续求解了
首先,从(2)式子能够解得A
A H A = R Z ( P ) = ( L ( P ) ) H ( L ( P ) ) A^H A = R_Z^{(P)} = (L^{(P)})^H (L^{(P)}) AHA=RZ(P)=(L(P))H(L(P))
A = L ( P ) ( A ) A = L^{(P)} \quad\quad (A) A=L(P)(A)
从(3)式子能够解得b
b = ( A H ) − 1 ∗ ( r 1 ∗ . . . r P ∗ ) = ( ( L ( P ) ) H ) − 1 ∗ ( r 1 ∗ . . . r P ∗ ) ( B ) b = (A^H)^{-1}*\begin{pmatrix} r_{1}^* \\ ... \\ r_{P}^* \end{pmatrix} \\ = ((L^{(P)})^H)^{-1} *\begin{pmatrix} r_{1}^* \\ ... \\ r_{P}^* \end{pmatrix} \quad\quad (B) b=(AH)−1∗⎝⎛r1∗...rP∗⎠⎞=((L(P))H)−1∗⎝⎛r1∗...rP∗⎠⎞(B)
从式子(1) 可以求得c
c = ( r 0 − b H b ) 1 2 ( C ) c = (r_0 - b^H b)^{\frac{1}{2}} \quad\quad (C) c=(r0−bHb)21(C)
求得了Abc之后,我们就可以进行公式的递推了
L ( P + 1 ) = ( ( r 0 − b H b ) 1 2 0 b L ( P ) ) L^{(P+1)} = \begin{pmatrix} (r_0 - b^H b)^{\frac{1}{2}} & 0 \\ b & L^{(P)} \end{pmatrix} L(P+1)=((r0−bHb)21b0L(P))
其中
b = ( ( L ( P ) ) H ) − 1 ∗ ( r 1 ∗ . . . r P ∗ ) b = ((L^{(P)})^H)^{-1} *\begin{pmatrix} r_{1}^* \\ ... \\ r_{P}^* \end{pmatrix} b=((L(P))H)−1∗⎝⎛r1∗...rP∗⎠⎞
4.3.2 误差模的平方与LU分解的关系
4.3.2.1 误差的表示
得到了相关矩阵LU分解的递推关系之后,我们可以回去考察我们的预测问题了
我们表示一下误差
ϵ P = Z ( k ) − ∑ m = 1 P α m Z ( k − m ) \epsilon_P = Z(k) - \sum_{m=1}^P \alpha_m Z(k-m) ϵP=Z(k)−m=1∑PαmZ(k−m)
由正交性原理可知,如果线性估计是最优的,误差会与所有的原材料正交
E ( ϵ P Z ∗ ( l ) ) = 0 l = { k − 1 , . . . , k − P } E(\epsilon_P Z^*(l)) = 0 \\ l = \{k-1,...,k-P \} E(ϵPZ∗(l))=0l={k−1,...,k−P}
对正交性原理式子进行转化可得
E ( ϵ P Z ∗ ( l ) ) = E ( ( Z ( k ) − ∑ m = 1 P α m Z ( k − m ) ( Z ∗ ( l ) ) ) = E ( Z ( k ) Z ∗ ( l ) ) − ∑ m = 1 P α m E ( Z ( k − m ) Z ∗ ( l ) ) = r ( k − l ) − ∑ m = 1 P α m r ( k − m − l ) = 0 E(\epsilon_P Z^*(l)) = E((Z(k) - \sum_{m=1}^P \alpha_m Z(k-m)(Z^*(l))) \\ = E(Z(k)Z^*(l)) - \sum_{m=1}^P \alpha_m E(Z(k-m)Z^*(l)) \\ = r_{(k-l)} - \sum_{m=1}^P \alpha_m r_{(k-m-l)} = 0 E(ϵPZ∗(l))=E((Z(k)−m=1∑PαmZ(k−m)(Z∗(l)))=E(Z(k)Z∗(l))−m=1∑PαmE(Z(k−m)Z∗(l))=r(k−l)−m=1∑Pαmr(k−m−l)=0
可得
r ( k − l ) = ∑ m = 1 P α m r ( k − m − l ) l = { k − 1 , . . . , k − P } r_{(k-l)} = \sum_{m=1}^P \alpha_m r_{(k-m-l)} \\ l = \{k-1,...,k-P \} r(k−l)=m=1∑Pαmr(k−m−l)l={k−1,...,k−P}
我们可以得到一个方程
( r 1 . . . r p ) = ( r 0 r − 1 . . . r ( − P + 1 ) r 1 r 0 . . . r ( − P + 2 ) r ( P − 1 ) r ( P − 2 ) . . . r 0 ) ( α 1 . . . α p ) ( i ) \begin{pmatrix} r_1 \\ ... \\ r_p \end{pmatrix} = \begin{pmatrix} r_0 & r_{-1} &... & r_{(-P+1)} \\ r_1 & r_0 & ... & r_{(-P+2)}\\ r_{(P-1)} & r_{(P-2)} & ... & r_0 \end{pmatrix} \begin{pmatrix} \alpha_1 \\ ... \\ \alpha_p \end{pmatrix} \quad\quad (i) ⎝⎛r1...rp⎠⎞=⎝⎛r0r1r(P−1)r−1r0r(P−2).........r(−P+1)r(−P+2)r0⎠⎞⎝⎛α1...αp⎠⎞(i)
我们发现正交性原理得到的矩阵与随机矢量的自相关矩阵相差一个转置
R Z ( P ) = ( r 0 r 1 . . . r ( P − 1 ) r − 1 r 0 . . . r ( P − 2 ) . . . . . . . . . . . . r ( − P + 1 ) r ( − P + 2 ) . . . r 0 ) ( i i ) R_Z^{(P)} = \begin{pmatrix} r_0 & r_1 &...& r_{(P-1)} \\ r_{-1} & r_0 &...& r_{(P-2)} \\ ... & ... &...& ... \\ r_{(-P+1)}& r_{(-P+2)}&...&r_0 \end{pmatrix} \quad\quad (ii) RZ(P)=⎝⎜⎜⎛r0r−1...r(−P+1)r1r0...r(−P+2)............r(P−1)r(P−2)...r0⎠⎟⎟⎞(ii)
我们把(ii)代入(i)中可以得到
( r 1 . . . r p ) = ( R Z ( P ) ) T ( α 1 . . . α p ) ( i i i ) \begin{pmatrix} r_1 \\ ... \\ r_p \end{pmatrix} = (R_Z^{(P)})^T \begin{pmatrix} \alpha_1 \\ ... \\ \alpha_p \end{pmatrix} \quad\quad (iii) ⎝⎛r1...rp⎠⎞=(RZ(P))T⎝⎛α1...αp⎠⎞(iii)
4.3.2.2 误差模的期望的表示
我们再来考察误差模的期望
E ∣ ϵ P ∣ 2 = E ( ϵ P ϵ P ‾ ) = E ( ϵ P ( Z ( k ) − ∑ m = 1 P α m Z ( k − m ) ‾ ) E| \epsilon_P|^2 = E(\epsilon_P \overline{\epsilon_P}) \\ = E(\epsilon_P \overline{(Z(k)-\sum_{m=1}^P \alpha_m Z(k-m)}) E∣ϵP∣2=E(ϵPϵP)=E(ϵP(Z(k)−m=1∑PαmZ(k−m))
这里有点不太懂,虽然误差与估计的原材料正交,与原材料的共轭也正交吗?
这是最优的估计,误差与原材料是正交的,可以得到
E ∣ ϵ P ∣ 2 = E ( ϵ P ( Z ( k ) − ∑ m = 1 P α m Z ( k − m ) ‾ ) = E ( ϵ P Z ∗ ( k ) ) = E ( Z ( k ) Z ∗ ( k ) − ∑ m = 1 P α m Z ( k − m ) Z ∗ ( k ) ) = r 0 − ∑ m = 1 P α m r ( − m ) = r 0 − ∑ m = 1 P α m r m ∗ = r 0 − ( r 1 ∗ . . . r P ∗ ) ∗ ( α 1 . . . α P ) ( i v ) E| \epsilon_P|^2 =E(\epsilon_P \overline{(Z(k)-\sum_{m=1}^P \alpha_m Z(k-m)}) \\ = E(\epsilon_P Z^*(k)) = E(Z(k)Z^*(k)-\sum_{m=1}^P \alpha_m Z(k-m)Z^*(k)) \\ = r_0 - \sum_{m=1}^P \alpha_m r_{(-m)} =r_0 - \sum_{m=1}^P \alpha_m r_m^* \\ = r_0 - \begin{pmatrix} r_1^* & ... &r_P^* \\ \end{pmatrix} *\begin{pmatrix} \alpha_1 \\ ... \\ \alpha_P \\ \end{pmatrix} \quad\quad (iv) E∣ϵP∣2=E(ϵP(Z(k)−m=1∑PαmZ(k−m))=E(ϵPZ∗(k))=E(Z(k)Z∗(k)−m=1∑PαmZ(k−m)Z∗(k))=r0−m=1∑Pαmr(−m)=r0−m=1∑Pαmrm∗=r0−(r1∗...rP∗)∗⎝⎛α1...αP⎠⎞(iv)
把(iii)代入(iv)中得
E ∣ ϵ P ∣ 2 = r 0 − ( r 1 ∗ . . . r P ∗ ) ∗ ( α 1 . . . α P ) = r 0 − ( r 1 ∗ . . . r P ∗ ) ∗ ( ( R Z ( P ) ) T ) ( − 1 ) ∗ ( r 1 . . . r p ) ( v ) E| \epsilon_P|^2 = r_0 - \begin{pmatrix} r_1^* & ... &r_P^* \\ \end{pmatrix} *\begin{pmatrix} \alpha_1 \\ ... \\ \alpha_P \\ \end{pmatrix} \quad\quad \\ = r_0 - \begin{pmatrix} r_1^* & ... &r_P^* \\ \end{pmatrix}*((R_Z^{(P)})^T)^{(-1)}*\begin{pmatrix} r_1 \\ ... \\ r_p \end{pmatrix} \quad\quad (v) E∣ϵP∣2=r0−(r1∗...rP∗)∗⎝⎛α1...αP⎠⎞=r0−(r1∗...rP∗)∗((RZ(P))T)(−1)∗⎝⎛r1...rp⎠⎞(v)
由(a)式可知
R Z ( P ) = ( L ( p ) ) H ( L ( p ) ) ( a ) R_Z^{(P)} = (L^{(p)})^H (L^{(p)}) \quad\quad(a) RZ(P)=(L(p))H(L(p))(a)
可得
( R Z ( P ) ) T = ( L ( p ) ) T ( L ( p ) ) ∗ ( ( R Z ( P ) ) T ) − 1 = ( ( L ( p ) ) ∗ ) − 1 ( ( L ( p ) ) T ) − 1 ( v i ) (R_Z^{(P)})^T = (L^{(p)})^T (L^{(p)})^* \\ ((R_Z^{(P)})^T)^{-1} = ((L^{(p)})^*)^{-1} ((L^{(p)})^T) ^{-1} \quad\quad (vi) (RZ(P))T=(L(p))T(L(p))∗((RZ(P))T)−1=((L(p))∗)−1((L(p))T)−1(vi)
(v)式可化为
r 0 − ( r 1 ∗ . . . r P ∗ ) ∗ ( ( R Z ( P ) ) T ) ( − 1 ) ∗ ( r 1 . . . r p ) = r 0 − ( r 1 ∗ . . . r P ∗ ) ∗ ( ( L ( P ) ) ∗ ) − 1 ( ( L ( P ) ) T ) − 1 ∗ ( r 1 . . . r p ) ( v i i ) r_0 - \begin{pmatrix} r_1^* & ... &r_P^* \\ \end{pmatrix}*((R_Z^{(P)})^T)^{(-1)}*\begin{pmatrix} r_1 \\ ... \\ r_p \end{pmatrix} \\ = r_0 - \begin{pmatrix} r_1^* & ... &r_P^* \\ \end{pmatrix}*((L^{(P)})^*)^{-1} ((L^{(P)})^T) ^{-1}*\begin{pmatrix} r_1 \\ ... \\ r_p \end{pmatrix} \quad\quad (vii) r0−(r1∗...rP∗)∗((RZ(P))T)(−1)∗⎝⎛r1...rp⎠⎞=r0−(r1∗...rP∗)∗((L(P))∗)−1((L(P))T)−1∗⎝⎛r1...rp⎠⎞(vii)
我们在4.3.2中得到了三个相关矩阵的递推公式,我们把(B)代入,可以发现误差模的平方与他们的关系
b = ( ( L ( P ) ) H ) − 1 ∗ ( r 1 ∗ . . . r P ∗ ) ( B ) b = ((L^{(P)})^H)^{-1} *\begin{pmatrix} r_{1}^* \\ ... \\ r_{P}^* \end{pmatrix} \quad\quad (B) b=((L(P))H)−1∗⎝⎛r1∗...rP∗⎠⎞(B)
对该式子做变换
b ‾ = ( ( L ( P ) ) T ) − 1 ∗ ( r 1 . . . r P ) ( v i i i ) b T = ( r 1 ∗ . . . r P ∗ ) ∗ ( ( L ( P ) ) ∗ ) − 1 ( i x ) \overline {b} = ((L^{(P)})^T)^{-1} *\begin{pmatrix} r_{1} \\ ... \\ r_{P} \end{pmatrix} \quad\quad (viii)\\ b^T = \begin{pmatrix} r_1^* & ... &r_P^* \\ \end{pmatrix}*((L^{(P)})^*)^{-1} \quad\quad (ix) b=((L(P))T)−1∗⎝⎛r1...rP⎠⎞(viii)bT=(r1∗...rP∗)∗((L(P))∗)−1(ix)
(vii)(ix)代入(vii)可得
E ∣ ϵ P ∣ 2 = r 0 − b T b ‾ = r 0 − b H b ( x ) E| \epsilon_P|^2 = r_0 - b^T \overline{b} \\ = r_0 - b^H b \quad\quad (x) E∣ϵP∣2=r0−bTb=r0−bHb(x)
把c的关系式©代入(x)可得
c = ( r 0 − b H b ) 1 2 ( C ) c = (r_0 - b^H b)^{\frac{1}{2}} \quad\quad (C) c=(r0−bHb)21(C)
E ∣ ϵ P ∣ 2 = c 2 = ∣ L 00 ( P + 1 ) ∣ 2 E| \epsilon_P|^2 = c^2 = |L^{(P+1)}_{00}|^2 E∣ϵP∣2=c2=∣L00(P+1)∣2
我们发现,随机过程线性预测误差与LU分解最左上角的元素(0,0)的平方有关系
4.5 复变函数分析法与谱分解
前面,我们通过引入相关矩阵的LU分解,发现随机过程线性预测误差与LU分解最左上角的元素有关,这里,我们会通过复变函数的分析方法,继续探讨预测。首先,我们会回到谱分解的问题上,求解B(ω)函数。然后把复分析得到的预测误差与线性代数得到的预测误差进行统一。
4.5.1 谱分解问题的求解
4.5.1.1 B(Z)的表示式
下面我们来分析一下谱分解的结果
因为解析函数B能够做如下的表示
B ( Z ) = b 0 + ∑ k = 1 ∞ b k Z k [ 1 ] B ∗ ( Z ) = b 0 + ∑ k = 1 ∞ b k ∗ ( Z ∗ ) k [ 2 ] B(Z) = b_0 + \sum_{k=1}^\infty b_k Z^k \quad\quad [1] \\ B^*(Z) = b_0 + \sum_{k=1}^\infty b_k^* (Z^*)^k \quad\quad [2] B(Z)=b0+k=1∑∞bkZk[1]B∗(Z)=b0+k=1∑∞bk∗(Z∗)k[2]
功率谱的表示
S Z ( ω ) = B ( ω ) B ∗ ( ω ) l o g S ( ω ) = l o g ( B ( ω ) ) + l o g ( B ∗ ( ω ) ) S_Z(\omega) = B(\omega) B^*(\omega) \\ logS(\omega) = log(B(\omega)) + log(B^*(\omega)) SZ(ω)=B(ω)B∗(ω)logS(ω)=log(B(ω))+log(B∗(ω))
B(ω)就是B(Z)在单位圆上的取值
我们假设
l o g B ( Z ) = d 0 + ∑ k = 1 ∞ d k Z k [ 3 ] log B(Z) = d_0 + \sum_{k=1}^\infty d_k Z^k \quad\quad [3] logB(Z)=d0+k=1∑∞dkZk[3]
则
B ( Z ) = e x p ( d 0 + ∑ k = 1 ∞ d k Z k ) = e x p ( d 0 ) e x p ( ∑ k = 1 ∞ d k Z k ) [ 4 ] B(Z) = exp(d_0 + \sum_{k=1}^\infty d_k Z^k) \\ = exp(d_0) exp(\sum_{k=1}^\infty d_k Z^k) \quad\quad [4] B(Z)=exp(d0+k=1∑∞dkZk)=exp(d0)exp(k=1∑∞dkZk)[4]
对后面式子进行talor展开
B ( Z ) = e x p ( d 0 ) ( 1 + ∑ k = 1 ∞ h k Z k ) = e x p ( d 0 ) + ∑ k = 1 ∞ h k ~ Z k [ 5 ] B(Z) = exp(d_0) (1+\sum_{k=1}^\infty h_k Z^k) \\ = exp(d_0) + \sum_{k=1}^\infty \widetilde{h_k} Z^k \quad\quad [5] B(Z)=exp(d0)(1+k=1∑∞hkZk)=exp(d0)+k=1∑∞hk Zk[5]
比较[1]和[5]得到
$$
\begin{cases}
B(Z) = b_0 + \sum_{k=1}^\infty b_k Z^k \quad\quad \
\
B(Z) = exp(d_0) + \sum_{k=1}^\infty \widetilde{h_k} Z^k
\end{cases} \
\Rightarrow b_0 = exp(d_0) \quad\quad [6]
$$
4.5.1.2 B(Z)的求解
下面,我们只要把B(Z)求解出来,其实就完成了谱分解
求解一些积分
1 2 π ∫ − π π l o g S ( ω ) d ω = 1 2 π ∫ − π π l o g ( B ( ω ) ) d ω + 1 2 π ∫ − π π l o g ( B ∗ ( ω ) ) d ω [ 7 ] \frac{1}{2\pi} \int_{-\pi}^{\pi} log S(\omega) d\omega \\ = \frac{1}{2\pi} \int_{-\pi}^{\pi} log(B(\omega)) d\omega + \frac{1}{2\pi} \int_{-\pi}^{\pi}log(B^*(\omega)) d\omega \quad\quad [7] 2π1∫−ππlogS(ω)dω=2π1∫−ππlog(B(ω))dω+2π1∫−ππlog(B∗(ω))dω[7]
[3]代入[7]可得
1 2 π ∫ − π π l o g S ( ω ) d ω = 1 2 π ∫ − π π ( d 0 + ∑ k = 1 ∞ d k e x p ( j ω k ) ) d ω + 1 2 π ∫ − π π ( d 0 + ∑ k = 1 ∞ d k ∗ e x p ( − j ω k ) ) d ω [ 8 ] \frac{1}{2\pi} \int_{-\pi}^{\pi} log S(\omega) d\omega \\ = \frac{1}{2\pi} \int_{-\pi}^{\pi} (d_0 + \sum_{k=1}^\infty d_k exp(j \omega k) ) d\omega+ \frac{1}{2\pi} \int_{-\pi}^{\pi}(d_0 + \sum_{k=1}^\infty d_k^* exp(-j \omega k)) d\omega \quad\quad[8] 2π1∫−ππlogS(ω)dω=2π1∫−ππ(d0+k=1∑∞dkexp(jωk))dω+2π1∫−ππ(d0+k=1∑∞dk∗exp(−jωk))dω[8]
由于复指函数在[-pi,pi]积分区间积分为0,因此,可以得到
1 2 π ∫ − π π l o g S ( ω ) d ω = 2 d 0 [ 9 ] \frac{1}{2\pi} \int_{-\pi}^{\pi} log S(\omega) d\omega = 2d_0 \quad\quad[9] 2π1∫−ππlogS(ω)dω=2d0[9]
可得
d 0 = 1 4 π ∫ − π π l o g S ( ω ) d ω d_0 =\frac{1}{4\pi} \int_{-\pi}^{\pi} log S(\omega) d\omega d0=4π1∫−ππlogS(ω)dω
由[6]式
$$
b_0 = exp(d_0) \
= exp(\frac{1}{4\pi} \int_{-\pi}^{\pi} log S(\omega) d\omega ) \quad\quad[10]
$$
我们其实就完成了对谱分解的解析
B ( Z ) = b 0 + ∑ k = 1 ∞ b k Z k B ∗ ( Z ) = b 0 + ∑ k = 1 ∞ b k ∗ ( Z ∗ ) k b 0 = e x p ( 1 4 π ∫ − π π l o g S ( ω ) d ω ) B(Z) = b_0 + \sum_{k=1}^\infty b_k Z^k \\ B^*(Z) = b_0 + \sum_{k=1}^\infty b_k^* (Z^*)^k \\ b_0 = exp(\frac{1}{4\pi} \int_{-\pi}^{\pi} log S(\omega) d\omega ) B(Z)=b0+k=1∑∞bkZkB∗(Z)=b0+k=1∑∞bk∗(Z∗)kb0=exp(4π1∫−ππlogS(ω)dω)
4.6 复变函数分析法与线性代数分析法的统一
4.6.1 两个结果
我们之前通过线性代数分析法,得到了线性预测的最优误差的计算结果。通过,基于复变函数的分析法,得到了谱分解的计算结果
- 线性代数分析结果
E ∣ ϵ P ∣ 2 = r 0 − b H b = c 2 = ∣ L 00 ( P + 1 ) ∣ 2 E| \epsilon_P|^2 = r_0 - b^H b=c^2 = |L^{(P+1)}_{00}|^2 E∣ϵP∣2=r0−bHb=c2=∣L00(P+1)∣2
- 复变函数分析结果
B ( Z ) = b 0 + ∑ k = 1 ∞ b k Z k B ∗ ( Z ) = b 0 + ∑ k = 1 ∞ b k ∗ ( Z ∗ ) k b 0 = e x p ( 1 4 π ∫ − π π l o g S ( ω ) d ω ) B(Z) = b_0 + \sum_{k=1}^\infty b_k Z^k \\ B^*(Z) = b_0 + \sum_{k=1}^\infty b_k^* (Z^*)^k \\ b_0 = exp(\frac{1}{4\pi} \int_{-\pi}^{\pi} log S(\omega) d\omega ) B(Z)=b0+k=1∑∞bkZkB∗(Z)=b0+k=1∑∞bk∗(Z∗)kb0=exp(4π1∫−ππlogS(ω)dω)
4.6.2 Kolmogorov 等式
在这里kolmogorov做出了一项伟大的工作,统一了复变函数分析与线性代数分析的结果。这个结论不进行证明
kolmogorov l i m P → ∞ ∣ L 00 ( P ) ∣ 2 = b 0 2 l i m P → ∞ ∣ L k m ( P ) ∣ 2 = b k − m 2 \text{kolmogorov} lim_{P \rightarrow \infty} |L^{(P)}_{00}|^2 = b_0^2 \\ lim_{P \rightarrow \infty} |L^{(P)}_{km}|^2 = b_{k-m}^2 kolmogorovlimP→∞∣L00(P)∣2=b02limP→∞∣Lkm(P)∣2=bk−m2
这个结论说明,当阶数P趋近于无穷大的时候,相关矩阵的LU分解的L矩阵的第(k,m)个元素,对应了谱分解泰勒展开的系数。
这个结论统一了线性代数分析结果和谱分析结果。既提供了一种分析随机过程线性预测问题最优误差的方法,又提供了一种对功率谱做谱分解的方法。
4.6.3 Kolmogorov-Szago 等式
同时,我们也就证明了随机过程的线性预测误差的极限,就是Kolmogorov-Szago等式
E ∣ ϵ P ∣ 2 = ∣ L 00 ( P + 1 ) ∣ 2 = b 0 2 = [ e x p ( 1 4 π ∫ − π π l o g S ( ω ) d ω ) ] 2 = e x p ( 1 2 π ∫ − π π l o g S ( ω ) d ω ) E| \epsilon_P|^2 = |L^{(P+1)}_{00}|^2 = b_0^2 = [exp(\frac{1}{4\pi} \int_{-\pi}^{\pi} log S(\omega) d\omega )]^2 = exp(\frac{1}{2\pi} \int_{-\pi}^{\pi} log S(\omega) d\omega ) E∣ϵP∣2=∣L00(P+1)∣2=b02=[exp(4π1∫−ππlogS(ω)dω)]2=exp(2π1∫−ππlogS(ω)dω)
5. 总结
这一节中,我们针对随机过程线性预测问题到底能做到什么程度,展开了讨论。
我们首先讨论了相关矩阵的正定性问题。并从相关矩阵的正定性问题中,引入了完美预测与非完美预测问题。因为相关矩阵是非负定的,可能存在奇异的情况。一旦相关矩阵是奇异的,得到的功率谱一定是线谱。也就是说明信号是谐波信号。因此,一定能够得到信号的完美预测。
( 1 ) R Z ( P ) ≥ 0 ⇔ Line Spectrum ⇔ Harmonic Signal ⇒ Perfect Prediction ⇒ ϵ P = 0 (1) \quad R_Z^{(P)} \geq 0 \Leftrightarrow \text{Line Spectrum} \Leftrightarrow \text{Harmonic Signal} \\ \Rightarrow \text{Perfect Prediction} \Rightarrow \epsilon_P = 0 (1)RZ(P)≥0⇔Line Spectrum⇔Harmonic Signal⇒Perfect Prediction⇒ϵP=0
然后,我们基于相关矩阵是非奇异的情况,研究了非完美预测问题的最优误差。非完美预测需要满足paley-wiener条件。如果功率谱足够宽,功率谱一定能够表示为谱分解的形式
( 2 ) R Z ( P ) > 0 ⇒ ( P a l e y − W i e n e r ) ⇒ 1 2 π ∫ − π + π l o g S Z ( ω ) d ω > − ∞ ⇒ S Z ( ω ) = ∣ B ( ω ) ∣ 2 (2) \quad R_Z^{(P)} > 0 \Rightarrow (Paley-Wiener) \Rightarrow \frac{1}{2 \pi} \int _{-\pi}^{+\pi} log S_Z(\omega) d\omega > -\infty \\ \Rightarrow S_Z(\omega) = |B(\omega)|^2 (2)RZ(P)>0⇒(Paley−Wiener)⇒2π1∫−π+πlogSZ(ω)dω>−∞⇒SZ(ω)=∣B(ω)∣2
然后,我们用复变函数的分析方法,解析了谱分解得到的B函数的具体表示方法
( 3 ) B ( ω ) = b 0 + ∑ k = 1 ∞ b k e x p ( j ω ) ⇒ b 0 = e x p ( 1 4 π ∫ − π π l o g S ( ω ) d ω ) (3) \quad B(\omega) = b_0 + \sum_{k=1}^\infty b_k exp(j \omega) \\ \Rightarrow b_0 = exp(\frac{1}{4\pi} \int_{-\pi}^{\pi} log S(\omega) d\omega ) (3)B(ω)=b0+k=1∑∞bkexp(jω)⇒b0=exp(4π1∫−ππlogS(ω)dω)
随后,我们又用线性代数分析法,相关矩阵的分解式以及三个跟递推有关的等式
( 4 ) R Z ( P ) = ( L ( P ) ) H ( L ( P ) ) L ( P ) = ( c 0 b A ) { A = L ( P ) b = ( ( L ( P ) ) H ) − 1 ∗ ( r 1 ∗ . . . r P ∗ ) c = ( r 0 − b H b ) 1 2 (4) \quad R_Z^{(P)} = (L^{(P)})^H (L^{(P)}) \\ L^{(P)} = \begin{pmatrix} c & 0 \\ b & A \end{pmatrix} \\ \begin{cases} A = L^{(P)} \\ b= ((L^{(P)})^H)^{-1} *\begin{pmatrix} r_{1}^* \\ ... \\ r_{P}^* \end{pmatrix} \\ c = (r_0 - b^H b)^{\frac{1}{2}} \end{cases} (4)RZ(P)=(L(P))H(L(P))L(P)=(cb0A)⎩⎪⎪⎪⎪⎨⎪⎪⎪⎪⎧A=L(P)b=((L(P))H)−1∗⎝⎛r1∗...rP∗⎠⎞c=(r0−bHb)21
而后我们发现,随机过程线性预测问题的误差均方值的极限,就是LU分解最左上角的值
( 5 ) E ∣ ϵ P ∣ 2 = ∣ L 00 ( P + 1 ) ∣ 2 (5) \quad E| \epsilon_P|^2 = |L^{(P+1)}_{00}|^2 (5)E∣ϵP∣2=∣L00(P+1)∣2
而后,我们通过kolmogorov等式,不但统一了复变函数分析得到的关系式和线性代数分析得到的关系式。而且为功率谱密度的谱分解提供了一种解决方案。
( 6 ) l i m P → ∞ ∣ L 00 ( P ) ∣ 2 = b 0 2 l i m P → ∞ ∣ L k m ( P ) ∣ 2 = b k − m 2 (6) \quad lim_{P \rightarrow \infty} |L^{(P)}_{00}|^2 = b_0^2 \\ lim_{P \rightarrow \infty} |L^{(P)}_{km}|^2 = b_{k-m}^2 (6)limP→∞∣L00(P)∣2=b02limP→∞∣Lkm(P)∣2=bk−m2
最后,我们证明了,Kolmogorov-Szago等式,得到了非完美预测问题误差均方值的极限
( 7 ) E ∣ ϵ P ∣ 2 = e x p ( 1 2 π ∫ − π π l o g S ( ω ) d ω ) (7) \quad E| \epsilon_P|^2 = exp(\frac{1}{2\pi} \int_{-\pi}^{\pi} log S(\omega) d\omega ) (7)E∣ϵP∣2=exp(2π1∫−ππlogS(ω)dω)
就此,我们就完成了宽平稳随机过程线性预测能力的评判