快速傅里叶变换(FFT)算法【详解】[通俗易懂]

快速傅里叶变换(FFT)算法【详解】[通俗易懂]快速傅里叶变换(FastFourierTransform)是信号处理与数据分析领域里最重要的算法之一。我打开一本老旧的算法书,欣赏了JWCooley和JohnTukey在1965年的文章

大家好,又见面了,我是你们的朋友全栈君。

快速傅里叶变换(Fast Fourier Transform)是信号处理与数据分析领域里最重要的算法之一。我打开一本老旧的算法书,欣赏了JW Cooley 和 John Tukey 在1965年的文章中,以看似简单的计算技巧来讲解这个东西。

本文的目标是,深入Cooley-Tukey  FFT 算法,解释作为其根源的“对称性”,并以一些直观的python代码将其理论转变为实际。我希望这次研究能对这个算法的背景原理有更全面的认识。

FFT(快速傅里叶变换)本身就是离散傅里叶变换(Discrete Fourier Transform)的快速算法,使算法复杂度由原本的O(N^2) 变为 O(NlogN),离散傅里叶变换DFT,如同更为人熟悉的连续傅里叶变换,有如下的正、逆定义形式:

快速傅里叶变换(FFT)算法【详解】[通俗易懂]

xn 到 Xk 的转化就是空域到频域的转换,这个转换有助于研究信号的功率谱,和使某些问题的计算更有效率。事实上,你还可以查看一下我们即将推出的天文学和统计学的图书的第十章(这里有一些图示和python代码)。作为一个例子,你可以查看下我的文章《用python求解薛定谔方程》,是如何利用FFT将原本复杂的微分方程简化。

正因为FFT在那么多领域里如此有用,python提供了很多标准工具和封装来计算它。NumPy 和 SciPy 都有经过充分测试的封装好的FFT库,分别位于子模块 numpy.fft 和 scipy.fftpack 。我所知的最快的FFT是在 FFTW包中 ,而你也可以在python的pyFFTW 包中使用它。

虽然说了这么远,但还是暂时先将这些库放一边,考虑一下怎样使用原始的python从头开始计算FFT。

计算离散傅里叶变换

简单起见,我们只关心正变换,因为逆变换也只是以很相似的方式就能做到。看一下上面的DFT表达式,它只是一个直观的线性运算:向量x的矩阵乘法,

快速傅里叶变换(FFT)算法【详解】[通俗易懂]

矩阵M可以表示为

快速傅里叶变换(FFT)算法【详解】[通俗易懂]

这么想的话,我们可以简单地利用矩阵乘法计算DFT:

1 import numpy as np
2 def DFT_slow(x):
3     """Compute the discrete Fourier Transform of the 1D array x"""
4     x = np.asarray(x, dtype=float)
5     N = x.shape[0]
6     n = np.arange(N)
7     k = n.reshape((N, 1))
8     M = np.exp(-2j * np.pi * k * n / N)
9     return np.dot(M, x)

对比numpy的内置FFT函数,我们来对结果进行仔细检查

x = np.random.random(1024)
np.allclose(DFT_slow(x), np.fft.fft(x))

 输出:

True

 现在为了验证我们的算法有多慢,对比下两者的执行时间

%timeit DFT_slow(x)
 
%timeit np.fft.fft(x)

 输出:

10 loops, best of 3: 75.4 ms per loop
10000 loops, best of 3: 25.5 µs per loop

使用这种简化的实现方法,正如所料,我们慢了一千多倍。但问题不是这么简单。对于长度为N的输入矢量,FFT是O(N logN)级的,而我们的慢算法是O(N^2)级的。这就意味着,FFT用50毫秒能干完的活,对于我们的慢算法来说,要差不多20小时! 那么FFT是怎么提速完事的呢?答案就在于他利用了对称性。

离散傅里叶变换中的对称性

算法设计者所掌握的最重要手段之一,就是利用问题的对称性。如果你能清晰地展示问题的某一部分与另一部分相关,那么你就只需计算子结果一次,从而节省了计算成本。

Cooley 和 Tukey 正是使用这种方法导出FFT。 首先我们来看下快速傅里叶变换(FFT)算法【详解】[通俗易懂]的值。根据上面的表达式,推出:

快速傅里叶变换(FFT)算法【详解】[通俗易懂]

对于所有的整数n,exp[2π i n]=1。

最后一行展示了DFT很好的对称性:

快速傅里叶变换(FFT)算法【详解】[通俗易懂]

简单地拓展一下:

快速傅里叶变换(FFT)算法【详解】[通俗易懂]

同理对于所有整数 i 。正如下面即将看到的,这个对称性能被利用于更快地计算DFT。

DFT 到 FFT:

利用对称性 Cooley 和 Tukey 证明了,DFT的计算可以分为两部分。从DFT的定义得:

快速傅里叶变换(FFT)算法【详解】[通俗易懂]

我们将单个DFT分成了看起来相似两个更小的DFT。一个奇,一个偶。目前为止,我们还没有节省计算开销,每一部分都包含(N/2)∗N的计算量,总的来说,就是N^2 。

技巧就是对每一部分利用对称性。因为 k 的范围是 0≤k<N , 而 n 的范围是 0≤n<M≡N/2 , 从上面的对称性特点来看,我们只需对每个子问题作一半的计算量。O(N^2) 变成了 O(M^2) 。

但我们不是到这步就停下来,只要我们的小傅里叶变换是偶倍数,就可以再作分治,直到分解出来的子问题小到无法通过分治提高效率,接近极限时,这个递归是 O(n logn) 级的。

这个递归算法能在python里快速实现,当子问题被分解到合适大小时,再用回原本那种“慢方法”。

 1 def FFT(x):
 2     """A recursive implementation of the 1D Cooley-Tukey FFT"""
 3     x = np.asarray(x, dtype=float)
 4     N = x.shape[0]
 5  
 6     if N % 2 > 0:
 7         raise ValueError("size of x must be a power of 2")
 8     elif N <= 32:  # this cutoff should be optimized
 9         return DFT_slow(x)
10     else:
11         X_even = FFT(x[::2])
12         X_odd = FFT(x[1::2])
13         factor = np.exp(-2j * np.pi * np.arange(N) / N)
14         return np.concatenate([X_even + factor[:N / 2] * X_odd,
15                                X_even + factor[N / 2:] * X_odd])

现在我们做个快速的检查,看结果是否正确:

x = np.random.random(1024)
np.allclose(FFT(x), np.fft.fft(x))

 输出:

True

 然后与“慢方法”的运行时间对比下:

%timeit DFT_slow(x)
 
%timeit FFT(x)
 
%timeit np.fft.fft(x)

输出:

10 loops, best of 3: 77.6 ms per loop
100 loops, best of 3: 4.07 ms per loop
10000 loops, best of 3: 24.7 µs per loop

现在的算法比之前的快了一个数量级。而且,我们的递归算法渐近于 O(n logn) 。我们实现了FFT 。

需要注意的是,我们还没做到numpy的内置FFT算法,这是意料之中的。numpy 的 fft 背后的FFTPACK算法 是以 Fortran 实现的,经过了多年的调优。此外,我们的NumPy的解决方案,同时涉及的Python堆栈递归和许多临时数组的分配,这显著地增加了计算时间。

还想加快速度的话,一个好的方法是使用Python/ NumPy的工作时,尽可能将重复计算向量化。我们是可以做到的,在计算过程中消除递归,使我们的python FFT更有效率。

向量化的NumPy

注意上面的递归FFT实现,在最底层的递归,我们做了N/32次的矩阵向量乘积。我们的算法会得益于将这些矩阵向量乘积化为一次性计算的矩阵-矩阵乘积。在每一层的递归,重复的计算也可以被向量化。因为NumPy很擅长这类操作,我们可以利用这一点来实现向量化的FFT

 1 def FFT_vectorized(x):
 2     """A vectorized, non-recursive version of the Cooley-Tukey FFT"""
 3     x = np.asarray(x, dtype=float)
 4     N = x.shape[0]
 5  
 6     if np.log2(N) % 1 > 0:
 7         raise ValueError("size of x must be a power of 2")
 8  
 9     # N_min here is equivalent to the stopping condition above,
10     # and should be a power of 2
11     N_min = min(N, 32)
12  
13     # Perform an O[N^2] DFT on all length-N_min sub-problems at once
14     n = np.arange(N_min)
15     k = n[:, None]
16     M = np.exp(-2j * np.pi * n * k / N_min)
17     X = np.dot(M, x.reshape((N_min, -1)))
18  
19     # build-up each level of the recursive calculation all at once
20     while X.shape[0] < N:
21         X_even = X[:, :X.shape[1] / 2]
22         X_odd = X[:, X.shape[1] / 2:]
23         factor = np.exp(-1j * np.pi * np.arange(X.shape[0])
24                         / X.shape[0])[:, None]
25         X = np.vstack([X_even + factor * X_odd,
26                        X_even - factor * X_odd])
27  
28     return X.ravel()

 

x = np.random.random(1024)
np.allclose(FFT_vectorized(x), np.fft.fft(x))

 输出:

True

因为我们的算法效率更大幅地提升了,所以来做个更大的测试(不包括DFT_slow)

 

x = np.random.random(1024 * 16)
%timeit FFT(x)
 
%timeit FFT_vectorized(x)
 
