思考

天文学是一门观测主导的科学,数据处理方法是所有天文工作者的必修课,而傅立叶变换作为最常见的数据处理方法,是处理天文观测中处理时序序列的基本方法。本文由我上学期《天文数据处理》这门课中傅立叶变换部分的笔记整理而来,动态分析了一个“傅立叶变换及其逆变换”过程中我们有意无意中进行的各种运算,算是我的一点学习心得吧。

零、全局符号意义说明

$N$:采样点数目

$\Delta_t$:时域采样间隔

$T_0=N\cdot \Delta_t$:时域采样长度

$W_0=\frac{2\pi}{\Delta_t}$:频域采样长度

$\Delta_\omega=\frac{W_0}{N}=\frac{2\pi}{T_0}$:频域采样间隔

其具体意义会在下文中说明

一、傅立叶变换(单变量)

傅立叶变换是一种积分变换,如果使用时间作为某个单变量函数的自变量时,傅立叶变换把它“分解”为不同振幅的一系列简协运动的叠加,公式表现为:

(需要说明的是,傅立叶变换有好几种定义方式,本文讨论以上这种形式,但本文中的思适用于所有的形式)

满足这样关系的一对函数称为一个傅立叶变换对,通常记做

如果$f(t)$是一个周期函数,那么有

离散傅里叶变换:

各个符号如第零章的意义,n,k分别表示时间的第n个采样点和频率的第k个采样点。

\begin{equation*} X(k)=\text{DFS}(x(n))=\sum_{n=0}^{N-1}x(n)\text{e}^{-\text{i}\frac{2\pi}{N}nk} \end{equation*}

\begin{equation} x(n)=\text{IDFS}(X(k))=\frac{1}{N}\sum_{k=0}^{N-1}X(k)\text{e}^{\text{i}\frac{2\pi}{N}nk} \end{equation}

二、傅里叶变换的几个性质

如果:

则有:

对偶性质:

时移性质:

频移性质:

卷积性质(见下文):

三、常见函数的傅立叶变换对

1、rect函数:

它的变种矩形窗函数:

由傅立叶变换的性质$f(t-a)\Longleftrightarrow {\text{e}^{-\text{i}a\omega}F(\omega)}和f(t)\cdot \text{e}^{\text{i}\omega_0 t}\Longleftrightarrow {F(\omega-\omega_0)}$

2、$\delta函数$(单位冲击函数):

满足:

3、冲击串函数:

满足:

四、卷积

卷积的定义通常写作

对于两个函数的卷积可以(从物理上)这样理解:

我们把一个函数$f(t)$理解为一种输入,这种输入必须具有线性可加性,比如往瓶子中倒水,我给出一个$f(t)$为t时刻的水流流量,则流量具有线性可加性,$f(t)$对t的积分就是倒入水的体积。再比如背单词,其输入就是背单词的个数,背一个单词就多会一个单词,具有线性可加性。

我们把卷积理解为$f(t)$在t上的一个积累,另外一个函数$g(t)$理解为“响应密度函数”,类似概率密度函数。

还说背单词,假如一个人背单词的速率是每分钟5个,那么另$f(t)=5$为其背单词的输入函数。我们定义$g(t)$为其遗忘函数,就是他在看了一个单词后t时刻还能记住的百分比,我们称之为其记忆对t时刻之前输入的“响应”,或者说能记住的t时刻之前的单词百分比,艾宾浩斯记忆曲线就是遗忘函数的一个示例,由于找不到这个曲线的具体形式,下面的例子中用作为遗忘函数,这样的函数有一个特点,就是$t<0$是其值为0,这很好理解,在记忆对“负的时刻之前(其实就是未来时刻)”的响应必须为零,他还没记忆这个东西,那么他就不会记住,我们称这样的相应密度函数为“因果的”,就是说现在t=0时刻的记忆状态,只由过去的输入决定,和将来无关。

怎么计算现在记住的单词个数呢?只需要把所有的输入,乘以对那个输入的记忆百分比,累加起来就是。离散的看来,我把输入分为好多小区间,近似认为在每个区间上的记忆率为定值,对时间轴上的所有区间乘以他们的记忆率累加就和就是现在的记忆量,如果区间取的无限小,求和变为积分,遗忘函数就变成了密度函数。回头看卷积的形式$h(\tau)=f(t) * g(t) = \int_{-\infty}^{\infty}f(t)g(\tau-t) \text{d}t$,令$\tau$为现在的时刻,$f(t)$为t时刻的输入,$g(\tau-t)$为$\tau-t$时间之前的响应密度,这里注意$\tau-t>0$为对过去的响应,$\tau-t<0$为对未来的响应,然后对时间积分,就是现在的总响应。

