快速傅里叶变换（FFT）原理及其加速实现

快速傅里叶变换（FFT）原理详解

快速傅里叶变换（Fast Fourier Transform, FFT）是计算离散傅里叶变换（Discrete Fourier Transform, DFT）的高效算法。FFT广泛应用于信号处理、图像处理、音频处理等领域。其核心思想是通过分治法（Divide and Conquer）将计算复杂度从 $O(N^2)$ 降低到 $O(N \log N)$。

离散傅里叶变换（DFT）

DFT的目的是将时域信号转换为频域信号。对于长度为 $N$ 的序列 $x[n]$，其DFT定义为：

$X[k] = \sum_{n=0}^{N-1} x[n] e^{-j \frac{2 \pi}{N} kn}, \quad k = 0, 1, \ldots, N-1$

其中，$X[k]$ 是频域信号，$x[n]$ 是时域信号，$j$ 是虚数单位，$\frac{2 \pi}{N}$ 是频率分辨率。

逆离散傅里叶变换（IDFT）

IDFT用于将频域信号转换回时域信号，其定义为：

$x[n] = \frac{1}{N} \sum_{k=0}^{N-1} X[k] e^{j \frac{2 \pi}{N} kn}, \quad n = 0, 1, \ldots, N-1$

快速傅里叶变换（FFT）

FFT利用DFT的对称性和周期性，通过分治法将计算复杂度从 $O(N^2)$ 降低到 $O(N \log N)$。

分治法

FFT 使用分治法 (Divide and Conquer) 将 $N$ 点 DFT 分解为两个 $\frac{N}{2}$ 点 DFT。这些 $\frac{N}{2}$ 点的 DFT 可以继续分解，直到分解到最小的 2 点 DFT。

假设 $N$ 是 2 的幂（如果不是可以用零填充补足），可以将序列 $x[n]$ 分成偶数项 $x[2m]$ 和奇数项 $x[2m+1]$，于是原始的 DFT 可以表示为：

$X[k] = \sum_{m=0}^{N/2-1} x[2m] \cdot e^{-j \frac{2\pi}{N} (2m)k} + \sum_{m=0}^{N/2-1} x[2m+1] \cdot e^{-j \frac{2\pi}{N} (2m+1)k}$

其中，偶数项和奇数项分别是 $\frac{N}{2}$ 点 DFT。

递归计算

通过递归地计算这些 $\frac{N}{2}$ 点 DFT，我们可以将 $N$ 点 DFT 的计算复杂度降到 $O(N\log N)$。

蝶形结构

在合并这些小的 DFT 结果时，我们利用了旋转因子的对称性和周期性，减少了重复计算：

$W_N^k = e^{-j \frac{2\pi}{N} k}$

通过利用这些性质，我们可以使用蝶形结构 (Butterfly Structure) 高效地合并结果，从而显著减少了计算量。

补充：选看——Cooley-Tukey算法

关于上述算法，有一个正式的名字：Cooley-Tukey 算法。Cooley-Tukey 算法是最常用的快速傅里叶变换（FFT）算法，通过递归地将离散傅里叶变换（DFT）分解为较小的 DFT，从而显著减少计算复杂度。让我们更详细地看看 Cooley-Tukey 算法的具体实现和原理。

Cooley-Tukey FFT 算法

Cooley-Tukey 算法的核心思想是将一个 $N$ 点的 DFT 分解成多个较小的 DFT，具体步骤如下：

1. 分解：将输入序列 $x[n]$ 分解为偶数序列和奇数序列。
2. 递归计算：递归地计算这些较小序列的 DFT。
3. 合并：利用蝶形结构合并较小 DFT 的结果。

算法详细步骤

1. 输入：长度为 $N$ 的信号 $x[n]$，其中 $N$ 是 2 的幂。

2. 分解： 将输入信号分成偶数和奇数部分：

   $x_{\text{even}}[m] = x[2m], \quad x_{\text{odd}}[m] = x[2m + 1], \quad m = 0, 1, \ldots, \frac{N}{2} - 1$

3. 递归计算： 分别计算偶数和奇数部分的 DFT：

   $X_{\text{even}}[k] = \text{DFT}(x_{\text{even}}), \quad X_{\text{odd}}[k] = \text{DFT}(x_{\text{odd}}), \quad k = 0, 1, \ldots, \frac{N}{2} - 1$

4. 合并： 利用旋转因子 $W_N^k = e^{-j \frac{2\pi}{N} k}$，合并结果：

   $X[k] = X_{\text{even}}[k] + W_N^k \cdot X_{\text{odd}}[k]$

   $X[k + N/2] = X_{\text{even}}[k] - W_N^k \cdot X_{\text{odd}}[k]$

   其中，$k = 0, 1, \ldots, \frac{N}{2} - 1$。

计算复杂度分析

• 分解：每次将问题规模减半，进行两次 $\frac{N}{2}$ 点的 DFT 计算。
• 递归：每次递归调用会进一步分解问题，直到最终的 2 点 DFT。
• 合并：每层的合并操作需要 $O(N)$ 次计算。

