使用python实现基-2FFT、基-4FFT快速傅里叶变换算法

FFT算法

快速傅立叶变换（FFT）是信号处理和数据分析中最重要的算法之一，很多人只是调用现成的库如FFTW，但为了知其所以然，加深对算法的理解，我们有必要搞懂FFT算法是怎么计算的，这里不讨论傅里叶变换的理论和推导，只讨论实际工程中怎样计算，由于python代码的可读性以及计算的方便性，使用python代码展示FFT计算过程

傅里叶变换

傅里叶变换FT（fourier transform）用于将时域信号x(t)和频域信号X(f)之间变换，公式如下所示：

$X(f)={\int }_{-\infty }^{\infty }x(t)e^{-j2\pi ft}dt={\int }_{-\infty }^{\infty }X(f)e^{j2\pi ft}df$

离散傅里叶变换

由于计算机只能处理有限长的离散信号，因此
因此必须建立对应的离散傅里叶变换DFT（Discrete Fourier Transform）：

$X_k=\sum _{n=0}^{N-1}{x}_n\cdot e^{-i2\pi k\frac{n}{N}}$

如果我们定义一个矩阵$M$
$M_{kn}=e^{-i2\pi k\frac{n}{N}}$

则很明显DFT的公式只是一个简单的线性变换：
$X=x\cdot M$
因此简单的使用矩阵乘法就能计算出DFT的结果，我们可以很容易的写出DFT的python代码

import numpy as np

def DFT(x):
    N = len(x)
    k,n = np.meshgrid(np.arange(N),np.arange(N))
    W = np.exp(-1j*2*np.pi*k*n/N)   
    return np.dot(x,W)

我们以2048点DFT为例，与numpy中内置的FFT做对比，看看速度相差多少

x=np.random.random(2048)
X=np.fft.fft(x)
X=DFT(x)

compute 2048 points dft using np.fft cost: 0.001677 ms
compute 2048 points dft using DFT cost: 0.321912 ms

可以看到，速度相差了差不多2000倍，对于每个值$X(k)$的计算需要N个复数乘法(4N个乘法和2N个加法)和N-1个复数加法(2N-2个加法),因此DFT的总计算量需要$N^2$个复数乘法和$N^2-N$个复数加法，复杂度是$\mathcal{O}[N^2]$，是不利于计算机进行实时信号处理的，因此为了优化DFT的计算量，便有了相关FFT算法，下面介绍快速傅里叶变换算法，对于快速傅里叶逆变换其优化方式非常相似，因此不做介绍

快速傅里叶变换

Cooley-Tukey快速傅里叶算法是常见的FFT算法，其思想是利用了DFT变换中的对称性和周期性来简化计算

首先我们定义
$W_N=e^{-i\frac{2\pi }{N}}$
$W_N$满足如下性质
周期性:$W_N^{k+N}=W_N^k$
对称性:$W_N^{k+\frac{N}{2}}=-W_N^k$
若$m$是$N$的约数:$W_N^{mkn}=W_{\frac{N}{m}}^{kn}$

我们只需几行代码就可验证上述特性

def Wn(k,N):
        return np.exp(-1j*2*np.pi*k/N)

定义如下一些变量：

N = 8
k = 3
m = 2
n = 2

验证周期性：

print(np.allclose(Wn(k,N),-Wn(k+N,,N)))

验证对称性：

