傅立叶变换

傅立叶变换(FT)的标准形式是积分,不过实际应用中,数据都是离散的,所以这里只研究他的离散形式:

\[X_k = \sum_{n=0}^{N-1} x_n\cdot e^{-i2\pi kn/N }\]

\[x_n = \frac{1}{N} \sum_{k=0}^{N-1} X_k\cdot e^{i2\pi kn/N}\]

离散傅立叶变换(DFT)的第一个公式的输入是时间维度的量,第二个是第一个的逆变换。

数学公式总是令人苦恼,更别说如此复杂的。

不过不要紧,理解傅立叶变换并不意味着一定要理解他的公式,从实际的事例去理解他反而更加直观。所以我们直接用他去处理数据,获取输出,通过各种直观的数据可视化方式去理解傅立叶变换到底干了什么。

实现

开源库 FFTW3 提供了快速成熟的傅立叶变换实现,后面我会直接用它来进行各种实验。

下面是 FFTW3 官方文档的一个一维傅立叶变换的示例代码:

#include <fftw3.h>
...
{
    fftw_complex *in, *out;
    fftw_plan p;
    ...
    in = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * N);
    out = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * N);
    p = fftw_plan_dft_1d(N, in, out, FFTW_FORWARD, FFTW_ESTIMATE);
    ...
    fftw_execute(p); /* repeat as needed */
    ...
    fftw_destroy_plan(p);
    fftw_free(in);
    fftw_free(out);
}

这个过程并不复杂:首先准备好输入输出缓冲区,再准备好输入数据,之后创建一个 FFTW 的 plan,可以把这个 plan 理解为上面的数学公式,execute 就相当于是把数据带入这个公式进行实际的计算,最后输出缓冲区内就得到傅立叶变换的结果了。

创建 plan 的时候第三个参数有点意思,有 FORWARD 的话,是不是就有 BACKWARD 的呢?确实如此。事实是,傅立叶变换是可逆的,大致观察一下最上面的公式就能发现他们两个非常相似。

可以注意到,DFT 的输入输出都是复数。回忆以前的数学知识,复数实际就是两个数组成的一种特殊数:一个是实数,一个是虚数。简单的,我们可以把他看成两个实数的组合,只是两个的意义不一样罢了。所以,后面会对这两个数单独进行观察、分析。

输入数据

接下来,得构造数据来进行实验。很简单,我们可以使用一个简单的单个参数的函数来填充傅立叶变换的输入的实数部分,因为实际中一般不会使用复数, 现在我们也不关心他,所以输入数据的虚数部分都是 0。

第一次使用的是一个正弦函数:

\[f(t)=sin(4t \times 2\pi)\]

这个函数的输入是时间,我们以 1/100 秒的间隔,从 0.0 秒开始,调用这个函数 100 次,得到 100 个数据。

一般的,把调用函数的这个过程称为「采样」,函数的结果称为「样本」,采样间隔一定是固定的,所以,1 秒内得到的样本数量称为「采样频率」。而调用的那个函数称之为「信号」。

然后把这些数据交给 FFTW,看看输出是个什么样子。

这样还没完,既然傅立叶变换还有个逆,那不妨把输出‘逆’过去,看看是什么效果。

输出

如果输出只是一堆数字,那就太无聊了,而且也看不出什么规律。所以可以借助可视化工具,把数据绘制成图表,就有趣多了。

下面是输入数据的折线图:

一个常见的正弦曲线。可以看到,这个曲线在 1 秒的间隔内,重复了 4 次,也就是它的频率为 4Hz, 当然了,曲线的最高最低处分别为 1、-1,也就是振幅为 1。这和输入函数非常吻合

输出,分别是实部和虚部:

几乎都是 0
大部分都是 0, 但是两个尖尖着实吸引眼球:(0.04, -50),(0.96, 50)

实数部分都非常接近 0, 所以暂且可以忽略他。虚数部分整条折线几乎非常平坦,只有两个尖尖有点‘不合时宜’。 而且这两个折线似乎是对称的,实数从中间开始左右对称,虚数也是从中间开始左右奇对称(一下一上)。不过右边差了一个数据,这可能是因为输入数据并没有在 t = 1.0 的时候调用(100 个数据是从 0 开始的)。

所以,就目前的情况看,DFT 的输出只有左边半部分有实际参考意义,右边是重复多余的。

最后是逆变换的结果:

居然和原始输入一模一样!
也几乎都是 0

逆变换的结果居然和输入一模一样(误差非常小,可以忽略)。

对比

但看一组数据很难看出什么名堂,所以接下来要修改一下那个输入函数,对比着看看输出是如何变化的。而且通过前面的第一次实验,输出的虚数部分是最吸引注意的,所以后面会专注于输出的虚数。

频率

首先,单独修改输入函数的频率,然后单独修改振幅, 对比修改之前的结果。

频率设为 3:

基本未变化,只有尖尖的位置:(0.03, -50)

频率设为 5:

基本也未变化,只有尖尖的位置:(0.05, -50)

现在可以确定,输入函数的频率是多少,尖尖就会出现在 频率/100 这个位置。频率为 4,那就是 0.04 秒。也就是说,从第一个位置开始,频率从 0 开始一一对应。如果频率是 0,是一种什么情况呢?频率 0 意味着输入函数输出的是一个常数,一个不会变化的数。

再次修改输入函数,让他始终输出 1:

虚数部分

奇怪的是,第一个位置并没有出现尖尖。相反,这次实数部分出现了:

实数部分

看来虚数部分不能反应出输入函数的常数特征,相反,实数部分能够得到反映。现在得补充一下之前的结论:频率为 0 时,第一个位置的虚数不会有影响,相反第一个实数的值会有变化。

振幅

接下来,单独修改振幅对比结果。

振幅设为 0.5:

基本未变化,只有尖尖的高度:(0.04, -25)

振幅设为 2:

基本未变化,只有尖尖的高度:(0.04, -100)

显然,振幅只影响尖尖的高度,位置没有影响。目前的实验来看,高度为振幅的 50 倍,这个 50 恰好是采样次数的一半。为了验证这一点,接下来再单独修改采样次数,也就是采样频率,对比结果。

采样频率

还是使用最初的那个输入函数,不做如何修改。因为这次样本数量非常多,所以图表里只记录了前一部分。

采样 200 次:

基本未变化,只有尖尖的高度:(0.04, -100)

采样 400 次:

基本未变化,只有尖尖的高度:(0.04, -200)

所以可以下结论:虚数部分的尖尖的高度与对应频率信号的振幅的关系即:\(振幅 \times 采样次数的一半 = 高度\)

结论

有了上面所有的实验结果,可以发现,DFT 可以提取一个信号的频率、振幅信息,而且每个位置对应一个频率,绝对值对应该频率的正弦曲线的振幅。

然而单纯就一个正弦曲线来说(正弦曲线在实际中的一个典型表现就是「波」),一般会用三个特性去描述他:频率、振幅和相位。

还是拿最开始的那个输入函数来看,但 t = 0 时,弧度即为 0, 那如果弧度不是 0 呢?简单的,可以给他加上一个常量:

\[f(t)=sin(4t \times 2\pi + \pi/2)\]

这样,时间为 0 时,弧度从 \(\pi/2\) 开始。所以说,起始弧度就是「相位」。

现在我们知道 DFT 给了我们一个正弦曲线的频率和振幅信息,那相位哪去了呢?