【随机过程】19 - 随机过程的线性预测问题

随机过程的线性预测问题

1. 随机过程的估计问题概述

 我们假设X是随机变量构成的概率空间,Xs是X的子空间,Xs是已知的。而Y是存在于X中的一个未知的随机变量,估计的目的就是为了利用Xs中的信息对Y做出统计推断,一般来说有这样的情况

1.1 预测问题

  利用过去和现在的已知信息,对未来进行估计,这类问题叫做预测问题,表达为

Prediction \text{Prediction} Prediction

X = { X ( t ) , t ∈ R } X s = { X ( u ) , u ≤ t } Y = ( t + τ ) X s ⇒ Y X = \{ X(t),t\in \R\}\\ X_s = \{ X(u),u\leq t\} \\ Y = (t+\tau) \\ X_s \Rightarrow Y X={X(t),tR}Xs={X(u),ut}Y=(t+τ)XsY

1.2 内插问题

  而已知随机过程两端的信息,来估计中间时刻的未知信息,叫做内插问题

Interpolation \text{Interpolation} Interpolation

X = { X ( t ) , t ∈ R } X s = { X ( u ) , u = A ∪ B , A ∈ ( − ∞ , t ) , B ∈ ( t , + ∞ ) } Y = X ( t ) X = \{ X(t),t\in \R\} \\ X_s = \{X(u),u=A \cup B,A\in(-\infty,t),B\in(t,+\infty) \} \\ Y = X(t) X={X(t),tR}Xs={X(u),u=AB,A(,t),B(t,+)}Y=X(t)

1.3 滤波问题

  一个原本的随机过程在观测的时候会引入噪声,利用带有噪声的观测去估计原本的随机过程叫做滤波问题

model { X ( t ) , t ∈ R } { W ( t ) , t ∈ R } X ⊃ { X ( t ) , t ∈ R } ∪ { W ( t ) , t ∈ R } Z ( t ) = X ( t ) + W ( t ) subspace X s = { Z ( u ) , u ≤ t } Y = X ( t ) \text{model} \{ X(t),t\in \R\} \quad \{ W(t),t\in \R\} \\ X \supset \{ X(t),t\in \R\} \cup \{ W(t),t\in \R\} \\ Z(t) = X(t)+W(t) \text{subspace} \\ X_s = \{ Z(u),u\leq t\} \\ Y = X(t) model{X(t),tR}{W(t),tR}X{X(t),tR}{W(t),tR}Z(t)=X(t)+W(t)subspaceXs={Z(u),ut}Y=X(t)

  这里我们希望观测的是X,但是实际测量的是Z,Z是X叠加一个加性噪声的结果。我们希望利用测量数据估计出原本的状态X

2. 随机过程的可预测性

  首先,我们需要搞清楚的第一个问题是,什么样的随机过程是可预测的

2.1 新息过程

2.1.1 信息过程的定义

  我们定义这样的参数

  假设有一个离散的随机过程X

{ X k , k = 1 , . . . , n } ∼ Discrete Stochastic Processes \{X_k,k=1,...,n\} \sim \text{Discrete Stochastic Processes} {Xk,k=1,...,n}Discrete Stochastic Processes

  并且定义{X1,…,Xk}构成的子空间是LkX

L k X = s p a n { X 1 , . . . , X k } ⇒ L n − 1 X = s p a n { X k , k ≤ k − 1 } L_k^X = span\{X_1,...,X_k\} \\ \Rightarrow L^X_{n-1} = span\{X_k,k\leq k-1\} LkX=span{X1,...,Xk}Ln1X=span{Xk,kk1}

  事实上,随机过程X的信息可能并不是一次获得的,而是随着时间的变化,逐步获取的。也就是这些随机变量Xk构成的子空间LXk也是在不断扩充的。因此,我们每次采样得到的一个数据Xn是可以分解为旧空间信息和新空间信息的

X n = ( X n ∣ L n − 1 X ) + ( X n − X n ∣ L n − 1 X ) X_n = (X_n|L_{n-1}^X) + (X_n - X_n|L_{n-1}^X) Xn=(XnLn1X)+(XnXnLn1X)

  其中前一部分数据表示新数据在旧空间上的投影,而后一部分数据在旧空间LXn-1的正交补空间的投影。如果后面一部分不为0的话,说明新的数据含有新的信息,而后数据子空间会得到扩展。

  我们把数据Xn在旧空间的正交补中的投影定义为In,这个变量表示新数据是否含有新的子空间信息

I n = X n − X n ∣ L n − 1 X I_n = X_n - X_n|L_{n-1}^X In=XnXnLn1X

  因此,新的采样数据可以表示为两个空间投影的和

X n = X n ∣ L n − 1 X + I n X_n = X_n|L_{n-1}^X +I_n Xn=XnLn1X+In

  由于In是真正含有新信息的部分,因此,我们称In是Xn的新息过程

innovations process \text{innovations process} innovations process

2.1.2 估计的子空间分解

  因此,对于我们的估计Y,也可以有新的表示形式。如果我们要用X1,…Xn去估计Y,得到的其实就是Y在子空间LXn上的投影

Y ^ = Y ∣ L n X \hat Y = Y | L^X_n Y^=YLnX

  而我们知道投影具有这样的性质,如果子空间能够分解为两个正交的子空间,那么在原来子空间上的投影,就等于在两个新正交子空间上的投影的和

Subspace  W , U , V W = U + V U ⊥ V ⇒ Y ∣ W = Y ∣ U + Y ∣ V \text{Subspace } W,U,V \\ W = U+V \quad U \perp V \\ \Rightarrow Y|W = Y|U+Y|V Subspace W,U,VW=U+VUVYW=YU+YV

  因为LXn-1与In具有正交补的关系,因此Y的估计也可以做子空间分解

Y ^ = Y ∣ L n X = Y ∣ L n − 1 X + Y ∣ I n \hat Y = Y | L^X_n = Y | L^X_{n-1} +Y|I_n Y^=YLnX=YLn1X+YIn

2.1.3 新息过程的性质

  下面介绍信息过程的性质,最重要的是,新息过程是原来采样数据所构成空间的正交基,新息过程和采样数据张成同一个子空间。

s p a n { I k , k ≤ n } = s p a n { X k , k ≤ n } n ∈ Z < I k , I m > = E ( I k I m ‾ ) = 0 k , m ∈ Z k = m span\{I_k,k\leq n\} = span\{X_k,k\leq n\} \quad n\in \Z \\ <I_k,I_m> = E(I_k \overline{I_m}) = 0 \quad k,m\in \Z \quad k \cancel=m span{Ik,kn}=span{Xk,kn}nZ<Ik,Im>=E(IkIm)=0k,mZk= m

  我们可以证明一下这个事情

L n X = L n − 1 X + X n = L n − 1 X + X n ∣ L n − 1 X + I n = L n − 1 X ⨁ I n = L n − 2 X ⨁ I n − 1 ⨁ I n = ⨁ k = − ∞ + ∞ I k L^X_{n} = L^{X}_{n-1} +X_n \\ = L^{X}_{n-1} +X_n|L^{X}_{n-1} +I_n \\ = L^{X}_{n-1}\bigoplus I_n \\ = L^{X}_{n-2}\bigoplus I_{n-1} \bigoplus I_n \\ = \bigoplus_{k=-\infty}^{+\infty}I_k LnX=Ln1X+Xn=Ln1X+XnLn1X+In=Ln1XIn=Ln2XIn1In=k=+Ik

2.2 随机过程的正则性与奇异性

2.2.1 正则性和奇异性概述

  我们对一个随机过程进行预测的时候,可能会考虑这样的问题,对其进行预测的时候,误差有可能达到0吗?如果误差不能达到0,最小可以达到什么程度呢?

  对于这个问题,我们先从直观角度进行考虑,如果一个随机过程每次采样得到的新数据都能完全投影到旧的子空间上,正交补空间没有任何新的信息,那么说明随着时间的推移,这个随机过程不会引入新的信息。那么就说明对这个随机过程的预测会完全没有误差,因为其数据空间在开始的时候就已经确定了。

  而如果每一个新数据都不能投影到原来的数据子空间上,而是全部都在正交补空间上,那么就说明这个数据是不可能准确预测的,因为数据之间是完全没有相关性的,其预测误差就取决于数据的方差大小了。

  因此,如果随机过程的随机性不随着时间发生变化,就是完全可以预测的,就叫做奇异过程。而如果随机过程的随机性完全随机时间变化,数据之前完全不相关,那么这个随机过程就是不可能准确预测的,这样的随机过程就叫做正则过程

Singular Purely deterministic \text{Singular} \quad \text{Purely deterministic} \\ SingularPurely deterministic

在这里插入图片描述

Regular Purely non-deterministic \text{Regular} \quad \text{Purely non-deterministic} RegularPurely non-deterministic

在这里插入图片描述

  然后我们给这个奇异性和正则性进行数学上的定义。

  我们知道采样数据构成的子空间具有包含关系,后面的子空间一定是包含前面的子空间的,也就是

L n X ⊃ L n − 1 X ⊃ L n − 2 X ⊃ . . . L_{n}^X \supset L_{n-1}^X \supset L_{n-2}^X \supset ... LnXLn1XLn2X...

  这个子空间是有下界和上界的

l i m n → − ∞ L n X = L − ∞ X l i m n → + ∞ L n X = L + ∞ X lim_{n\rightarrow -\infty}L^X_{n} = L^X_{-\infty} \\ lim_{n\rightarrow +\infty}L^X_{n} = L^X_{ +\infty} limnLnX=LXlimn+LnX=L+X

  并且显然,下界是所有子空间的交集,而上界是所有子空间的并集

L − ∞ X = ⋂ k = − ∞ k = + ∞ L k X L + ∞ X = ⋃ k = − ∞ k = + ∞ L k X L^X_{ -\infty} = \bigcap_{k=-\infty}^{k=+\infty} L^X_k \\ L^X_{ +\infty} = \bigcup_{k=-\infty}^{k=+\infty} L^X_k \\ LX=k=k=+LkXL+X=k=k=+LkX

  当这个子空间是上界和下界相等的时候,说明随机性不随着时间改变,就是奇异的

Singular ⇒ L − ∞ X = L + ∞ X \text{Singular} \Rightarrow L^X_{ -\infty} = L^X_{ +\infty} SingularLX=L+X

  而如果下界是0,说明随机性是完全随机时间改变的,就是正则的

