谱表示与长球面波函数PSWF
文章目录
1. 问题引入
这一部分会引入PSWF (Prolate Spherical Wave Function)的方法进行谱估计。我们这里会使用到特征函数这个重要概念
Prolate Spherical Wave Function \text{Prolate Spherical Wave Function} Prolate Spherical Wave Function
首先,我们对随机信号的谱分析进行一下回顾。
与确定信号相比,由于随机信号的傅里叶变换存在不收敛的问题,因此,我们没有办法把确定性信号的傅里叶分析方法使用在随机信号上。
但是我们前面也介绍过,为了解决这个问题,我们可以使用两个思路
第一个思路是使用功率谱密度的方法,但是相对于频谱,使用功率谱密度损失也会非常大。首先,功率谱密度是个二阶量,这样就意味着二阶量没有线性形式。其次, 这个变换是不保真的,损失了相位信息,变换是不可逆的。不过功率谱密度最终可以得到一个良好的性质
第二个思路就是使用谱表示的方法。谱表示更加接近随机过程的表示方法,我们会在这里对谱表示进行一个适当的介绍,并且基于谱表示介绍基于PSWF的方法进行谱估计。
谱表示(Spectral Representation)一般认为是瑞典学者Cramer首先研究的。
Spectral Representation \text{Spectral Representation} Spectral Representation
2. 随机矢量去相关化与KL展开
我们首先从一个随机矢量的去相关化问题开始介绍。
假设我们有一个随机矢量Z,Z的各个维度都是随机变量
Z = ( z 1 , . . . , z n ) T ∈ R n z k r . v . E ( z i z j ) = r i j = 0 Z = (z_1,...,z_n)^T \in R^n \quad z_k \quad r.v. \quad \\ E(z_iz_j) = r_{ij} \cancel = 0 Z=(z1,...,zn)T∈Rnzkr.v.E(zizj)=rij= 0
一般来说随机矢量都是有相关性的,也就是zk不但是随机变量,并且彼此之间是相关的。
但是现在我们希望能够找到一个函数,一个从n维空间映射到n维空间的函数,使得随机矢量Z映射成为另外一个随机矢量,并且这个随机矢量的各个分量之间没有相关性
g : R n → R n Y = g ( Z ) → E ( Y i Y j ) = r i δ i j = { r i i = j 0 i = j g:R^n \rightarrow R^n \\ Y = g(Z) \rightarrow E(Y_iY_j) = r_i \delta_{ij} =\begin{cases} r_i & i =j \\ 0 & i \cancel = j \end{cases} g:Rn→RnY=g(Z)→E(YiYj)=riδij={ri0i=ji= j
我们做的这个映射叫做白化,也叫做正交化
Whitening Orthogonal \text{Whitening} \\ \text{Orthogonal} WhiteningOrthogonal
这里我们把映射g简化为线性变换
g ( Z ) = A Z A ∈ R n ∗ n Y = A Z E ( Y i Y j ) = r i δ i j g(Z) = AZ \quad A \in R^{n*n} \\ Y = AZ \\ E(Y_iY_j) = r_i \delta_{ij} g(Z)=AZA∈Rn∗nY=AZE(YiYj)=riδij
这看起来像是个欠定问题,因为方程A中有n*n个未知数,而方程只有n*(n-1)/2个
我们求一下Y的相关矩阵
Y = A Z ⇒ R Y = E ( Y Y T ) = E ( A Z Z T A ) = A E ( Z Z T ) A T = A R Z A T Y = AZ \Rightarrow R_Y = E(YY^T)=E(AZZ^TA) = AE(ZZ^T)A^T =A R_Z A^T Y=AZ⇒RY=E(YYT)=E(AZZTA)=AE(ZZT)AT=ARZAT
如果不相关,那么Y的相关矩阵应该是一个对角阵
R Y = A R Z A T = Λ Λ = d i a g ( i ) R_Y=A R_Z A^T = \Lambda \quad \Lambda = diag \quad\quad(i) RY=ARZAT=ΛΛ=diag(i)
由于RZ是一个可三角化的对称矩阵,可以做特征分解(谱分解),得到一个正交的特征矩阵
R Z = U Λ U T R_Z= U \Lambda U^T RZ=UΛUT
该式子可以化为
U T R Z U = Λ = d i a g ( λ 1 , . . . , λ n ) ( i i ) U^T R_Z U= \Lambda = diag( \lambda_1,...,\lambda_n) \quad\quad(ii) UTRZU=Λ=diag(λ1,...,λn)(ii)
对比(i)和(ii),我们可以求得这个线性变换A
A = U T U = ( u 1 , . . . , u n ) A = U^T \\ U = (u_1,...,u_n) A=UTU=(u1,...,un)
基于这种方法,我们就得到了随机矢量的去相关化
Y = A Z ⇒ Y = U T Z ⇒ Z = U Y ⇒ Z = ∑ k = 1 n y k u k Y = AZ \Rightarrow Y = U^T Z \Rightarrow Z = UY \Rightarrow Z = \sum_{k=1}^ny_ku_k Y=AZ⇒Y=UTZ⇒Z=UY⇒Z=k=1∑nykuk
我们发现,我们可以使用一组正交规范基去对Z进行展开。其中U是正交规范基
y是个标量,也是个随机变量,y在随机变量空间中,同样也是正交的
因此这里对Z的展开是一个正交展开,我们把这个展开叫做Karhumen-Loeve expansion
Karhumen-Loeve expansion \text{Karhumen-Loeve expansion } Karhumen-Loeve expansion
3. 随机过程的KL展开与谱表示
3.1 泛函分析概述
在了解了随机矢量的KL展开之后,我们现在希望把研究的问题从随机矢量扩展到随机过程上去
Z ( t ) = ∑ k = − ∞ + ∞ α k ϕ k ( t ) Z(t) = \sum_{k=-\infty}^{+\infty} \alpha_k \phi_k(t) Z(t)=k=−∞∑+∞αkϕk(t)
这个问题要在函数空间中进行,其实我们仍然可以把这个问题当做n维度的线性空间去看待。这算是个泛函分析问题,就是线性代数在函数空间中的延伸。
3.2 特征函数与积分方程
我们现在研究个这个随机过程的问题,更刚才研究的随机矢量是一回事。唯一去区别在于,矢量展开现在变成了函数展开问题。但是这个问题并不本质。因为函数也看做是无穷长的矢量。现在我们希望这个展开也具有双正交的性质
Z ( t ) = ∑ k = − ∞ + ∞ α k ϕ k ( t ) Z(t) = \sum_{k=-\infty}^{+\infty} \alpha_k \phi_k(t) Z(t)=k=−∞∑+∞αkϕk(t)
{ E ( α i α j ) = σ i δ i j ∫ I ϕ i ( t ) ϕ j ( t ) d t = β i δ i j \begin{cases} E(\alpha_i \alpha_j) = \sigma_i \delta_{ij} \\ \\ \int_I \phi_i(t) \phi_j(t) dt = \beta_i \delta_{ij} \end{cases} ⎩⎪⎨⎪⎧E(αiαj)=σiδij∫Iϕi(t)ϕj(t)dt=βiδij
这里其实有个问题需要留意,就是积分区间。我们把这个问题可以暂时一放。
如果φ(t)是个正交函数,那么系数将非常好求
α k = 1 β i ∫ I Z ( t ) ϕ i ( t ) d t \alpha_k =\frac{1}{\beta_i}\int_I Z(t) \phi_i(t) dt αk=βi1∫IZ(t)ϕi(t)dt
为了方便计算,我们假设β和σ全部是1
E ( α i α j ) = E ( 1 β i β j ∫ I Z ( t ) ϕ i ( t ) d t ∫ I Z ( s ) ϕ j ( s ) d s ) = E ( ∫ I Z ( t ) ϕ i ( t ) d t ∫ I Z ( s ) ϕ j ( s ) d s ) = ∫ I ∫ I E ( Z ( t ) Z ( s ) ) ϕ i ( t ) ϕ j ( t ) d t d s ( 1 ) E(\alpha_i \alpha_j) = E(\frac{1}{\beta_i \beta_j}\int_I Z(t) \phi_i(t) dt\int_I Z(s) \phi_j(s) ds) \\ =E(\int_I Z(t) \phi_i(t) dt\int_I Z(s) \phi_j(s) ds) \\ = \int_I \int_IE(Z(t)Z(s)) \phi_i(t) \phi_j(t)dtds \quad\quad (1) E(αiαj)=E(βiβj1∫IZ(t)ϕi(t)dt∫IZ(s)ϕj(s)ds)=E(∫IZ(t)ϕi(t)dt∫IZ(s)ϕj(s)ds)=∫I∫IE(Z(t)Z(s))ϕi(t)ϕj(t)dtds(1)
这里我们先不假设随机过程是个宽平稳的。但是随机过程的相关仍然是时间的函数
R ( t , s ) = E ( Z ( t ) Z ( s ) ) ( 2 ) R(t,s) = E(Z(t)Z(s)) \quad\quad (2) R(t,s)=E(Z(t)Z(s))(2)
把(2)代入(1)中去
E ( α i α j ) = ∫ I ∫ I E ( Z ( t ) Z ( s ) ) ϕ i ( t ) ϕ j ( t ) d t d s = ∫ I ( ∫ I R ( t , s ) ϕ i ( s ) d s ) ϕ j ( t ) d t = δ i j ( 3 ) E(\alpha_i \alpha_j) = \int_I \int_IE(Z(t)Z(s)) \phi_i(t) \phi_j(t)dtds \\ = \int_I (\int_I R(t,s) \phi_i(s)ds) \phi_j(t)dt = \delta_{ij} \quad\quad (3) E(αiαj)=∫I∫IE(Z(t)Z(s))ϕi(t)ϕj(t)dtds=∫I(∫IR(t,s)ϕi(s)ds)ϕj(t)dt=δij(3)
如果我们能够找到一个函数φ(t)满足这样的式子
∫ I R ( t , s ) ϕ i ( s ) d s = λ j ϕ j ( t ) ( 4 ) \int_I R(t,s) \phi_i(s)ds = \lambda_j \phi_j(t) \quad\quad (4) ∫IR(t,s)ϕi(s)ds=λjϕj(t)(4)
那样的话(3)式子就成立
在随机矢量的KL展开的问题中,我们是通过特征分解来解决问题的。我们发现(4)式子其实与特征向量非常像
R Z ϕ j = λ j ϕ j R_Z \phi_j = \lambda_j \phi_j RZϕj=λjϕj
其实在(4)式子中,φ(t)被称为特征函数,与线性代数的特征矢量相对。特征函数求解是个积分方程求解问题。
R(t,s)在积分方程中被称为kernal。
如果我们能够找到这样的函数φ,α就是正交的了。同时,我们还必须保证φ本身是个正交基函数。与特征分解相同,如果R(t,s)能够保证自己是个对称函数,那么得到的特征函数就是正交基函数
R Z ( t , s ) = R Z ( s , t ) R_Z(t,s) = R_Z(s,t) RZ(t,s)=RZ(s,t)
我们要把这个事情证明一下。首先我们从特征分解开始证明。如果矩阵R的对称的,R的特征向量之间都是正交的
{ R ϕ i = λ i ϕ i R ϕ j = λ j ϕ j λ i = λ j \begin{cases} R \phi_i = \lambda_i \phi_i \\ R \phi_j = \lambda_j \phi_j \\ \lambda_i \cancel = \lambda_j \end{cases} ⎩⎪⎨⎪⎧Rϕi=λiϕiRϕj=λjϕjλi= λj
⇒ ϕ i T ϕ j = ( ϕ i T R ϕ j ) λ j = ( ϕ i T R ϕ j ) T λ j = ( ϕ j T R ϕ i ) λ j = ( ϕ j T λ i ϕ i ) λ j = λ i λ j ( ϕ j T ϕ i ) \Rightarrow \phi_i^T \phi_j =\frac{ (\phi_i^T R \phi_j) }{\lambda_j} = \frac{ (\phi_i^T R \phi_j)^T }{\lambda_j} \\ =\frac{ (\phi_j^T R \phi_i)}{\lambda_j} = \frac{ (\phi_j^T \lambda_i \phi_i)}{\lambda_j} = \frac{\lambda_i}{\lambda_j}(\phi_j^T \phi_i) ⇒ϕiTϕj=λj(ϕiTRϕj)=λj(ϕiTRϕj)T=λj(ϕjTRϕi)=λj(ϕjTλiϕi)=λjλi(ϕjTϕi)
因为正交矩阵是个对称矩阵,左右应该相当,但是λi与λi不相等,只能说明左右都是0,也就是
ϕ i T ϕ j = 0 \phi_i^T \phi_j = 0 ϕiTϕj=0
证明对称矩阵的特征矢量是正交的。
然后我们再来证明一下,对称函数的特征函数也是正交的
∫ I ϕ i ( t ) ϕ j ( t ) d t = 1 λ j ∫ I ϕ i ( t ) ( ∫ I R ( t , s ) ϕ j ( s ) d s ) d t = 1 λ j ∫ I ∫ I R ( t , s ) ϕ i ( t ) ϕ j ( s ) d t d s = 1 λ j ∫ I ( ∫ I R Z ( s , t ) ϕ i ( t ) d t ) ϕ j ( s ) d s = λ i λ j ∫ I ϕ i ( s ) ϕ j ( s ) d s \int_I \phi_i(t) \phi_j(t) dt =\frac{1}{ \lambda_j} \int_I \phi_i(t)(\int_I R(t,s) \phi_j(s)ds )dt \\ = \frac{1}{ \lambda_j} \int_I \int_I R(t,s) \phi_i(t) \phi_j(s) dtds \\ = \frac{1}{ \lambda_j} \int_I (\int_I R_Z(s,t) \phi_i(t)dt)\phi_j(s) ds \\ = \frac{\lambda_i}{ \lambda_j} \int_I \phi_i(s)\phi_j(s) ds ∫Iϕi(t)ϕj(t)dt=λj1∫Iϕi(t)(∫IR(t,s)ϕj(s)ds)dt=λj1∫I∫IR(t,s)ϕi(t)ϕj(s)dtds=λj1∫I(∫IRZ(s,t)ϕi(t)dt)ϕj(s)ds=λjλi∫Iϕi(s)ϕj(s)ds
同理可得,左右都是0,即对称函数的特征函数是正交的。
∫ I ϕ i ( s ) ϕ j ( s ) d s = 0 \int_I \phi_i(s)\phi_j(s) ds = 0 ∫Iϕi(s)ϕj(s)ds=0
现在我们就有条件做KL展开了
Z ( t ) = ∑ k = − ∞ + ∞ α k ϕ k ( t ) Z(t) = \sum_{k=-\infty}^{+\infty} \alpha_k \phi_k(t) Z(t)=k=−∞∑+∞αkϕk(t)
3.3 宽平稳周期性随机过程的KL展开
我们现在假设这个随机过程是宽平稳的,我们就可以很好的求这个特征函数了。这个特征函数其实就是复指函数。我们可以验证一下。这样的话,KL展开就和傅里叶展开是一样的了
Z ( t ) = ∑ k = − ∞ + ∞ α k ϕ k ( t ) ⇒ ϕ k ( t ) = e x p ( j ω k t ) Z(t) = \sum_{k=-\infty}^{+\infty} \alpha_k \phi_k(t) \Rightarrow \phi_k(t) = exp(j \omega_k t) Z(t)=k=−∞∑+∞αkϕk(t)⇒ϕk(t)=exp(jωkt)
我们只需要验证一下,这个基函数是否是满足特征函数的即可
∫ I R ( t , s ) ϕ i ( s ) d s = λ j ϕ j ( t ) \int_I R(t,s) \phi_i(s)ds = \lambda_j \phi_j(t) \quad\quad ∫IR(t,s)ϕi(s)ds=λjϕj(t)
∫ I R Z ( t , s ) e x p ( j ω k s ) d s = ∫ I R Z ( t − s ) e x p ( j ω k s ) d s \int_I R_Z(t,s)exp(j\omega_ks)ds = \int_I R_Z(t-s)exp(j\omega_ks)ds ∫IRZ(t,s)exp(jωks)ds=∫IRZ(t−s)exp(jωks)ds
换元
Let t − s = s ′ Then s = t − s ′ ∫ I R Z ( t − s ) e x p ( j ω k s ) d s = ∫ I ′ ( t ) R Z ( s ′ ) e x p ( − j ω k s ′ ) d s ′ e x p ( j ω k t ) \text{Let } t-s = s' \\ \text{Then } s = t-s' \\ \int_I R_Z(t-s)exp(j\omega_ks)ds \\ = \int_{I'(t)}R_Z(s')exp(-j \omega_ks')ds' exp(j\omega_kt) Let t−s=s′Then s=t−s′∫IRZ(t−s)exp(jωks)ds=∫I′(t)RZ(s′)exp(−jωks′)ds′exp(jωkt)
验证的过程中,我们发现了一个问题,积分区间在换元之后是与时间有关的,这样积分一般不是个常数。
因此,我们在假设宽平稳随机过程的同时,引入周期性的条件,假设Z的周期是T
E ∣ Z ( t ) − Z ( t + T ) ∣ 2 = 0 ∀ t E|Z(t) - Z(t+T)|^2 = 0 \quad \forall t E∣Z(t)−Z(t+T)∣2=0∀t
接下来我们证明一下随机过程是周期的与相关函数是周期的是充要条件
E ∣ Z ( t ) − Z ( t + T ) ∣ 2 = 0 ↔ R Z ( τ ) = R Z ( τ + T ) E|Z(t) - Z(t+T)|^2 = 0 \leftrightarrow R_Z(\tau) = R_Z(\tau+T) E∣Z(t)−Z(t+T)∣2=0↔RZ(τ)=RZ(τ+T)
从左边往右边证明
E ∣ Z ( t ) − Z ( t + T ) ∣ 2 = E ( Z ( t ) 2 + Z ( t + T ) 2 − 2 Z ( t ) Z ( t + T ) ) = E ( Z ( t ) 2 ) + E ( Z ( t + T ) 2 ) − 2 E ( Z ( t ) Z ( t + T ) ) = 2 R Z ( 0 ) − 2 R Z ( T ) = 0 E|Z(t) - Z(t+T)|^2 = E(Z(t)^2 +Z(t+T)^2 - 2Z(t)Z(t+T)) \\ = E(Z(t)^2) +E(Z(t+T)^2) - 2E(Z(t)Z(t+T)) \\ = 2R_Z(0) - 2R_Z(T) = 0 E∣Z(t)−Z(t+T)∣2=E(Z(t)2+Z(t+T)2−2Z(t)Z(t+T))=E(Z(t)2)+E(Z(t+T)2)−2E(Z(t)Z(t+T))=2RZ(0)−2RZ(T)=0
从右边往左边证明
R Z ( τ ) = R Z ( τ ) ∣ R Z ( τ ) − R Z ( τ ) ∣ = ∣ E ( Z ( 0 ) ( Z ( τ + T ) − Z ( τ ) ) ∣ ≤ E ( ∣ Z ( 0 ) ∣ ∗ ∣ Z ( τ + T ) − Z ( τ ) ∣ ) = ( ∣ E ( Z ( 0 ) ) ∣ 2 ∗ ∣ E ( Z ( τ + T ) − Z ( τ ) ) ∣ 2 ) 1 2 = 0 R_Z(\tau) = R_Z(\tau) \\ |R_Z(\tau) - R_Z(\tau)| = |E(Z(0)(Z(\tau+T)-Z(\tau))| \\ \leq E(|Z(0)|*|Z(\tau+T)-Z(\tau)|) = (|E(Z(0))|^2*|E(Z(\tau+T)-Z(\tau))|^2)^{\frac{1}{2}} = 0 RZ(τ)=RZ(τ)∣RZ(τ)−RZ(τ)∣=∣E(Z(0)(Z(τ+T)−Z(τ))∣≤E(∣Z(0)∣∗∣Z(τ+T)−Z(τ)∣)=(∣E(Z(0))∣2∗∣E(Z(τ+T)−Z(τ))∣2)21=0
⇒ ( E ∣ Z ( τ + T ) − Z ( τ ) ∣ ) 2 = 0 \Rightarrow (E|Z(\tau+T)-Z(\tau)|)^2 = 0 ⇒(E∣Z(τ+T)−Z(τ)∣)2=0
我们这次在一个周期中进行展开,证明特征函数
∫
−
T
2
T
2
R
(
t
,
s
)
ϕ
i
(
s
)
d
s
=
λ
j
ϕ
j
(
t
)
\int_{-\frac{T}{2}}^{\frac{T}{2}} R(t,s) \phi_i(s)ds = \lambda_j \phi_j(t) \quad\quad
∫−2T2TR(t,s)ϕi(s)ds=λjϕj(t)
∫ − T 2 T 2 R Z ( t , s ) e x p ( j ω k s ) d s = ∫ − T 2 T 2 R Z ( t − s ) e x p ( j ω k s ) d s \int_{-\frac{T}{2}}^{\frac{T}{2}} R_Z(t,s)exp(j\omega_ks)ds = \int_{-\frac{T}{2}}^{\frac{T}{2}} R_Z(t-s)exp(j\omega_ks)ds ∫−2T2TRZ(t,s)exp(jωks)ds=∫−2T2TRZ(t−s)exp(jωks)ds
换元
KaTeX parse error: Undefined control sequence: \ at position 284: …xp(j\omega_kt) \̲ ̲= \lambda_k exp…
这里利用周期函数积分区域与起始位置无关的性质。同时我们就证明了复指函数确实是我们需要的特征函数。
因此我们就完成了宽平稳周期性随机过程的KL展开
Z ( t ) = ∑ k = − ∞ + ∞ α k e x p ( j 2 k π T t ) t ∈ [ − T 2 , T 2 ] E ( α i α j ) = 0 ( i = j ) Z(t) = \sum_{k=-\infty}^{+\infty} \alpha_k exp(j \frac{2k\pi}{T}t) \\ t \in [-\frac{T}{2},\frac{T}{2}] \\ E(\alpha_i \alpha_j) = 0(i \cancel = j) Z(t)=k=−∞∑+∞αkexp(jT2kπt)t∈[−2T,2T]E(αiαj)=0(i= j)
3.4 宽平稳非周期性随机过程的KL展开与谱表示
接下来,我们要继续把KL展开推广到宽平稳,非周期的情况,并且由此引出谱表示的概念
推广到非周期就意味着周期无穷大,求和就要变成求积分了
在求αk的时候,我们又会面临积分发散不收敛的问题了。
之前我们是通过定义二阶矩的办法解决的,现在我们要换一种思路
Stieltjes integral \text{Stieltjes integral} Stieltjes integral
Z ( t ) = ∫ − ∞ + ∞ e x p ( j ω t ) d F Z ( ω ) Z(t) = \int_{-\infty}^{+\infty} exp(j \omega t) dF_Z(\omega) Z(t)=∫−∞+∞exp(jωt)dFZ(ω)
我们就把这个系数表示为一个函数的积分。如果能求导我们就求,不能求导就放到这里。
我们对这个积分做两点注解
- 这是个什么
这是个增量形式变换。这里相当于把增量定义为了函数的形式
∑ k f ( k k ) ( x k − x k − 1 ) → ∫ f ( x ) d x ∑ k f ( k k ) ( g ( x k ) − g ( x k − 1 ) ) → ∫ f ( x ) d g ( x ) \sum_k f(k_k)(x_k-x_{k-1}) \rightarrow \int f(x)dx \\ \sum_k f(k_k)(g(x_k)-g(x_{k-1})) \rightarrow \int f(x)dg(x) k∑f(kk)(xk−xk−1)→∫f(x)dxk∑f(kk)(g(xk)−g(xk−1))→∫f(x)dg(x)
- 为什么这么写
我们相当于使用一个符号来回避了积分发散的困难,如果F可导就导出来,如果F有奇异点,就不求导,我们只保留增量的形式。
其中从随机矢量的KL展开,写到周期性随机过程的KL展开,再写到非周期性随机过程的KL展开,只是为了类比一件事情。这里面的F是个正交的。因为F的位置实际上对应的是系数α
E ( d F Z ( ω ) d F Z ( ω ′ ) ‾ ) = 0 ω = ω ′ E(dF_Z(\omega)\overline{dF_Z(\omega')}) = 0 \\ \omega \cancel = \omega' E(dFZ(ω)dFZ(ω′))=0ω= ω′
我们用增量的形式代替微元形式来看这个正交性
d F Z ( ω ) = F Z ( ω + Δ ) − F Z ( ω ) E [ ( F Z ( ω + Δ ) − F Z ( ω ) ) ( F Z ( ω ′ + Δ ) − F Z ( ω ′ ) ) ] = 0 ω = ω ′ dF_Z(\omega) = F_Z(\omega + \Delta) - F_Z(\omega) \\ E[(F_Z(\omega + \Delta) - F_Z(\omega))(F_Z(\omega' + \Delta) - F_Z(\omega'))] = 0 \\ \omega \cancel = \omega' dFZ(ω)=FZ(ω+Δ)−FZ(ω)E[(FZ(ω+Δ)−FZ(ω))(FZ(ω′+Δ)−FZ(ω′))]=0ω= ω′
这里面的FZ(ω)是随机过程Z的谱过程
F Z ( ω ) Spectrum Process F_Z(\omega) \text{ Spectrum Process} FZ(ω) Spectrum Process
随机过程就是一组随机变量,用某种标记,连在一起。这个标记可以是时间,但是也不一定非得是时间,如果是空间的话,就叫做随机场,如果是频率的话,就是谱过程。含有Z(t)的所有谱信息
这样我们就完成了非周期的宽平稳随机过程的谱表示
连续时间的谱表示为
Z ( t ) = ∫ − ∞ + ∞ e x p ( j ω t ) d F Z ( ω ) Z(t) = \int_{-\infty}^{+\infty} exp(j \omega t) dF_Z(\omega) Z(t)=∫−∞+∞exp(jωt)dFZ(ω)
离散时间的谱表示为
Z ( k ) = ∫ − π π e x p ( j ω k ) d F Z ( ω ) Z(k) = \int_{-\pi}^{\pi} exp(j \omega k) dF_Z(\omega) Z(k)=∫−ππexp(jωk)dFZ(ω)
因此谱表示,从形式上来看,就是个反傅里叶变换的形式。
4. 谱表示与功率谱密度的联系
4.1 从谱表示到功率谱密度
现在,我们想知道,谱表示和功率谱之间是否是兼容的
我们做如下规定
E ∣ d F Z ( ω ) ∣ 2 = 1 2 π S Z ( ω ) d ω E|dF_Z(\omega)|^2 = \frac{1}{2 \pi} S_Z(\omega) d \omega E∣dFZ(ω)∣2=2π1SZ(ω)dω
如果不做特殊说明,我们的信号一般是实信号,但是这里可能出现复数,因此,我们假定我们的信号是复信号。因此,我们计算宽平稳随机过程的相关为
R Z ( t − s ) = E ( Z ( t ) Z ( s ) ‾ ) = E ( ∫ − ∞ + ∞ e x p ( j ω t ) d F Z ( ω ) ∗ ∫ − ∞ + ∞ e x p ( j ω ′ s ) d F Z ( ω ′ ) ‾ ) = ∫ − ∞ + ∞ ∫ − ∞ + ∞ e x p ( j ω ( t − s ) ) E ( d F Z ( ω ) d F Z ( ω ′ ) ) R_Z(t-s) = E(Z(t) \overline{Z(s)}) \\ = E(\int_{-\infty}^{+\infty} exp(j \omega t) dF_Z(\omega)*\overline {\int_{-\infty}^{+\infty} exp(j \omega' s) dF_Z(\omega')}) \\ = \int_{-\infty}^{+\infty} \int_{-\infty}^{+\infty} exp(j \omega(t-s)) E(dF_Z(\omega) dF_Z(\omega')) RZ(t−s)=E(Z(t)Z(s))=E(∫−∞+∞exp(jωt)dFZ(ω)∗∫−∞+∞exp(jω′s)dFZ(ω′))=∫−∞+∞∫−∞+∞exp(jω(t−s))E(dFZ(ω)dFZ(ω′))
由于dF_Z(ω) 是正交的
E
(
d
F
Z
(
ω
)
d
F
Z
(
ω
′
)
‾
)
=
0
ω
=
ω
′
E(dF_Z(\omega)\overline{dF_Z(\omega')}) = 0 \\ \omega \cancel = \omega'
E(dFZ(ω)dFZ(ω′))=0ω=
ω′
因此上式实质是是一个一重积分。
R Z ( t − s ) = ∫ − ∞ + ∞ e x p ( j ω ( t − s ) ) E ∣ d F Z ( ω ) ∣ 2 = 1 2 π ∫ − ∞ + ∞ e x p ( j ω ( t − s ) ) S Z ( ω ) d ω R_Z(t-s)= \int_{-\infty}^{+\infty} exp(j \omega(t-s)) E|dF_Z(\omega)|^2 \\ = \frac{1}{2 \pi} \int_{-\infty}^{+\infty} exp(j \omega(t-s)) S_Z(\omega) d \omega RZ(t−s)=∫−∞+∞exp(jω(t−s))E∣dFZ(ω)∣2=2π1∫−∞+∞exp(jω(t−s))SZ(ω)dω
最终,我们可以看到,谱表示和功率谱密度能够得到相同的表达式,对谱表示的Z求相关,最后能够得到功率谱密度
4.2 从谱表示到功率谱密度的线性响应
我们也可以从另外一个角度看到谱表示与功率谱密度之间的关系
我们假设Y(t)是Z(t)通过线性系统变换得到的
Y ( t ) = ∫ − ∞ + ∞ h ( t − τ ) Z ( τ ) d τ = ∫ − ∞ + ∞ h ( t − τ ) ∫ − ∞ + ∞ e x p ( j ω τ ) d F Z ( ω ) d τ = ∫ − ∞ + ∞ ∫ − ∞ + ∞ h ( t − τ ) e x p ( j ω τ ) d τ d F Z ( ω ) Y(t) = \int_{-\infty}^{+\infty} h(t-\tau) Z(\tau) d\tau \\ = \int_{-\infty}^{+\infty} h(t-\tau) \int_{-\infty}^{+\infty} exp(j \omega \tau) dF_Z(\omega) d\tau \\ = \int_{-\infty}^{+\infty} \int_{-\infty}^{+\infty} h(t-\tau) exp(j \omega \tau) d\tau dF_Z(\omega) Y(t)=∫−∞+∞h(t−τ)Z(τ)dτ=∫−∞+∞h(t−τ)∫−∞+∞exp(jωτ)dFZ(ω)dτ=∫−∞+∞∫−∞+∞h(t−τ)exp(jωτ)dτdFZ(ω)
交换积分顺序后换元
Let τ ′ = t − τ \text{Let } \tau' = t - \tau Let τ′=t−τ
Y ( t ) = ∫ − ∞ + ∞ ∫ − ∞ + ∞ h ( τ ′ ) e x p ( − j ω τ ′ ) d τ ′ e x p ( j ω t ) d F Z ( ω ) = ∫ − ∞ + ∞ H ( j ω ) e x p ( j ω t ) d F Z ( ω ) [ 1 ] Y(t) = \int_{-\infty}^{+\infty} \int_{-\infty}^{+\infty} h(\tau') exp(-j \omega \tau') d\tau' exp(j \omega t) dF_Z(\omega) \\ = \int_{-\infty}^{+\infty} H(j \omega) exp(j \omega t) dF_Z(\omega) \quad\quad [1] Y(t)=∫−∞+∞∫−∞+∞h(τ′)exp(−jωτ′)dτ′exp(jωt)dFZ(ω)=∫−∞+∞H(jω)exp(jωt)dFZ(ω)[1]
对Y(t)做谱表示
Y ( t ) = ∫ − ∞ + ∞ e x p ( j ω t ) d F Y ( ω ) [ 2 ] Y(t) = \int_{-\infty}^{+\infty} exp(j \omega t) dF_Y(\omega) \quad\quad [2] Y(t)=∫−∞+∞exp(jωt)dFY(ω)[2]
[1]和[2]式子相等,因此我们可以得到
∫ − ∞ + ∞ H ( j ω ) e x p ( j ω t ) d F Z ( ω ) = ∫ − ∞ + ∞ e x p ( j ω t ) d F Y ( ω ) \int_{-\infty}^{+\infty} H(j \omega) exp(j \omega t) dF_Z(\omega) = \int_{-\infty}^{+\infty} exp(j \omega t) dF_Y(\omega) ∫−∞+∞H(jω)exp(jωt)dFZ(ω)=∫−∞+∞exp(jωt)dFY(ω)
即
H ( j ω ) d F Z ( ω ) = d F Y ( ω ) H(j \omega) dF_Z(\omega) = dF_Y(\omega) H(jω)dFZ(ω)=dFY(ω)
已知
E ∣ d F Y ( ω ) ∣ 2 = 1 2 π S Y ( ω ) d ω E|dF_Y(\omega)|^2 = \frac{1}{2 \pi} S_Y(\omega) d \omega E∣dFY(ω)∣2=2π1SY(ω)dω
可得
1 2 π S Y ( ω ) d ω = E ∣ d F Y ( ω ) ∣ 2 = E ∣ H ( j ω ) d F Z ( ω ) ∣ 2 = ∣ H ( j ω ) ∣ 2 E ∣ d F Z ( ω ) ∣ 2 = ∣ H ( j ω ) ∣ 2 1 2 π S Z ( ω ) d ω \frac{1}{2 \pi} S_Y(\omega) d \omega = E|dF_Y(\omega)|^2 \\ = E|H(j \omega) dF_Z(\omega)|^2 = |H(j \omega)|^2 E |dF_Z(\omega)|^2 \\ = |H(j \omega)|^2 \frac{1}{2 \pi} S_Z(\omega) d \omega 2π1SY(ω)dω=E∣dFY(ω)∣2=E∣H(jω)dFZ(ω)∣2=∣H(jω)∣2E∣dFZ(ω)∣2=∣H(jω)∣22π1SZ(ω)dω
可得
S Y ( ω ) = ∣ H ( j ω ) ∣ 2 S Z ( ω ) S_Y(\omega) = |H(j \omega)|^2 S_Z(\omega) SY(ω)=∣H(jω)∣2SZ(ω)
也可以得到与功率谱密度相同的表达
5. 基于PSWF与谱表示进行谱估计
5.1 基于谱表示进行谱估计
我们之前学习了基于周期图的谱估计,现在我们想利用谱表示技术实现谱估计。同样的,对于实际的信号,是离散的,长度有限的
Z ( k ) = { Z ( 1 ) , . . . , Z ( N ) } Z(k) = \{Z(1),...,Z(N) \} Z(k)={Z(1),...,Z(N)}
我们最Z(k)做离散傅里叶变换,得到频谱
Z ^ ( ω ) = ∑ k = 1 N Z ( k ) e x p ( − j ω k ) DFT \hat Z(\omega) = \sum_{k=1}^N Z(k) exp(-j \omega k) \\ \text{DFT} Z^(ω)=k=1∑NZ(k)exp(−jωk)DFT
- 基于功率谱密度和周期图的谱估计
如果我们用周期图进行谱估计,我们做的是这样的事情
S ^ Z ( ω ) = ∣ Z ^ ( ω ) ∣ 2 \hat S_Z(\omega) = | \hat Z(\omega) |^2 S^Z(ω)=∣Z^(ω)∣2
但是对数据进行截断之后,周期图估计误差比较大,功率谱的估计并不是非常好,我们希望能够借助谱表示获取更多的谱信息
- 基于谱表示和PSWF的谱估计
Z ^ ( ω ) = ∑ k = 1 N Z ( k ) e x p ( − j ω k ) = ∑ k = 1 N ∫ − π π e x p ( j ω ′ k ) d F Z ( ω ′ ) e x p ( − j ω k ) = ∫ − π π ( ∑ k = 1 N e x p ( − j k ( ω − ω ′ ) ) ) d F Z ( ω ′ ) \hat Z(\omega) = \sum_{k=1}^N Z(k) exp(-j \omega k) \\ = \sum_{k=1}^N \int_{-\pi}^{\pi} exp(j \omega' k) dF_Z(\omega') exp(-j \omega k) \\ = \int_{-\pi}^{\pi} (\sum_{k=1}^N exp(-jk(\omega - \omega'))) dF_Z(\omega') Z^(ω)=k=1∑NZ(k)exp(−jωk)=k=1∑N∫−ππexp(jω′k)dFZ(ω′)exp(−jωk)=∫−ππ(k=1∑Nexp(−jk(ω−ω′)))dFZ(ω′)
令
Let D ( ω ) = ∑ k = 1 N e x p ( − j k ( ω ) ) Then D ( ω − ω ′ ) = ∑ k = 1 N e x p ( − j k ( ω − ω ′ ) ) \text{Let } D(\omega) = \sum_{k=1}^N exp(-jk(\omega )) \\ \text{Then } D(\omega- \omega') = \sum_{k=1}^N exp(-jk(\omega - \omega')) Let D(ω)=k=1∑Nexp(−jk(ω))Then D(ω−ω′)=k=1∑Nexp(−jk(ω−ω′))
则,原式为
Z ^ ( ω ) = ∫ − π π D ( ω − ω ′ ) d F Z ( ω ′ ) Fundamental Equation \hat Z(\omega) = \int_{-\pi}^{\pi} D(\omega - \omega') dF_Z(\omega') \\ \text{Fundamental Equation} Z^(ω)=∫−ππD(ω−ω′)dFZ(ω′)Fundamental Equation
这个等式是个非常关键的等式,因为方程的左边是用截断数据做的傅里叶变换,方程右边的F是真实的谱信息。左边是能够得到的谱信息。因此,我们希望能够通过截断数据的谱,使他尽可能得到真实的谱信息。
5.2 基本方程的求解
如果我们希望能够得到真实的谱信息,我们就需要求解这个基本方程。这其实是个线性积分方程。如果我们想要求解这个方程,就需要泛函分析的手段。线性积分方程的求解其实是通过特征函数求解的。
我们与线性方程的求解做个类比
5.2.1 线性方程的求解
如果我们直接求解
Z = D ∗ F F = D − 1 Z Z = D*F \\ F = D^{-1}Z Z=D∗FF=D−1Z
但是积分方程不能求逆。我们得想其他方法。 我们可以使用特征向量的方法进行求解
- (1) 计算特征向量
D U k = λ k U k k = 1 , 2 , . . . D U_k = \lambda_k U_k \\ k = 1,2,... DUk=λkUkk=1,2,...
- (2)以特征向量为基表示F
F = ∑ i α i U i Z = D F = ∑ i α i D U i = ∑ i α i λ i U i F = \sum_i \alpha_i U_i \\ Z=DF = \sum_i \alpha_i D U_i = \sum_i \alpha_i \lambda_i U_i F=i∑αiUiZ=DF=i∑αiDUi=i∑αiλiUi
- (3) 求解系数
U k T Z = U k T ∗ ∑ i α i λ i U i = α k λ k ⇒ α k = U k T Z λ k U_k ^TZ = U_k^T *\sum_i \alpha_i \lambda_i U_i = \alpha_k \lambda_k \\ \Rightarrow \alpha_k = \frac{U_k ^TZ}{\lambda_k} UkTZ=UkT∗i∑αiλiUi=αkλk⇒αk=λkUkTZ
- (4) 求解F
F = ∑ i U i T Z λ i U i F = \sum_i \frac{U_i ^TZ}{\lambda_i} U_i F=i∑λiUiTZUi
5.2.2 线性积分方程的求解–PSWF
我们用同样的方法求解线性积分方程
Z ^ ( ω ) = ∫ − π π D ( ω − ω ′ ) d F Z ( ω ′ ) \hat Z(\omega) = \int_{-\pi}^{\pi} D(\omega - \omega') dF_Z(\omega') Z^(ω)=∫−ππD(ω−ω′)dFZ(ω′)
- (1) 求解特征函数
∫ − π π D ( ω − ω ′ ) U k ( ω ′ ) d ω = λ k U k ( ω ) \int_{-\pi}^{\pi} D(\omega - \omega') U_k(\omega') d\omega = \lambda_k U_k(\omega) ∫−ππD(ω−ω′)Uk(ω′)dω=λkUk(ω)
这个Uk叫做长球面波函数,长球面波函数就是以复指函数为核的积分方程的特征函数,因此特别适合用于截断以后的傅里叶变换。
PSWF最成功的应用就是在信号处理中。但是用的很不普遍,因为这个方程数值解不好求,而且数值稳定性很差
U k ( ω ) ∼ P S W F U_k(\omega) \sim PSWF Uk(ω)∼PSWF
- (2) 用PSWF表示真实谱
d F Z ( ω ′ ) = ∑ i α i U i ( ω ) d ω dF_Z(\omega') = \sum_i \alpha_i U_i(\omega) d \omega dFZ(ω′)=i∑αiUi(ω)dω
Z ^ ( ω ) = ∫ − π π D ( ω − ω ′ ) d F Z ( ω ′ ) = ∑ i α i ∫ − π π D ( ω − ω ′ ) U i ( ω ) d ω = ∑ i α i λ i U i ( ω ) \hat Z(\omega) = \int_{-\pi}^{\pi} D(\omega - \omega') dF_Z(\omega') \\ = \sum_i \alpha_i \int_{-\pi}^{\pi} D(\omega - \omega')U_i(\omega) d \omega \\ = \sum_i \alpha_i \lambda_i U_i(\omega) Z^(ω)=∫−ππD(ω−ω′)dFZ(ω′)=i∑αi∫−ππD(ω−ω′)Ui(ω)dω=i∑αiλiUi(ω)
- (3) 求解系数
∫ − π π Z ^ ( ω ) U k ( ω ) d ω = < Z ^ ( ω ) , U k ( ω ) > = ∑ i α i λ i < U i ( ω ) , U k ( ω ) > = α k λ k ⇒ α k = < Z ^ ( ω ) , U k ( ω ) > λ k \int_{-\pi}^{\pi} \hat Z(\omega) U_k(\omega) d \omega = <\hat Z(\omega),U_k(\omega)> = \sum_i \alpha_i \lambda_i <U_i(\omega),U_k(\omega)> = \alpha_k \lambda_k \\ \Rightarrow \alpha_k = \frac{<\hat Z(\omega),U_k(\omega)> }{\lambda_k} ∫−ππZ^(ω)Uk(ω)dω=<Z^(ω),Uk(ω)>=i∑αiλi<Ui(ω),Uk(ω)>=αkλk⇒αk=λk<Z^(ω),Uk(ω)>
- (4) 求解真实谱
Z ^ n e w ( ω ) = ∑ i α i U i ( ω ) = ∑ i < Z ^ ( ω ) , U i ( ω ) > λ i U i ( ω ) \hat Z_{new}(\omega) = \sum_i \alpha_i U_i(\omega) \\ = \sum_i \frac{<\hat Z(\omega),U_i(\omega)> }{\lambda_i} U_i(\omega) Z^new(ω)=i∑αiUi(ω)=i∑λi<Z^(ω),Ui(ω)>Ui(ω)
估计的谱和精准的谱通过基本方程的方法联系在了一起,并且通过PSWF,可以让估计的谱更加精确的去表示精准谱
积分方程逆不了,特征函数是求解积分方程的重要工具,特征函数是可以offline计算的,不影响谱估计的运算效果。计算中只有第三项才是真正的计算复杂度。
由于特征函数的求解很不容易的,所以这种方法使用的不是非常广泛
6. 小结
这一节我们从正交化开始介绍,先以矢量的KL展开为例进行引入
Z ∈ R n Z = U Y = ∑ k Y k U k Z \in R^n \\ Z = UY = \sum_k Y_k U_k Z∈RnZ=UY=k∑YkUk
然后我们推广到随机矢量的去相关化,介绍了对随机矢量进行双正交展开(KL展开)
Z ( t ) = ∑ k α k ϕ ( t ) t ∈ I Z(t) = \sum_k \alpha_k \phi(t) \\ t \in I Z(t)=k∑αkϕ(t)t∈I
然后又从随机矢量推广到了随机过程,先从周期性的宽平稳随机过程开始
Z ( t ) = ∑ k α k e x p ( j 2 k π T t ) t ∈ [ − T 2 , T 2 ] Z(t) = \sum_k \alpha_k exp(j \frac{2k\pi}{T}t) \\ t\in[-\frac{T}{2},\frac{T}{2}] Z(t)=k∑αkexp(jT2kπt)t∈[−2T,2T]
最后引入谱表示,推广到了更加普通的非周期性的宽平稳随机过程
Z ( t ) = ∫ − ∞ + ∞ e x p ( j ω ) d F Z ( ω ) Z(t) = \int_{-\infty}^{+\infty} exp(j \omega) dF_Z (\omega) Z(t)=∫−∞+∞exp(jω)dFZ(ω)
其中F是个正交分量,是随机过程的谱过程