【现代信号处理】 15 - 谱分析基础和周期图谱分析

谱分析基础和周期图谱估计

1. 什么是谱

  所谓的谱，就是把一个复杂的系统，按照某种指标进行分解。常见的谱就比如

• 频谱:分解的指标是频率，信号在每一个频率上强度的大小
• 光谱：分解的指标是光的波长，光在每一个波长上分量的大小
• 抗菌谱：广谱/窄谱，是指能够对抗的细菌种类

  我们这里的谱针对的是随机信号/随机过程

  接下来，我们要解决两个问题

• 确定信号中已经有了谱的概念，随机信号的概念是否可以直接搬过来呢?
• 如果不能直接搬过来，要做哪些工作呢?

2. 确定信号的频谱

  我们先对确定信号的谱进行一个回顾。

2.1 周期信号的频谱

  我们先假设确定信号具有周期性

$$\begin{gathered} \text{Period T} \\ X(t+T)=X(t)\forall t \end{gathered}$$

  周期信号可以做傅里叶级数展开。但是注意，一旦谈到傅里叶级数展开，就一定要规定区间。不谈区间的傅里叶级数展开都是耍流氓

$$\begin{gathered} X(t)=\sum _{k=-\infty }^{+\infty }{\alpha }_{k}exp(j\frac{2k\pi }{T}t)\text{ (Fourier Series)} \\ t\in [-\frac{T}{2},\frac{T}{2}] \end{gathered}$$
  我们对傅里叶级数展开具有两个认识

  第一，傅里叶级数展开是一种正交展开，因为傅里叶级数的基是有正交性的。这种正交性一定是体现在某一个区间上的。在这个区间上，基函数正交。基函数的内积如下

$$\frac{1}{T}{\int }_{-\frac{T}{2}}^{\frac{T}{2}}exp(j\frac{2k\pi }{T}t)exp(-j\frac{2m\pi }{T}t)dt={\delta }_{km}$$

  因为是正交的，这个展开的系数也非常好求，只要对傅里叶级数做内积，就能得到系数

$${\alpha }_k=\frac{1}{T}{\int }_{-\frac{T}{2}}^{\frac{T}{2}}X(t)exp(-j\frac{2k\pi }{T}t)dt$$

  第二个认识，这里1/T有专门的说法，叫做基频。周期函数并非在每一个位置都有能量，只有在基频整数倍的位置上，才有能量，因此，周期函数谱，被称为离散谱

$$\begin{gathered} \frac{2\pi }{T}or\frac{1}{T}\text{ BaseBand Frequency} \\ \text{Discrete Spectrum} \end{gathered}$$

2.2 非周期信号的频谱

  我们在周期信号的基础上，对我们的认识继续延伸

  也就是从周期信号延伸到非周期信号。就是令T趋近于无穷大。一旦周期函数的周期趋近于无穷大，也就变成了非周期函数

  我们来看一下，从周期函数变化到非周期函数，到底有哪些特性发生了变化

  我们把周期信号的系数用积分的形式进行代入

$$\begin{gathered} Periodic\to Non-Periodic \\ T\to \infty X(t)=\sum _{k=-\infty }^{+\infty }(\frac{1}{T}{\int }_{-\frac{T}{2}}^{+\frac{T}{2}}X(s)exp(-j\frac{2k\pi }{T}s)ds)exp(j\frac{2k\pi }{T}t)[-\frac{T}{2},-\frac{T}{2}] \end{gathered}$$

  我们可以对非周期函数的傅里叶级数展开进行变换

$$X(t)=\frac{1}{2\pi }\sum _{k=-\infty }^{+\infty }({\int }_{-\frac{T}{2}}^{+\frac{T}{2}}X(s)exp(-j\frac{2k\pi }{T}s)ds)exp(j\frac{2k\pi }{T}t)\frac{2\pi }{T}[-\frac{T}{2},-\frac{T}{2}]$$

  在这个式子上，我们让T趋近于无穷大，也就是T的倒数趋近于0

  这个式子，很容易让我们联想到积分和，事实上，此时，我们已经可以把式子写成积分的形式了,并且我们的积分函数和积分变量都有了一个新的表达

$$\begin{gathered} T\to +\infty \\ X(t)=\frac{1}{2\pi }{\int }_{-\infty }^{+\infty }Z(\omega )exp(j\omega t)d\omega \end{gathered}$$

  这就是我们大名鼎鼎的傅里叶变换