Regular ⇒ L − ∞ X = { 0 } \text{Regular} \Rightarrow L^X_{ -\infty} = \{0 \} RegularLX={0}

2.2.2 奇异性的充要条件

  下面从误差的角度,给出奇异性和正则性的充要条件。不过这个定义比较paperwork,不具有解析意义

  首先给出随机过程线性预测的估计误差,我们假设当前时刻为n,我们要预测第n+m个时刻的数据Xn+m,最优线性预测就是Xn+m在子空间LXn的投影

ϵ ( m ) = ∣ ∣ X n + m − X n + m ∣ L n X ∣ ∣ = ( E ( ∣ X n + m − X n + m ∣ L n X ∣ 2 ) ) 1 2 \epsilon(m) = ||X_{n+m} - X_{n+m}|L_n^X|| = (E(|X_{n+m} - X_{n+m}|L_n^X|^2))^{\frac{1}{2}} ϵ(m)=Xn+mXn+mLnX=(E(Xn+mXn+mLnX2))21

  利用这个误差,可以给出奇异过程的充要条件。

{ X k , k ∈ Z } X ∼ W . S . S . X is Singular ⇔ ∃ m 0 ϵ ( m 0 ) = 0 \{X_k,k\in \Z \} \quad X \sim W.S.S. \text{X is Singular} \Leftrightarrow \exist m_0 \quad \epsilon(m_0) = 0 {Xk,kZ}XW.S.S.X is Singularm0ϵ(m0)=0

  首先证明必要性

X is Singular ⇒ X n + m = X n + m ∣ L n X ⇒ ϵ ( m ) = 0 ∀ m ∈ Z \text{X is Singular} \Rightarrow X_{n+m} = X_{n+m}|L_n^X \\ \Rightarrow \epsilon(m) = 0 \quad \forall m \in \Z X is SingularXn+m=Xn+mLnXϵ(m)=0mZ

  然后证明充分性

∃ m 0 ϵ ( m 0 ) = 0 ⇒ ϵ ( m ) = 0 m ≤ m 0 X n + m = X n + m ∣ L n X ∈ L n X ∀ n ∈ Z , m ≤ m 0 ⇒ ∀ m ≤ n X n ∈ L m X ⇒ L n X ∈ L m X ⇒ L − ∞ X = L + ∞ X \exist m_0\quad \epsilon(m_0) = 0 \\ \Rightarrow \epsilon(m) = 0 \quad m \leq m_0 \\ X_{n+m} = X_{n+m}|L_n^X \in L_n^X \quad \forall n \in \Z,m \leq m_0 \\ \Rightarrow \forall m\leq n \quad X_n \in L_m^X \Rightarrow L_n^X \in L_m^X \Rightarrow L^X_{-\infty}=L^X_{+\infty} m0ϵ(m0)=0ϵ(m)=0mm0Xn+m=Xn+mLnXLnXnZ,mm0mnXnLmXLnXLmXLX=L+X

2.2.3 正则性的充要条件

  然后也可以利用估计误差来定义正则性的充要条件

{ X k , k ∈ Z } X ∼ W . S . S . X is Regular ⇔ ϵ ( m ) = R X ( 0 ) , m → ∞ \{X_k,k\in \Z \} \quad X \sim W.S.S. \text{X is Regular} \Leftrightarrow \epsilon(m) = R_X(0),m\rightarrow \infty {Xk,kZ}XW.S.S.X is Regularϵ(m)=RX(0),m

  证明

E ( X n 2 ) = E ( X n ∣ L m X + X n − X n ∣ L m X ) 2 = E ( X n ∣ L m X ) 2 + E ( X n − X n ∣ L m X ) 2 + 0 = E ( X n ∣ L m X ) 2 + ϵ ( n − m ) m → − ∞ ⇒ R X ( 0 ) = E ( X n ∣ L − ∞ X ) 2 + ϵ ( ∞ ) ϵ ( ∞ ) = R X ( 0 ) ⇔ X n ∣ L − ∞ = { 0 } , ∀ n ∈ Z ⇔ L − ∞ = { 0 } E(X_n^2) = E(X_n|L_m^X+X_n-X_n|L_m^X)^2 = E(X_n|L_m^X)^2+E(X_n-X_n|L_m^X)^2 +0 \\ = E(X_n|L_m^X)^2+ \epsilon(n-m) \\ m \rightarrow -\infty \Rightarrow R_X(0) = E(X_n|L_{-\infty}^X)^2+\epsilon(\infty)\\ \epsilon(\infty) = R_X(0) \Leftrightarrow X_n|L_{-\infty} = \{0\}, \forall n \in \Z \Leftrightarrow L_{-\infty} = \{0\} E(Xn2)=E(XnLmX+XnXnLmX)2=E(XnLmX)2+E(XnXnLmX)2+0=E(XnLmX)2+ϵ(nm)mRX(0)=E(XnLX)2+ϵ()ϵ()=RX(0)XnL={0},nZL={0}

  因此,如果一个随机过程是奇异的,一定有某种预测方法,可以使得预测误差为0。而如果一个随机过程是正则的,不论采取何种预测方法,预测误差都与随机过程自身的波动起伏有关,没有办法给出更加清晰的预测。

2.3 正则性的解析判据–Paley-Wiener条件

Paley-Wiener \text{Paley-Wiener} Paley-Wiener

  刚才给出的正则性和奇异性的定义在实际解析计算中用处不大,这里给出一个正则过程的充分必要条件

  对于离散的宽平稳随机过程Xn,判断其为正则过程的判据为

∫ − ∞ + ∞ l o g S ( ω ) d ω > − ∞ \int_{-\infty}^{+\infty} log S(\omega) d\omega >-\infty +logS(ω)dω>

  这个判据的证明在信号处理一的笔记中有。

  这里加入注解,我们知道,如果随机过程归一化了,其功率谱的积分是1。如果对功率谱的log积分是负无穷,比如说明大部分频点能量很低,也就是说功率谱能量会集中在零点。功率谱集中,必然说明相关函数是很分散,各个点之间都有比较强的相关性,就有办法预测。如果这个积分不是负无穷,就说明功率谱比较分散,相关函数就比较集中,各个点之间没有相关性,就难以预测,就越接近正则。

在这里插入图片描述

  对于连续的宽平稳随机过程Xn,判断其为正则过程的判据为

∫ − ∞ + ∞ ∣ l o g S ( ω ) ∣ 1 + ω 2 d ω < ∞ \int_{-\infty}^{+\infty} \frac{|log S(\omega)|}{1+\omega^2} d \omega <\infty +1+ω2logS(ω)dω<

2.4 随机过程的分解

  奇异过程和正则过程是事物的两个极端,一个代表完全可以预测,另外一个代表几乎无法预测。大部分宽平稳随机过程介于这两者之间。一般随机过程中既包含奇异分量,也包含正则的分量。因此,随机过程是可以做分解的。

2.4.1 Wold分解
  • 离散形式

  假设X(t)是零均值二阶矩随机过程,X可以做如下分解

{ X n , n ∈ Z } X n = Y n + Z n Y n ∼ Regular Process Z n ∼ Singular Process \{X_n,n\in \Z\} \\ X_n = Y_n + Z_n \\ Y_n \sim \text{Regular Process} \\ Z_n \sim \text{Singular Process} {Xn,nZ}Xn=Yn+ZnYnRegular ProcessZnSingular Process

  对于任意的m,n,Ym和Zn正交

  不过Wold分解只是证明了分解的存在性,并没有给出求解随机过程结构参数的方法。而且是在最优线性预测意义下定义的

  其中正则部分可以做如下表示

X n = ∑ k = 0 ∞ h k U n − k X_n = \sum_{k=0}^{\infty}h_k U_{n-k} Xn=k=0hkUnk

  • 连续形式

X ( t ) = ξ ( t ) + η ( t ) X(t) = \xi(t) + \eta(t) X(t)=ξ(t)+η(t)

  其中分解得到的两部分是互相正交的随机过程,一个是奇异过程,另外一个是正则过程

  其中正则部分可以做如下表示

η ( t ) = ∫ − ∞ + ∞ h ( t − τ ) d U ( τ ) \eta(t) = \int_{-\infty}^{+\infty}h(t-\tau) dU(\tau) η(t)=+h(tτ)dU(τ)

2.4.2 Doob-Meyer分解

  设随机过程X

{ X n , n ∈ Z + } E ∣ X n ∣ < ∞ ⇒ X n = Y n + Z n \{X_n,n\in Z_+ \} \quad E|X_n| <\infty \\ \Rightarrow X_n = Y_n+Z_n {Xn,nZ+}EXn<Xn=Yn+Zn

  其中Y和Z满足这样的关系

E ( Δ Y n ∣ X n , X n − 1 , . . . , X 0 ) = Δ Y n E ( Δ Z n ∣ X n , X n − 1 , . . . , X 0 ) = 0 E(\Delta Y_n|X_n,X_{n-1},...,X_0) = \Delta Y_n \\ E(\Delta Z_n | X_n,X_{n-1},...,X_0) = 0 E(ΔYnXn,Xn1,...,X0)=ΔYnE(ΔZnXn,Xn1,...,X0)=0

  也就是Y的增量可以被X精确的预测出来,Z的增量无法进行预测。

  这个分解更加一般,是在非线性最优估计条件下的分解,不过得到的条件期望会丧失线性估计的各种优势,并且分解的两个部分不再具有正交关系。由于这些原因,人们往往更加关注线性预测。

3. 随机过程的白化与误差分析

  通过刚才的介绍,我们了解到,正则过程是几乎无法预测的,也就一定是有误差的,而奇异过程是可以做到无误差预测的。并且随机过程是可以基于Wold分解,分解成奇异部分和正则部分的。因此,对随机过程的预测误差主要来自于其中的正则部分,如果我们有办法得到正则部分的线性预测误差,我们就可以得到这个随机过程的线性预测误差。

  因此,接下来,我们会介绍如何分析一个正则过程的预测误差。

3.2 谱分布函数

  引入谱分布函数的概念。首先,把一个随机过程进行谱表示,我们假设随机过程Xn是一个离散的宽平稳随机过程,其谱表示具有这样的关系

