从傅里叶变换到加窗傅里叶变换

最近结束了一门课,叫《基于波动理论的目标探测》。名字如此高深,实际上还是还是在讲信号处理的一些知识。其中的核心内容是加窗傅里叶变换(或者叫短时傅里叶变换)。这里也想说一说自己的一些理解。

从傅里叶变换到加窗傅里叶变换

傅里叶变换是我们所熟悉的,它把我们所要分析的信号从时间域变换到了频率域,这样最大的好处是能让我们看的更清楚信号是由哪些些基本的“原子”所组成的。这些“原子”其实就是三角级数,或者说是我们所熟悉的: ejωt ,这样的变换能够告诉我们原始信号的频谱,这无疑是非常好的工具。但我们现在来看个例子:
在matlab的命令窗口中输入xpsound我们观察其中的任意一个信号都会发现,尽管我们能够从时间上一眼看处声音的频率其实是在一个局部有一个固定的频率而在另外的一个局部又有另外一个固定的频率。然而在其中我们看每一个信号的幅度谱(或者能量谱)都无法看处这一点。
我们再来说的更形象一些,如果一个美妙的音乐片段,它是在不同的时间由不同的音符所组成的,我们清楚地知道在不同位置(时间)有不同的音符(频率)。然而我们却不能够通过傅里叶变换的手段来轻松获得这一分析。原因正是因为我们取得的时间片段太长了,以至于这些不同频率的内容都混杂在了一起所以不能够区分开来。那么解决的办法也很简单,就是把原来的片段“切成”每个音符的大小(假定每个音符持续的时间都一样)的小片段,再在每个小片段上做傅里叶变换。这样我们能够直接定位到频率随时间的变化。这种想法和思路实际上就是加窗傅里叶变换。从傅里叶变换到加窗傅里叶变换就是这样过度过来的。

如何实现加窗傅里叶变换(WFT)

后面为了书写方便,我们暂时都把加窗傅里叶变换叫做是WFT。有了上面的一些基础的加窗傅里叶变化的由来,其实还是很容易想到实现加窗傅里叶变换的方法的。
我们还是从傅里叶变换(FT)出发来导出WFT的表达式,在FT中我们熟悉它的表达式是

F(jω)=f(t)ejωtdt

根据上面所说的WFT的由来,其实我们是不难推出它的表达式的:
首先,待分析的信号是 f(t)L2(R) 它的意思就是说 f(t) 是一个能量有限的信号,窗函数 g(t)L2(R) 也是能量有限的信号。我们记 g¯(t) 是窗函数 g(t) 的共轭(这里我们可以先把窗理解成是一个实函数,它的共轭和它本身是一回事儿,写成共轭形式对于一会儿另外一个概念有用)。
那么WFT的思路就是,先用窗的不同的时移去截取信号,即
ft(u)=g¯(ut)f(u)

再对截取的信号取傅里叶变换,我们把变换的结果记为 f˜(ω,t) ,则变换结果为(这里的频率单位是Hz,并且我们还是记为 ω ,下同)
f˜(ω,t)=g¯(ut)f(u)e2πjωudu

这样我们可以想象得来此时的 f˜(ω,t) 是一个二维信号,并且它能够从频率和时间两个维度去刻画一个信号。而且一个直观的感受是,我们所说的“窗”如果取得越细密,那么好像可以刻画出信号的刻画的更加完美,但实际上一个极端特例就是“窗”细的不能更细,变成了一个冲击,那么不就变成了连续不断的采样?这么来说在t这个维度上看,每个t都是一个白谱,那么WFT也就失去了意义。这告诉我们,显然不可以把这个时间的“窗”开的很细,否则频率域去看会是很大一片(分辨率极低),这也从另一个角度刻画了时间和频率只能折衷,而不能二者都取得性能极佳!

时频局域化与不确定性定理

实际上从数学的表达式上,也可以去刻画上面我们所描述的时频分辨率这一特性。
我们先定义两个概念:时间的“重心”和频率的“重心”。
对于窗函数 g(t) 以及它对应的傅里叶变换 g^(ω) 分别有如下的定义:

t0=1g(t)2t|g(t)|2dt

ω0=1g^(ω)2ω|g^(ω)|2dω

窗嘛一般来说都是有一个能量聚集中心的,不可能是像噪声那样不规则分布的,所以有一个指标就很重要了,它能够从某种角度说明“分辨率”这个概念,它就是标准差。直观的想如果标准差越大,那么在这个维度上它的能量是分散,所以分辨率自然也不会高,反之亦然。
时间和频率的标准差分别如下
T2=(tt0)2|g(t)|2dt|g(t)|2dt

