FFT算法实现，python，Java

FFT算法实践报告

FFT基本原理

代码链接: 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分后对应的点。

单位根的性质复平面上的单位根有这样一个性质

于是正负对出现了。
我们的程序就可以按照这个框架来写了。

python实现

DFT

#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;
	}
	
	/**
	 * 加法运算
	 * @param w
	 * @return
	 */
	public Complex add(Complex w) {
   
		if (w ==null) {
   
			return new Complex();
		}
		return new Complex(this.realPart+w.getRealPart(),this.imaginaryPart+w.getImaginaryPart());
	}
	
	
	/**
	 * 减法运算
	 * @param w
	 * @return
	 */
	public Complex subt(Complex w) {
   
		if (w ==null) {
   
			return new Complex();
		}
		return new Complex(this.realPart-w.getRealPart(),this.imaginaryPart-w.getImaginaryPart());
	}
	
	 /**
	  *乘法 
	  * @param w
	  * @return
	  */
	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());
	}


	/**
	 * 除法
	 * @param w
	 * @return
	 */
	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(