{ X n , n ∈ Z } X n = ∫ − π + π e x p ( j ω n ) d Z X ( ω ) E ∣ Z X ∣ 2 = 1 2 π S X ( ω ) d ω \{ X_n,n \in \Z \} X_n = \int_{-\pi}^{+\pi} exp(j\omega n) d Z_X(\omega) E|Z_X|^2 = \frac{1}{2\pi}S_X(\omega) d\omega {Xn,nZ}Xn=π+πexp(jωn)dZX(ω)EZX2=2π1SX(ω)dω

  而谱分布定义为

d F X ( ω ) = E ∣ Z X ∣ 2 ⇒ d F X ( ω ) = 1 2 π S X ( ω ) d ω dF_X(\omega) = E|Z_X|^2 \Rightarrow dF_X(\omega) = \frac{1}{2\pi}S_X(\omega) d\omega dFX(ω)=EZX2dFX(ω)=2π1SX(ω)dω

3.1 滑动平均表示

  滑动平均表示是一种白化的手段。也就是说,随机过程是可以表示成为白噪声通过线性系统的响应的,一旦变成了白噪声,各个分量之间就是不相关的了,这样其实也实现了随机过程各个分量的解耦。

  如果一个离散的宽平稳随机过程X的谱分布函数足够光滑,这个宽平稳随机过程就可以表示为白噪声的滑动平均表示。

d F X ( ω ) = 1 2 π S X ( ω ) d ω dF_X(\omega) = \frac{1}{2\pi}S_X(\omega) d\omega dFX(ω)=2π1SX(ω)dω

  即上式成立,也就是SX(ω)存在的时候。离散宽平稳随机过程Xn具有如下表示

X n = ∑ k = − ∞ + ∞ h k U n − k X_n = \sum_{k=-\infty}^{+\infty}h_k U_{n-k} Xn=k=+hkUnk

  其中

U k ∼ White Noise ∑ k ∣ h k ∣ 2 < ∞ U_k \sim \text{White Noise} \\ \sum_k |h_k|^2 < \infty UkWhite Noisekhk2<

  可以证明一下这个事情

X n = ∫ − π + π e x p ( j ω n ) d Z X ( ω ) ⇔ X n = ∑ k = − ∞ + ∞ h k U n − k X_n = \int_{-\pi}^{+\pi} exp(j\omega n) d Z_X(\omega) \Leftrightarrow X_n = \sum_{k=-\infty}^{+\infty}h_k U_{n-k} Xn=π+πexp(jωn)dZX(ω)Xn=k=+hkUnk

  首先写一个白噪声的谱表示

U n = ∫ − π + π e x p ( j ω n ) d Z U ( ω ) E ( Z U ( ω ) ) 2 = 1 2 π S U ( ω ) d ω = 1 2 π d ω U_n = \int_{-\pi}^{+\pi} exp(j\omega n) d Z_U(\omega) \\ E(Z_U(\omega))^2 = \frac{1}{2\pi}S_U(\omega) d\omega = \frac{1}{2\pi} d\omega Un=π+πexp(jωn)dZU(ω)E(ZU(ω))2=2π1SU(ω)dω=2π1dω

  先证明一下必要性

X n = ∑ k = − ∞ + ∞ h k U n − k = ∑ k = − ∞ + ∞ h k ∫ − π + π e x p ( j ω ( n − k ) ) d Z U ( ω ) = ∫ − π + π ∑ k = − ∞ + ∞ h k e x p ( − j ω k ) e x p ( j ω n ) d Z U ( ω ) = ∫ − π + π H ( ω ) e x p ( j ω n ) d Z U ( ω ) = ∫ − π + π e x p ( j ω n ) d H ( ω ) Z U ( ω ) = ∫ − π + π e x p ( j ω n ) d Z X ( ω ) X_n = \sum_{k=-\infty}^{+\infty}h_k U_{n-k} \\ =\sum_{k=-\infty}^{+\infty}h_k \int_{-\pi}^{+\pi} exp(j\omega (n-k)) d Z_U(\omega) \\ = \int_{-\pi}^{+\pi} \sum_{k=-\infty}^{+\infty}h_k exp(-j\omega k) exp(j\omega n) d Z_U(\omega) \\ = \int_{-\pi}^{+\pi} H(\omega)exp(j\omega n) d Z_U(\omega) \\ = \int_{-\pi}^{+\pi} exp(j\omega n) d H(\omega)Z_U(\omega) = \int_{-\pi}^{+\pi} exp(j\omega n) d Z_X(\omega) Xn=k=+hkUnk=k=+hkπ+πexp(jω(nk))dZU(ω)=π+πk=+hkexp(jωk)exp(jωn)dZU(ω)=π+πH(ω)exp(jωn)dZU(ω)=π+πexp(jωn)dH(ω)ZU(ω)=π+πexp(jωn)dZX(ω)

  其中X的谱分布函数可以表示为

E ∣ Z X ( ω ) ∣ 2 = E ∣ H ( ω ) Z U ( ω ) ∣ 2 = ∣ H ( ω ) ∣ 2 1 2 π d ω E|Z_X(\omega)|^2 = E|H(\omega)Z_U(\omega)|^2 = |H(\omega)|^2 \frac{1}{2\pi} d\omega EZX(ω)2=EH(ω)ZU(ω)2=H(ω)22π1dω

  再证明充分性

X n = ∫ − π + π e x p ( j ω n ) d Z X ( ω ) = ∫ − π + π e x p ( j ω n ) S X ( ω ) d Z X ( ω ) S X ( ω ) = ∫ − π + π e x p ( j ω n ) S X ( ω ) d Z X ( ω ) ~ E ( d Z X ( ω ) ~ ) 2 = 1 2 π d ω X_n = \int_{-\pi}^{+\pi} exp(j\omega n) d Z_X(\omega) \\ = \int_{-\pi}^{+\pi} exp(j\omega n) \sqrt{S_X(\omega)}d \frac{Z_X(\omega)}{\sqrt{S_X(\omega)}} \\ = \int_{-\pi}^{+\pi} exp(j\omega n) \sqrt{S_X(\omega)}d \widetilde{Z_X(\omega)} \\ E(d \widetilde{Z_X(\omega)} )^2 = \frac{1}{2\pi} d\omega Xn=π+πexp(jωn)dZX(ω)=π+πexp(jωn)SX(ω) dSX(ω) ZX(ω)=π+πexp(jωn)SX(ω) dZX(ω) E(dZX(ω) )2=2π1dω

  做傅里叶展开

S X ( ω ) = ∑ k = − ∞ + ∞ c k e x p ( − j ω k ) \sqrt{S_X(\omega)} = \sum_{k=-\infty}^{+\infty}c_k exp(-j\omega k) SX(ω) =k=+ckexp(jωk)

  代入

X n = ∫ − π + π e x p ( j ω n ) ∑ k = − ∞ + ∞ c k e x p ( − j ω k ) d Z X ( ω ) ~ = ∑ k = − ∞ + ∞ c k ∫ − π + π e x p ( j ω ( n − k ) ) d Z U ( ω ) = ∑ k = − ∞ + ∞ c k U n − k X_n = \int_{-\pi}^{+\pi} exp(j\omega n) \sum_{k=-\infty}^{+\infty}c_k exp(-j\omega k) d \widetilde{Z_X(\omega)} \\ = \sum_{k=-\infty}^{+\infty}c_k \int_{-\pi}^{+\pi} exp(j\omega(n-k)) d Z_U(\omega) \\ = \sum_{k=-\infty}^{+\infty}c_k U_{n-k} Xn=π+πexp(jωn)k=+ckexp(jωk)dZX(ω) =k=+ckπ+πexp(jω(nk))dZU(ω)=k=+ckUnk

3.2 具有因果性的滑动平均表示

  不过刚才的那个滑动平均表示不具备因果性,我们无法获取未来的时刻的数据。因此,如果宽平稳随机过程满足一定条件的话,可以表示为具有因果性的滑动平均表示

X n = ∑ k = 0 ∞ h k U n − k X_n = \sum_{k=0}^{\infty}h_k U_{n-k} Xn=k=0hkUnk

  其中

U k ∼ White Noise ∑ k ∣ h k ∣ 2 < ∞ U_k \sim \text{White Noise} \\ \sum_k |h_k|^2 < \infty UkWhite Noisekhk2<

  这个分解要求宽平稳随机过程是一个正则过程。也就是要求Xn功率谱满足paley-wiener条件

X n  Regular ⇒ ∫ − ∞ + ∞ l o g S ( ω ) d ω > − ∞ X_n \text{ Regular} \Rightarrow \int_{-\infty}^{+\infty} log S(\omega) d\omega >-\infty Xn Regular+logS(ω)dω>

  如果X是一个连续的宽平稳随机过程,满足paley-wiener条件也可以变成连续的滑动平均表示

X ( t )  Regular ⇒ ∫ − ∞ + ∞ ∣ l o g S ( ω ) ∣ 1 + ω 2 d ω < ∞ X(t) \text{ Regular} \Rightarrow \int_{-\infty}^{+\infty} \frac{|log S(\omega)|}{1+\omega^2} d \omega <\infty X(t) Regular+1+ω2logS(ω)dω<

X ( t ) = ∫ − ∞ + ∞ h ( t − τ ) d U ( τ ) X(t) = \int_{-\infty}^{+\infty}h(t-\tau) dU(\tau) X(t)=+h(tτ)dU(τ)

3.2 正则过程的误差分析

  然后我们可以基于滑动平滑表示分析正则过程的线性预测误差了。

  假设Xn是一个离散时间宽平稳随机过程,并且具有正则性。则对X的线性预测表示为

X ^ n + m = ∑ k = 0 ∞ d k X n − k \hat X_{n+m} = \sum_{k=0}^{\infty} d_kX_{n-k} X^n+m=k=0dkXnk

  满足

E ∣ X n + m − X ^ n + m ∣ 2 = m i n Y ∈ s p a n { X k , k ≤ n } E ∣ X n + m − Y ∣ 2 E| X_{n+m}-\hat X_{n+m}|^2 = min_{Y \in span \{ X_k,k\leq n\} } E|X_{n+m} - Y|^2 EXn+mX^n+m2=minYspan{Xk,kn}EXn+mY2

  由于X可以进行滑动平均表示,因此

s p a n { X k , k ≤ n } = s p a n { U k , k ≤ n } span \{ X_k,k\leq n\} = span \{ U_k,k\leq n\} span{Xk,kn}=span{Uk,kn}

  最优误差为