Ω2=(ωω0)2|g^(ω)|2dω|g^(ω)|2dω

下面这个图叫时频盒子图,它形象地用图像的形式给出了描述了上述的时间和频率分辨率的方法
Time-Frequency boxes(Heisenberg rectangles)
显然我们都希望在时间和频率的分辨率都能够越高越好,但是瓶颈在于我们已经知道不可能达到两者的分辨率都很好,所以我们试图去折衷。
如果能够让 TΩ 尽可能的小,那么我们就认为达到了目的。已经证明了,当窗函数 g(t)=et2 中所谓的高斯窗的形式时,能够让 TΩ 尽可能的小,譬如 g(t)=(2a)14eπat2 他满足能量的归一化,即 g=1 。在这种情况下 g^(ω)=(2a)2eπω2/a ,并且 TΩ=14π 对于其他形式的窗函数, TΩ>14π 。这也正是我们窗函数为什么要选择高斯窗的原因。

从“原子”的角度来理解WFT

在上面我们说道傅里叶变换是要看看一个信号是由哪些基本的“原子”所构成的。这也从另一个角度说明了傅里叶变换的本质——做投影。把信号投影在一个个“相互正交”的原子上,就获得了信号是由每个多大的“原子”所构成的。那么WFT既然是从FT演化而来的,我们能不能以相同的角度来理解呢?答案是肯定的。
再回到WFT的表达式:

f˜(ω,t)=g¯(ut)f(u)e2πjωudu

对比傅里叶变换的表达式:
F(jω)=f(t)ejωtdt

对于FT它的原子是 ejωt ,那么很自然的WFT的原子就是:
gω,t(u)=e2πjωug(ut)

同样的,其实求解WFT的过程也是一个“投影”(内积)的过程,在 L2(R) 空间中,内积的表达式为:
<f,g>=f(t)g¯(t)dt

他有三个比较重要的性质,分别是三角不等式:
|<f,g>|fg

Plancherel定理(其实就是能量守恒定理):
f=f^

两个域中的能量是一定的。
以及它的一个推广Parseval恒等式(时域内积等于频域内积):
<f,g>=<f^,g^>

有了以上的数学铺垫,我们可以轻松写出WFT的投影形式:
f˜(ω,t)=<f(u),gω,t(u)>

并且由Parseval定理,我们还可以得到从频率域来计算WFT
f˜(ω,t)=<f˜,g˜ω,t>

其中的 f˜ g˜ω,t 分别是 f(u) gω,t(u) 对应的傅里叶变换。
那么,从频率域得到WFT可以表述为:
<f˜,g˜ω,t>=f˜(v)g˜¯ω,t(v)dv

=f˜(v)e2πjvtg˜¯ω,t(vω)dv

看这个形式是不是非常像有一个频率的窗叫 g˜ω,t ,它在频率域与原始信号相乘(加窗),然后再做一个傅里叶变换呢(此时的积分变量是频率,此处我们用v来表示)?
以上内容说明了一个事实,那就是WFT仍然是一个分解原子的过程,只不过原子不再是之前FT中简单的原子了,而变成了窗函数的时移再乘以一个复数调制 e2πjωt

关于高斯窗

最后关于WFT,我们再来说说为什么要选用高斯窗。上面已经给过了一个最有力的条件,但是其实他还有几个好处。

  1. 它是一种平滑的函数,这相比于一个我们熟悉的矩形窗就要好很多,它能够避免频率域的震荡。
  2. 它是一种从零开始平滑到极值再平滑地衰落到0的函数,这符合一个我们常见信号的特点——平滑到增长到极大值然后再平滑地衰落。
  3. 它的FT也是一个高斯函数。
  4. 它是偶对称的,有很多非常良好的性质。

其实高斯窗的应用非常广泛,图像处理中也经常使用到它。我个人认为,它是一种最接近实际真实信号的一个有力模型,它的很多性质比如平滑,以及增加和衰落的过程都能比其他信号更好地描述真实情况。这是它应用广泛的原因。 (注:我所说的常见信号,比如像音乐,或者是图像的亮度信号,他们都有一定的持续期,但是由于惯性,它们都不可能立刻跳到极值(或者稳定值),都有一个充能以及在衰减时候释放能量的过程。选择高斯窗应该说会更好地贴合这个特点。)