$$\text{Fourier Transform}\{\begin{array}{l}X(t)=\frac{1}{2\pi }{\int }_{-\infty }^{+\infty }Z(\omega )exp(j\omega t)d\omega \\ Z(\omega )={\int }_{-\infty }^{+\infty }X(t)exp(-j\omega t)dt\end{array}$$

  因此，对于非周期函数来说，仍然可以做某种展开，这种展开是连续的，也就是我们的积分。傅里叶变换可以看做是傅里叶级数展开的连续版本。因为我们仍然是以复指函数作为基函数进行展开的。差别在于，傅里叶变换中，已经不局限于某个区间进行展开了，因为区间就是整个实数域

  傅里叶变换一定是在实数轴上进行展开的。同时，傅里叶变换的系数也不局限于某个基频的整数倍了，而是充满了整个实数轴

  所以，这里形成的谱叫做连续谱，每一个频点上都有能量。这就是非周期函数，或者是一般函数的傅里叶分析

$$\text{Continuous Spectrum}$$

3. 随机信号的谱

3.1 随机信号做傅里叶分析的困难

  在了解了确定信号的频谱分析之后，我们就要转移到随机信号的分析上了。假设，我们想把确定信号的谱分析工具，平行的搬移到随机信号的谱分析中去。

  现在我们假设X(t)是随机信号。与确定信号相比，随机信号还是非常不同的。假如，我们想在随机信号上，重复确定性信号的傅里叶变换，这确实有一定的困难。主要困难在积分区间上，因为积分区间从有限区间到无限区间的变化，就需要考虑积分是否是收敛的。

  如果傅里叶积分是收敛的，我们就要求他绝对可积

$$\int |X(t)|dt<\infty$$

  绝对可积的这个条件，对于确定性信号很多时候是可以满足的，因为确定性信号在无穷远处会有衰减的趋势。

  但是随机信号往往不满足这个条件。比如有一大类随机信号，就是宽平稳随机过程。随机信号普遍具有这样的性质，就是不可能在无穷远处衰减，而是在不停的反复震荡。

  就比如，宽平稳意味着两个时刻的相关，只和两个时刻的距离有关系，而和时间轴上的位置没有关系。如果在无穷远处信号有衰减，衰减前和衰减后的相关，是不可能一样的

  我们看到的很多信号，比如电报信号，很多都是在时间轴上反复震荡的

  因此宽平稳随机过程很多时候都不是绝对可积的

$$\begin{gathered} Z(t)\text{ Stochastic} \\ \text{Wide-Sense Stationary}\cancel{\Rightarrow }{\int }_{-\infty }^{+\infty }|Z(t)|dt<\infty \end{gathered}$$

  这就使得我们在做傅里叶变换的时候会遇到很多的困难。人们证明，宽平稳随机过程，在很多信号的频点上都是发散的

  因此，对于随机信号，如果我们想做傅里叶变换，我们的理论就不够牢靠。我们无法阻止积分的发散。我们也就无法定义和使用傅里叶变换

  我们就一定要用新的方法来解决这个问题，通常来说解决的方法有两个

• 走向功率谱
• 走向谱表示

3.2 功率谱密度的物理和数学表示形式

3.2.1 功率谱密度的物理表示

  下面我们要用功率谱的手段去对随机信号的谱进行表示。在随机信号上，我们需要重新定义谱分析的工具

$${\int }_{-\frac{T}{2}}^{+\frac{T}{2}}Z(t)exp(-j\omega t)dt$$

  如果我们直接让T趋近于无穷大，会导致积分发散，我们需要做三步操作

  第一步，求模取平方
$$|{\int }_{-\frac{T}{2}}^{+\frac{T}{2}}Z(t)exp(-j\omega t)dt{|}^2$$

  这样做承受的损失是很大的，因为相位信息没有了

  这里使用的谱分析工具就已经和傅里叶变换有区别了。求模取平方就变不回来了,因此变换不具有可逆性

  第二步，我们求期望

  求期望是因为这个积分具有随机性，我们希望消除这种随机性。

$$E(|{\int }_{-\frac{T}{2}}^{+\frac{T}{2}}Z(t)exp(-j\omega t)dt{|}^2)$$

  第三步，我们还要除以T，实现在时间上的归一化

$$\frac{1}{T}E(|{\int }_{-\frac{T}{2}}^{+\frac{T}{2}}Z(t)exp(-j\omega t)dt{|}^2)$$

  然后我们再让T趋近于无穷，这样就定义出来了一个新的量，我们称之为功率谱密度

$$\text{Power Spectral Density}li{m}_{T\to \infty }\frac{1}{T}E(|{\int }_{-\frac{T}{2}}^{+\frac{T}{2}}Z(t)exp(-j\omega t)dt{|}^2)=S_Z(\omega )$$

  之所以叫谱，是因为这个公式与傅里叶级数和傅里叶变换具有一定的相似性。

  为什么叫做功率，是因为我们做的是求模取平方的操作。

  这样定义的功率谱只是一种定义方法，这是物理的定义方法

$$\text{Physical}{S}_Z(\omega )=li{m}_{T\to \infty }\frac{1}{T}E(|{\int }_{-\frac{T}{2}}^{+\frac{T}{2}}Z(t)exp(-j\omega t)dt{|}^2)$$

3.2.2 功率谱密度的数学表示

  下面，我们想对物理形式的功率谱密度进行计算和变形，最终得到数学形式的功率谱密度

  这个功率谱密度中有一个平方，对于复数来说，求平方就是复数与其共轭的积。

$$\begin{gathered} E|{\int }_{-\frac{T}{2}}^{+\frac{T}{2}}Z(t)exp(-j\omega t)dt{|}^2 \\ =E({\int }_{-\frac{T}{2}}^{+\frac{T}{2}}Z(t)exp(-j\omega t)dt)({\int }_{-\frac{T}{2}}^{+\frac{T}{2}}Z(s)exp(j\omega t)ds) \\ ={\int }_{-\frac{T}{2}}^{+\frac{T}{2}}{\int }_{-\frac{T}{2}}^{+\frac{T}{2}}E(Z(t)Z(s))exp(-j\omega (t-s))dtds \\ ={\int }_{-\frac{T}{2}}^{+\frac{T}{2}}{\int }_{-\frac{T}{2}}^{+\frac{T}{2}}{R}_Z(t-s)exp(-j\omega (t-s))dtds \end{gathered}$$
  因为我们假设他是宽平稳的，这个里面的期望是可以写成相关函数的。宽平稳的功率谱具有明确定义，非宽平稳的还不好说

  下面我们要做积分换元了。积分换元一般有三个步骤

• 旧元变新元
• 把雅克比行列式计算清楚
• 更换积分域

  第一步，换元

$$\begin{gathered} \text{Let }u=t-s \\ \text{Let }v=t+s \end{gathered}$$

  第二步，计算雅克比行列式

$$\begin{gathered} dtds\to dudv \\ dtds=|det(\frac{\partial (t,s)}{\partial (u,v)})|dudv \end{gathered}$$

  因为我们现在函数的表示形式是用t和s表示u和v，因此计算u和v对t和s的偏导数更加方便，因此我们要对雅克比行列式进行变形

$$(\frac{\partial (t,s)}{\partial (u,v)})=(\frac{\partial (u,v)}{\partial (t,s))}{)}^{-1}$$

  在计算行列式的时候，结果就是倒数关系