俗话说的好,没图你说个XX,现在上图。(下图为8M的GIF)

图:其中蓝色虚线表示LearnCurve,底部的黑色虚线为ForgetCurve,绿色部分为\(\text{learnCurve}\times \text{forgetCurve}(\tau-t)\)

红色实线为卷积结果,注意到红线在一点的取值的取值为那时绿色区域的面积,也就是过去对过去输入强度按照(关于Y轴反射的)响应密度函数进行加权的加权和

此图展示了我们的上面的记忆模型。输入为图中LearnCurve所示曲线

就是前20分钟每分钟输入5个单词,后来不专心了,效率越来越低,到40分钟就背不下去了……

他的遗忘函数为分段指数函数,当然实际中我们的记忆力不会这么差…

图:遗忘函数ForgetCurve

卷积的结果作为的函数,表示时刻的记忆总量。

值得注意的是,如果交换f和g,则这个运算就没有直观的物理意义,而最终结果却也为记忆总量,看来正确恰当的理解方式,对于分析问题至关重要。

函数的卷积:

可以这样理解:函数其实是“即时响应函数”,在时刻只对时刻本身的信号有“单位冲击响应”,与其他时间的信号无关。

与冲击串函数的卷积:

冲击串函数其实就是一串冲击函数(汗…),其卷积效果其实就是“无限等间距平移相加”,如图

图:一个函数,令他为\(f(t)\)

图:\(f(t)\)以15为间隔无限平移

图:红色实线是f(t)函数的无限平移相加,其实就是\(f(t) * \delta_{15}(t)\),“边界”处出现的“误差”称为混叠效应

五、实例

一般以处理数据为目而学习傅立叶变换的学习时间也就几个课时:几个定义,几个性质,稍作练习后就转到了其应用:使用FFT(快速傅立叶变换)处理各种数据,画频谱,一看:“哇,图真漂亮!”,但我们的注意力此时却是感性的,看到了图的相对形状,而搞不清他的绝对坐标的意义。

本节要解决的问题描述如下:

假如我们对一个物理过程感兴趣,我们要观察和分析这个过程中的某个量。

现实世界中的物理过程一般都为连续的,有一个函数可以精确的完全的去描述他,但是的具体形式也许只有拉普拉斯妖才知道,我们称为这个过程的“原始形式”或“上帝形式”。我们对其观察和采样,得到实验数据:一堆数据点。这里不考虑实验误差,我们得到的数据其实就是在若干t时刻的精确值,称这个数据为过程的“采样形式”。然后把数据扔给FFT程序,得到一个“频谱”,把这个“频谱”所对应的真实物理过程称为“逆变换形式”,我们要探讨的就是“上帝形式”、“采样形式”和“逆变换形式”的区别和联系。

(说明,由于频谱值一般为复数,下面的频谱图Y轴均为\(\text{Abs}(F(\omega)))\)

现在我是上帝,我知道某个物理过程的信号如图所示,然后我还知道他的“真实”频谱

图:\(f(t)\Longleftrightarrow F(\omega)\)

啊,愚蠢的人类,你们没有能力去观察整个时间轴上的数据,你们只能获得这有限时间上的数据,正所谓开窗取样

等效结果就是乘以了一个如图所示的窗函数

图:\(\text{timeWindow}(t,T_0)\Longleftrightarrow {\text{ItimeWindow}(\omega)}\)

我们来复习一下傅立叶变换的卷积性质:,时域的相乘对应频域的卷积。加过窗后的函数频谱等于原频谱与窗函数频谱的卷积。

如果我们不加窗,或者认为加了一个“单位窗”,则窗函数的频谱为结果不变

现在看看我们的窗函数的样子:其对附近的相应比较大,而对的地方也有相应,相应不是因果的,卷积的效果就使我们的频谱失真了,某个频率的值受到其附近频率的影响,我们称之为频谱泄漏,窗开的越大,窗函数就越接近函数,频谱泄漏现象就越小,但是开无限大的窗是不可能的,这就是我们和“上帝频谱”之间的第一个差别。

然后“加过窗”的新函数及其傅立叶频谱如图

图:\(f(t)\cdot{\text{timeWindow}(t,T_0)}\Longleftrightarrow F(\omega) * (\frac{1}{2\pi}\cdot \text{ItiemWindow}(\omega))\)