σ m 2 = E ∣ X m + n − X m + n ∣ s p a n ( X k , k ≤ n ) ∣ 2 = E ∣ X m + n − X m + n ∣ s p a n ( U k , k ≤ n ) ∣ 2 = E ∣ h 0 U n + m + h 1 U n + m − 1 + . . . + h m − 1 U n + 1 ∣ 2 = ∣ h 0 ∣ 2 + . . . + ∣ h m − 1 ∣ 2 \sigma_m^2 = E|X_{m+n} - X_{m+n}|span(X_k,k\leq n)|^2 \\ = E|X_{m+n} - X_{m+n}|span(U_k,k\leq n)|^2 \\ = E|h_0 U_{n+m}+h_1U_{n+m-1}+...+h_{m-1}U_{n+1}|^2 \\ = |h_0|^2 + ...+ |h_{m-1}|^2 σm2=EXm+nXm+nspan(Xk,kn)2=EXm+nXm+nspan(Uk,kn)2=Eh0Un+m+h1Un+m1+...+hm1Un+12=h02+...+hm12

  其中h0是傅里叶展开的首项。

  因此,一步预测误差可以表示为

σ 1 = h 0 = e x p ( ∫ − ∞ + ∞ l o g S ( ω ) d ω ) \sigma_1 = h_0 = exp(\int_{-\infty}^{+\infty} log S(\omega) d\omega) σ1=h0=exp(+logS(ω)dω)

  后面exp中的部分就是正则过程的判据,如果这个积分是负无穷,误差就是0,而正则过程要求这个积分是可积的,得到的就是一个具体的数字。证明了正则过程是不可预测的,一定会有预测误差的。

4. 随机过程的谱因式分解

  在前面,我们讨论了,随机过程是可以转换为白噪声通过线性系统的响应的。随机过程一旦转换成了白噪声,由于其正交性的一些优点,能够使得我们处理各种问题更加的简便。这一部分,我们来讨论这个线性系统是怎么求解出来的。

4.1 正则过程的谱因式分解

  从白噪声的谱表示开始。

U n = ∫ − ∞ + ∞ e x p ( j ω n ) d Z U ( ω ) U_n = \int_{-\infty}^{+\infty}exp(j\omega n) dZ_U(\omega) Un=+exp(jωn)dZU(ω)

  其中

E ( Z U ( ω ) ) 2 = d ω 2 π E(Z_U(\omega))^2 = \frac{d\omega}{2\pi} E(ZU(ω))2=2πdω

  如果Xn为正则过程,Xn可以表示为

X n = ∑ k = 0 ∞ h k U n − k = ∑ k = 0 ∞ h k ∫ − π + π e x p ( j ω ( n − k ) ) d Z U ( ω ) = ∫ − π + π ∑ k = 0 ∞ h k e x p ( − j ω k ) e x p ( j ω n ) d Z U ( ω ) = ∫ − π + π e x p ( j ω n ) d H ( ω ) Z U ( ω ) = ∫ − π + π e x p ( j ω n ) d Z X ( ω ) X_n = \sum_{k=0}^{\infty} h_k U_{n-k} \\ = \sum_{k=0}^{\infty}h_k \int_{-\pi}^{+\pi}exp(j\omega(n-k))dZ_U(\omega) \\ = \int_{-\pi}^{+\pi} \sum_{k=0}^{\infty}h_k exp(-j\omega k) exp(j\omega n)dZ_U(\omega) \\ = \int_{-\pi}^{+\pi} exp(j\omega n)dH(\omega) Z_U(\omega) \\ = \int_{-\pi}^{+\pi} exp(j\omega n)dZ_X(\omega) Xn=k=0hkUnk=k=0hkπ+πexp(jω(nk))dZU(ω)=π+πk=0hkexp(jωk)exp(jωn)dZU(ω)=π+πexp(jωn)dH(ω)ZU(ω)=π+πexp(jωn)dZX(ω)

  代入功率谱和谱表示之间的关系可以得到

E ( Z X ( ω ) ) 2 = ∣ H ( ω ) ∣ 2 1 2 π d ω E ( Z X ( ω ) ) 2 = 1 2 π S X ( ω ) d ω E(Z_X(\omega))^2 = |H(\omega)|^2 \frac{1}{2 \pi} d\omega \\ E(Z_X(\omega))^2 = \frac{1}{2\pi}S_X(\omega) d\omega \\ E(ZX(ω))2=H(ω)22π1dωE(ZX(ω))2=2π1SX(ω)dω

S X ( ω ) = ∣ H ( ω ) ∣ 2 S_X(\omega)=|H(\omega)|^2 SX(ω)=H(ω)2

  我们可以发现,这个传递函数就是对功率谱进行分解得到的

∣ H ( ω ) ∣ 2 = S X ( ω ) |H(\omega)|^2 = S_X(\omega) H(ω)2=SX(ω)

4.2 因果谱因式分解的多样性

  这里我们说明一点,因为传递函数的模的平方等于功率谱,相当于传递函数是损失了相位信息的。因此,对功率谱做分解得到传递函数必定不是唯一的,因此,假设滤波器hat H(ω)是因果的,且满足

∣ H ^ ( ω ) ∣ 2 = S X ( ω ) = ∣ H ( ω ) ∣ 2 |\hat H(\omega)|^2 = S_X(\omega) = |H(\omega)|^2 H^(ω)2=SX(ω)=H(ω)2

  则必定存在一个因果的传递函数θ,满足

θ ( ω ) ∣ θ ( ω ) ∣ = 1 ⇒ H ^ ( ω ) = θ ( ω ) H ( ω ) \theta(\omega) \quad |\theta(\omega)| = 1 \\ \Rightarrow \hat H(\omega) = \theta(\omega) H(\omega) θ(ω)θ(ω)=1H^(ω)=θ(ω)H(ω)

  如果要在这些幅度响应相同的因果滤波器中做出选择,还需要进一步附加条件。最小相位是最为常见的附加条件。如果H(ω)和H-1(ω)都是因果的,就称系统H(ω)是最小相位系统。满足最小相位条件的系统在所有给定幅度响应的系统中是唯一确定的。

  最小相位系统的传递函数,其极点和零点都在单位圆内

H ( Z ) = ∑ k = 0 ∞ h k Z − k H(Z) = \sum_{k=0}^{\infty} h_k Z^{-k} H(Z)=k=0hkZk

  对于最小相位系统,不但可以用Un对Xn进行滑动平均表示,也可以用Xn对Un进行滑动平均表示

4.3 谱因式分解过程

  如果随机过程的功率谱密度满足一定条件,通过谱因式分解是可以求得最小相位系统的

  如果logS(ω)可以做傅里叶级数展开

l o g S ( ω ) = ∑ k = − ∞ k = + ∞ b k e x p ( − j ω k ) log S(\omega) = \sum_{k=-\infty}^{k=+\infty} b_k exp(-j\omega k) logS(ω)=k=k=+bkexp(jωk)

  那么满足|H(ω)|2 = S(ω)的最小相位因果系统为

H ( ω ) = e x p ( b 0 2 + ∑ k = 1 ∞ b k e x p ( − j ω k ) ) H(\omega) = exp(\frac{b_0}{2}+\sum_{k=1}^{\infty}b_k exp(-j\omega k )) H(ω)=exp(2b0+k=1bkexp(jωk))

  我们可以证明一下这个事情,首先证明传递函数的平方等于功率谱

S ( ω ) = H ( ω ) H ( ω ) ‾ = e x p ( b 0 2 + ∑ k = 1 ∞ b k e x p ( − j ω k ) ) e x p ( b 0 2 + ∑ k = 1 ∞ b k e x p ( − j ω k ) ) ‾ = e x p ( b 0 2 + ∑ k = 1 ∞ b k e x p ( − j ω k ) ) e x p ( b 0 2 + ∑ k = 1 ∞ b − k e x p ( j ω k ) ) = e x p ( b 0 + ∑ k = 1 ∞ b k e x p ( − j ω k ) + ∑ − ∞ k = − 1 b k e x p ( − j ω k ) ) = e x p ( ∑ k = − ∞ + ∞ b k e x p ( − j ω k ) ) = e x p ( l o g ( S ( ω ) ) ) = S ( ω ) S(\omega) = H(\omega) \overline{H(\omega)} \\ = exp(\frac{b_0}{2}+\sum_{k=1}^{\infty}b_k exp(-j\omega k )) \overline{exp(\frac{b_0}{2}+\sum_{k=1}^{\infty}b_k exp(-j\omega k ))} \\ = exp(\frac{b_0}{2}+\sum_{k=1}^{\infty}b_k exp(-j\omega k ))exp(\frac{b_0}{2}+\sum_{k=1}^{\infty}b_{-k} exp(j\omega k )) \\ = exp(b_0 + \sum_{k=1}^{\infty}b_k exp(-j\omega k) +\sum^{k=-1}_{-\infty}b_k exp(-j\omega k )) \\ = exp(\sum_{k=-\infty}^{+\infty}b_k exp(-j\omega k )) = exp(log(S(\omega))) = S(ω) S(ω)=H(ω)H(ω)=exp(2b0+k=1bkexp(jωk))exp(2b0+k=1bkexp(jωk))=exp(2b0+k=1bkexp(jωk))exp(2b0+k=1bkexp(jωk))=exp(b0+k=1bkexp(jωk)+k=1bkexp(jωk))=exp(k=+bkexp(jωk))=exp(log(S(ω)))=S(ω)

  然后再证明因果性条件

  考虑到,指数的幂级数展开

e x p ( g ( z ) ) = 1 + g ( z ) + g 2 ( z ) 2 ! + . . . + g k ( z ) k ! + . . . exp(g(z)) = 1 + g(z)+\frac{g^2(z)}{2!} +...+\frac{g^k(z)}{k!}+... exp(g(z))=1+g(z)+2!g2(z)+...+k!gk(z)+...

  如果g(z)的幂级数展开只有非负幂次项,则exp(g(z))的展开中也只有非负幂次项。因此,如果g(z)是因果的,exp(g(z))也是因果的

  由于

g ( ω ) = b 0 2 + ∑ k = 1 ∞ b k e x p ( − j ω k ) g(\omega) = \frac{b_0}{2}+\sum_{k=1}^{\infty}b_k exp(-j\omega k ) g(ω)=2b0+k=1bkexp(jωk)

  因此该传递函数一定也是个因果的

  第三点,还需要验证其最小相位的问题。也就是验证其零点都在单位圆里面,根据指数函数不等于0的特点,这个传递函数没有零点,也满足最小相位的要求