$$(\frac{\partial (u,v)}{\partial (t,s))})=\left(\begin{array}{cc}\frac{\partial u}{\partial t}& \frac{\partial u}{\partial s}\\ & \\ \frac{\partial v}{\partial t}& \frac{\partial v}{\partial s}\end{array}\right)=\left(\begin{array}{cc}1& -1\\ 1& 1\end{array}\right)det|(\frac{\partial (u,v)}{\partial (t,s))})|=2$$

  因此

$$|det(\frac{\partial (t,s)}{\partial (u,v)})|=\frac{1}{2}$$

$$dtds=\frac{1}{2}dudv$$

$$R_Z(t-s)exp(-j\omega (t-s))dtds\Rightarrow R_Z(u)exp(-j\omega u)\frac{1}{2}dudv$$

  第三步，积分域变换

  因为我们换元做的是一个线性变换，因此积分域大致形状还是四边形，我们只要找到四个顶点就可以

$$(0,T),(-T,0),(0,-T),(T,0)$$

$$\begin{gathered} E|{\int }_{-\frac{T}{2}}^{+\frac{T}{2}}Z(t)exp(-j\omega t)dt{|}^2 \\ ={\int }_{-\frac{T}{2}}^{+\frac{T}{2}}{\int }_{-\frac{T}{2}}^{+\frac{T}{2}}{R}_Z(t-s)exp(-j\omega (t-s))dtds \\ =({\int }_{-T}^0{\int }_{-T-u}^{T+u}+{\int }_0^T{\int }_{u-T}^{u+T})R_Z(u)exp(-ju\omega )\frac{1}{2}dudv \\ ={\int }_{-T}^T{\int }_{-T+|u|}^{-|u|+T}{R}_Z(u)exp(-ju\omega )\frac{1}{2}dudv \\ ={\int }_{-T}^T(T-|u|)R_Z(u)exp(-j\omega u)du \end{gathered}$$

  因此，我们可以得到一个相关函数傅里叶变换的形式

$$\begin{gathered} S_Z(\omega )=li{m}_{T\to \infty }\frac{1}{T}E(|{\int }_{-\frac{T}{2}}^{+\frac{T}{2}}Z(t)exp(-j\omega t)dt{|}^2) \\ =li{m}_{T\to \infty }\frac{1}{T}{\int }_{-T}^T(T-|u|)R_Z(u)exp(-j\omega u)du \\ =li{m}_{T\to \infty }{\int }_{-T}^T(1-\frac{|u|}{T})R_Z(u)exp(-j\omega u)du \\ ={\int }_{-\infty }^{+\infty }{R}_Z(u)exp(-j\omega u)du \end{gathered}$$

  然后我们就可以得到功率谱密度和相关函数是一个傅里叶变换对。这个关系叫做

$$wiener-Khinchin$$
  功率谱密度的这种描述方式，相比于物理的描述，更加数学化。因此这个称为功率谱密度的数学描述

$$\text{Mathematical}$$

  至此，我们得到了随机信号谱分析的有效工具。就是功率谱密度与相关函数，他们是傅里叶变换对

$$\begin{gathered} S_Z(\omega )={\int }_{-\infty }^{+\infty }{R}_Z(\tau )exp(-j\omega \tau )d\tau \\ {R}_Z(\tau )=\frac{1}{2\pi }{\int }_{-\infty }^{+\infty }{S}_Z(\omega )exp(j\omega \tau )d\omega \end{gathered}$$

  这里要解释一下它为什么要叫做功率谱密度。功率和谱的说明前面已经提到了，这里说一下密度的概念

  对于反傅里叶变换的形式，如果我们取相关函数在0点的值，可以得到

$$\begin{gathered} 2\pi R_Z(0)={\int }_{-\infty }^{+\infty }{S}_Z(\omega )d\omega \\ 1=\frac{{\int }_{-\infty }^{+\infty }{S}_Z(\omega )d\omega }{2\pi R_Z(0)} \end{gathered}$$

  我们可以看到功率谱密度的积分其实是一个常数，并且得到归一化以后，积分值就是1。同时功率谱密度是一个大于0的数字，这个性质非常像概率密度，因此把这个就叫做功率谱密度

3.3 功率谱密度的量纲

  我们接着讨论一下功率谱密度的量纲，功率谱密度的量纲的焦耳。可以有两种解读的方法

  一方面，我们从物理描述出发