(下文中令\(f(t)\cdot{timeWindow(t,T_0)}=f_1(t)\),\(F(\omega) * (\frac{1}{2\pi}\cdot ItimeWindow(\omega))=F_1(\omega))\)

红色虚线是“上帝频谱”,蓝色实线加过窗的频谱,其中的虚假波动是由频谱泄露引起的

而后呢,即使是在有限的时间区域内,我们也不能得到连续的,只能对他进行采样,所谓的采样,就是乘以一个单位冲击串函数

“等等….”,你会说:“我的采样明明是在一个时间序列上的每个点的一个数值啊,不是一串奇怪的函数啊?”

“额…这是个好问题,不过先..别管他吧,我们就看看我们这样做能得到什么吧…”

“好吧….”

再复习一下,单位冲击串函数的傅立叶变换是另外一个非单位的冲击串函数

那么时域上乘以一个冲击串函数,在频域上对应的就是与另外一个冲击串函数的卷积(好像还有其他奇怪的东西,比如W_0这个系数),如图

图:\(f_2(t)\),黑色箭头示意采样点

上图的蓝色实线和下图的黑线分别表示:\(f_1(t)\cdot \delta_{\Delta_t}(t)\Longleftrightarrow F_1(\omega) * \frac{1}{2\pi}\cdot W_0\delta_{W_0}(w)\)

下图中绿线为\(F_1(\omega)\)的无限平移,红线为\frac{1}{2\pi}\cdot W_0 F_1(\omega)\)的无限平移,黑线为卷积结果

进一步计算得

(下文中令\(f_1(t)\cdot \delta_{\Delta_t}(t)=f_2(t),F_1(\omega) * \frac{1}{\Delta_t}\delta_{W_0}(w)=F_2(\omega))\)

左边与冲击串的卷积相当于“无限平移相加”

A:“哈哈,加窗、采样后我们进行傅立叶变换的到的频谱,除了稍微有点失真、变成了一个奇怪的周期函数、还有乘以了一个奇怪的系数以外,和上帝模式的形状长得差不多么,啦啦啦~~”

B:“咳咳…好了,不管你这是什么乱七八糟的东西,赶紧给我频谱数据,老板还等着我的结果呢,你这个数据好像还得乘以吧,不然和上帝模式的差好多倍呢”

A:“现在这个有点奇怪的频谱结果如果直接逆变换回去,会得到一个冲击函数序列,每个点的冲击函数和那个点你的采样值倍的函数一样,并不是你原来的采样序列,恩,是有点奇怪,不管了,就把数据给你吧,不过….愚蠢的人类,这是一个无限长的周期频谱啊,你只能取其中的一段”

B:“既然这个一个周期函数,我就取它的一个周期吧,我要这段的数据”

A:“恩,不错,你相当与在频域上加了一个长度为的窗,我们令你这个加过窗的频谱为

但是你胡乱加这个窗后再逆变回去后的结果是个啥啊?

你看我还记得,这频域里加窗等效于时域中的卷积,但频域中timeWindow窗的傅里叶逆变换甚至都不是实数,它与时域中冲击串的卷积相当于’加权无限平移相加’,卷积结果甚至都不一定是实数,你加过窗的频谱还能用么?”

B:“啥?…..那我们就先看看频域中timeWindow窗逆傅里叶变换以后的形式ItimeWindow吧..如图”

图:ItimeWindow(w)

我们注意到,这个函数除了在处不为0外,每隔一个间隔就取值为0,而这个间隔正好就是(证明留作习题)!而我们这个函数是要和一个等间隔的(非单位)冲击函数卷积的,我们已经知道和一个非单位的冲串函数卷积的结果是其本身乘以冲击函数前面的系数,而和一个非单位的冲击串函数卷积则是重复上面的过程,平移相加,另外计算可知(通过取极限得到),这样一来,卷积的结果就是在的冲击值处取值,我们的终于从一个奇怪的冲击串函数变成了一个可以正常取值的常数

图:上下蓝色实线分别画的是\(F_2(\omega)\cdot timeWindow(w,W_0)\)和\(f_{20}(t)=\frac{2\pi}{W_0}*f_2(t) *ItimeWindow(t)\)

\(f_{20}(t)\)最前面的系数为了使其形状和\(f_2(t)\)比较容易比较,而这个傅里叶变换对的真正形式应该如下

图中红线为\(f_1(t)\),黑色箭头为\(f_2(t)\)

绿色实线为ItimeWindow(t)与某个采样点处\(\delta\)函数的卷积,

