新手学FFT

# 小白学 FFT 来源：https://entropicthoughts.com/fft fft.jpg 我从来没做过任何时频变换。我理解大致的概念，但从未深入接触过细节。今天我终于有了一个使用它的理由，但首先得学习基础知识。 ## 纯正弦波（https://entropicthoughts.com/fft#pure-sine-wave） 下面是一个频率为 5 Hz 的正弦波，时长三秒，采样率为 100 Hz。¹¹嗯，它是合成生成的，但这就是它所代表的含义。 fft-tutorial-01.svg 如果我们用 R 中的`fft`函数处理这个信号，会得到一堆复数——数量与样本数一样多。本例中有 300 个。这些数值的大小取决于信号的振幅，因此对于这个信号，它们应该最大为 1。然而，从 R 的`fft`函数输出后，这些值还被乘以样本数进行了缩放。在下面的图中，我除以了样本数以恢复振幅。 fft-tutorial-02.svg 每个点对应信号在特定频率上的振幅。图中 300 个点中有 298 个位于原点，因为数据中不存在这些频率的信号。 远离原点的两个点让人觉得我们的数据有*两个*频率，但我们明明只有一个！我得承认这是傅里叶变换的一个怪癖，我并未完全理解，但显然这个变换把我们的振幅为 1 的 5 Hz 信号看作两个信号： - 一个振幅为 0.5 的 5 Hz 信号，以及 - 一个振幅为 −0.5 的 −5 Hz 信号。 它们加在一起就是振幅为 1 的 5 Hz 信号。大概是这么回事。 总之，在最新的这张图中，X 轴（`fft`输出的实部）对应该频点上的余弦信号，而 Y 轴（虚部）对应该频点上的正弦信号。由于我们只有正弦波，所以所有频率的实部都为零。 然而，这样看有点误导。每个频率输出是复数的原因在于，正弦加余弦可以表示任意相位的波。因此，通过复数的大小和角度来观察它们更为直观。角度是最无趣的部分。我们先从它开始。 如前所述，`fft`的输出是与原始数据等长的向量。该向量中的每个值对应一个频率，第一个值是“零频”，即数据的均值。之后的值对应的频率线性排列，最后一个频率刚好低于采样率。在下面的图中，我把 X 轴从向量索引转换为了相应的频率。图中最高频率几乎等于每秒 100 的采样率。 fft-tutorial-03.svg 相位从 \\(-\pi\\) 到 \\(\pi\\)，在我们的例子中它到处都是。但这并不奇怪！在大多数频率上我们没有任何信号，所以不关心那些地方的相位。我们关心的是每个频率上信号的振幅，而那正是我们期望的形状。这里我们绘制`fft`输出复数的大小，而不是它们的角度。 fft-tutorial-04.svg 这两根柱子对应 5 Hz 和 −5 Hz。我们可以通过把镜像的振幅相加并只绘制频谱的前半部分来简化此图。毕竟，奈奎斯特告诉我们，我们无法恢复频率高于采样率一半的信号。 fft-tutorial-05.svg 成了！我们恢复了原始 5 Hz 纯正弦波信号的全部振幅。原来已经有这么多细节是我之前不知道的。 ## 截断信号（https://entropicthoughts.com/fft#truncated-signal） 在上述例子中，我们展示了一个不自然地干净的信号。数据的构造甚至保证了采样窗口正好覆盖了整数个周期。如果我们提前截断采样，就会出现频谱泄漏。 下面是同一信号经过`fft`后的原始复数值，只不过我们不是在 3 秒时停止记录，而是在 2.9 秒时。 fft-tutorial-06.svg 这看起来和之前的结果完全不同。从各个频率的振幅来看，会更容易理解发生了什么。 fft-tutorial-07.svg 截断会导致频谱泄漏，把干净的 5 Hz 正弦波涂抹成一组紧密相关的信号。就实际而言，这没什么大不了的：我们仍然可以完美重建原始信号。泄漏只是让频谱看起来令人困惑，实际上并不会破坏任何数据。 处理截断的一种常见方法是对数据施加一个包络函数，使其在端点附近平滑地过渡到零。下面的图显示的是与之前完全相同的信号，只是在两端有所衰减。 fft-tutorial-08.svg 这个包络给频谱增加了一些模糊性，但确实有效地消除了频谱泄漏。 ## 含噪的复合信号（https://entropicthoughts.com/fft#a-noisy--composite-signal） 如果我们给信号加一点噪声，`fft`仍然能从中得到一个有代表性的频谱。 fft-tutorial-09.svg 我们可以看到 5 Hz 附近的峰值，但其他所有频率上都有噪声底。在下一组图中，我向信号中添加了另一个正弦波。在这个噪声水平下很难分辨，但借助傅里叶变换的魔力，我相信你能找出那个频率。 fft-tutorial-10.svg 没错，5 Hz 的波还在，现在又多了一个 23 Hz 的波。 ## 奇特的信号（https://entropicthoughts.com/fft#odd-signals） 下面是一些奇特的信号类型，我们可能想知道它们在频域中长什么样。首先是脉冲。 fft-tutorial-11.svg 虽然这里不太明显，但脉冲抬高了频谱的底噪——它对所有频率都有贡献。 相比之下，阶跃函数对低频的贡献呈衰减趋势。 fft-tutorial-12.svg 方波贡献的是衰减的奇次谐波。 fft-tutorial-13.svg 自回归/衰减脉冲看起来与阶跃函数相似，但更平滑。 fft-tutorial-14.svg 不过玩得差不多了。有一些真实数据让我感兴趣。 ## 计算真实增长（https://entropicthoughts.com/fft#computing-real-growth） 我不能说这些数据到底意味着什么，但幸好我准备了一个老套的掩护故事：在我孩子的卧室外面，他们正在建造一栋楼。我们假装我有一个测量设备，每小时记录一次建筑高度的增长速度，而且已经记录了整整一年。 由于建筑工人主要在白天的工作日施工，我预期会出现日周期（建筑高度逐日增长）和周周期（周末增长停滞）。在下面的图中，频率谱被改成了周期谱。因此，频谱中 7 处的峰值对应的是一个周周期，而不是每天 7 次的频率。 fft-tutorial-15.svg 频谱中最高的柱子对应周期为 365 天的循环。显然，仅凭一年的数据，我们无法知道这是否是一个重复出现的循环；这根柱子实际上只代表年末总体上的增长趋势，而这很可能是真的：建筑增长速度似乎提高了。 除了我们预期的一天和七天周期之外，还有一个显著的周期为 3.5 天的循环。这是因为工作日模式实际上是一个占空比为 5 开 2 关的方波。我们预期该方波突出的谐波是周期为 3.5 天、2.4 天、1.4 天和 1.2 天的循环。频谱中确实有对应所有这些周期的峰值！ 虽然在这张小图上不可见，频谱中还有一个周期为两小时的循环。看看是否存在小时周期会很有趣，但奈奎斯特极限告诉我们，从小时数据中能提取的最短周期是两小时。 但是，如果我们把目光从频谱移回到时间序列上，会看到数据中有一些尖锐的峰值。这些通常不是真实的增长。继续用我们老套的例子，我们可以假装这些是有人撞到放置测量设备的桌子造成的伪信号，导致了错误的读数。我想看一整年的真实增长，所以我需要以某种方式*减去这些碰撞*带来的影响。 回想我们之前的实验，我们可能会意识到这些碰撞就是脉冲！也许我们可以假设它们对所有频率的贡献相等。然后我们应该可以从所有频率中减去一个恒定幅度，从而消除最严重的桌子碰撞影响。我本来没指望它能奏效，但它确实有点用。 fft-tutorial-16.svg 看上面滤波后的时间序列，在滤除最严重的伪信号之后，周周期变得明显多了。看起来越来越像是一周中有某一天他们进行了大部分实际建筑施工。²²我注意到我的掩护故事开始站不住脚了。请您暂且放下怀疑。 我们可以从周期谱看出，滤波操作的实际效果是衰减了大量高频噪声，即周期短于一天的噪声。结果有些类似于滚动中位数，但基于`fft`的方法能更好地保留逐日变化。滚动中位数要达到足够的碰撞滤波强度，其窗口大小就必须足够大，而大到那个程度时也会抹平掉大部分日周期。 如果我们将这些峰值视为高频分量而不是单独的脉冲，我们可能会尝试用低通滤波器将它们滤除。这也会把频谱中的短周期分量归零，但对长周期分量没有任何影响。我也试过这种方法，但效果不如前者。 我怀疑我们不能把这些峰值视为高频分量的原因是它们分布得有些随机。它们并不构成一个统一的高频信号。也许如果我们对时间序列的某个放大局部做频谱，就可以把它们视为高频分量了。让我们把这个想法稍微形式化一下。 ## 制作频谱图（https://entropicthoughts.com/fft#making-a-spectrogram） 在上述分析中，我们对一整年的小时增长数据生成了一个频谱。我们可以设想把这一年切成更小的片段，然后对每个片段计算频谱。这样做的好处是，如果频谱在一年中发生变化，我们可以为每个片段得到定制的频谱，而不是一个过于笼统而无法使用的东西。 选择多大的片段并不显然。较小的片段让我们能对频谱的较小变化做出反应，但另一方面，片段的长度也限制了我们能检测到的循环周期长度。以下是一些备选方案： - 我们可以立即排除日片段，因为信号中比日循环更快的基本只有噪声。 - 我们可以尝试周片段，因为这样可以同时发现日循环和周循环。然而，尝试过周片段后，实际中`fft`每个片段得到的数据太少，根本提取不出任何有意义的循环。 - 双周片段有与周片段类似的问题。 - 月片段有点噪，但效果还可以。 - 季度片段也有效，但前三个季度完全无趣。 我们从中得到的叫做*频谱图*。下面使用的是月片段，但无论片段大小如何，把频谱图绘制在时间序列下方的妙处在于，我们可以立即看出它们之间的关系。 fft-tutorial-17.svg 频谱图中较暗的列对应时间序列中活动较少的区域。明亮的列有更多活动。这一点一目了然！ 再仔细看，在一年中的前 200 天左右，周循环实际上并不可见。日循环很微弱，但几乎在所有月份中都能看到。在倒数第三个月，周循环似乎被双周循环取代了！周循环的泛音在最后两个月非常突出。 像这样制作频谱图的缺点是，如果有用的信息恰好跨越了月份的边界，那么任一个月的频谱都可能捕捉不到它。我们可以通过在时间序列上滑动窗口来解决这个问题。下图仍然使用月窗口，只不过不是逐月跳跃，而是每次滑动 12 小时。 fft-tutorial-18.svg 我还不完全确定我理解出现的所有模式，但看着它思考很有趣！ 但现在，也许我们可以把峰值视为高频分量了？我们可以取每个窗口的频谱，通过低通滤波器处理，然后从这些连续的频谱重建原始时间序列。想想就令人兴奋！不过这写代码听起来有点繁琐，我在这篇文章上已经花了太长时间，而且还有一个想法我想先探索一下。 这标志着我目前阶段`fft`实验的结束。我确实对它熟悉了很多，现在认为它是我的工具箱中以前没有的工具。感谢一路相伴！