$$\text{Physical}{S}_Z(\omega )=li{m}_{T\to \infty }\frac{1}{T}E(|{\int }_{-\frac{T}{2}}^{+\frac{T}{2}}Z(t)exp(-j\omega t)dt{|}^2)$$

  因为信号的量纲不是电流就是电压，我们假设是电流。复指函数没有量纲，积分的量纲是时间，积分就是(I*t)²，然后1/T去掉一个t，就得到了焦耳的单位I2T

  另外一方面，我们从数学描述出发

$$S_Z(\omega )={\int }_{-\infty }^{+\infty }{R}_Z(\tau )exp(-j\omega \tau )d\tau$$
$$R_Z(t)=E(|Z(t){|}^2)$$

  相关R是随机过程的模的平方取期望，这是随机过程功率的定义。因此积分完了以后，得到的是总的功率，因此功率谱密度积分完了以后得到的是总的功率，换句话说，就是功率在频率轴上进行分解的结果，因此量纲应该是功率除以频率，也就是功率乘以时间，也就是能量

3.4 功率谱密度与线性变换

  假设我们有一个随机信号，我们现在想知道，这个随机信号通过线性系统之后，得到的信号的功率谱密度是怎么样的。

  如果Z是确定信号，Z的频谱和Y的频谱之间的关系就是相差一个传递函数

$$\begin{gathered} Z(t)deterministic\underrightarrow{Z(t)}\boxed{H}\underrightarrow{Y(t)} \\ Z(\omega )H(\omega )=Y(\omega ) \end{gathered}$$

  我们知道冲激响应和传递函数之间的傅里叶变换对的关系

$$h(t)\underleftrightarrow{F}H(\omega )$$

  如果Z(t)是随机信号，要计算功率谱的话，首先要求相关

$$\begin{gathered} R_Y(t,s)=E(Y(t)Y(s))=E({\int }_{-\infty }^{+\infty }h(t-\tau )Z(\tau )d\tau {\int }_{-\infty }^{+\infty }h(s-r)Z(r)dr) \\ ={\int }_{-\infty }^{+\infty }{\int }_{-\infty }^{+\infty }{R}_Z(\tau -r)h(t-\tau )h(s-r)d\tau dr \end{gathered}$$

  对待卷积有两个技巧

• 被积函数的变量加在一起，应该刚好能够把积分变量消掉
• 卷积出来的结果一定是个函数，被积函数的变量加在一起是什么，就在那个地方进行取值

  我们发现现在这个式子中被积函数中的变量相加并不能消掉积分变量，但是变个符号可以

$$Let\tilde{h}(t)=h(-t)$$

$$\begin{gathered} R_Y(t,s)={\int }_{-\infty }^{+\infty }{\int }_{-\infty }^{+\infty }{R}_Z(\tau -r)h(t-\tau )h(s-r)d\tau dr \\ ={\int }_{-\infty }^{+\infty }{\int }_{-\infty }^{+\infty }{R}_Z(\tau -r)h(t-\tau )\tilde{h}(r-s)d\tau dr \end{gathered}$$

  我们发现

$$\tau -r+t-\tau +r-s=t-s$$

  因此得到的是一个卷积的关系

$$\begin{gathered} R_Y(t,s)={\int }_{-\infty }^{+\infty }{\int }_{-\infty }^{+\infty }{R}_Z(\tau -r)h(t-\tau )\tilde{h}(r-s)d\tau dr \\ =(R_Z⊛h⊛\tilde{h})(t-s) \end{gathered}$$

  如果我们对相关函数做傅里叶变换，得到的就是功率谱密度

$$\begin{gathered} S_Y(\omega )=\mathcal{F}(R_Y(\tau )) \\ =S_Z(\omega )H(\omega )\tilde{H}(\omega ) \end{gathered}$$

  两个传递函数之间是共轭关系

$$\begin{gathered} \tilde{H}(\omega )={\int }_{-\infty }^{+\infty }\tilde{h}(t)exp(-j\omega t)dt \\ ={\int }_{-\infty }^{+\infty }h(-t)exp(-j\omega t)dt \\ ={\int }_{-\infty }^{+\infty }h(t)exp(j\omega t)dt \\ =\overline{{\int }_{-\infty }^{+\infty }h(t)exp(-j\omega t)dt}=\overline{H(\omega )} \end{gathered}$$

  因此，我们可以得到

$$\begin{gathered} S_Y(\omega )=\mathcal{F}(R_Y(\tau )) \\ =S_Z(\omega )H(\omega )\tilde{H}(\omega ) \\ =S_Z(\omega )|H(\omega ){|}^2 \end{gathered}$$

  我们就得到了输入和输出函数的功率谱密度在频域上的关系

3.5 小结

  以上就是随机过程谱分析的基本理论，我们要理解两件事情

• 功率谱密度与频谱不一样。一个是一阶量，一个是二阶量；一个是可逆变换，一个是不可逆变换；一个是线性变换，一个是非线性变换。两个信号和的功率谱密度，不等于两个信号功率谱密度的和。
• 功率谱密度和相关函数有关系，是相关函数的傅里叶变换，而不是频谱所规定的频谱本身的傅里叶变换

4. 随机信号的谱估计

  现在，我们要把我们的知识从天上落实到地面上

  这里我们有两个问题。首先，一般随机信号经过采样之后，得到的数据是离散的，而不是连续的。因此我们现在要考察的是离散时间信号/离散时间样本。在离散时间信号上，我们仍然可以使用相关进行研究