可以看出其在所有其他采样点处的模为零,不影响这些采样点处的值

B:“好了,我们这一对傅里叶变换总算稍微正常了一点,两边长的都像是正常的函数了,不过比我们理想中的结果还多了个系数,我们是不是考虑把右边的这个频谱再乘个后再给我?正如上图中的处理一样?”

A:“别急,你不是要数据么,右边的有限区间的连续函数还是没法给你,我们还得采样。”

B:“额,既然也要采样,那么就和之前T的采样一样多得了,也取N个值,采样区间就是吧”

A:“好的,我们给频域乘以一个单位冲击串函数

B:“什么?好不容易两端都是正常的函数了,你又把右边弄成了奇怪的冲击串?”

A:“好问题+1….我们不管他,变过去看看什么样再说吧….”

B:“好….吧….别给我说你科研都是这么搞的”

A:“(躺枪….)”

上图之前先简单推导点东西:注意到的各个值是从上帝模式的精确采样过来的(然后又乘以了一个单位冲击函数…)

为了比较上式左边和上帝模式的形状,我们进一步演算

于是令

图:红色虚线和黑色箭头同上,蓝色实线表示\(f_{40}(t)\)

图:这两个图是在周期开始的0和周期结束的15附近的特写,注意到下图最右的采样点似乎不对,但那个箭头其实是第N+1和采样点,而其值应该与第一个采样点,既\(t=0\)那一点相等,事实确实如此

A:“大功告成,你看我右边的频谱终于可以逆变换得到和上帝模式长的一模一样的函数了~~额,不对,是在你的采样点处和上帝模式长的一模一样,但在其他点处跟上帝模式差别就大了…..还有,我们的图是逆变换过去后又除以得到的,恩,你怎么看?”

B:“尼玛我要数据!我要频谱!你这左边又变成了一个奇怪的有限长的冲击串,你给我个啥啊?”

A:“…..我们看看是怎么变回左边的吧

\begin{eqnarray} f_4(t)&=&\frac{1}{2\pi}\int_{-\infty}^{\infty}F_4(\omega)\cdot \text{e}^{\text{i}wt}\text{d}\omega\ &=&\frac{1}{2\pi}\sum_{k=0}^{N-1}F_3(k\Delta_\omega)\cdot \text{e}^{\text{i}k\Delta_\omega t}\ &=&\frac{N}{2\pi}f_2(t)\ &t=n\cdot \Delta_t , n=0,1,2…N-1&
\end{eqnarray}

你看这样吧,我把在采样点处的值(而不是其对应的冲击函数)给你,你回去使用

这种逆变换方法就能变回上帝模式在各个t的采样处的值了,咦,这好像就是吧?擦….”

B:“你这半天又是卷积又是加窗,还用什么奇怪的冲击函数去所谓的采样,早告诉我用matlab中的fft和ifft不就得了,你看我扔给他数据就直接出来了和你这折腾半天一样的结果,直接plot,这图长的多漂亮啊,这就是频谱啦?拿数据走人,多谢不送了~(闪人)”

A:“喂,你到底弄清楚了你得到的是什么东西了么?人呢?…”

A:“算了,自己总结个吧….一般搞实验的除了自己谁都得相信你的结果,我还是自己再检查一下吧”

总结:

一般来说我们“进行傅里叶变换”指的是对一信号进行采样、得到离散数据点,这个过程中我们相当于对原信号进行了:1、加窗;2、采样两个运算。这个“采样”其实是指乘以一个单位的冲激串函数,这样过后,我们得到的这个采样数据所对应的傅里叶频谱,与我们求之不得的“上帝频谱”相比,存在三点不同:1、有加窗引入的频谱泄露;2、由采样间隔不能无限小造成的频谱混叠;3、得到的频谱在以上两个误差的基础上又多乘了一个系数

为了从我们得到的还不错的频谱中还原出原来的采样值,我们同样对得到的频谱进行了加窗和采样,但由于加窗长度和采样频率的特殊性,我们逆变换回去的值并没有出现类似频谱混叠和频谱泄露的现象,只是比原来的采样值多了一个系数,如果我们在逆变换的时候不使用连续傅里叶逆变换的公式,而是使用离散逆变换的公式,那么按照公式(??)可以直接逆变换为上帝模式的采样值。

现在回答最后一个问题:为什么用奇怪的冲激串函数进行采样

答案是:采完样后进行傅里叶或逆傅里叶变换,一积分,就变成了和离散傅里叶变换同样的形式,相当于在采样点处取值,相加。