傅立叶变换。 速度与激情

通常，在开发算法时，我们会遇到计算复杂性的限制，这似乎是无法克服的。 傅立叶变换具有复杂性 $$$O（n ^ 2）$ 和一个快速版本，由众议院¹在1805年左右提出（并于1965年由James Cooley和John Tukey发明） $$$O（nlog（n））$ 。 在本文中，我想向您展示可以在线性时间内获得转换结果 $$$O（n）$ 甚至达到持续的困难 $$$O（1）$ 在实际问题中发现的某些条件下。


当我面对编写实时分析声音系统传递函数的程序的任务时，我和其他所有人一样，首先转向了快速转换。 一切都很好，但是随着时间窗口的增大，CPU负载变得很大，必须执行某些操作。 决定暂停并再次研究转换，同时寻找解决问题的方法。 让我们回到约瑟夫·傅里叶²的原始转换：

$$

$$f（x）= \和\极限_ {-\ infty} ^ {+ \ infty} c_ke ^ {2 \ pi ikx / T} \\ c_k = \ frac {1} {T} \ int \ limit_0 ^ Tf（ x）e ^ {-2 \ pi ikx / T} dx$$

我们将仔细研究这里的情况。 频域中的每个输出值 $$$c_k$ 是信号所有输入值的总和 $$$f（x）$ 乘以 $$$e ^ {-2 \ pi ikx / T}$ 。 为了进行计算，我们需要遍历每个输出值的所有输入数据，即 满足那些 $$$n ^ 2$ 操作。

摆脱n

让我提醒您，最初的任务是实时分析声音数据。 为此，将大小为N的所选时间窗口（基本上是缓冲区）填充有频率f _d对应于采样频率的数据。 在周期T中，输入数据从时间窗口转换为频率窗口。 如果查看实数，则N在2 ¹⁴ （16 384）到2 ¹⁶ （65 536）个样本之间变化（这些值是从FFT继承的，其中窗口大小必须是2的幂）。 时间T = 80ms（每秒12.5次），这使您可以非常方便地查看更改，并且不会使CPU和GPU过载。 采样频率f _d为标准频率，为48 kHz。 让我们计算一下时间窗口中各维度之间有多少数据变化。 在时间T期间，它进入缓冲区 $$$48000 * 0.08 = $ 384$ 样品。 因此，窗口中只有5％到23％的数据被更新。 在最坏的情况下，尽管事实上在先前的迭代中已经对它们进行了处理，但仍有95％（最好是73％，这也很多！）将一次又一次地进入转换。

那时，细心的读者会举手说：“等等，但是系数呢？ $$$e ^ {-2 \ pi ikx / T}$ ？ 毕竟，对于每个新的变换，相同的数据将位于序列的新位置，因此系数不同？” 对于每五个需要照顾的人，让我们回忆一下通常被遗忘的转换中的一个重要细节。 在函数值的研究中 $$$f（t）$ 在从0到t的时间间隔内，该函数被认为是周期性的，这使您可以轻松地在时间上向左或向右移动函数。 结果，我们有权不在末尾插入新值并从头开始删除旧值，而是循环替换缓冲区中的数据。

为了清楚起见，您可以以表格形式编写缓冲区的变化方式：

t = 0	f（0）	f（1）	f（2）	f（3）	f（4）	f（5）	（（6）	f（7）	（（8）	f（9）
t = 1	（（10）	f（1）	f（2）	f（3）	f（4）	f（5）	（（6）	f（7）	（（8）	f（9）
t = 2	（（10）	（（11）	f（2）	f（3）	f（4）	f（5）	（（6）	f（7）	（（8）	f（9）
t = 3	（（10）	（（11）	（（12）	f（3）	f（4）	f（5）	（（6）	f（7）	（（8）	f（9）
t = 4	（（10）	（（11）	（（12）	（（13）	f（4）	f（5）	（（6）	f（7）	（（8）	f（9）

您可以写出时间变换如何从t ₁变为t ₂ ：

$$

$$F_t = F_ {t-1} + \ Delta F \\ \ Delta F：\ Delta c_k = \ frac {1} {T} \ int \ limits_ {t_1} ^ {t_2}（f_t（x）-f_ { t-1}（x））e ^ {-2 \ pi ikx / T} dx$$

价值 $$$F_ {t-1}（x）$ 是上一次转换的结果，以及计算的复杂性 $$$\ Delta f（x）$ 不依赖于时间窗口的大小，因此是恒定的。 结果，转换的复杂度将是 $$$O（n）$ ^*因为 对于我们来说，剩下的一切就是一次通过频率窗口，并对随时间变化的T个样本应用更改。 我还想提请您注意以下事实： $$$e ^ {-2 \ pi ikx / T}$ 可以预先进行计算，从而提高了生产率，并且循环中仅剩两个操作：将实数相减并将实数乘以一个复数，实际上，这两个操作都很简单且便宜。