{ Z ( k ) } k = − ∞ + ∞  Discrete-Time Signal(Samples) Stochastic \{Z(k)\}_{k=-\infty} ^{+\infty} \text{ Discrete-Time Signal(Samples) Stochastic} \\ {Z(k)}k=−∞+∞​ Discrete-Time Signal(Samples) Stochastic

  我们假设它宽平稳。这与连续的时间信号的相关没有什么区别

  我们可以在离散时间上定义功率谱了。与连续时间形成了直接对应。

• 第一，用傅里叶级数求和取代了积分
• 第二，用离散时间的相关，取代了连续时间的相关

  我们定义一下离散随机信号的相关以及功率谱

$$r_Z(l)=E(Z(k+l)Z(k))\Rightarrow S_Z(\omega )=\sum _{l=-\infty }^{+\infty }{r}_Z(l)exp(-j\omega l)$$

  第二个问题，我们拿到的数据，只有有限长度。

{ Z ( 1 ) , . . . , Z ( N ) }  Finite Samples \{Z(1),...,Z(N) \} \text{ Finite Samples} {Z(1),...,Z(N)} Finite Samples

  把谱分析问题转换为实际问题以后就是，在离散的，有限长的样本的条件下，如何对随机信号做谱分析?

$$\text{Spectral Analysis}$$

  有限样本下，我们能做到的仅仅是功率谱密度的一个估计。我们要让这个估计尽可能的接近实际的功率谱密度。这就是我们的谱估计的问题。在上个世纪，这是信号处理的核心问题，人们进行了大量的深入研究。我们这里只介绍谱估计问题中最经典的事情。

5 周期图估计

5.1 周期图的定义

  现在我们来研究谱估计问题

  我们从功率谱的物理定义出发。首先是连续时间的功率谱物理定义为

$$\text{Physical}{S}_Z(\omega )=li{m}_{T\to \infty }\frac{1}{T}E(|{\int }_{-\frac{T}{2}}^{+\frac{T}{2}}Z(t)exp(-j\omega t)dt{|}^2)$$

  如果现在变成离散时间，并且有限长度，物理定义就变成了

$${\hat{S}}_Z(\omega )=\frac{1}{N}|\sum _{k=1}^N(Z(k)exp(-j\omega k)){|}^2$$

  我们记得对于连续时间的随机变量来说，与傅里叶变换相比，功率谱估计的物理定义做了三件事情

• 求模取平方
• 取期望
• 做时间的归一化

  而这里，我们没有求期望。因为对于实际情况来说，期望没法求，所以干脆就不求了，但是时间的归一化还是可以做的。

  因此，我们就得到了物理形式的谱估计，这个叫做周期图

$$Periodogram$$

5.2 周期图的求解

  就像我们把功率谱从物理形式变成数学形式一样，这里我们也要做同样的事情,就是把物理形式的周期图估计进行求解

$$\begin{gathered} {\hat{S}}_Z(\omega )=\frac{1}{N}|\sum _{k=1}^N(Z(k)exp(-j\omega k)){|}^2 \\ =\frac{1}{N}[\sum _{k=1}^N\sum _{i=1}^{N}Z(k)Z(i)exp(-j\omega (k-i))] \end{gathered}$$

  累加和其实也是可以换元的，并且相比于积分来说，累加和不需要计算雅克比行列式

• 旧元换新元
• 更换求和区间

$$\begin{gathered} Letl=k-i \\ Letn=i \end{gathered}$$

  跟积分换元一样，求和区间变换也只要确定边界点就行

$$\begin{gathered} {\hat{S}}_Z(\omega )=\frac{1}{N}[\sum _{k=1}^N\sum _{i=1}^{N}Z(k)Z(i)exp(-j\omega (k-i))] \\ =\frac{1}{N}(\sum _{l=-N+1}^0\sum _{n=-l+1}^N+\sum _{l=0}^{N-1}\sum _{n=1}^{-l+N})Z(l+n)Z(n)exp(-j\omega l) \\ =\sum _{l=-N+1}^{N-1}{\hat{r}}_Z(l)exp(-j\omega l) \end{gathered}$$

  其中