H ( ω ) = e x p ( b 0 2 + ∑ k = 1 ∞ b k e x p ( − j ω k ) ) = 0 H(\omega) = exp(\frac{b_0}{2}+\sum_{k=1}^{\infty}b_k exp(-j\omega k )) \cancel = 0 H(ω)=exp(2b0+k=1bkexp(jωk))= 0

  由功率谱出发计算幅度已知的最小相因果滤波器的过程称为随机过程的谱因式分解。谱因式分解是一种正交化步骤,是一种白化。白噪声过程是随机过程X对应的信息过程

4.4 有理谱因式分解过程

  一般的谱因式分解需要用到傅里叶级数展开,非常的麻烦,如果功率谱有其他特点,谱因式分解会得到简化。

  这里介绍有理谱

S ( ω ) = S ( Z ) ∣ Z = e x p ( j ω ) = ∣ P ( Z ) Q ( Z ) ∣ Z = e x p ( j ω ) ∣ 2 S(\omega) = S(Z)|_{Z=exp(j\omega)} = |\frac{P(Z)}{Q(Z)}|_{Z=exp(j\omega)}|^2 S(ω)=S(Z)Z=exp(jω)=Q(Z)P(Z)Z=exp(jω)2

  假设差分方程定义的离散系统

∑ k = 0 N b k X n − k = ∑ l = 0 M a l U n − l H ( Z ) = a 0 + a 1 z − 1 + . . . + a M Z − M b 0 + b 1 z − 1 + . . . + b N Z − N S X ( ω ) = σ U 2 ∣ a 0 + a 1 z − 1 + . . . + a M Z − M b 0 + b 1 z − 1 + . . . + b N Z − N ∣ 2 = σ U 2 ∣ P ( e x p ( − j ω ) ) ∣ 2 ∣ Q ( e x p ( − j ω ) ) ∣ 2 = σ U 2 ∣ P ( Z ) ∣ 2 ∣ Q ( Z ) ∣ 2 ∣ Z = e x p ( j ω ) \sum_{k=0}^N b_k X_{n-k} = \sum_{l=0}^M a_l U_{n-l} \\ H(Z) = \frac{a_0+a_1z^{-1}+...+a_M Z^{-M}}{b_0+b_1z^{-1}+...+b_N Z^{-N}} \\ S_X(\omega) = \sigma_U^2|\frac{a_0+a_1z^{-1}+...+a_M Z^{-M}}{b_0+b_1z^{-1}+...+b_N Z^{-N}} |^2 \\ = \sigma_U^2 \frac{|P(exp(-j\omega))|^2}{|Q(exp(-j\omega))|^2} \\ = \sigma_U^2 \frac{|P(Z)|^2}{|Q(Z)|^2}|_{Z=exp(j\omega)} k=0NbkXnk=l=0MalUnlH(Z)=b0+b1z1+...+bNZNa0+a1z1+...+aMZMSX(ω)=σU2b0+b1z1+...+bNZNa0+a1z1+...+aMZM2=σU2Q(exp(jω))2P(exp(jω))2=σU2Q(Z)2P(Z)2Z=exp(jω)

  上式可以转换为

S ( Z ) = ( Z − α 1 ) . . . ( Z − α M ) . . . ( Z − 1 − α 1 ) ( Z − 1 − α M ) ( Z − β 1 ) . . . ( Z − β N ) . . . ( Z − 1 − β 1 ) ( Z − 1 − β N ) = N ( Z ) N ( Z − 1 ) D ( Z ) D ( Z − 1 ) S + ( Z ) = ( Z − α 1 ) . . . ( Z − α M ) ( Z − β 1 ) . . . ( Z − β N ) = N ( Z ) D ( Z ) S(Z)= \frac{(Z-\alpha_1)...(Z-\alpha_M)...(Z^{-1}-\alpha_1)(Z^{-1}-\alpha_M)}{(Z-\beta_1)...(Z-\beta_N)...(Z^{-1}-\beta_1)(Z^{-1}-\beta_N)} \\ = \frac{N(Z)N(Z^{-1})}{D(Z)D(Z^{-1})} \\ S^+ (Z) = \frac{(Z-\alpha_1)...(Z-\alpha_M)}{(Z-\beta_1)...(Z-\beta_N)} =\frac{N(Z)}{D(Z)} S(Z)=(Zβ1)...(ZβN)...(Z1β1)(Z1βN)(Zα1)...(ZαM)...(Z1α1)(Z1αM)=D(Z)D(Z1)N(Z)N(Z1)S+(Z)=(Zβ1)...(ZβN)(Zα1)...(ZαM)=D(Z)N(Z)

  差分方程也有类似的谱因式分解,不过对于离散的情况衡量的是零点和极点在单位圆内外,到了连续情况衡量的就是零点和极点在虚轴左边和右边。

5. 线性预测滤波器的具体形式

  前面介绍了很多随机过程是否能够被完美预测,以及预测的误差一些问题,下面介绍两种具体的线性预测方法,代表了线性预测的两种思路。

5.1 Wiener滤波器

