由于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()