$${\hat{r}}_Z(l)=\{\begin{array}{ll}\frac{1}{N}{\sum }_{n=-l+1}^{N}Z(l+n)Z(n)& l<0\\ & \\ \frac{1}{N}{\sum }_{n=1}^{-l+N}Z(l+n)Z(n)& l\ge 0\end{array}$$

  该函数具有这样的性质

$${\hat{r}}_Z(l)={\hat{r}}_Z(-l)$$
  证明

$$\begin{gathered} \text{assume }l>0 \\ \text{then }-l<0 \\ {\hat{r}}_Z(-l)=\frac{1}{N}\sum _{n=l+1}^{N}Z(n-l)Z(n) \\ \text{Let }{n}^{\prime }=n-l \\ \text{then }{\hat{r}}_Z(-l)=\frac{1}{N}\sum _{n^{\prime }=1}^{N-l}Z(n^{\prime })Z(n^{\prime }+l) \\ =\frac{1}{N}\sum _{n=1}^{N-l}Z(n)Z(n+l)={\hat{r}}_Z(l) \end{gathered}$$

  所以相关矩阵的估计是个偶函数

  最终，我们得到了周期图的另外一种表示形式

$${\hat{S}}_Z(\omega )=\sum _{l=-N+1}^{N-1}{\hat{r}}_Z(l)exp(-j\omega l)$$

5.3 周期图估计与离散时间的功率谱之间的比较

5.3.1 概述

  现在我们来比较一下周期图谱估计与离散时间功率谱之间的差异

$${\hat{S}}_Z(\omega )=\sum _{l=-N+1}^{N-1}{\hat{r}}_Z(l)exp(-j\omega l)$$

$$r_Z(l)=E(Z(k+l)Z(k))\Rightarrow S_Z(\omega )=\sum _{l=-\infty }^{+\infty }{r}_Z(l)exp(-j\omega l)$$

  首先，由于使用周期图进行谱估计的时候，我们没有求期望，因为相关函数一定是一个随机变量，因此周期图也一定是个随机变量。如果我们想对周期图的估计方法进行评估，最好是计算一下期望和方差。

5.3.2 周期图估计的期望问题

(1) 周期图估计期望的表示

$$\begin{gathered} \text{assume }l\ge 0 \\ E({\hat{S}}_Z(\omega ))=\sum _{l=-N+1}^{N-1}E({\hat{r}}_Z(l))exp(-j\omega l)(a) \end{gathered}$$

  其中

$$\begin{gathered} E({\hat{r}}_Z(l))=E(\frac{1}{N}\sum _{n=1}^{-l+N}Z(l+n)Z(n))) \\ =\frac{1}{N}\sum _{n=1}^{-l+N}E(Z(l+n)Z(n)) \\ =\frac{1}{N}\sum _{n=1}^{-l+N}{r}_Z(l)=\frac{1}{N}(N-l)r_Z(l) \\ =(1-\frac{l}{N})r_Z(l) \end{gathered}$$

  因为

$${\hat{r}}_Z(l)={\hat{r}}_Z(-l)$$

  所以，对于全部的l

$$E({\hat{r}}_Z(l))=(1-\frac{|l|}{N})r_Z(l)(b)$$

  在此，我们首先得到一个结论，就是对相关函数的估计是个有偏估计。

  将(b)代入(a)可得

$$E({\hat{S}}_Z(\omega ))=\sum _{l=-N+1}^{N-1}(1-\frac{|l|}{N})r_Z(l)exp(-j\omega l)$$

  如果N趋近于无穷，得到的结果就是离散形式的功率谱密度了

$$S_Z(\omega )=\sum _{l=-\infty }^{+\infty }{r}_Z(l)exp(-j\omega l)$$

  在这个意义上，周期图对功率谱的估计是渐进无偏的。

$$\text{Asymptotically Unbias}$$

(2) 期望的计算与窗函数

  现在我们想让谱估计的期望的表示形式与功率谱密度接近，于是我们引入一个窗函数