%timeit np.fft.fft(x)

 输出:

10 loops, best of 3: 72.8 ms per loop
100 loops, best of 3: 4.11 ms per loop
1000 loops, best of 3: 505 µs per loop

我们的实现又提升了一个级别。这里我们是以 FFTPACK中大约10以内的因数基准,用了仅仅几十行 Python + NumPy代码。虽然没有相应的计算来证明, Python版本是远优于 FFTPACK源,这个你可以从这里浏览到。

那么 FFTPACK是怎么获得这个最后一点的加速的呢?也许它只是一个详细的记录簿, FFTPACK花了大量时间来保证任何的子计算能够被复用。我们这里的numpy版本涉及到额外的内存的分配和复制,对于如Fortran的一些低级语言就能够很容易的控制和最小化内存的使用。并且Cooley-Tukey算法还能够使其分成超过两部分(正如我们这里用到的Cooley-Tukey FFT基2算法),而且,其它更为先进的FFT算法或许也可以能够得到应用,包括基于卷积的从根本上不同的方法(例如Bluestein的算法和Rader的算法)。结合以上的思路延伸和方法,就可使阵列大小即使不满足2的幂,FFT也能快速执行。

版权声明:本文内容由互联网用户自发贡献,该文观点仅代表作者本人。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如发现本站有涉嫌侵权/违法违规的内容, 请发送邮件至 举报,一经查实,本站将立刻删除。

发布者:全栈程序员-用户IM,转载请注明出处:https://javaforall.cn/155732.html原文链接:https://javaforall.cn

【正版授权,激活自己账号】: Jetbrains全家桶Ide使用,1年售后保障,每天仅需1毛

【官方授权 正版激活】: 官方授权 正版激活 支持Jetbrains家族下所有IDE 使用个人JB账号...

(0)
blank

相关推荐

  • Shiro框架:缓存、session会话、自定义FormAuthenticationFilter、RemenberMe

    Shiro框架:缓存、session会话、自定义FormAuthenticationFilter、RemenberMe

  • python中的pop函数和append函数

    python中的pop函数和append函数pop()函数1、描述pop()函数用于移除列表中的一个元素(默认最后一个元素),并且返回该元素的值。语法pop()方法语法:list.pop(obj=list[-1])2、参数obj–可选参数,要移除列表元素的对象。3、返回值该方法返回从列表中移除的元素对象。4、实例以下实例展示了pop()函数的使用方法:#!/usr/bin/pythonaList=[123,’xyz’,’z

  • 图标变成了一张白纸_电脑图标是白纸形式但能打开

    图标变成了一张白纸_电脑图标是白纸形式但能打开大家好,我是波导终结者。WIN10到现在也有些年头了,虽然好的不学学坏的,天天搞强制升级有些烦人,有时候新版本的BUG也比较致命,但是整体的性能,功能和稳定性上还是有飞跃性的提升。其实WIN10自带不少实用工具,而且兼容性和稳定性没得说,很多时候我们可以不必寻找第三方工具,省了很多钱,也提高了不少效率。今天就一起来看看吧。WIN10自带截图工具WIN10自带截图工具,虽然肯定不比专业的截图工具强,…

    2022年10月19日
  • java 论坛_5 个最好用的 Java 开源论坛系统

    java 论坛_5 个最好用的 Java 开源论坛系统大家好!我是Guide哥,Java后端开发。一个会一点前端,喜欢烹饪的自由少年。最近有点小忙。但是,由于前几天答应了一位读者自己会推荐一些开源的论坛系统,所以,昨晚就简单地熬了个夜,对比了很多个开源论坛系统之后,总结成了这篇文章。这篇文章我一共推荐了5个论坛类开源项目,除了有1个是基于PHP开发之外,其他都是基于Java,并且大部分都是基于SpringBoot这个主流框…

  • RestFul风格「建议收藏」

    RestFul风格「建议收藏」RestFul风格概念Restful就是一个资源定位及资源操作的风格。不是标准也不是协议,只是一种风格。基于这个风格设计的软件可以更简洁,更有层次,更易于实现缓存等机制。功能资源:互联网所有的事物都可以被抽象为资源资源操作:使用POST、DELETE、PUT、GET,使用不同方法对资源进行操作。分别对应添加、删除、修改、查询。传统方式操作资源:通过不同的参数来实现不同的效果!方法单一,post和get​ http://127.0.0.1/item/queryItem.actio

  • 2019版idea激活码(破解版激活)[通俗易懂]

    2019版idea激活码(破解版激活),https://javaforall.cn/100143.html。详细ieda激活码不妨到全栈程序员必看教程网一起来了解一下吧!

发表回复

您的电子邮箱地址不会被公开。

关注全栈程序员社区公众号