∫ − ∞ + ∞ h ( t ′ − τ ′ ) R Z ( τ ′ ) d τ ′ = R Z Y ( t ′ ) \int_{-\infty}^{+\infty} h(t' - \tau') R_Z(\tau')d\tau' = R_{ZY}(t') +h(tτ)RZ(τ)dτ=RZY(t)

H ( ω ) ∗ S Z ( ω ) = S Z Y ( ω ) H(\omega)*S_Z(\omega) = S_{ZY}(\omega) H(ω)SZ(ω)=SZY(ω)

  wiener是满足wiener-hopf关系的一种线性预测滤波器。其思路是基于正交性原理进行预测,把要预测的变量投影到数据空间进行估计。

  不过这个滤波器是非因果的维纳滤波器,因为我们不可能得到时间为负的数据,我们又不能直接截断,因为直接阶段的话得到的并不是最优估计。

  但是可以把信号分解为白噪声,白噪声是彼此正交的,可以做截断,就可以得到因果的wiener滤波器

Z ( t ) → h 1 → U ( t ) → h 2 → U ^ ( t ) Z(t)\rightarrow \boxed{h_1} \rightarrow U(t) \rightarrow \boxed{h_2} \rightarrow \hat U(t) Z(t)h1U(t)h2U^(t)

Z ( t ) → Y ( t ) = > [ S Z Y ( ω ) S Z − ( ω ) ] + ∗ 1 S Z + ( ω ) Z(t) \rightarrow Y(t) =>[ \frac{S_{ZY}(\omega)}{ S^{-} _Z (\omega) }]_+* \frac{1}{S^{+} _Z (\omega)} Z(t)Y(t)=>[SZ(ω)SZY(ω)]+SZ+(ω)1

  具体过程见正交化与维纳滤波

5.2 Kalman 滤波器

  模型
{ Z n + 1 = F n Z n + V n Y n = H n Z n + W n \begin{cases} Z_{n+1} = F_n Z_n +V_n\\ Y_n = H_n Z_n + W_n \end{cases} {Zn+1=FnZn+VnYn=HnZn+Wn

E ( V n ) = E ( W n ) = 0 E ( V n ∗ W n T ) = 0 E ( V n V m T ) = Q n δ m n E ( W n W m T ) = R n δ m n E(V_n) = E(W_n) = 0 \quad E(V_n*W_n^T) = 0 \\ E(V_n V_m^T) = Q_n \delta_{mn} \quad E(W_n W_m^T) = R_n \delta_{mn} E(Vn)=E(Wn)=0E(VnWnT)=0E(VnVmT)=QnδmnE(WnWmT)=Rnδmn

  我们得到了卡尔曼滤波的完整方程

  • 预测:在原有估计的基础上,通过状态方程进行预测
  • 校正:在预测的基础上,通过卡尔曼增益,实现校正
  • 卡尔曼增益:也算出来结果了
  • 协方差矩阵的递推

Z ^ n + 1 ∣ n = F n ∗ Z ^ n ∣ n ( 1 ) Z ^ n + 1 ∣ n + 1 = Z ^ n + 1 ∣ n + K n + 1 ∗ ( Y n + 1 − H n + 1 Z ^ n + 1 ∣ n ) ( 2 ) K n + 1 = P n + 1 ∣ n ∗ H n + 1 T ∗ [ H n + 1 ∗ P n + 1 ∣ n ∗ H n + 1 T + R n + 1 ] − 1 ( 3 ) P n + 1 ∣ n = F n P n ∣ n F n T + Q n ( 4 ) P n + 1 ∣ n + 1 = ( I − K n + 1 H n + 1 ) ∗ P n + 1 ∣ n ( 5 ) \hat Z_{n+1|n} = F_n * \hat Z_{n|n} \quad\quad(1) \\ \hat Z_{n+1|n+1} = \hat Z_{n+1|n} + K_{n+1}*(Y_{n+1}-H_{n+1} \hat Z_{n+1|n}) \quad\quad(2) \\ K_{n+1} = P_{n+1|n} * H_{n+1}^T *[H_{n+1}*P_{n+1|n} * H_{n+1}^T +R_{n+1}]^{-1} \quad\quad(3) \\ P_{n+1|n} = F_n P_{n|n} F_n^T + Q_n \quad\quad(4) \\ P_{n+1|n+1} =(I -K_{n+1}H_{n+1} )*P_{n+1|n} \quad\quad(5) Z^n+1n=FnZ^nn(1)Z^n+1n+1=Z^n+1n+Kn+1(Yn+1Hn+1Z^n+1n)(2)Kn+1=Pn+1nHn+1T[Hn+1Pn+1nHn+1T+Rn+1]1(3)Pn+1n=FnPnnFnT+Qn(4)Pn+1n+1=(IKn+1Hn+1)Pn+1n(5)

  kalman filter通过建立状态空间模型进行预测,需要有先验知识,具体推导见卡尔曼滤波

6. 噪声中直流分量的估计

  下面我们举一个具体的例子,就是对噪声的直流分量进行估计,我们分别使用wiener滤波器和kalman滤波器进行

7.1 基于Wiener滤波器分析

  在这个问题中,我们有已知的采样数据Z

Z = ( Z 1 , . . . , Z n ) T Z = (Z_1,...,Z_n)^T Z=(Z1,...,Zn)T

   我们要用采样数据Z去估计A,我们令Y=A,这样看起来似乎奇怪,因为别的线性估计的两个参数Z和Y都是已知的,需要确定的是组合系数,但是在这个问题中A是未知的,组合系数也是未知的。但是我们尽可这样做,因为后面我们可以发现,我们求解的系数与A的关系并没有那么密切。

Y = A Y = A Y=A

  其实按照我们以往的经验,对采样数据求平均值,基本就是个很好的估计,也就是

A = 1 n ∑ k = 1 n Z k A = \frac{1}{n}\sum_{k=1}^{n}Z_k A=n1k=1nZk

  我们来计算看看是不是这个样子的。

  采样数据满足如下关系

Z 1 = A + N 1 Z 2 = A + N 2 . . . Z n = A + N n Z_1 = A + N_1 \\ Z_2 = A + N_2 \\ ...\\ Z_n = A + N_n Z1=A+N1Z2=A+N2...Zn=A+Nn

  Z的模型可以写做
Z = A ∗ e + N Z = A*e +N Z=Ae+N

  其中
Z = ( Z 1 , . . . , Z n ) T e = ( 1 , . . . , 1 ) T N = ( N 1 , . . . , N n ) T Z = (Z_1,...,Z_n)^T \\ e = (1,...,1)^T \\ N = (N_1,...,N_n)^T Z=(Z1,...,Zn)Te=(1,...,1)TN=(N1,...,Nn)T

  根据Wiener-Hopf方程,我们需要计算一下两个变量的自相关和互相关

R Z = E ( Z ∗ Z T ) r Z Y = E ( Z A ) R_Z = E(Z*Z^T) \\ r_{ZY} = E(ZA) RZ=E(ZZT)rZY=E(ZA)

  计算自相关
R Z = E ( ( A ∗ e + N ) ∗ ( A ∗ e + N ) T ) = A 2 ∗ E ( e ∗ e T ) + E ( N ∗ N T ) R_Z = E((A*e+N)*(A*e+N)^T) = A^2*E(e*e^T)+E(N*N^T) RZ=E((Ae+N)(Ae+N)T)=A2E(eeT)+E(NNT)

  我们发现式子中,只有N是随机变量,并且N的自相关是对角线为σ2的对角阵,因此上式可以化为

R Z = A 2 e e T + σ 2 I R_Z = A^2 e e^T + \sigma^2 I RZ=A2eeT+σ2I

  计算互相关

r Z Y = E ( Z Y ) = E ( ( A ∗ e + N ) ∗ A ) = A 2 ∗ e r_{ZY} = E(ZY) =E((A*e+N)*A) = A^2*e rZY=E(ZY)=E((Ae+N)A)=A2e

  根据Wiener-Hopf方程,我们可以得到系数向量

α = R z − 1 ∗ r Z Y = ( A 2 e e T + σ 2 I ) − 1 ∗ ( A 2 ∗ e ) \alpha = R_z^{-1}*r_{ZY} = (A^2 e e^T + \sigma^2 I)^{-1} * ( A^2*e ) α=Rz1rZY=(A2eeT+σ2I)1(A2e)


  我们需要求组合矩阵的逆矩阵,我们在这里补充一点知识

形如  ( I + B b T ) − 1  的逆矩阵一定可以写成下面的形式 ( I + b b T ) − 1 = I − β ∗ b ∗ b T \text{形如 }(I+ Bb^T)^{-1} \text{ 的逆矩阵一定可以写成下面的形式}\\ (I+ bb^T)^{-1} = I - \beta*b*b^T 形如 (I+BbT)1 的逆矩阵一定可以写成下面的形式(I+bbT)1=IβbbT

  我们求解一下这个β。

  我们记住一点,这里面的b是个列向量,所以bTb是个常数,并且是可以提到乘法外面的。

( I + b b T ) − 1 = I − β ∗ b ∗ b T = > ( I + b b T ) ∗ ( I − β ∗ b ∗ b T ) = I (I+ bb^T)^{-1} = I - \beta*b*b^T => (I+ bb^T)*( I - \beta*b*b^T) = I (I+bbT)1=IβbbT=>(I+bbT)(IβbbT)=I
  把式子展开

I − β ∗ b ∗ b T + b ∗ b T − β ∗ b ∗ ( b T ∗ b ) ∗ b T = I I + ( 1 − β − β ∗ b T ∗ b ) ∗ ( b ∗ b T ) = I I - \beta*b*b^T +b*b^T-\beta*b*(b^T*b)*b^T = I\\ I + (1-\beta -\beta*b^T*b)*(b*b^T) = I IβbbT+bbTβb(bTb)bT=II+(1ββbTb)(bbT)=I

  因此

1 − β − β ∗ b T ∗ b = 0 1-\beta -\beta*b^T*b = 0 1ββbTb=0

β = 1 1 + b T ∗ b \beta = \frac{1}{1+b^T*b} β=1+bTb1

I + b ∗ b T = I − 1 1 + b T ∗ b ∗ b ∗ b T I + b*b^T = I - \frac{1}{1+b^T*b}*b*b^T I+bbT=I1+bTb1bbT


  我们回到我们原来的式子中

α = R z − 1 ∗ r Z Y = ( A 2 e e T + σ 2 I ) − 1 ∗ ( A 2 ∗ e ) = 1 σ 2 ∗ ( A 2 σ 2 ∗ e ∗ e T + I ) − 1 ∗ ( A 2 ∗ e ) ( 1 ) \alpha = R_z^{-1}*r_{ZY} = (A^2 e e^T + \sigma^2 I)^{-1} * ( A^2*e ) \\ = \frac{1}{\sigma^2}*(\frac{A^2}{\sigma^2}*e*e^T+I)^{-1}*( A^2*e ) \quad\quad \quad\quad(1) α=Rz1rZY=(A2eeT+σ2I)1(A2e)=σ21(σ2A2eeT+I)1(A2e)(1)

  我们令

A σ ∗ e = b ( 2 ) \frac{A}{\sigma}*e = b \quad\quad \quad\quad(2) σAe=b(2)

  可得

( A 2 σ 2 ∗ e ∗ e T + I ) − 1 = ( I + b ∗ b T ) − 1 = I − 1 1 + b T ∗ b ∗ b ∗ b T ( 3 ) (\frac{A^2}{\sigma^2}*e*e^T+I)^{-1} = (I + b*b^T)^{-1} = I - \frac{1}{1+b^T*b}*b*b^T \quad\quad \quad\quad(3) (σ2A2eeT+I)1=(I+bbT)1=I1+bTb1bbT(3)

  (2)代入(3)可得

( A 2 σ 2 e e T + I ) − 1 = I − A 2 σ 2 1 + A 2 σ 2 ∗ e T e ∗ e e T (\frac{A^2}{\sigma^2} e e^T+I)^{-1} = I - \frac{\frac{A^2}{\sigma^2}}{1+\frac{A^2}{\sigma^2}*e^Te}*e e^T (σ2A2eeT+I)1=I1+σ2A2eTeσ2A2eeT
  由于

e T e = n e^Te = n eTe=n

  所以上式等于
( A 2 σ 2 e e T + I ) − 1 = I − A 2 σ 2 1 + A 2 σ 2 ∗ n ∗ e e T ( 4 ) (\frac{A^2}{\sigma^2} e e^T+I)^{-1} = I - \frac{\frac{A^2}{\sigma^2}}{1+\frac{A^2}{\sigma^2}*n}*e e^T \quad\quad \quad\quad(4) (σ2A2eeT+I)1=I1+σ2A2nσ2A2eeT(4)

  (4)代入(1)

α = 1 σ 2 ∗ ( I − A 2 σ 2 1 + A 2 σ 2 ∗ n ∗ e e T ) ∗ ( A 2 ∗ e ) = 1 σ 2 ∗ ( I ∗ e − A 2 σ 2 1 + A 2 σ 2 ∗ n e ( e T ∗ e ) ) ∗ ( A 2 ) = 1 σ 2 ∗ ( I ∗ e − A 2 σ 2 1 + A 2 σ 2 ∗ n e n ) ∗ ( A 2 ) = 1 σ 2 ∗ ( I − A 2 σ 2 1 + A 2 σ 2 ∗ n n ∗ I ) ∗ ( A 2 e ) = A 2 σ 2 ∗ ( 1 1 + A 2 σ 2 ∗ n ) ∗ e = ( A 2 σ 2 1 + A 2 σ 2 ∗ n ) ∗ e = 1 n ∗ ( A 2 σ 2 1 n + A 2 σ 2 ) ∗ e \alpha = \frac{1}{\sigma^2} * ( I - \frac{\frac{A^2}{\sigma^2}}{1+\frac{A^2}{\sigma^2}*n}*e e^T ) *( A^2*e ) \\ = \frac{1}{\sigma^2} * ( I*e - \frac{\frac{A^2}{\sigma^2}}{1+\frac{A^2}{\sigma^2}*n}e (e^T*e) ) *( A^2) \\ = \frac{1}{\sigma^2} * ( I*e - \frac{\frac{A^2}{\sigma^2}}{1+\frac{A^2}{\sigma^2}*n}e n ) *( A^2) \\ =\frac{1}{\sigma^2} * ( I - \frac{\frac{A^2}{\sigma^2}}{1+\frac{A^2}{\sigma^2}*n} n*I ) *( A^2 e) \\ = \frac{A^2}{\sigma^2}*(\frac{1}{1+\frac{A^2}{\sigma^2}*n})*e \\ = (\frac{\frac{A^2}{\sigma^2}}{1+\frac{A^2}{\sigma^2}*n})*e \\ = \frac{1}{n}* (\frac{\frac{A^2}{\sigma^2}}{\frac{1}{n}+\frac{A^2}{\sigma^2}})*e α=σ21(I1+σ2A2nσ2A2eeT)(A2e)=σ21(Ie1+σ2A2nσ2A2e(eTe))(A2)=σ21(Ie1+σ2A2nσ2A2en)(A2)=σ21(I1+σ2A2nσ2A2nI)(A2e)=σ2A2(1+σ2A2n1)e=(1+σ2A2nσ2A2)e=n1(n1+σ2A2σ2A2)e

  我们观测这个估计器的结构,发现e前面是个常数。因此,每个系数都是相等的,可以看出来,本身最优的估计就是对所有的Z等权重相加,然后取平均。

  但是我们希望后面的东西是1,因为我们不知道A,甚至也不知道σ。这样我们有两种途径

  • n充分大,后面的常数就可以忽略不计
  • 信噪比充分高(信噪比就是信号的平均功率与噪声的平均功率之比),后面的常数就约等于1

7.2 基于Kalman滤波器分析

  我们也是建立一样的模型,假设要观测的信号是

Y k = A + N k Y_k = A+N_k Yk=A+Nk

  首先要建立状态空间模型

X k + 1 = X k , X 0 = A Y k = X k + N k , Y 0 = 0 X_{k+1} = X_k, X_0 = A \\ Y_k = X_k + N_k,Y_0 = 0 Xk+1=Xk,X0=AYk=Xk+Nk,Y0=0

E ( V n ) = E ( W n ) = 0 E ( V n ∗ W n T ) = 0 E ( V n V m T ) = 0 E ( W n W m T ) = σ 2 δ m n E(V_n) = E(W_n) = 0 \quad E(V_n*W_n^T) = 0 \\ E(V_n V_m^T) = 0 \quad E(W_n W_m^T) = \sigma^2 \delta_{mn} E(Vn)=E(Wn)=0E(VnWnT)=0E(VnVmT)=0E(WnWmT)=σ2δmn

  递推关系

X ^ k + 1 ∣ k + 1 = X ^ k + 1 ∣ k + K n + 1 ( Y k + 1 − X ^ k + 1 ∣ k ) K k + 1 = P k + 1 ∣ k P k + 1 ∣ k + σ 2 P k + 1 ∣ k = P k ∣ k P k + 1 ∣ k + 1 = P k + 1 ∣ k − P k + 1 ∣ k 2 P k + 1 ∣ k + σ 2 \hat X_{k+1|k+1} = \hat X_{k+1|k} + K_{n+1}(Y_{k+1}- \hat X_{k+1|k}) K_{k+1} = \frac{P_{k+1|k}}{P_{k+1|k}+\sigma^2} \\ P_{k+1|k} = P_{k|k}\\ P_{k+1|k+1} = P_{k+1|k} - \frac{P_{k+1|k}^2}{P_{k+1|k}+\sigma^2} X^k+1k+1=X^k+1k+Kn+1(Yk+1X^k+1k)Kk+1=Pk+1k+σ2Pk+1kPk+1k=PkkPk+1k+1=Pk+1kPk+1k+σ2Pk+1k2

  设定初始值

X 1 ∣ 0 = 0 P 1 ∣ 0 = σ 2 Y 0 = 0 X_{1|0} = 0 \\ P_{1|0} = \sigma^2 \\ Y_0 = 0 X10=0P10=σ2Y0=0

  可得

X ^ 2 ∣ 1 = 1 2 Y 1 = 1 2 ( Y 1 + Y 0 ) X ^ 3 ∣ 2 = X ^ 2 ∣ 1 + 1 3 ( Y 2 − X ^ 2 ∣ 1 ) = 1 3 ( Y 2 + Y 1 + Y 0 ) . . . X ^ k ∣ k − 1 = 1 k ( Y k − 1 + . . . + Y 1 + Y 0 ) \hat X_{2|1} = \frac{1}{2} Y_1 = \frac{1}{2} (Y_1+ Y_0) \\ \hat X_{3|2} = \hat X_{2|1}+\frac{1}{3} (Y_2 - \hat X_{2|1}) = \frac{1}{3} (Y_2+Y_1+ Y_0) \\ ... \\ \hat X_{k|k-1} = \frac{1}{k} (Y_{k-1}+...+Y_1+ Y_0) X^21=21Y1=21(Y1+Y0)X^32=X^21+31(Y2X^21)=31(Y2+Y1+Y0)...X^kk1=k1(Yk1+...+Y1+Y0)

7. 匹配滤波器

  维纳滤波器和卡尔曼滤波器是均方意义下的线性估计,最优准则都是相同的。判断线性滤波器的最优准则可以不唯一。

  下面介绍一种以输出信噪比最大为最优准则的滤波器,叫做匹配滤波器。

  假设S(t)是确定性信号,频谱为Ss(ω),N(t)是背景噪声功率。功率谱密度SN(ω),实际上观测数据是

X ( t ) = s ( t ) + N ( t ) X(t) = s(t) +N(t) X(t)=s(t)+N(t)

  设计一个以X(t)为输入的滤波器,使得在t0时刻输出的信噪比最大。令滤波器的冲激响应为h(t),传递函数为H(ω),则输出信号和噪声分别是

Y s ( t 0 ) = ∫ − ∞ + ∞ s ( τ ) h ( t 0 − τ ) d τ Y N ( t 0 ) = ∫ − ∞ + ∞ N ( τ ) h ( t 0 − τ ) d τ Y_s(t_0) = \int_{-\infty}^{+\infty}s(\tau) h(t_0 - \tau) d\tau \\ Y_N(t_0) = \int_{-\infty}^{+\infty}N(\tau) h(t_0 - \tau) d\tau Ys(t0)=+s(τ)h(t0τ)dτYN(t0)=+N(τ)h(t0τ)dτ

  信噪比

S N R ( t 0 ) = ∣ Y s ( t 0 ) ∣ 2 E ( Y N 2 ( t 0 ) ) = ∣ 1 2 π ∫ − ∞ + ∞ s ( τ ) h ( t 0 − τ ) d τ ∣ 2 1 2 π R Y N ( 0 ) = ∣ 1 2 π ∫ − ∞ + ∞ s ( τ ) h ( t 0 − τ ) d τ ∣ 2 1 2 π ∫ − ∞ + ∞ S Y N ( ω ) d ω = ∣ 1 2 π ∫ − ∞ + ∞ s ( τ ) h ( t 0 − τ ) d τ ∣ 2 1 2 π ∫ − ∞ + ∞ ∣ H ( ω ) ∣ 2 S N ( ω ) d ω = ∣ ∫ − ∞ + ∞ H ( ω ) S s ( ω ) e x p ( j ω t 0 ) d ω ∣ 2 2 π ∫ − ∞ + ∞ ∣ H ( ω ) ∣ 2 S N ( ω ) d ω = ∣ ∫ − ∞ + ∞ H ( ω ) S N ( ω ) 1 2 S s ( ω ) S N ( ω ) 1 2 e x p ( j ω t 0 ) d ω ∣ 2 2 π ∫ − ∞ + ∞ ∣ H ( ω ) ∣ 2 S N ( ω ) d ω SNR(t_0) = \frac{|Y_s(t_0)|^2}{E(Y_N^2(t_0))} \\ = \frac{ |\frac{1}{2\pi}\int_{-\infty}^{+\infty}s(\tau) h(t_0 - \tau) d\tau |^2 }{\frac{1}{2\pi}R_{Y_N}(0)} \\ =\frac{ |\frac{1}{2\pi}\int_{-\infty}^{+\infty}s(\tau) h(t_0 - \tau) d\tau|^2 }{\frac{1}{2\pi}\int_{-\infty}^{+\infty}S_{Y_N}(\omega)d\omega} \\ = \frac{ |\frac{1}{2\pi}\int_{-\infty}^{+\infty}s(\tau) h(t_0 - \tau) d\tau|^2}{\frac{1}{2\pi}\int_{-\infty}^{+\infty}|H(\omega)|^2S_{N}(\omega)d\omega} \\ =\frac{ |\int_{-\infty}^{+\infty} H(\omega) S_s(\omega) exp(j\omega t_0) d\omega|^2}{2\pi\int_{-\infty}^{+\infty}|H(\omega)|^2S_{N}(\omega)d\omega} =\frac{ |\int_{-\infty}^{+\infty} H(\omega)S_N(\omega)^{\frac{1}{2}} \frac{S_s(\omega)}{S_N(\omega)^{\frac{1}{2}}} exp(j\omega t_0) d\omega|^2}{2\pi\int_{-\infty}^{+\infty}|H(\omega)|^2S_{N}(\omega)d\omega} SNR(t0)=E(YN2(t0))Ys(t0)2=2π1RYN(0)2π1+s(τ)h(t0τ)dτ2=2π1+SYN(ω)dω2π1+s(τ)h(t0τ)dτ2=2π1+H(ω)2SN(ω)dω2π1+s(τ)h(t0τ)dτ2=2π+H(ω)2SN(ω)dω+H(ω)Ss(ω)exp(jωt0)dω2=2π+H(ω)2SN(ω)dω+H(ω)SN(ω)21SN(ω)21Ss(ω)exp(jωt0)dω2

  由柯西不等式

( ∫ − ∞ + ∞ f ( x ) g ( x ) d x ) 2 ≤ ∫ − ∞ + ∞ f ( x ) 2 d x ∫ − ∞ + ∞ g ( x ) 2 d x ( ∫ − ∞ + ∞ f ( x ) g ( x ) d x ) 2 ∫ − ∞ + ∞ f ( x ) 2 d x ≤ ∫ − ∞ + ∞ g ( x ) 2 d x (\int_{-\infty}^{+\infty} f(x)g(x)dx)^2 \leq \int_{-\infty}^{+\infty} f(x)^2dx\int_{-\infty}^{+\infty} g(x)^2dx \\ \frac{(\int_{-\infty}^{+\infty} f(x)g(x)dx)^2}{\int_{-\infty}^{+\infty} f(x)^2dx} \leq \int_{-\infty}^{+\infty} g(x)^2dx (+f(x)g(x)dx)2+f(x)2dx+g(x)2dx+f(x)2dx(+f(x)g(x)dx)2+g(x)2dx

  因此,上式可以表示为

S N R ≤ 1 2 π ∫ − ∞ + ∞ S s ( ω ) 2 S N ( ω ) ∣ e x p ( j ω t 0 ) ∣ 2 d ω = 1 2 π ∫ − ∞ + ∞ S s ( ω ) 2 S N ( ω ) d ω SNR \leq \frac{1}{2\pi} \int_{-\infty}^{+\infty} \frac{S_s(\omega)^2}{S_N(\omega)} |exp(j\omega t_0)|^2 d \omega = \frac{1}{2\pi} \int_{-\infty}^{+\infty} \frac{S_s(\omega)^2}{S_N(\omega)} d \omega SNR2π1+SN(ω)Ss(ω)2exp(jωt0)2dω=2π1+SN(ω)Ss(ω)2dω

  等号成立条件为

H ( ω ) S N ( ω ) 1 2 = k S s ( ω ) S N ( ω ) 1 2 e x p ( j ω t 0 ) ‾ H ( ω ) = k S s ( ω ) ‾ ∣ S N ( ω ) ∣ 2 e x p ( − j ω t 0 ) S N R m a x = 1 2 π ∫ − ∞ + ∞ S s ( ω ) 2 S N ( ω ) d ω H(\omega)S_N(\omega)^{\frac{1}{2}} = k \overline{\frac{S_s(\omega)}{S_N(\omega)^{\frac{1}{2}}} exp(j\omega t_0)} \\ H(ω) = k \frac{\overline{S_s(\omega)}}{|S_N(\omega)|^2} exp(-j\omega t_0) \\ SNR_{max} = \frac{1}{2\pi} \int_{-\infty}^{+\infty} \frac{S_s(\omega)^2}{S_N(\omega)} d \omega H(ω)SN(ω)21=kSN(ω)21Ss(ω)exp(jωt0)H(ω)=kSN(ω)2Ss(ω)exp(jωt0)SNRmax=2π1+SN(ω)Ss(ω)2dω
  如果背景噪声是白噪声可以表示为

H ( ω ) = k S s ( ω ) ‾ N 0 e x p ( − j ω t 0 ) H(ω) = k \frac{\overline{S_s(\omega)}}{N_0} exp(-j\omega t_0) H(ω)=kN0Ss(ω)exp(jωt0)

  如果不计常数和相位,传递函数与确定信号的频谱是共轭关系。

h ( t ) = k N 0 ∫ − ∞ + ∞ S s ( ω ) ‾ e x p ( − j ω t 0 ) e x p ( j ω t ) d ω = k N 0 ∫ − ∞ + ∞ S s ( − ω ) e x p ( − j ω ( t 0 − t ) ) d ω = k N 0 ∫ − ∞ + ∞ S s ( ω ) e x p ( j ω ( t 0 − t ) ) d ω = k N 0 s ( t − t 0 ) h(t) = \frac{k}{N_0} \int_{-\infty}^{+\infty}\overline{S_s(\omega)} exp(-j\omega t_0) exp(j\omega t)d\omega \\ =\frac{k}{N_0} \int_{-\infty}^{+\infty}S_s(-\omega) exp(-j\omega (t_0-t))d\omega \\ =\frac{k}{N_0} \int_{-\infty}^{+\infty}S_s(\omega) exp(j\omega (t_0-t))d\omega \\ = \frac{k}{N_0} s(t-t_0) h(t)=N0k+Ss(ω)exp(jωt0)exp(jωt)dω=N0k+Ss(ω)exp(jω(t0t))dω=N0k+Ss(ω)exp(jω(t0t))dω=N0ks(tt0)

  滤波器的冲激响应就是信号的时域波形经过反向和延时得到的。这就是匹配滤波器名字的来源。匹配率滤波器可以使得信号频谱中有的频率全部过滤,而信号中没有的频率会全部被抑制掉。

8. 小结

  这一部分介绍了随机过程的线性预测问题。

  • 首先介绍了什么样的随机过程可以被完全预测,什么样的随机过程不能被完美预测,会有误差。能被完美预测的是奇异过程,完全不能被预测的是正则过程
  • 然后介绍了如何通过分解手段把不能完美预测的部分分离出来研究
  • 然后研究了正则过程是满足paley-wiener条件的,可以被白化为白噪声的冲激响应,并且基于该手段分析了正则过程的预测误差
  • 然后利用谱因式分解,了解了如果用白噪声的冲激响应来表示正则过程,这个冲激响应的系数是怎么计算出来的。
  • 继而介绍了维纳滤波和卡尔曼滤波两种线性预测的手段,并举例说明了怎么用来估计实际问题
  • 最后从信噪比最大的角度介绍了匹配滤波器这种方法。
  • 0
    点赞
  • 4
    收藏
    觉得还不错? 一键收藏
  • 3
    评论
波士顿房价数据集是一个经典的机器学习数据集,其中包含了波士顿市不同地区房屋的价格及其相关特征,如犯罪率、房屋所在位置的贫困指数、房屋的平均房间数等。 在本项目中,我们将使用随机森林和线性回归两种模型来预测波士顿房价。 ## 1. 数据探索 我们首先加载数据,查看数据的基本信息: ```python import pandas as pd import numpy as np from sklearn.datasets import load_boston boston = load_boston() df = pd.DataFrame(boston.data, columns=boston.feature_names) df['PRICE'] = boston.target print(df.head()) print(df.describe()) ``` 输出结果如下: ``` CRIM ZN INDUS CHAS NOX RM AGE DIS RAD TAX PTRATIO B LSTAT PRICE 0 0.00632 18.0 2.31 0.0 0.538 6.575 65.2 4.0900 1.0 296.0 15.3 396.90 4.98 24.0 1 0.02731 0.0 7.07 0.0 0.469 6.421 78.9 4.9671 2.0 242.0 17.8 396.90 9.14 21.6 2 0.02729 0.0 7.07 0.0 0.469 7.185 61.1 4.9671 2.0 242.0 17.8 392.83 4.03 34.7 3 0.03237 0.0 2.18 0.0 0.458 6.998 45.8 6.0622 3.0 222.0 18.7 394.63 2.94 33.4 4 0.06905 0.0 2.18 0.0 0.458 7.147 54.2 6.0622 3.0 222.0 18.7 396.90 5.33 36.2 CRIM ZN INDUS CHAS NOX RM AGE DIS RAD TAX PTRATIO B LSTAT PRICE count 506.000000 506.000000 506.000000 506.000000 506.000000 506.000000 506.000000 506.000000 506.000000 506.000000 506.000000 506.000000 506.000000 506.000000 mean 3.613524 11.363636 11.136779 0.069170 0.554695 6.284634 68.574901 3.795043 9.549407 408.237154 18.455534 356.674032 12.653063 22.532806 std 8.601545 23.322453 6.860353 0.253994 0.115878 0.702617 28.148861 2.105710 8.707259 168.537116 2.164946 91.294864 7.141062 9.197104 min 0.006320 0.000000 0.460000 0.000000 0.385000 3.561000 2.900000 1.129600 1.000000 187.000000 12.600000 0.320000 1.730000 5.000000 25% 0.082045 0.000000 5.190000 0.000000 0.449000 5.885500 45.025000 2.100175 4.000000 279.000000 17.400000 375.377500 6.950000 17.025000 50% 0.256510 0.000000 9.690000 0.000000 0.538000 6.208500 77.500000 3.207450 5.000000 330.000000 19.050000 391.440000 11.360000 21.200000 75% 3.677082 12.500000 18.100000 0.000000 0.624000 6.623500 94.075000 5.188425 24.000000 666.000000 20.200000 396.225000 16.955000 25.000000 max 88.976200 100.000000 27.740000 1.000000 0.871000 8.780000 100.000000 12.126500 24.000000 711.000000 22.000000 396.900000 37.970000 50.000000 ``` 从输出结果可以看出,数据集共有506个样本,每个样本包含13个特征和1个目标变量(即房价)。特征包括: - CRIM:城镇人均犯罪率 - ZN:住宅用地所占比例 - INDUS:非零售商业用地所占比例 - CHAS:是否靠近查尔斯河(1表示靠近,0表示不靠近) - NOX:一氧化氮浓度 - RM:每个住宅的平均房间数 - AGE:1940年以前建成的自住单位的比例 - DIS:距离五个波士顿就业中心的加权距离 - RAD:距离高速公路的便利指数 - TAX:每万美元的不动产税率 - PTRATIO:城镇师生比例 - B:黑人比例 - LSTAT:低收入人群比例 ## 2. 数据预处理 我们需要将数据集分成训练集和测试集,其中训练集用于训练模型,测试集用于评估模型性能。我们将训练集和测试集的比例设置为7:3。 ```python from sklearn.model_selection import train_test_split X = df.drop('PRICE', axis=1) y = df['PRICE'] X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.3, random_state=42) ``` 我们还需要对数据进行标准化处理,以便于模型更好地学习数据集。 ```python from sklearn.preprocessing import StandardScaler scaler = StandardScaler() X_train_scaled = scaler.fit_transform(X_train) X_test_scaled = scaler.transform(X_test) ``` ## 3. 线性回归模型 首先,我们使用线性回归模型来预测波士顿房价。 ```python from sklearn.linear_model import LinearRegression lr = LinearRegression() lr.fit(X_train_scaled, y_train) print(lr.coef_) print(lr.intercept_) ``` 输出结果如下: ``` [-0.96124802 1.04288045 0.13821887 0.65948615 -1.92937869 2.60118157 0.14167752 -3.08023279 2.54724553 -1.81543372 -2.33302112 0.85664044 -3.74868433] 22.74580335752233 ``` 接下来,我们使用测试集对模型进行评估。 ```python from sklearn.metrics import r2_score, mean_squared_error y_pred = lr.predict(X_test_scaled) r2 = r2_score(y_test, y_pred) mse = mean_squared_error(y_test, y_pred) rmse = np.sqrt(mse) print('R2 score:', r2) print('MSE:', mse) print('RMSE:', rmse) ``` 输出结果如下: ``` R2 score: 0.711226005748496 MSE: 21.517444231177425 RMSE: 4.637781108552027 ``` ## 4. 随机森林模型 接下来,我们使用随机森林模型来预测波士顿房价。 ```python from sklearn.ensemble import RandomForestRegressor rf = RandomForestRegressor(n_estimators=100, random_state=42) rf.fit(X_train_scaled, y_train) ``` 接下来,我们使用测试集对模型进行评估。 ```python y_pred = rf.predict(X_test_scaled) r2 = r2_score(y_test, y_pred) mse = mean_squared_error(y_test, y_pred) rmse = np.sqrt(mse) print('R2 score:', r2) print('MSE:', mse) print('RMSE:', rmse) ``` 输出结果如下: ``` R2 score: 0.8463672544193478 MSE: 11.153723157894741 RMSE: 3.339417217974791 ``` ## 5. 结论 从以上结果可以看出,使用随机森林模型预测波士顿房价的性能比线性回归模型更好,其R2得分为0.85,平均方差为11.15,平均根方差为3.34。因此,我们可以选择随机森林模型作为预测波士顿房价的模型。
评论 3
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值