可恶啊,又让我想起这玩意儿了,然后又忘记怎么推了,只能回去查一查了。其实我困扰的是,CNN的卷积为啥叫卷积啊,卷积不是吗,那个卷积核,分明可以直接对应元素相乘吗,网上有些图居然也是直接乘了。我以为卷积也有什么快速算法啊,可是,普通的卷积又不具有傅里叶变换的性质。
其实傅里叶变换也忘得差不多了,不过,就先这样推着吧。
参考:The Fast Fourier Transform and its Applications
Notation
性质
总而言之,有个周期
FFT推导
让我们来看看这第一次分解所消耗的计算量,计算一个消耗r次,那么计算全部的就是次。在全部知道的情况下,计算一个是次,计算所有的需要次,故本次分解消耗次计算。如果不分解的话,大概是级别。
更加恐怖的是这是第一次分解,可以看到,成了周期为的傅里叶变换,成了
周期为的傅里叶变换,所以,如果,还能够分解的话,计算量还能进一步减少。
假设:,那么,计算总共消耗次,这时变为周期为的傅里叶变换,可以预见,将分解为,计算应当消耗次,但第二次需要计算个这种情况(因为总共有N个,分成了r个周期为s的傅里叶变换),所以第二次消耗的计算量也为,以此类推,最后结果为:
又
论文里的推导过程
代码
import numpy as np
import time
from scipy.fftpack import fft, ifft
def number_fc(N): #因式分解 比较蠢的方法 我在数论里好像看到过更好的就这样吧
for i in range(2, int(np.sqrt(N)) + 1):
if N % i == 0:
s = i
r = N / i
return int(s), int(r)
return 0, 0
def conv(x, k): #普通的运算
N = len(x)
w = np.array([np.exp(-2 * i * k * np.pi * 1j / N) for i in range(N)]) #论文中是+的 不是-的 但是scipy库里的是-的所以我这里也取-
return x @ w
def w_s_N(j0, j1, s, N): #求 怎么说呢 看懂式子就懂这个了
w_s_j1 = np.array([np.exp(-2 * n0 * j1 * np.pi * 1j / s) for n0 in range(s)])
w_N_n0 = np.array([np.exp(-2 * n0 * j0 * np.pi * 1j / N) for n0 in range(s)])
return w_s_j1 * w_N_n0
def FFT(x): #FFT
N = len(x)
s, r = number_fc(N)
if not s: #当s或r=0也就是N不能再分解了求直接返回傅里叶变换后的
return np.array([conv(x, k) for k in range(N)])
else:
A0 = np.zeros(N, dtype=complex) #相当于X_j1_j0
A1 = np.array([FFT(x[n0::s]) for n0 in range(s)]) #计算A1_j0_n0 输出的是一个矩阵,s*r的,每个元素是一个A1_J0_N0
for j1 in range(s):
for j0 in range(r):
A0[j1 * r + j0] = A1[:, j0] @ w_s_N(j0, j1, s, N)
return A0
测试代码:
N = 4096
x = np.arange(N)
t1 = time.time()
y1 = FFT(x)
t2 = time.time()
print(t2 - t1) #0.5311074256896973
t1 = time.time()
y2 = [conv(x, k) for k in range(N)]
t2 = time.time()
print(t2-t1) #26.705339193344116
t1 = time.time()
y3 = fft(x)
t2 = time.time()
print(t2 - t1) #0.0
从上面可以看到,当N = 4096的时候,二者的差距已经十分明显了。第三个方案是,scipy库里面的,即便N都这么大了,依然动不了他分毫。大概是我写代码的水平还是太low了吧。
在写代码的时候,对于时间复杂的计算也有了新的认识,不是那么想当然的,果然实践才是检验真理的唯一标准。