为了完整显示图片，它仅表示初始状态，但是这里的一切都很简单：

$$

$$F_0（x）= 0$$

*-当然，整个转换的最终复杂性将保持不变 $$$O（n ^ 2）$ ，但是它将在更新缓冲区时逐步执行n次迭代。 $$$O（n）$ -这是更新数据的复杂性，但这正是我们所需要的（使用FFT时，每次转换的复杂性 $$$O（nlog（n））$ ）

但是，如果您深入研究该怎么办。 或摆脱第二个

我想立即保留一个保留意见，即只有在您不打算对结果执行逆变换（以校正信号或获得脉冲响应）时，才适用后续步骤。 首先，我想提醒您，转换后，我们得到了一个对称的数据数组，该数组立即使我们可以将转换次数减少一半。

现在，在给出问题条件的情况下，让我们分析结果数据集。 我们有一组复数，每个复数描述特定频率下的振荡幅度和相位。 频率可以通过以下公式确定： $$$f [j] = j \ frac {fd} {N}$ 为 $$$j <\ frac {N} {2}$ 。 让我们根据数据评估频率窗口的阶跃： $$$\ Delta f = \ frac {fd} {N}$ 对于N = 2 ¹⁴ ：2.93 Hz（对于2 = ¹⁶ ：0.73 Hz）。 因此，在1 kHz至2 kHz的范围内，我们可获得341个结果。 尝试独立评估N = 65536时，从8 kHz到16 kHz范围内将有多少数据。 好多！ 我们需要这么多数据吗？ 当然，在显示声音系统的频率特性的问题上，答案是否定的。 另一方面，对于低频区域的分析，小步长非常有用。 别忘了还有一些时间表可以解决这些问题（ $$$\ frac {N} {2}$ ）转换为易于理解的形式（取平均值，样条曲线或平滑度）并将其显示在屏幕上。 在高频率下，即使具有4K屏幕并以对数频率轴在全屏模式下显示图形，步长也会很快小于1个像素。

根据经验，您发现每个八度音阶只有48个点就足够了，并且让数据更平滑，更平均，我建议停在96。在20 Hz至20 kHz的音频频率范围内，仅计算10个八度音阶很容易：20、40、80 ，160、320、640、1280、2560、5120、10240、20480，它们中的每一个都可以划分为给定数量的子带（不要忘记应该以几何方式而不是算术方式进行分区），因此，仅对获得960个频率 Performan是比原来的版本更小的16 ... 65倍。

因此，结合两种方法，我们可以获得数据更新算法的恒定复杂度 $$$O（1）$ 。

平方蜂蜜和一匙焦油

现在我们可以放心地说 $$$O（n ^ 2）$ 我们变得复杂 $$$O（1）$ 使用两个简单的生活技巧：

• 在分析了问题之后，我们注意到数据是逐渐添加的，并且时间窗口的完全更新的周期比变换的周期要高得多，然后继续进行傅立叶变换的差的计算。
• 从频率窗口中的算术步骤变为仅受指定值限制，这可以大大减少转换次数。

但是，当然，生活不是童话，而是童话。 这两种方法的应用使我们能够真正卸载CPU，以便猜测它可以计算傅里叶变换并在屏幕上显示结果。 $$$N = 2 ^ {20}$ 很难。 但是，当现实中的信号不是周期性的时（这种情况对于获得正确的转换结果是必要的）并且无法选择合适的窗口大小，惩罚并没有让自己等待时间，因此有必要使用各种窗口功能，这不再允许您完全使用第一步。 实践表明，使用窗函数对于研究频率小于2的信号至关重要。 $$$0,1f_d$ 。 在高频率下，落入时间窗口的周期数会显着衰减原始函数中由于存在一阶间隙（介于f（0）和f（N-1）之间）而引起的失真。

最后，我也拒绝了第二步，回到了FFT，因为 这个任务的收益已经很小了。

总结

如果您的数据具有明显的周期性并且需要使用较大的时间窗口对时间进行分析，则可以采用第一种方法，我记得这不必是2级，即 任何自然数。
如果仅在数据中分析特定的一小部分频率集，则第二种方法适用（即使考虑到窗口功能）。

las，对于我这个问题，这仍然只是一点数学上的娱乐，但我希望它能激发您在假期中根据输入数据随时间的变化研究其他算法:)

文学作品

1. Cooley – Tukey FFT算法
2. E.A. Vlasova Rows。 MSTU出版社以鲍曼（N.E. Bauman）的名字命名。 莫斯科 2002年

图片由涩谷道雄（Michio Shibuya）漫画改编而成，“令人兴奋的数学。 傅立叶分析