FFT

可恶啊，又让我想起这玩意儿了，然后又忘记怎么推了，只能回去查一查了。其实我困扰的是，CNN的卷积为啥叫卷积啊，卷积不是 $h(t)=\mathop{\sum}\limits_kf(k)g(t-k)$ 吗，那个卷积核，分明可以直接对应元素相乘吗，网上有些图居然也是直接乘了。我以为卷积也有什么快速算法啊，可是，普通的卷积又不具有傅里叶变换的性质。

其实傅里叶变换也忘得差不多了，不过，就先这样推着吧。
参考：The Fast Fourier Transform and its Applications

在这里插入图片描述

Notation

$W_N = exp\{\frac{2\pi i}{N}\}$
$X(j)=\mathop{\sum}\limits_{n=0}^{N-1}A(n)W_N^{jn}$
$A(n)=\frac{1}{N}\mathop{\sum}\limits_{n=0}^{N-1}X(j)W_N^{-jn}$

性质

$W_N^N=1, \quad W_N^{j+N}=W_N^j$
$W_N^j=W_N^{j \:mod\:N}$
$X(j)=X(j\:mod\:N)$
$A(n) = A(n \: mod\: N)$
总而言之，有个周期 $N$
$\mathop{\sum}\limits_{n=0}^{N-1}W_N^{nj}W_N^{mj}=N \quad if \: n=m\: mod\:N否则0$

FFT推导

$假设N= r\times s，那么存在j_1,j_0，对于任意的j可用下列式子表示:$
$j = j_1r+j_0 \quad j_1=0,1,\ldots,s-1,\quad j_0=0,1,\ldots,r-1$
$同样，对于n来说，存在n_1, n_0:$
$n = n_1s + n_0 \quad n_1=0,1,\ldots,r-1,\quad n_0=0,1,\ldots,s-1$
$W_N^{jn}=W_N^{j_1n_1rs}W_N^{j_1rn_0}W_N^{j_0n_1s}W_N^{j_0n_0}=W_s^{j_1n_0}W_r^{j_0n_1}W_N^{j_0n_0}$
$X(j)=X(j_1,j_0)=\mathop{\sum}\limits_{n_0=0}^{s-1}\mathop{\sum}\limits_{n_1=0}^{r-1}A(n_1,n_0)W_r^{j_0n_1}W_N^{j_0n_0}W_s^{j_1n_0}$
$令：A_1(j_0, n_0)=W_N^{j_0n_0}\mathop{\sum}\limits_{n_1=0}^{r-1}A(n_1,n_0)W_r^{j_0n_1}$
$X(j_1,j_0)=\mathop{\sum}\limits_{n_0=0}^{s-1}A_1(n_1,n_0)W_s^{j_1n_0}$

让我们来看看这第一次分解所消耗的计算量，计算一个 $A_1(j_0,n_0)$ 消耗r次，那么计算全部的 $A_1$ 就是 $rsr=Nr$ 次。在 $A_1$ 全部知道的情况下，计算一个 $X(j_1,j_0)$ 是 $s$ 次，计算所有的 $X$ 需要 $rss=Ns$ 次，故本次分解消耗 $N(r+s)$ 次计算。如果不分解的话，大概是 $N^2$ 级别。
更加恐怖的是这是第一次分解，可以看到， $X(j_1,j_0)$ 成了周期为 $s$ 的傅里叶变换， $A_1(j_0,n_0)$ 成了
周期为 $r$ 的傅里叶变换，所以，如果， $r,s$ 还能够分解的话，计算量还能进一步减少。
假设： $N=r^m,s=r^{m-1}$ ,那么，计算 $A_1$ 总共消耗 $Nr=r^{m+1}$ 次，这时 $X$ 变为周期为 $s=r^{m-1}$ 的傅里叶变换，可以预见，将 $s$ 分解为 $r^{m-2}\times r$ ，计算 $A_2$ 应当消耗 $r^{m}$ 次，但第二次需要计算 $r$ 个这种情况（因为 $X(j_1,j_0)$ 总共有N个，分成了r个周期为s的傅里叶变换），所以第二次消耗的计算量也为 $r^{m+1}$ ,以此类推，最后结果为:
$m \times r^{m+1}=mNr$
又 $N=r^m,所以m=log_rN，多么牛啊，原本二次的计算量近似线性了！$
$一般情况下，N=r_1\times r_2\times \cdots r_m,最后的计算量为N(r_1+r_2+\ldots +r_m)$

论文里的推导过程

在这里插入图片描述

代码

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了吧。

在写代码的时候，对于时间复杂的计算也有了新的认识，不是那么想当然的，果然实践才是检验真理的唯一标准。

Notation

性质

FFT推导

论文里的推导过程

代码

推荐阅读更多精彩内容