时间序列的傅里叶变换:使用numpy解释快速卷积
时间序列的傅里叶变换:使用numpy解释快速卷积
从头实现 vs 使用numpy
傅里叶变换算法被认为是数学领域中最重要的发现之一。法国数学家让-巴普蒂斯特·约瑟夫·傅里叶在他的书《热分析理论》中于1822年奠定了谐波分析的基础。如今,傅里叶变换及其各种变体构成了我们现代世界的基础,驱动着诸如压缩、通信、图像处理等技术。
![法国数学家让-巴普蒂斯特·约瑟夫·傅里叶(Jean Baptiste Joseph Fourier)的雕刻肖像(1768-1830年),19世纪早期。[来源:维基百科,图片来自公共领域]](https://miro.medium.com/v2/resize:fit:640/format:webp/1*U5yW8KvFsE_9uIIDI8UD5A.png)
这个奇妙的框架还为分析时间序列提供了很好的工具……这也是我们在这里的原因!
本文是傅里叶变换系列文章的一部分。今天我们将讨论卷积以及傅里叶变换提供的最快的实现方式。
所有的图形和方程都由作者制作。
离散傅里叶变换(DFT)的定义
让我们从基本定义开始。对于一个具有N个元素的离散时间序列x,其离散傅里叶变换为:

其中k表示x频谱的第k个频率。注意,一些作者在该定义中添加了一个缩放因子1/N,但对于本文来说并不重要——总的来说,这只是一个定义的问题和坚持使用的问题。
反傅里叶变换则是(根据正向傅里叶变换的定义):

话虽如此,傅里叶变换的一个最重要的定理是,一个空间中的卷积等价于另一个空间中的乘法。换句话说,两个信号的乘积的傅里叶变换是它们各自傅里叶谱的卷积,而两个信号的卷积的傅里叶变换是它们各自傅里叶谱的乘积。

和

其中点表示标准乘积(乘法),圆圈中的星号表示循环卷积。
两个重要的注意事项:
- 周期信号:傅里叶分析框架认为我们处理的信号是周期性的。换句话说,它们从负无穷重复到正无穷。然而,用有限内存的计算机处理这样的信号并不总是实际可行的,所以我们只“玩弄”一个周期,正如我们将在下文中看到的。
- 循环卷积:卷积定理规定乘法等效于循环卷积,这与我们更熟悉的线性卷积有点不同。正如我们将看到的,它并没有那么不同,也没有那么复杂。
循环卷积与线性卷积
如果你熟悉线性卷积,通常简称为“卷积”,你就不会对循环卷积感到困惑。基本上,循环卷积只是卷积周期信号的方式。你可以猜到,线性卷积只对有限长度的信号有意义,这些信号不会从负无穷延伸到正无穷。在我们的情况下,在傅里叶分析的背景下,我们的信号是周期性的,因此不满足这个条件。我们无法讨论(线性)卷积。
然而,我们仍然可以在我们的周期信号上直观地进行类似于线性卷积的操作:只需将周期信号卷积一次周期长度。这就是循环卷积所做的:它将两个长度相同的周期信号沿一个周期跨度进行卷积。
为了进一步说明区别,比较离散线性卷积和离散循环卷积的公式:


注意以下差异:
– 边界:线性卷积使用由负无穷到正无穷的样本,正如前面所述,在这种情况下,x和y具有有限的能量,求和是有意义的。对于循环卷积,我们只需要在一个周期范围内发生的情况,因此求和只在一个周期内进行。
– 循环索引:在循环卷积中,我们使用长度为N的取模运算对y的索引进行“包装”。这只是一种确保y被视为具有周期N的方法:当我们想要知道位置k处的y的值时,我们只需使用位置k%N处的y的值 – 由于y是N周期性的,我们得到正确的值。同样,这只是处理周期性无限长度样本序列的数学方法。
在numpy中的实现
Numpy为有限长度信号提供了很好的工具:这是一个好消息,因为正如我们所看到的,我们的无限长度周期信号可以用一个周期表示。
让我们创建一个简单的类来表示这样的信号。我们添加了一个快速绘制数组的方法,以及“基本”数组之前和之后的额外周期,以记住我们正在处理周期序列。
让我们看两个例子:首先是采样的正弦序列,然后是线性序列。这两个序列都被认为是N周期性的(在本例中为N = 10)。

循环卷积,慢速方式
现在让我们实现上述循环卷积方程。使用索引和模运算符,这是非常直观的:

这太棒了,我们现在可以看到两个信号之间的循环卷积是什么样子的。把所有内容放在一个图中:

现在这个解决方案的效果非常好,但有一个主要缺点:它很慢。正如你所看到的,我们必须通过2个嵌套循环来计算结果:一个用于结果数组中的每个位置,另一个用于计算该位置的结果:我们称该算法为O(N²),随着N的增加,操作的数量将随N的平方增加。
对于像示例中那样的小数组,这不是问题,但是当数组增长时,它将成为一个主要问题。
此外,对数值数据进行循环大多数情况下被认为是Python中的一种不良实践。一定会有更好的方法…
循环卷积,傅里叶的方法
这就是福里叶变换和卷积定理发挥作用的地方。由于离散傅里叶变换的实现方式,使用快速傅里叶变换(FFT)以非常快速和优化的方式实现,所以该操作非常快速(我们说FFT的时间复杂度是O(N log N),比O(N²)好得多)。
使用卷积定理,我们可以利用两个序列的离散傅里叶变换的乘积,当通过逆离散傅里叶变换再次转换回时域时,得到输入时间序列的卷积。换句话说,我们有:

其中DFT表示离散傅里叶变换,IDFT表示逆操作。
然后我们可以使用numpy轻松地实现此算法来计算x和y的卷积:
数值和时间比较
最后,让我们验证这两种方法产生的结果是否相同,并比较使用两种技术计算循环卷积所需的时间:

这是一个完美的匹配!两者在数值上严格等价。
现在来看时间比较:
结果是:
- 对于N=10个样本,DFT快了6倍
- 对于N=1000个样本,DFT快了大约10000倍
这太大了!现在想象一下当你分析成千上万个样本的时间序列时,它能带给你什么!
总结
我们在本文中看到,傅里叶变换是一个强大的工具,特别是由于允许我们以非常高效的方式计算卷积的卷积定理。我们已经看到线性和循环并不完全相同,但都基于卷积。
订阅以便将来直接获取有关傅里叶变换的帖子到您的订阅源中!
此外,请查看我的其他帖子,如果您喜欢其中任何一个,请订阅,这对我实现100个订阅者的目标非常有帮助:
使用NumPy将有限差分方法的分辨率提高300倍 | 作者:Yoann Mocquin | Towards Data Science (VoAGI.com)
PCA/LDA/ICA:成分分析算法比较 | 作者:Yoann Mocquin | Towards Data Science (VoAGI.com)
使用容器方法封装NumPy的数组 | 作者:Yoann Mocquin | Towards Data Science (VoAGI.com)
深入研究Seaborn:调色板 | 作者:Yoann Mocquin | Analytics Vidhya | VoAGI