FFT算法原理及其代码实现

0. FFT简介

FFT(Fast Fourier Transform, FFT)，是实现快速计算序列的离散傅里叶变换(DFT)的方法。它将DFT的复杂度由$\mathrm{O}(n^2)$ 降低到$\mathrm{O}(n \log n)$.

1. DFT的定义及计算

DFT的定义如下：

$$ X[k] = \sum_{n=0}^{N-1} x[n]\mathrm{e}^{-\mathrm{j2\pi} \frac{kn}{N}}, \quad n = 0, \dots, N-1. $$

其中$x[n],n=0,\dots,N-1$，为时域下的信号序列，$\mathrm{j}$为虚数单位。

Octave代码实现如下(实际上不需要做循环，直接矩阵相乘就可得到)：

function X = mydft(x)
  NF = length(x);
  X = zeros(NF,1);
  for k=0:NF-1
    for n=0:NF-1
      X(k+1) = X(k+1)+x(n+1)*exp(-j*2*pi*k*n/NF);
    endfor
  endfor
endfunction

由于计算DFT，需要进行两层循环，故时间复杂度为$\mathrm{O}(n^2)$。当信号序列很长时，使用此算法会消耗很长的时间，故需要FFT算法。

2. FFT原理

FFT算法有多种实现，这里介绍一种常见的FFT算法——库利-图基算法。在此算法中，如果不加其它的处理，只能进行长度为$N=2^n$的信号序列的变换，其中$n \in \mathbb{N}$.

FFT实际上是一种分治算法，将长度为$N$的信号分解为两个长度为$\frac{N}{2}$的信号进行处理。

下面，将进行FFT算法的数学推导：

为了表达简洁，令$W_N = e^{-\mathrm{j}\frac{2\mathrm{\pi}}{N}}$，则由DFT的定义可知，$X[k] = \sum_{n=0}^{N-1} x[n]W_N^{kn}$.

$$ \begin{aligned} X[k] &= \sum_{n=0}^{N-1} x[n]W_N^{kn} \\ &= \sum_{n=0}^{\frac{N}{2}-1}x[2n]W_N^{k(2n)} + \sum_{n=0}^{\frac{N}{2}-1}x[2n+1]W_N^{k(2n+1)}\\ &=\sum_{n=0}^{\frac{N}{2}-1}x[2n]W_N^{k(2n)} + W_N^{k} \sum_{n=0}^{\frac{N}{2}-1}x[2n+1]W_N^{k(2n)}\\ &=Ev[k] + W_N^k Od[k] \end{aligned} $$

其中：${Ev}[k] =\sum_{n=0}^{\frac{N}{2}-1}x[2n]W_N^{k(2n)} $，代表偶数样本点的DFT；${Od}[k] =\sum_{n=0}^{\frac{N}{2}-1}x[2n+1]W_N^{k(2n)} $，代表奇数样本点的DFT。

这样，就将长度为N的信号分解为了两个长度为$\frac{N}{2}$的信号。

另外，在计算中，还可以利用周期性，来简化计算。

其中：

$W_N^{kn} = W_N^{k(n+N)} = W_N^{n(k+N)}$.

$Ev[k] = Ev[k+\frac{N}{2}],\quad Od[k] = Od[k+\frac{N}{2}]$.

3. FFT软件代码实现

Octave代码实现：

function X = myfft(x)
  N = length(x);
  if(N == 1)
    X = [x];
  else

    Xev = myfft(x(1:2:N)); %偶数样本点的DFT
    Xod = myfft(x(2:2:N)); %奇数样本点的DFT

    X = zeros(1,N);
    for k=0:N/2-1
      p = Xev(k+1);
      q = exp(-2*pi*i/N*k) * Xod(k+1);
      X(k+1) = p + q;
      X(k+1+N/2) = p - q;
    endfor
  endif
endfunction

对于代码中的X(k+1) = p + q;，是因为：

$$ \boxed{W[k] = Ev[k] + W_N^{k}Od[k]} $$

k+1是因为Octave的数组索引从1开始，而非从0开始。

对于代码中的X(k+1+N/2) = p - q;，则是利用了$Od[k]$和$Ev[k]$的周期性：

$$ \begin{aligned} X[k+\frac{N}{2}] &=Ev[k+\frac{N}{2}]+W_N^{k+\frac{N}{2}}Od[k+\frac{N}{2}]\\ &=Ev[k] + W_N^{k+\frac{N}{2}}Od[k] \end{aligned} $$

又因为：

$$ \begin{aligned} W_N^{k+\frac{N}{2}} &= \mathrm{e}^{-\mathrm{j}\frac{\mathrm{2\pi}}{N}(k+\frac{N}{2})}\\ &= \mathrm{e}^{-\mathrm{j}(\frac{\mathrm{2\pi}}{N}k-\mathrm{\pi})}\\ &= -\mathrm{e}^{-\mathrm{j}\frac{\mathrm{2\pi}}{N}k}\\ &= -W_N^{k} \end{aligned} $$

所以：

$$ \boxed{W[k+\frac{N}{2}] = Ev[k] - W_N^{k}Od[k]} $$

C代码实现：

此处的C代码只能处理实信号。

double complex* fft(double *x, size_t length)
{
    double complex* X = (double complex*)malloc(sizeof(double complex)*length);
    memset(X,0,sizeof(double complex)*length);
    ditfft2(x,X,length,1);
    return X;
}

// x:信号序列，X:输出的傅里叶序列，N：信号长度，s：数组地址偏移量
void ditfft2(double *x, double complex* X, size_t N, size_t s)
{
    if(N==1)
    {
        X[0] = x[0] + 0*I;
        return;
    }
    
    ditfft2(x,X,N/2,2*s);        //奇
    ditfft2(x+s,X+N/2,N/2,2*s);  //偶
    
    for(size_t k=0; k<N/2; k++)
    {
        double complex p = X[k];
        double complex q = cexp(-2*PI*I/N*k) * X[k+N/2];
        X[k] = p + q;
        X[k + N/2] = p - q;
    }
    return;
}

参考

维基百科.库利-图基快速傅里叶变换算法.

杨宇翔.如何理解和掌握快速傅里叶变换的计算和概念.知乎

Wikipedia.Cooley–Tukey FFT algorithm. (维基百科中文页面和英文页面的内容并不相同)