本文章将介绍如何用python一行代码实现基二时间抽选FFT函数的定义。

在我们进入正题之前,先来热个身,用python实现一行快速排序,这个是相对轻松的,列表推导式是一个很方便的东西,因此我们只需要:

quick_sort = lambda x :quick_sort([i for i in x if i<x[0]])+[i for i in x if i==x[0]]+quick_sort([i for i in x if i>x[0]]) if len(x) else []

非常简单一个递归就可以做到。好了,那我们现在来尝试一下用一行代码实现FFT算法好了。

虽说是一行,但是出于代码的简洁,我们先将π与e两个数学极限定义(楼主标题党!),但是毕竟只是常数赋值,并且只是为了代码简洁,因此还请各位谅解。 

import math
e = math.exp(1)
pi = math.pi

 首先我们遇到的第一个问题,基二时间FFT算法要求信号长度必须为2的n次幂,不然的话我们就得在信号的后面补零,直到信号长度变为大于它原有长度的2的最小n次幂。那么如何只用一行实现这个功能呢?

如果我们要判断一个正整数是否为2的n次方,首先我们也许会想到将它反复除以2,直到有余数出现或者它等于1。但是这显然并不是一行代码可以做到的。这时候我们可以想到python的位运算,众所周知数字在python里其实是通过二进制保存的,因此可以利用位运算判断一个数是否为2的n次幂。这里直接给出结论,如果x&(x-1)为0,那么它一定是2的n次幂(值得注意的是x=0时也满足,因此要单独判断)。

那么接下来,我们应该如何判断大于一个数的最小2的n次幂呢?bin函数可以把一个数字转化为它的二进制数字的字符串(格式为“0bxxxxx”),接下来只需要判断这个字符串的长度就可以知道到底是多少啦,这里直接上写好的结果。

y = lambda x:x + [0] * (2 ** (len(bin(len(x))) - 2) - len(x)) if (len(x) & (len(x) - 1)) and (len(x) > 2) else x

这样我们就得到了把x自动补0的函数 ,但是这和FFT算法还差了十万八千里,因此我们接下来就要考虑当信号长度符合要求时,如何递归实现FFT。

FFT = lambda x :FFT(x + [0] * (2 ** (len(bin(len(x))) - 2) - len(x))) if (len(x) & (len(x) - 1)) and (len(x) > 2) else

首先是这样,判断如果信号长度大于2且长度不符合要求,那自动补0以后递归调用FFT。那么在这个else后面,就可以直接写我们对信号的处理了。

FFT = lambda x :FFT(x + [0] * (2 ** (len(bin(len(x))) - 2) - len(x))) if (len(x) & (len(x) - 1)) and (len(x) > 2) else [x1 + x2 * (len(x) / 2 - j - 0.5)/abs(len(x) / 2 - j - 0.5) for j ,[x1 ,x2] in enumerate([[x1 ,x2 * (e ** ((-1j) * 2 * pi * i / len(x)))] for i ,(x1 ,x2) in enumerate(zip(FFT([x[i] for i in range(len(x)) if not(i%2)]),FFT([x[i] for i in range(len(x)) if i%2])))] * 2)] if len(x)>2 else

接下来这一长段代码也许有点难理解,但是也是最核心的部分。

我们单独分析之前第一个else后面的部分,我们知道,基二时间抽选算法是把长信号的离散傅里叶变换转换为短信号的离散傅里叶变换,并且递归直到信号长度为2的算法。

因此我们判断当信号长度为2时,我们将信号的偶数序列和奇数序列分割,然后再次递归调用,并且为了防止每一层递归运算是重复计算下一层导致时间复杂度指数增长,我们只将结果计算一次以后复制需要使用的次数(也就是那里的*2的作用),并且结合enumerate与zip函数打包成列表,以便使用我们之前用过的列表推导式对下一层返回的奇偶序列傅里叶变换结果做出相应的处理。提一嘴中间那个(len(x) / 2 - j - 0.5)/abs(len(x) / 2 - j - 0.5)的作用,它负责将前半与后半部分分开,当i索引为序列的前一半时,它的值为+1,反之为-1,用这样的方式来代替if语句实现前后半数列的分别处理。(如果分别处理并且用+号连接,那么将会不得不计算两次下一层FFT导致时间复杂度提升非常多)

FFT = lambda x :FFT(x + [0] * (2 ** (len(bin(len(x))) - 2) - len(x))) if (len(x) & (len(x) - 1)) and (len(x) > 2) else [x1 + x2 * (len(x) / 2 - j - 0.5)/abs(len(x) / 2 - j - 0.5) for j ,[x1 ,x2] in enumerate([[x1 ,x2 * (e ** ((-1j) * 2 * pi * i / len(x)))] for i ,(x1 ,x2) in enumerate(zip(FFT([x[i] for i in range(len(x)) if not(i%2)]),FFT([x[i] for i in range(len(x)) if i%2])))] * 2)] if len(x)>2 else [x[0] + x[1],x[0] - x[1]] if len(x)==2 else x

最后就是大功告成的一步,当我们递归一层层深入,直到信号长度为2时,我们终于可以直接返回它的值!另外出于代码的健壮性,当x的信号长度为1时,我们也返回它的DFT结果——它本身。这样我们就完成了使用一行代码实现的FFT算法.!

使用matalb来进行验证:

X = [1,4,4,1,4,1,1,4]
print(FFT(X))

——————————————————————————

结果为:
[(20+0j), (1.2426406871192857-3j), 0j, (-7.242640687119286+2.9999999999999996j), 0j, (-7.242640687119286-3j), 0j, (1.2426406871192857+3.0000000000000004j)]

Python fft变换之后 python中fft函数_后端

 完美perfect!