昨天写了两篇关乎傅立叶变换的小短文。

第一篇以9副图展示了傅立叶变换把信号分解成若干个正弦波的能力，并且实现了两个For循环版本的DFT,最后分析了FFT蝶形算法的实现思路以及计算量。

第一篇文章地址：https://www.toutiao.com/a7047374338911224357

第二篇主要展示了卷积操作与傅立叶变换的等效性，都具有滤波能力。在某些场合，比如相关滤波做跟踪，我们就用FFT方法实现滤波。在深度神经网络中，我们直接用点集来实现滤波。其实，无论哪种滤波方法，都是可以相互转换的。

第二篇文章地址：https://www.toutiao.com/a7047478476533694987

今天我们就把蝶形算法的代码实现出来。

代码下载链接：https://gitee.com/anjiang2020_admin/fft/blob/master/fft_zmm.py

推导

先看一下DFT的公式：

为了简化写法：

要求解信号f(n)中含有频率为k的正选信号的幅值，使用公式（1）写代码：

F(k)=0
for n in range(0:N):
       F(k)+=f(n)*W(N,n,k)

可以看到，理论上非常简单的。我们再对k取0，1，2，，，，N值，即可得到F(0),F(1),....F(N)

由于W(N,n,k)与W(N,n,K+N/2)的关系，当n为欧时，两者相等。当n为奇数时，两者相反。

又由于：F(k) 的关于 W(N,n,k)的for循环，而F(k+N/2)是关于W(N,n,k+N/2)的for 循环，所以，我们用n的奇偶行将将F(k)的实现拆为两部分：

拆成n为偶数，n为奇数两部分后，关键就是求出每个部分的值了，我么将F(k)表示为：