print(np.allclose(Wn(k,N),-Wn(k+N//2,,N)))

验证可约性：

print(np.allclose(Wn(m*k*n,N),Wn(k*n,N//m)))

结果：

True
True
True

基2 FFT

根据上面的对称性，我们可以将DFT计算分为两个较小的部分

$X_k=\sum _{n=0}^{N-1}{x}_n\cdot W_N^{kn}$
$=\sum _{m=0}^{N/2-1}{x}_{2m}\cdot W_N^{2mk}+{\sum }_{m=0}^{N/2-1}{x}_{2m+1}\cdot W_N^{(2m+1)k}$
$=\sum _{m=0}^{N/2-1}{x}_{2m}\cdot W_{\frac{N}{2}}^{km}+W_N^k\sum _{m=0}^{N/2-1}{x}_{2m+1}\cdot W_{\frac{N}{2}}^{km}$
$=F_1(k)+W_N^k{F}_2(k)$
这样一个N点变换就分解为了两个N/2点变换，这里$F_1(k)$和$F_2(k)$分别是序列x中的奇数号和偶数号序列的$N/2$点DFT变换，根据以上公式我们也能很快写出python代码：

def R2FFT(x):
   N = len(x)
   N2 = N // 2
   k,n = np.meshgrid(np.arange(N2),np.arange(N2))
   W = np.exp(-1j*2*np.pi*k*n/N2)
   G = np.exp(-2j * np.pi * np.arange(N2) / N)
   X_even = np.dot(x[::2],W)
   X_odd = G*np.dot(x[1::2],W)
   return np.concatenate([X_even+X_odd,X_even-X_odd])

同样计算2048点DFT，速度如下：

compute 2048 points dft using R2FFT cost: 0.081140 ms

对于$N=2^r$,很显然两个N/2点的DFT变换还可以继续分解下去，分解为4个N/4的更短的序列，N/4的序列还可以将序列继续分解下去，直到分解为N/2个2点的DFT变换，2点的DFT变换只需要复数加法和减法就能实现，复数乘法计算量减小至$(N/2)logN$,复数加法计算量减小至$NlogN$，算法复杂度为$\mathcal{O}[NlogN]$，大大减少了DFT的计算量，这就是Cooley-Tukey快速傅里叶变换的基本原理，我们将一个DFT变换分解为两个较小的DFT变换，即基2FFT，我们可以通过递归来实现该算法：

def RecursiveR2FFT(x):
    N = x.shape[0]
    if N <= 2:
        return [x[0]+x[1],x[0]-x[1]]
    else:
        X_even = RecursiveR2FFT(x[::2])
        X_odd = RecursiveR2FFT(x[1::2])
        factor = np.exp(-2j * np.pi * np.arange(N//2) / N)
        return np.concatenate([X_even + factor * X_odd,
                               X_even - factor * X_odd])

计算2048点DFT速度如下：

compute 2048 points dft using RecursiveR2FFT cost: 0.081140 ms

相比前面的版本速度并没有提升,是因为python的递归版本并不高效,并且没有进行并行化的计算,因此,通过观察基2fft的规律我们可以将递归调用的向量乘法转换为并行计算的矩阵乘法以删除递归调用以及并行计算，python代码如下：

def NonRecursiveR2FFT(x):
    L = len(x)
    N_base = 2
    base = L//N_base
    X = np.reshape(x,(-1,base))
    X = np.vstack([X[0]+X[1],X[0]-X[1]]).T
    for n in range(int(np.log2(base))):
        N = X.shape[1]
        W = np.exp(-1j * np.pi * np.arange(N) / N)
        X_even = X[:X.shape[0]//2]
        X_odd = W*X[X.shape[0]//2:]
        X = np.concatenate([X_even+X_odd,X_even-X_odd],axis=-1)
    return X.ravel()

计算2048点DFT速度如下：

compute 2048 points dft using NonRecursiveR2FFT cost: 0.000327 ms

可以看到,速度又提高了一个数量级，相比numpy的fft只差了1倍

基4FFT

当DFT点数N为4的幂时，我们当然可以使用基2FFT算法进行计算，但对于这种情况使用基4FFT算法更为高效，基4FFT的原理与基2FFT类似，只不过是将N点DFT序列拆分成4个N/4的子序列：
$X_k=\sum _{n=0}^{N-1}{x}_n\cdot W_N^{kn}$
$=W_N^0{F}_0(k)+W_N^k{F}_1(k)+W_N^{2k}{F}_2(k)+W_N^{3k}{F}_3(k)$
在这里直接给出非递归的基4FFT代码：

def NonRecursiveR4FFT(x):
    L = len(x)
    N_base = 2
    base = L//N_base
    X = np.reshape(x,(-1,base))
    X = np.vstack([X[0]+X[1],X[0]-X[1]]).T
    butterfly_matrix = np.array([[1,1,1,1],[1,-1j,-1,1j],[1,-1,1,-1],[1,1j,-1,-1j]])
    for n in range(int(np.log(base)/np.log(4))):
        N4 = X.shape[1]
        N = N4*4
        G = np.exp(-2j * np.pi * (np.arange(N4).reshape(1,-1)*np.arange(4).reshape(-1,1)) / N)
        X = G*X.reshape((4,-1,N4)).transpose([1,0,2])
        X = np.dot(butterfly_matrix,X).transpose([1,0,2]).reshape((-1,N))
    return X.ravel()

注意，由于$N=2048=2\cdot 4^5$,因此我们最后分成了1024个2点FFT，如果$N$是4的幂例如N=1024,那么最后会得到2个512点的结果，并不满足基4FFT的条件，那么我们可以将这2个512点序列按照基2FFT原理进行计算，最终得到1024个FFT点的计算结果，这实际上是一个混合了基2FFt和基4FFT的混合基FFT算法。

结论

由于我们的python算法涉及较多的内存分配和复制等操作，因此内存利用率和效率方面还是较numpy的fft差了一些，在实际的FFT算法中会使用C或者汇编来编写高度优化的FFT算法，例如会使用蝶形运算结构来优化计算流程，以最小化内存的使用，对于其他长度的序列，有更复杂的混合基FFT以及一些其他算法。

使用python实现该算法的目的也并不是为了效率，而是因为python代码在可读性方面远优于使用C或者汇编。

尽管在工程中我们有许多开源的FFT库可以使用，但是我相信理解FFT这种基础的算法以及其工程实现是非常有必要的，在必要的情况下，我们也可以按照自己的需求实现特定的FFT算法，这对于我们的技术能力的提升以及知识的积累和沉淀是有非常大的帮助的。

欢迎感兴趣的同学加微信交流：xu19315519614