代码链接: link.
DFT在讨论FFT之前，我们需要先了解以下DFT。所谓的DFT其实就是两个矩阵做点乘。
多项式可以有两种表示方法，一种是系数表示法，另一种是点值表示法。
这两种表示法之间是可以转换的，系数表示法到点值表示法非常简单，就是随便取几个点带入求值，而点值表示法到系数表示法就需要用到插值法，关于插值法又不在本文的讨论范围之内了。
所谓的点值表示法，设A(x)是一个n阶的多项式，那么至少可以用n+1对(x0,A(x0)),(x1,A(x1))…这样的点对来表示，也就是说确定了这n+1个点，就可以唯一的确定一个多项式，形如a0+a1x+a2x²+…这种。证明的方法利用到了幼儿园级别的线性代数知识。
当两个多项式相乘时，我们可以很自然的想到将对应的点值相乘，这样就得到了结果的点值表示法。
假设现在有两个n阶的多项式相乘，那么利用点值表示法，我们需要将对应的n个点分别相乘，也就是要做n*n次乘法运算，这个时间复杂度是O(n²),与我们直接将系数表示下的多项式相乘复杂度一致。
这已经是离散傅里叶变化的思想了。
所谓离散傅里叶变换就是
其本质就是M矩阵和X向量做点乘。
FFT我们想要找到突破口就需要在将多项式转化为点值表示这一步做一点小操作。这一点操作就是利用对称性，减少计算量。
我们注意到一个n阶多项式A(x) = a0+a1x²+a2x³+…其奇数项是奇函数，偶数项是偶函数。而奇数项又可以分离成x乘一个偶函数。
因此，我们利用偶函数的性质可以减少很多计算量。
比方现在有一个多项式x+x²+x³，那么我们应该至少需要4对点才可以求得该多项式的点值表示法，我可以找x1,x2,-x1,-x2这样的正负对，所以只需要找2对即可，因为另一半的值就无非就是符号不同。
但是只是这样的一次对称只能减少一半的计算量，只能让原本O(n²)的算法变为O(n²/2)，其提升不是很大。如果可以一直利用这种对称性就好了，这种对称性建立的前提是正负对的存在，但是显然在实数域我们无法做到这一点，所以必须要在复数域讨论。
所以Fourier提出，对于一个n-1阶的多项式，我们不妨直接用n个复数w0,w1,w2…wn代替原来的x0,x1,x2…xn，这样得到一种特殊的点值表示法就是被成为离散傅里叶变化(DFT)。而这n个复数不是随便选取的，而是在复平面将单位圆等分n分后对应的点。
单位根的性质复平面上的单位根有这样一个性质
于是正负对出现了。
我们的程序就可以按照这个框架来写了。
#DFT，本质就是两个矩阵相乘
def DFT(x):
x = np.asarray(x,dtype=float)
N = x.shape[0]
M = [[j for i in range(N)]for j in range(N)]
M = np.asarray(M,dtype=complex)
w = np.exp(-2j*np.pi/N)
for i in range(N):
for j in range(N):
M[i][j] = np.power(w,i*j)
return np.dot(M,x)
FFT
递归
#递归FFT，利用分治思想的dft
def fft_recurrence(x):
x = np.asarray(x,dtype=float)
N = x.shape[0]
if N<2:
return DFT(x)
x_even = fft_recurrence(x[0::2])
x_odd = fft_recurrence(x[1::2])
factor = np.exp(-2j * np.pi * np.arange(N) / N)
return np.concatenate([x_even + factor[:int(N / 2)] * x_odd,
x_even + factor[int(N / 2):] * x_odd])
非递归
上面那种FFT是最基础的，同时因为递归的原因执行速度很慢，所以提出了非递归的优化方法。
观察 n =4时候的子序列位置变换情况
我们发现子序列最终的位置其实就是初始位置的二进制翻转后结果，如3(11)变为了3(11)，而1(01)变为了2(10)。
Java实现 Complex复数类首先创建一个复数类，定义一系列的加减乘除操作。
package fft;
public class Complex {
private double realPart;
private double imaginaryPart;
public Complex() {
this.realPart =0;
this.imaginaryPart =0;
}
public Complex(double realPart,double imaginaryPart) {
this.realPart = realPart;
this.imaginaryPart = imaginaryPart;
}
public Complex add(Complex w) {
if (w ==null) {
return new Complex();
}
return new Complex(this.realPart+w.getRealPart(),this.imaginaryPart+w.getImaginaryPart());
}
public Complex subt(Complex w) {
if (w ==null) {
return new Complex();
}
return new Complex(this.realPart-w.getRealPart(),this.imaginaryPart-w.getImaginaryPart());
}
public Complex mult(Complex w) {
if (w ==null) {
return new Complex();
}
return new Complex(this.realPart*w.getRealPart()-this.imaginaryPart*w.getImaginaryPart(),
this.realPart*w.getImaginaryPart()+this.imaginaryPart*w.getRealPart());
}
public Complex divide(Complex w) {
if (w ==null) {
return new Complex();
}
double W = w.getImaginaryPart()*w.getImaginaryPart() - w.getRealPart()*w.getRealPart();
return new Complex((this.realPart*w.getRealPart()+this.imaginaryPart*w.getImaginaryPart())/W,
(this.imaginaryPart*w.getRealPart()-this.realPart*w.getImaginaryPart())/W);
}
@Override
public String toString() {
// TODO Auto-generated method stub
String s = (this.imaginaryPart>=0)? "+":"";
return this.realPart+s+this.imaginaryPart+"i";
}
public double getRealPart() {
return realPart;
}
public void setRealPart(double realPart) {
this.realPart = realPart;
}
public double getImaginaryPart() {
return imaginaryPart;
}
public void setImaginaryPart(double imaginaryPart) {
this.imaginaryPart = imaginaryPart;
}
}
Matrix复数矩阵类
Matrix类是一个复数矩阵类，定义了一系列我们要用到的矩阵操作
package fft;
public class Matrix {
public static Complex[][] subMatrix(Complex[][] a,int x1,int x2){
int n = a.length;
int m = a[0].length;
int l = x2-x1;
Complex[][] b = new Complex[l][m];
for(int i =0;i=x1&&j=x1&&k=y1&&p=na)? b[i-na]:a[i];
return c;
}
public static void show(Complex[][] a) {
for (int i =0;i
FFT类
主要定义FFT递归版，以及DFT操作，因为我们要求的n必须是2的幂，所以如果给定的向量不满足，则会调用genArray()生成一个满足条件的向量将给定向量末尾补上0。
package fft;
public class FFT {
Complex omega(int n,int k) {
return new Complex((Math.cos(2*Math.PI*k/n)),(-Math.sin(2*Math.PI*k/n)));
}
int genN(int n) {
int s =2;
while(s
测试类
package fft;
public class Main {
public static void main(String[] args) {
int n = 1024;
FFT fft = new FFT();
//随机生成一个列向量，该向量的每一个元素都是随机的复数
Complex[][] x = new Complex[n][1];
for (int i =0;i
测试结果
python
运行代码
if __name__ == '__main__':
#随机生成 维度为1024的列向量
x = np.random.random(1024)
print('x矩阵为:')
print(x)
print('fft结果为：')
print(fft_recurrence(x))
print('')
#打印不同方法的fft耗时对比
start_1 = time.perf_counter()
fft_recurrence(x)
end_1 = time.perf_counter()
print("fft_recurrence cost：",(end_1-start_1)*1000,'毫秒')
start_2 = time.perf_counter()
np.fft.fft(x)
end_2 = time.perf_counter()
print('numpy.fft.fft() cost:', (end_2 - start_2)*1000,'毫秒')
start_3 = time.perf_counter()
DFT(x)
end_3 = time.perf_counter()
print('dft cost:', (end_3 - start_3)*1000,'毫秒')
# 用numpy.fft的结果和我的fft对比，返回True则两者相同
result = np.allclose(fft_recurrence(x),np.fft.fft(x))
print('my fft numpy.fft计算出的结果相同？',result)
OUTPUT
我们调用了numpy.fft用于验证我们的FFT是否正确，可以看到python实现的FFT得到的结果是正确的。
现在我们对代码进行一点修改来测验java的FFT是否正确。
替换了x向量为
x = []
for i in range(16):
x.append(i)
打印结果为
Java
同样生成16维的向量，[0,1,2…15]
int n = 16;
FFT fft = new FFT();
//随机生成一个列向量，该向量的每一个元素都是随机的复数
Complex[][] x = new Complex[n][1];
for (int i =0;i
OUTPUT
这里显示0毫秒是因为取的时间是long，长整型，而花费时间小于1毫秒所以显示为0
对比
我们不妨让python版和Java版同时计算一个大向量来看看结果。
现在我们将n变为1024，
java耗时
python耗时
可以看到同样的递归版本，Java耗时比python要少了一半多。这还是在没有进行优化的情况下。
