随机过程的线性预测问题
文章目录
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),t∈R}Xs={X(u),u≤t}Y=(t+τ)Xs⇒Y
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),t∈R}Xs={X(u),u=A∪B,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),t∈R}{W(t),t∈R}X⊃{X(t),t∈R}∪{W(t),t∈R}Z(t)=X(t)+W(t)subspaceXs={Z(u),u≤t}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}⇒Ln−1X=span{Xk,k≤k−1}
事实上,随机过程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=(Xn∣Ln−1X)+(Xn−Xn∣Ln−1X)
其中前一部分数据表示新数据在旧空间上的投影,而后一部分数据在旧空间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=Xn−Xn∣Ln−1X
因此,新的采样数据可以表示为两个空间投影的和
X n = X n ∣ L n − 1 X + I n X_n = X_n|L_{n-1}^X +I_n Xn=Xn∣Ln−1X+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^=Y∣LnX
而我们知道投影具有这样的性质,如果子空间能够分解为两个正交的子空间,那么在原来子空间上的投影,就等于在两个新正交子空间上的投影的和
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+VU⊥V⇒Y∣W=Y∣U+Y∣V
因为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^=Y∣LnX=Y∣Ln−1X+Y∣In
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,k≤n}=span{Xk,k≤n}n∈Z<Ik,Im>=E(IkIm)=0k,m∈Zk= 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=Ln−1X+Xn=Ln−1X+Xn∣Ln−1X+In=Ln−1X⨁In=Ln−2X⨁In−1⨁In=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 ... LnX⊃Ln−1X⊃Ln−2X⊃...
这个子空间是有下界和上界的
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} limn→−∞LnX=L−∞Xlimn→+∞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 \\ L−∞X=k=−∞⋂k=+∞LkXL+∞X=k=−∞⋃k=+∞LkX
当这个子空间是上界和下界相等的时候,说明随机性不随着时间改变,就是奇异的
Singular ⇒ L − ∞ X = L + ∞ X \text{Singular} \Rightarrow L^X_{ -\infty} = L^X_{ +\infty} Singular⇒L−∞X=L+∞X
而如果下界是0,说明随机性是完全随机时间改变的,就是正则的
Regular ⇒ L − ∞ X = { 0 } \text{Regular} \Rightarrow L^X_{ -\infty} = \{0 \} Regular⇒L−∞X={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+m−Xn+m∣LnX∣∣=(E(∣Xn+m−Xn+m∣LnX∣2))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,k∈Z}X∼W.S.S.X is Singular⇔∃m0ϵ(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 Singular⇒Xn+m=Xn+m∣LnX⇒ϵ(m)=0∀m∈Z
然后证明充分性
∃ 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)=0m≤m0Xn+m=Xn+m∣LnX∈LnX∀n∈Z,m≤m0⇒∀m≤nXn∈LmX⇒LnX∈LmX⇒L−∞X=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,k∈Z}X∼W.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(Xn∣LmX+Xn−Xn∣LmX)2=E(Xn∣LmX)2+E(Xn−Xn∣LmX)2+0=E(Xn∣LmX)2+ϵ(n−m)m→−∞⇒RX(0)=E(Xn∣L−∞X)2+ϵ(∞)ϵ(∞)=RX(0)⇔Xn∣L−∞={0},∀n∈Z⇔L−∞={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+ω2∣logS(ω)∣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,n∈Z}Xn=Yn+ZnYn∼Regular ProcessZn∼Singular 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=0∑∞hkUn−k
- 连续形式
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,n∈Z+}E∣Xn∣<∞⇒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(ΔYn∣Xn,Xn−1,...,X0)=ΔYnE(ΔZn∣Xn,Xn−1,...,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,n∈Z}Xn=∫−π+πexp(jωn)dZX(ω)E∣ZX∣2=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(ω)=E∣ZX∣2⇒dFX(ω)=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=−∞∑+∞hkUn−k
其中
U k ∼ White Noise ∑ k ∣ h k ∣ 2 < ∞ U_k \sim \text{White Noise} \\ \sum_k |h_k|^2 < \infty Uk∼White Noisek∑∣hk∣2<∞
可以证明一下这个事情
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=−∞∑+∞hkUn−k
首先写一个白噪声的谱表示
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=−∞∑+∞hkUn−k=k=−∞∑+∞hk∫−π+πexp(jω(n−k))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 E∣ZX(ω)∣2=E∣H(ω)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ω(n−k))dZU(ω)=k=−∞∑+∞ckUn−k
3.2 具有因果性的滑动平均表示
不过刚才的那个滑动平均表示不具备因果性,我们无法获取未来的时刻的数据。因此,如果宽平稳随机过程满足一定条件的话,可以表示为具有因果性的滑动平均表示
X n = ∑ k = 0 ∞ h k U n − k X_n = \sum_{k=0}^{\infty}h_k U_{n-k} Xn=k=0∑∞hkUn−k
其中
U k ∼ White Noise ∑ k ∣ h k ∣ 2 < ∞ U_k \sim \text{White Noise} \\ \sum_k |h_k|^2 < \infty Uk∼White Noisek∑∣hk∣2<∞
这个分解要求宽平稳随机过程是一个正则过程。也就是要求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+ω2∣logS(ω)∣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=0∑∞dkXn−k
满足
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 E∣Xn+m−X^n+m∣2=minY∈span{Xk,k≤n}E∣Xn+m−Y∣2
由于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,k≤n}=span{Uk,k≤n}
最优误差为
σ 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=E∣Xm+n−Xm+n∣span(Xk,k≤n)∣2=E∣Xm+n−Xm+n∣span(Uk,k≤n)∣2=E∣h0Un+m+h1Un+m−1+...+hm−1Un+1∣2=∣h0∣2+...+∣hm−1∣2
其中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=0∑∞hkUn−k=k=0∑∞hk∫−π+πexp(jω(n−k))dZU(ω)=∫−π+πk=0∑∞hkexp(−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) θ(ω)∣θ(ω)∣=1⇒H^(ω)=θ(ω)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=0∑∞hkZ−k
对于最小相位系统,不但可以用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=1∑∞bkexp(−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=1∑∞bkexp(−jωk))exp(2b0+k=1∑∞bkexp(−jωk))=exp(2b0+k=1∑∞bkexp(−jωk))exp(2b0+k=1∑∞b−kexp(jωk))=exp(b0+k=1∑∞bkexp(−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=1∑∞bkexp(−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=1∑∞bkexp(−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=0∑NbkXn−k=l=0∑MalUn−lH(Z)=b0+b1z−1+...+bNZ−Na0+a1z−1+...+aMZ−MSX(ω)=σU2∣b0+b1z−1+...+bNZ−Na0+a1z−1+...+aMZ−M∣2=σU2∣Q(exp(−jω))∣2∣P(exp(−jω))∣2=σU2∣Q(Z)∣2∣P(Z)∣2∣Z=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)...(Z−1−β1)(Z−1−βN)(Z−α1)...(Z−αM)...(Z−1−α1)(Z−1−αM)=D(Z)D(Z−1)N(Z)N(Z−1)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)→h1→U(t)→h2→U^(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(Vn∗WnT)=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+1∣n=Fn∗Z^n∣n(1)Z^n+1∣n+1=Z^n+1∣n+Kn+1∗(Yn+1−Hn+1Z^n+1∣n)(2)Kn+1=Pn+1∣n∗Hn+1T∗[Hn+1∗Pn+1∣n∗Hn+1T+Rn+1]−1(3)Pn+1∣n=FnPn∣nFnT+Qn(4)Pn+1∣n+1=(I−Kn+1Hn+1)∗Pn+1∣n(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=1∑nZk
我们来计算看看是不是这个样子的。
采样数据满足如下关系
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=A∗e+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(Z∗ZT)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((A∗e+N)∗(A∗e+N)T)=A2∗E(e∗eT)+E(N∗NT)
我们发现式子中,只有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((A∗e+N)∗A)=A2∗e
根据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 ) α=Rz−1∗rZY=(A2eeT+σ2I)−1∗(A2∗e)
我们需要求组合矩阵的逆矩阵,我们在这里补充一点知识
形如 ( 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−β∗b∗bT
我们求解一下这个β。
我们记住一点,这里面的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−β∗b∗bT=>(I+bbT)∗(I−β∗b∗bT)=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−β∗b∗bT+b∗bT−β∗b∗(bT∗b)∗bT=II+(1−β−β∗bT∗b)∗(b∗bT)=I
因此
1 − β − β ∗ b T ∗ b = 0 1-\beta -\beta*b^T*b = 0 1−β−β∗bT∗b=0
β = 1 1 + b T ∗ b \beta = \frac{1}{1+b^T*b} β=1+bT∗b1
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+b∗bT=I−1+bT∗b1∗b∗bT
我们回到我们原来的式子中
α = 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) α=Rz−1∗rZY=(A2eeT+σ2I)−1∗(A2∗e)=σ21∗(σ2A2∗e∗eT+I)−1∗(A2∗e)(1)
我们令
A σ ∗ e = b ( 2 ) \frac{A}{\sigma}*e = b \quad\quad \quad\quad(2) σA∗e=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) (σ2A2∗e∗eT+I)−1=(I+b∗bT)−1=I−1+bT∗b1∗b∗bT(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=I−1+σ2A2∗eTeσ2A2∗eeT
由于
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=I−1+σ2A2∗nσ2A2∗eeT(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∗(I−1+σ2A2∗nσ2A2∗eeT)∗(A2∗e)=σ21∗(I∗e−1+σ2A2∗nσ2A2e(eT∗e))∗(A2)=σ21∗(I∗e−1+σ2A2∗nσ2A2en)∗(A2)=σ21∗(I−1+σ2A2∗nσ2A2n∗I)∗(A2e)=σ2A2∗(1+σ2A2∗n1)∗e=(1+σ2A2∗nσ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(Vn∗WnT)=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+1∣k+1=X^k+1∣k+Kn+1(Yk+1−X^k+1∣k)Kk+1=Pk+1∣k+σ2Pk+1∣kPk+1∣k=Pk∣kPk+1∣k+1=Pk+1∣k−Pk+1∣k+σ2Pk+1∣k2
设定初始值
X 1 ∣ 0 = 0 P 1 ∣ 0 = σ 2 Y 0 = 0 X_{1|0} = 0 \\ P_{1|0} = \sigma^2 \\ Y_0 = 0 X1∣0=0P1∣0=σ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^2∣1=21Y1=21(Y1+Y0)X^3∣2=X^2∣1+31(Y2−X^2∣1)=31(Y2+Y1+Y0)...X^k∣k−1=k1(Yk−1+...+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 SNR≤2π1∫−∞+∞SN(ω)Ss(ω)2∣exp(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(ω)=k∣SN(ω)∣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ω(t0−t))dω=N0k∫−∞+∞Ss(ω)exp(jω(t0−t))dω=N0ks(t−t0)
滤波器的冲激响应就是信号的时域波形经过反向和延时得到的。这就是匹配滤波器名字的来源。匹配率滤波器可以使得信号频谱中有的频率全部过滤,而信号中没有的频率会全部被抑制掉。
8. 小结
这一部分介绍了随机过程的线性预测问题。
- 首先介绍了什么样的随机过程可以被完全预测,什么样的随机过程不能被完美预测,会有误差。能被完美预测的是奇异过程,完全不能被预测的是正则过程
- 然后介绍了如何通过分解手段把不能完美预测的部分分离出来研究
- 然后研究了正则过程是满足paley-wiener条件的,可以被白化为白噪声的冲激响应,并且基于该手段分析了正则过程的预测误差
- 然后利用谱因式分解,了解了如果用白噪声的冲激响应来表示正则过程,这个冲激响应的系数是怎么计算出来的。
- 继而介绍了维纳滤波和卡尔曼滤波两种线性预测的手段,并举例说明了怎么用来估计实际问题
- 最后从信噪比最大的角度介绍了匹配滤波器这种方法。