F(k)=E(k)+O(k),其中E(k）为偶数部分，O(k)为奇数部分。

E(k):

E(k)=0
for n in 0:N/2:
    E(k)+=f(2*n)*W(N,2*n,k)

O(k):

O(k)=0
for n in 0:N/2:
    O(k)+=f(2*n+1)*W(N,2*n+1,k)

只要分别求出E(k)，O(k)的值，F(k)=E(k)+O(k),F(k+N/2)=E(k)-O(k).

其实这里写代码也就可以了，写for循环就能求出相应的值。

但是，这个因为还是有for循环，所以虽然加速了，但是仍然慢，我们要想办法加速，看能不能把for循环消灭掉。

如何求E(k)的值？

将f(n)中偶数部分单独拿出来作为一个信号：fe(m), 则E(k)的公式可改写为：

E(k):

E(k)=0
for m in 0:N/2:
    E(k)+=fe(m)*W(N,2*m,k)

也可以写为：

E(k)=0
for m in 0:N/2:
    E(k)+=fe(m)*W(N/2,m,k)

我们再接着把N/2表示为：M，则也可以写为：

E(k)=0
for m in 0:M:
    E(k)+=f1(m)*W(M,m,k)

注意到这个写法，与我们求F(k)的写法一样：

F(k)=0
for n in range(0:N):
       F(k)+=f(n)*W(N,n,k)

观察E(k)与F(k)，发现结构一摸一样。

回想一下，我们求F(k),k=1,2,3,4,....N时，只用求前半部分F(1),F(2),,,,,F(N/2)即可，F(N/2+1)...F(N）的值，都可以由前半部分给出。

对于F(k),将其分解成了F(k) = E(k) + O(k),k=1,2,3,4,...N/2

这列E(k)的写法与F(k)一摸一样了，F(k)的性质， E(k）同样具有。

也就是要求E(k),k=1,2,3,...N/2，只用求E(1),E(2),...E(N/4)即可，从E(N/4+1),..到E(N/2)都可以用前半部分来求。

对E(k)做奇偶分解：E(k)=E_E(k)+E_O(k)

总的来说，就是将F(k)做奇偶分解得到的E和O ,其中E 又可以做奇偶分解。

如果O能够做奇偶分解的话，就完美了，此时，我们只用实现一个奇偶分解的函数，然后递归调用就能计算出F(k)的值了。

如何求O(k)的值？

观察O:

O(k):

O(k)=0
for n in 0:N/2:
    O(k)+=f(2*n+1)*W(N,2*n+1,k)

f(n)中的奇数部分表示为fo(m), 然后令：M=N/2，得：

O(k):

O(k)=0
for m in 0:M:
    O(k)+=fo(m)*W(2M,2*m+1,k)

而W(2M,2m+1,k) = W(2M,2m,k)*W(2M,1,k)

又因为：W(N,2n,k) = W(N/2,n,k)，

所以上式W(2M,2m,k)*W(2M,1,k)=W(2M/2,m,k)*W(2M,1,k)=W(M,m,k)*W(2M,1,k)

?：

O(k):

O(k)=0
for m in 0:M:
    O(k)+=fo(m)*W(M,m,k)*W(2M,1,k) 

注意到，在for循环中，m是变化的，所以W(2M,1,k)是个常量，我们可以将其提取出来，for循环结束后再乘，即：

O(k):

O(k)=0
for m in 0:M:
    O(k)+=fo(m)*W(M,m,k)
O(k)=O(k)*W(2M,1,k)

我们只看前三行，是不是形式与F(k)求取过程又一样了。这样形式是可以进行奇偶分解的：O(k) = [OE(k) + OO(k) ]* W(2M,1,k)

我们总结一下：

F(k)=E(k)+O(k)

E(k)=EE(k) + EO(k)

O(k)=W(2M,1,k)*[ OE(k) + OO(k) ]

也就是说，F(k)进行奇偶分解后得到的E和O,又都可以进行奇偶分解。

显然，EE(k)也是可以进行奇偶分解的，EO(k),OE(k),OO(k)都是可以进行奇偶分解的。

直到M=1时便不能再分解了。

所以，我们的伪代码为：

设输入信号f(n),n=1,2,3,4...N 要求其FFT为F(k),k=1,2,3,...N

则：

def getF[ f(n) ]
     E,O=奇偶分解[f(n)]
     F[k]=E+O,k=1,2,3,4....N/2
     F[k]=E-O,k=N/2+1,N/2+2,...N.
     return F
def 奇偶分解[f(n)]:  
     E = getF( f[0:2:n] ),偶数部分可以用得到F的方法得到，因为他们公式相同， 输入的信号变为f[0:2:n]了，只取偶数部分
    O = getF( f[1:2:n] ),偶数部分可以用得到F的方法得到，因为他们公式相同输入的信号变为f[0:2:n]了，只取偶数部分
    O=O*W(2M,1,k)
    retunr E,O

可以看到，我们写了两个函数： getF 和奇偶分解，他们组成递归调用，就可以递归的进行奇偶分解，递归调用的过程中，每进入一次奇偶分解函数，输入到getF 函数的序列长度就减半，当这个序列的长度为2时，就可以不再递归了：

N=len(f(n))

如果 N=2:

那么 E= f(0)*W(N,0,0)

O=f(0)*W(N,0,0)*W(2*M,1,0)

此时M=N/2=1

我们把这段写进奇偶分解函数里：

def 奇偶分解[f(n)]:
if N==2:
    E= f(0)*W(N,0,0)   
    M=N/2
    O=f(0)*W(N,0,0)*W(2*M,1,0) 
  
else:
    E = getF( f[0:2:n] ),偶数部分可以用得到F的方法得到，因为他们公式相同， 输入的信号变为f[0:2:n]了，只取偶数部分
    O = getF( f[1:2:n] ),偶数部分可以用得到F的方法得到，因为他们公式相同输入的信号变为f[0:2:n]了，只取偶数部分
    O=O*W(2M,1,k)
    retunr E,O

到这里，完成了连个函数getF,和奇偶分解函数，利用这两个函数递归的调用，完成了对奇偶分解后的值递归的进行奇偶分解的作用。

当 N==2时，就不需要再用for循环了，所以，就没有必要再调用getF函数了，直接写出E 和 O值即可。

python实现

奇偶分解函数的实现：

注意到，这里fft_add 及时getF函数

MultiWnk是一个函数，实现了复数W(N,n,k) 与复数 W(N,1,k)乘法的功能。

getF函数的实现：

使用这两个函数完成FFT的蝶形算法，并与scipy的fft进行了对比，代码为：

#coding:utf-8
import numpy as np
import matplotlib.pyplot as plt
from scipy.fftpack import fft,ifft

# input : ys 
# output: F


def Wnk(N,n,k):
   N=float(N)
   n=float(n)
   k=float(k)
   return np.array([np.cos(2*np.pi/N*n*k),np.sin(2*np.pi/N*n*k)]) 
#(a+bi)(c+di)=(ac-bd)+(bc+ad)i
def MultiWnk(wnk1,wnk2):
    a=wnk1[0]
    b=wnk1[1]
    c=wnk2[0]
    d=wnk2[1]
    return np.array([a*c-b*d,b*c+a*d])

def EO(xs):
   N=len(xs)
   if N==2:
       E = xs[0]*Wnk(N,0,0)
       #O = xs[1]*Wnk(N,0,0)*Wnk(N,1,0)
       O = xs[1]*Wnk(N,0+1,0)
       #O = xs[1]*MultiWnk(Wnk(N,0,0),Wnk(N,1,0))
       #O=[O[k,:]*Wnk(N,1,k) for k in range(0,N/2)]
       #O = O*Wnk(N,1,0)
   else:
       E=fft_add(xs[range(0,N,2)])
       O=fft_add(xs[range(1,N,2)])
       #O=[O[k,:]*Wnk(N,1,k) for k in range(0,N/2)]
       # 使用复数乘法
       O=[MultiWnk(O[k,:],Wnk(N,1,k)) for k in range(0,N/2)]
   return E,O
def EO1(xs):
    N = len(xs)
    E=np.zeros((N/2,2))
    O=np.zeros((N/2,2))
    for k in range(0,N/2):
        for i in range(0,N/2):
            E[k,:]+=xs[2*i]*Wnk(N,2*i,k) 
            #O[k,:]+=xs[2*i+1]*Wnk(N,2*i+1,k)
            O[k,:]+=xs[2*i+1]*MultiWnk(Wnk(N,2*i,k),Wnk(N,1,k))
    return E,O

def fft_add(xs):
   N = len(xs)
   F=np.zeros((N,2))
   E,O=EO(xs)
   F[range(N/2),:]=E+O;
   F[range(N/2,N),:]=E-O;
   print("F:%s"%(F))
   return F

def absF(F):
   Fabs=np.zeros((F.shape[0]))
   for i in range(len(F)):
       Fabs[i]=(F[i,0]**2+F[i,1]**2)**0.5
   return Fabs

if __name__=="__main__":
   print("hello")
   xs=np.linspace(0,1,32)
   print(xs)
   ys=5.3*np.sin(2*np.pi*5*xs)
   print(ys)
   plt.subplot(221);plt.plot(np.arange(len(xs)),ys);plt.title("ys")

   F=fft_add(ys)
   print(absF(F)) 
   plt.subplot(222);plt.plot(np.arange(len(xs)),absF(F)/len(xs));plt.title("ys's fft by zmm")
   plt.subplot(224);plt.plot(np.arange(len(xs)),abs(fft(ys))/len(xs));plt.title("ys' fft by scipy")
   plt.show()
   

可以看到，ys's fft by zmm是文本的FFT 实现的结果。ys'fft by scipy是scipy库中的fft,结果一致，说明代码理论上是 正确的。

代码中的采样频率为32，信号本身频率为5，可以看到，在5的位置有峰值出现。

读者可以更改信号本身频率，以及采样频率，接着去做实验。

不过，目前代码只支持ys的长度刚好为2的几次方，不是2的几次方，奇偶分解就进行不下去了。

解决这个的改进思路：如果遇到长度不是2的几次方的ys，就把ys用线性插值法，把长度插值成2的几次方，再进行蝶形FFT就可以了.

如果，不使用蝶形FFT,用for循环实现DFT, 就不用担心ys的长度是否刚好为2的几次方。

至此，从本科时期就想做的事情：用蝶形算法实现DFT的加速：FFT。

今天终于完成了