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

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

©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 216,470评论 6 501
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 92,393评论 3 392
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 162,577评论 0 353
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 58,176评论 1 292
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 67,189评论 6 388
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 51,155评论 1 299
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 40,041评论 3 418
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 38,903评论 0 274
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 45,319评论 1 310
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 37,539评论 2 332
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 39,703评论 1 348
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 35,417评论 5 343
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 41,013评论 3 325
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 31,664评论 0 22
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 32,818评论 1 269
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 47,711评论 2 368
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 44,601评论 2 353

推荐阅读更多精彩内容