$$Let\omega (l)=\{\begin{array}{ll}1-\frac{|l|}{n}& -N+1\le l\le N-1\\ 0& others\end{array}$$

$$\begin{gathered} E({\hat{S}}_Z(\omega ))=\sum _{l=-N+1}^{N-1}(1-\frac{|l|}{N})r_Z(l)exp(-j\omega l) \\ =\sum _{l=-\infty }^{+\infty }\omega (l)r_Z(l)exp(-j\omega l) \end{gathered}$$

  我们可以看出，这个地方是两个函数乘积的离散傅里叶变换。而两个信号乘积的离散傅里叶变换，就等于两个信号各自傅里叶变换的卷积。因为是离散的傅里叶变换，所以卷积的范围是[-pi,pi]

$$\sum _{l=-\infty }^{+\infty }\omega (l)r_Z(l)exp(-j\omega l)=\frac{1}{2\pi }{\int }_{-\pi }^{\pi }W(\omega -{\omega }^{\prime })S_Z({\omega }^{\prime })d{\omega }^{\prime }$$

  如果我们想把这个功率谱密度的期望计算出来，就必须要计算窗函数的傅里叶变换

$$\begin{gathered} W(\omega )=\sum _{l=-\infty }^{+\infty }\omega (l)exp(-j\omega l) \\ =\sum _{l=-N+1}^{N-1}(1-\frac{|l|}{n})exp(-j\omega l) \end{gathered}$$

  其实ω(l)是一个三角窗，三角窗是矩形窗函数的相关计算得到的。我们假设h是个门函数(矩形窗)

$$\begin{gathered} W(\omega )=\sum _{l=--N+1}^{N-1}(1-\frac{|l|}{n})exp(-j\omega l) \\ =\sum _{l=-N+1}^{N-1}(h\ast h)(l)exp(-j\omega l) \\ =\sum _{l=-\infty }^{+\infty }(h\ast h)(l)exp(-j\omega l) \\ =\sum _{l=-\infty }^{+\infty }\sum _{k=-\infty }^{+\infty }h(k)h(k-l)exp(-j\omega l) \end{gathered}$$

  因为卷积的傅里叶变换是频域相乘。然后相关又和卷积相差一个符号。因此

$$\text{Let }\tilde{h}(l)=h(-l)$$

  则

$$\begin{gathered} W(\omega )=\sum _{l=-\infty }^{+\infty }\sum _{k=-\infty }^{+\infty }h(k)\tilde{h}(l-k)exp(-j\omega l) \\ =\sum _{l=-\infty }^{+\infty }h(l)⊛\tilde{h}(l)exp(-j\omega l) \\ =H(\omega )\ast \tilde{H}(\omega )=|H(\omega ){|}^2 \end{gathered}$$

  证明

$$\begin{gathered} \tilde{H}(\omega )=\sum _{l=-\infty }^{+\infty }h(-l)exp(-j\omega l) \\ =\sum _{l=-\infty }^{+\infty }h(l)exp(j\omega l) \\ =\overline{\sum _{l=-\infty }^{+\infty }h(l)exp(-j\omega l)} \\ =\overline{H(\omega )} \end{gathered}$$

  因此

$$H(\omega )\ast \tilde{H}(\omega )=|H(\omega ){|}^2$$

  窗函数 h(l)的傅里叶变换如下

$$H(\omega )=\sum _{k=1}^{N}exp(-j\omega k)$$

  用等比数列求和公式

$$H(\omega )=exp(-j\omega )\frac{exp(-j\omega N)-1}{exp(-j\omega )-1}$$

$$\begin{gathered} |H(\omega ){|}^2=\frac{4si{n}^2\frac{N\omega }{2}}{4si{n}^2\frac{\omega }{2}} \\ =\frac{si{n}^2\frac{N\omega }{2}}{si{n}^2\frac{\omega }{2}} \end{gathered}$$

  因此我们得到了窗函数的傅里叶变换

$$W(\omega )=\frac{si{n}^2\frac{N\omega }{2}}{si{n}^2\frac{\omega }{2}}$$

  我们得到的这个结果叫做Feijer Kernal

  Feijer核的函数图像表示为

  这个核第一个零点是2pi/N

  从这个核上，我们也能得到渐进无偏估计的另外一个解释。因为如果N趋近于无穷大，相当于W(ω)就是个冲激函数。冲激函数和功率谱密度的卷积，就是功率谱密度。

  如果N不是无穷大，对功率谱密度的估计和功率谱密度的区别就体现在两个方面

• 主瓣有宽度。本来窗函数的谱应该是冲激，现在有宽度了，相当于发生了谱模糊。Spectral Smearing。
• 而旁瓣我们本来不希望有，现在有了，相当于发生了泄漏Spectral Leakage。
• 如果我们想把模糊和泄漏减轻，我们就必须增大数据长度

(3) Blackman-Turkey与周期图比较

  我们知道周期图估计是一个有偏估计

$$\{\begin{array}{ll}{\hat{S}}_Z(\omega )=\frac{1}{N}|{\sum }_{k=1}^N(Z(k)exp(-j\omega k)){|}^2& \text{Physical}\\ & \\ {\hat{S}}_Z(\omega )={\sum }_{l=-N+1}^{N-1}{\hat{r}}_Z(l)exp(-j\omega l)& \text{Mathematical}\end{array}$$

$$\begin{gathered} E({\hat{r}}_Z(l))=E(\frac{1}{N}\sum _{n=1}^{-l+N}Z(l+n)Z(n))) \\ =(1-\frac{|l|}{N})r_Z(l) \end{gathered}$$

  如果我们在做谱估计的时候引入权重，其实就能消除这个偏差，比如

$$\begin{gathered} {\hat{S}}_Z(\omega )=\sum _{l=-N+1}^{N-1}\omega (l){\hat{r}}_Z(l)exp(-j\omega l) \\ \Rightarrow {\hat{S}}_Z(\omega )=\sum _{l=-N+1}^{N-1}\frac{N}{N-l}{\hat{r}}_Z(l)exp(-j\omega l) \\ \text{Blackman-Turkey} \end{gathered}$$

  这种方法叫做Blackman-Turkey，其中Turkey是发明FFT的人。

  引入这个权重以后，最终得到的就不是一个三角窗了，而是一个矩形窗。当然，引入权重的方法还有很多，比如加汉明窗等

  但是人们其实并不是非常希望通过在周期图中加入权重，得到无偏估计。因为效率问题。

  没有权重的周期图物理的形式如下

$${\hat{S}}_Z(\omega )=\frac{1}{N}|\sum _{k=1}^N(Z(k)exp(-j\omega k)){|}^2$$

  这里是可以直接使用FFT做运算的。FFT时间复杂度是O(n²)因此，周期图的效率是最高的。而如果引入权重进行主动加权的话，要从数学形式去做，首先要计算 \hat r_Z(l),计算这个至少需要O(n²)的时间复杂度

$${\hat{r}}_Z(l)=\{\begin{array}{ll}\frac{1}{N}{\sum }_{n=-l+1}^{N}Z(l+n)Z(n)& l<0\\ & \\ \frac{1}{N}{\sum }_{n=1}^{-l+N}Z(l+n)Z(n)& l\ge 0\end{array}$$

  因此很多时候，人们宁愿选择渐进无偏的周期图估计，也不选择Blackman-Turkey这种无偏估计。

5.3.3 周期图估计的方差问题

(1) 周期图方差的求解

  算清楚均值之后，我们就该分析周期图的方差了

$$\begin{gathered} Var({\hat{S}}_Z(\omega ))=E[({\hat{S}}_Z(\omega )-E({\hat{S}}_Z(\omega ){)}^2] \\ =E[{\hat{S}}_Z(\omega {)}^2]-2E[{\hat{S}}_Z(\omega )\ast E({\hat{S}}_Z(\omega ))]+E[E({\hat{S}}_Z(\omega ){)}^2] \\ =E[{\hat{S}}_Z(\omega {)}^2]-2E({\hat{S}}_Z(\omega ))\ast E({\hat{S}}_Z(\omega ))+[E({\hat{S}}_Z(\omega )){]}^2 \\ =E[{\hat{S}}_Z^2(\omega )]-[E({\hat{S}}_Z(\omega )){]}^2 \end{gathered}$$

  我们前面分析均值的时候，已经做了谱估计的期望，然后平方即可得到第二项。现在只要分析谱估计平方的期望即可

  这里我们要做个假定,假定数据是高斯白噪声，并且规定二阶矩。

$$\begin{gathered} \text{Assume }Z\sim \text{ Gaussian White Noise } \\ E(Z(m)Z(n))={\sigma }^2{\delta }_{nm} \\ {\delta }_{nm}=\{\begin{array}{ll}1& m=n\\ 0& m\cancel{=}n\end{array} \end{gathered}$$

  则我们计算一下谱估计平方的期望

$$\begin{gathered} E[{\hat{S}}_Z^2(\omega )]=E(\frac{1}{N^2}|\sum _{k=1}^N(Z(k)exp(-j\omega k)){|}^4) \\ =\frac{1}{N^2}\sum _{k=1}^N\sum _{n=1}^N\sum _{i=1}^N\sum _{l=1}^{N}E(Z(k)Z(n)Z(i)Z(l))exp(-j\omega (k-n)-j\omega (i-l)) \end{gathered}$$

  高斯的四阶矩可以做如下的变换,变成三个乘积的和

$$\begin{gathered} E(Z(k)Z(n)Z(i)Z(l))=E(Z(k)Z(n))E(Z(i)Z(l))+E(Z(k)Z(i))E(Z(n)Z(l))+E(Z(k)Z(l))E(Z(i)Z(n)) \\ ={\sigma }^2{\delta }_{kn}\ast {\sigma }^2{\delta }_{il}+{\sigma }^2{\delta }_{ki}\ast {\sigma }^2{\delta }_{nl}+{\sigma }^2{\delta }_{kl}\ast {\sigma }^2{\delta }_{in} \end{gathered}$$

  因此我们原来的式子可以写成

$$\begin{gathered} E[{\hat{S}}_Z^2(\omega )]=\frac{1}{N^2}\sum _{k=1}^N\sum _{n=1}^N\sum _{i=1}^N\sum _{l=1}^N({\sigma }^2{\delta }_{kn}\ast {\sigma }^2{\delta }_{il}+{\sigma }^2{\delta }_{ki}\ast {\sigma }^2{\delta }_{nl}+{\sigma }^2{\delta }_{kl}\ast {\sigma }^2{\delta }_{in})exp(-j\omega (k-n)-j\omega (i-l)) \\ =\frac{1}{N^2}(\sum _{k=1}^N\sum _{i=1}^N{\sigma }^4+\sum _{k=1}^N\sum _{n=1}^N{\sigma }^{4}exp(-2j\omega (k-n))+\sum _{k=1}^N\sum _{i=1}^N{\sigma }^4) \\ =2{\sigma }^4+\sum _{k=1}^N\sum _{n=1}^N\frac{{\sigma }^4}{N^2}exp(-2j\omega (k-n)) \end{gathered}$$

  因为后面一部分是对复指函数求和，N趋近于无穷大的时候一定是有界的，所以上式当N趋近于无穷的时候极限是可求的

$$E[{\hat{S}}_Z^2(\omega )]=2{\sigma }^4+\frac{{\sigma }^4}{N^2}\ast C\underrightarrow{N\to \infty }2{\sigma }^4$$

(2) 周期图方差的弊端

  算到这里，我们使用周期图进行谱估计的缺点就暴露了出来。因为随着样本量的增加，方差并不会减小。也就是方差是不相合的。周期图谱估计的方差并不能趋近于0

  我们可以自己写程序去验证这个事情。用matlab实现周期图谱估计非常容易

N = 1024;
plot((abs(fft(randn(1,N))).^2) /N);

N = 16384;

plot((abs(fft(randn(1,N))).^2) /N);

  我们发现样本数量增加，并不能使得谱估计样本的方差减小。

(3) 周期图求方差的改进-Bartlett & Welch

  为了能够减小方差，我们唯一能做的事情，就是对采集到的数据分段，分别做周期图谱估计，然后求平均值

$$Var(\frac{S_1+...S_N}{n})=\frac{Var(S_1)}{n}$$

  但是这个时候我们的数据的长度就减小了，在均值上的损失就变大了，我们又遇到了bias variance tradeoff问题

  这是我们第三次遇到这个问题了

• 随机误差与系统误差
• 过拟合
• 数据分段进行谱估计

  我们这种降低谱估计方差的方法叫做Bartlett

$$Bartlett$$

  我们用matalb做仿真看看得到的结果。这是把数据分为8组合并求谱估计和分开求谱估计

$$Welch$$

  除此之外，还有一个方法叫做welch，兼顾了均值和方差，就是在分段取数据的时候，重叠着取

  其实这样并不好，因为重叠着采样，数据是有相关性的，求方差的时候达不到原来方差的1/n,而且数据重叠性越高，对方差的修正能力越差

5.4 周期图谱估计小结

  到目前为止，周期图是应用广泛的功率谱估计方法。电子系的用这个比较low，其他专业的用这个比较多

  周期图是成熟技术。所谓成熟技术，都是我们对其优点和缺点都有很详细的认识。

• 优点：性能好，应用广泛
• 缺点：均值上有谱模糊和谱泄漏问题，方差上有不相合的问题

---