通过对上述过程的分析，总的计算复杂度为：$T(N) = 2T(N/2) + O(N)$

根据递推关系，可以得到 FFT 的总复杂度为： $T(N) = O(N \log N)$

示例

以长度为8的序列为例：

$x = [x_0, x_1, x_2, x_3, x_4, x_5, x_6, x_7]$

将其分解为偶数项和奇数项：

$x_{\text{even}} = [x_0, x_2, x_4, x_6], \quad x_{\text{odd}} = [x_1, x_3, x_5, x_7]$

递归计算 $x_{\text{even}}$ 和 $x_{\text{odd}}$ 的DFT，最终合并得到原序列的DFT。

最后，通过蝶形运算合并结果：

$X[0]=E[0]+W_8^0​O[0], X[4]=E[0]−W_8^0​O[0]$

$X[1]=E[1]+W_8^1​O[1],X[5]=E[1]−W_8^1​O[1]$

$X[2]=E[2]+W_8^2​O[2],X[6]=E[2]−W_8^2​O[2]$

$X[3]=E[3]+W_8^3​O[3],X[7]=E[3]−W_8^3​O[3]$

其中，$W_8^k = e^{-j \frac{2 \pi}{8} k}$是旋转因子。

通过上述过程，可以看到FFT的高效性来源于利用了DFT的对称性和周期性，从而大幅减少了计算量。蝶形运算在这个过程中起到了关键作用，是FFT算法的基本操作单元。

实现

#ifndef FFT_H
#define FFT_H
#include <stdint.h>
#include "common.h"

typedef struct {
    double real;
    double imag;
} Complex;

Complex * fftTwiddleFactor(int size);
int ones_32(uint32_t n);
uint32_t floor_log2_32(uint32_t x);
void fftChange(Complex *A, int size);

Complex * fft(Complex *x, uint32_t N);
Complex * ifft(Complex *x, uint32_t N);

int fftshift(Complex *x, uint32_t a);

extern const Complex *fftTwiddleFactorTables[];
#endif

ones_32

ones_32函数计算一个32位整数的二进制表示中“1”的个数。这个功能通常用于位操作或其他需要知道某个数的二进制中“1”的个数的情况。

int ones_32(uint32_t n)
{
    unsigned int c = 0;
    for (c = 0; n; ++c)
    {
        n &= (n - 1);
    }
    return c;
}

这个函数的工作原理是利用一个小技巧：每次将当前数 n 与 n-1 进行按位与操作，会将 n 的最低位的一个“1”变成“0”。因此，这个循环会执行 c 次，其中 c 是 n 的二进制表示中“1”的个数。

floor_log2_32

floor_log2_32函数计算一个正整数的二进制对数的下取整值（floor of log base 2）。该函数在确定FFT的分解级数时非常有用。

uint32_t floor_log2_32(uint32_t x)
{
    x |= (x >> 1);
    x |= (x >> 2);
    x |= (x >> 4);
    x |= (x >> 8);
    x |= (x >> 16);

    return (ones_32(x >> 1));
}

这个函数的工作原理是通过一系列的按位或操作，将 x 的最高位“1”之后的所有位都变成“1”。然后，通过调用 ones_32 函数，计算出这些“1”的个数，这个个数即为 x 的二进制对数的下取整值。

fftTwiddleFactor

fftTwiddleFactor函数生成并返回FFT算法中所需的旋转因子（twiddle factors），这些旋转因子用于蝶形运算中的复数乘法。

Complex * fftTwiddleFactor(int size)
{
    int i;
    Complex * W = (Complex *)malloc(sizeof(Complex) * size);
    double temp = 2 * M_PI / size;
    for (i = 0; i < size; i++)
    {
        W[i].real = cos(temp * i);
        W[i].imag = -sin(temp * i);
    }
    return W;
}

这个函数的工作原理是：对于给定的FFT规模 size，计算每个旋转因子。每个旋转因子的实部和虚部分别是 cos(2 * π * i / size) 和 -sin(2 * π * i / size)。这些旋转因子在FFT的每个蝶形运算中都会用到。

fftChange

fftChange函数实现了位反转排列（bit-reversal permutation），这是一种数据重排，使得FFT蝶形运算能够正确进行。

void fftChange(Complex *A, int size)
{
    Complex temp;
    int32_t i = 0, j = 0, k = 0;
    int t;
    int t_temp = floor_log2_32(size);
    for (i = 0; i < size; i++)
    {
        k = i;
        j = 0;
        t = t_temp;
        while ((t--) > 0)
        {
            j = j << 1;
            j |= (k & 1);
            k = k >> 1;
        }
        if (j > i)
        {
            temp = A[i];
            A[i] = A[j];
            A[j] = temp;
        }
    }
}

这个函数的工作原理是：对于每个索引 i，计算出其位反转后的索引 j，如果 j 大于 i，则交换 A[i] 和 A[j]。这样，整个数组就按照位反转顺序重排了。

bit_reverse

bit_reverse函数也是实现位反转排列的另一种方法。作用与fftChange函数相同。

void bit_reverse(Complex* x, int n)
{
    Complex temp;
    for (int i = 0, j = 1; j < n - 1; j++)
    {
        for (int k = n >> 1; k > (i ^= k); k >>= 1);
        if (i < j)
        {
            temp = x[i];
            x[i] = x[j];
            x[j] = temp;
        }
    }
}

这个函数的工作原理是：通过一个双层循环，逐步计算出位反转的索引，并进行相应的数据交换，使得整个数组按照位反转顺序重排。

fftshift

fftshift函数实现了频谱搬移操作，将频域信号从零频点中心移到数组中心。这在频谱分析和信号处理的可视化中非常有用。

int fftshift(Complex *x, uint32_t a)
{
    Complex temp;
    uint32_t n = a / 2;
    for (uint32_t i = 0; i < n; i++)
    {
        temp = x[i];
        x[i] = x[i + n];
        x[i + n] = temp;
    }
    return 0;
}

这个函数的工作原理是：将数组前半部分与后半部分交换，从而实现频谱的搬移。这样零频点就从数组的开始位置移到了中间位置。

FFT

fft函数使用了Cooley-Tukey FFT算法，该算法通过分治法将DFT的计算复杂度从 $O(N^2)$ 降低到 $O(N \log N)$。函数的具体实现包括以下几个步骤：

1. 位反转排列（Bit-Reversal Permutation） ：

   • fftChange(x, N)：这个步骤通过fftChange函数将输入数据重新排列。具体来说，它将每个元素的位置根据其二进制表示进行反转排列。例如，原位置为0011的元素会被移动到位置1100。

2. 生成旋转因子（Twiddle Factors） ：

   • Complex *W = fftTwiddleFactor(N)：这个步骤生成一个长度为N的旋转因子表W。这些因子用于后续的蝶形运算，减少了每次运算时重新计算旋转因子的开销。

3. 蝶形运算（Butterfly Operation） ：

   • for (step = 1; step < N; step <<= 1)：外层循环控制蝶形运算的级数，从最底层的两个元素开始，逐步合并到整个序列。

   • for (group = 0; group < N; group += (step << 1))：中层循环控制每一层中不同组的蝶形运算。

   • for (pair = 0; pair < step; pair++)：内层循环控制每一组内的蝶形运算。

   • temp 变量用于存储蝶形运算的临时结果。蝶形运算的核心是将两个频谱元素进行加法和减法运算，并乘以相应的旋转因子，更新结果。

4. 释放旋转因子表的内存：

   • free(W)：在完成所有运算后，释放之前分配的旋转因子表的内存，避免内存泄漏。

Complex * fft(Complex *x, uint32_t N)
{
    fftChange(x, N);
    Complex *W = fftTwiddleFactor(N);
    int step, group, pair;
    Complex temp;
    for (step = 1; step < N; step <<= 1)
    {
        for (group = 0; group < N; group += (step << 1))
        {
            for (pair = 0; pair < step; pair++)
            {
                temp.real = W[pair * (N / (step << 1))].real * x[group + pair + step].real -
                            W[pair * (N / (step << 1))].imag * x[group + pair + step].imag;
                temp.imag = W[pair * (N / (step << 1))].imag * x[group + pair + step].real +
                            W[pair * (N / (step << 1))].real * x[group + pair + step].imag;
                x[group + pair + step].real = x[group + pair].real - temp.real;
                x[group + pair + step].imag = x[group + pair].imag - temp.imag;
                x[group + pair].real += temp.real;
                x[group + pair].imag += temp.imag;
            }
        }
    }
    free(W);
    return x;
}

fftTwiddleFactorTable.c

这个文件包含了预先计算好的旋转因子表，用于FFT计算

加速原理

1. 预先计算的旋转因子表：fftTwiddleFactorTable.c文件中包含了各种规模的预先计算好的旋转因子表。这些旋转因子在FFT算法的蝶形运算中被频繁使用。预先计算并存储这些旋转因子表可以减少FFT计算过程中的重复计算，从而提高计算效率。
2. 位反转排列：fftChange 和 bit_reverse 函数实现了对输入数据的位反转排列，这是一种预处理步骤，使得FFT的蝶形操作能够正确进行。这样可以将输入序列按照频率排序，从而减少计算的复杂度。
3. 生成旋转因子：fftTwiddleFactor 函数生成旋转因子表，用于后续的蝶形运算。这些旋转因子是FFT计算中必不可少的，它们根据FFT的规模预先计算好，减少了重复计算的开销。
4. 蝶形运算：fft 函数中，通过双重循环实现了FFT的核心蝶形运算。在每一层循环中，逐步合并部分频谱结果，直至完成整个频谱的计算。每次合并时，利用预先计算好的旋转因子进行复数乘法和加法运算。