C++ Cooley-Tukey FFT 算法

2025年3月22日 | 阅读 12 分钟

引言

Cooley-Tukey 快速傅里叶变换 (FFT) 算法是一种广泛使用且高效的计算复数序列或数组离散傅里叶变换 (DFT) 的方法。它由 J.W. Cooley 和 John Tukey 于 1965 年提出，自此成为许多信号处理和数值分析应用中的基本算法。

Cooley-Tukey FFT 算法的基本思想是将任何复合大小 N = N1 * N2 的 DFT 递归地分解为许多更小的、大小为 N1 和 N2 的 DFT，并结合周期性旋转因子。这种递归分解一直持续到大小为 2 的 DFT 的基本情况。关键的见解是，通过利用旋转因子的周期性并重复利用递归计算，可以显著减少整体计算量。

方法 1：使用递归实现

通过将计算分解为更小的子问题直到达到基本情况，以递归方式实现 Cooley-Tukey FFT 算法。复数表示输入和中间值，从而在递归过程中可以有效地处理实部和虚部。

程序

#include <iostream>
#include <complex>
#include <cmath>
#include <vector>
typedef std::complex<double> Complex;
void fft(std::vector<Complex>& a, bool invert) {
    int n = a.size();
    if (n <= 1) return;
    // Divide the array into even and odd parts
    std::vector<Complex> a0(n / 2), a1(n / 2);
    for (int i = 0, j = 0; i < n; i += 2, ++j) {
        a0[j] = a[i];
        a1[j] = a[i + 1];
    }
    // Recursively compute FFT on even and odd parts
    fft(a0, invert);
    fft(a1, invert);
    // Combine the results
    double angle = (invert ? 2 : -2) * M_PI / n;
    Complex w(1), wn(cos(angle), sin(angle));
    for (int i = 0; i < n / 2; ++i) {
        Complex t = w * a1[i];
        a[i] = a0[i] + t;
        a[i + n / 2] = a0[i] - t;
        if (invert) {
            a[i] /= 2;
            a[i + n / 2] /= 2;
        }
        w *= wn;
    }
}
// Wrapper function for convenience
void cooleyTukeyFFT(std::vector<Complex>& a) {
    fft(a, false);
}
// Wrapper function for convenience
void inverseCooleyTukeyFFT(std::vector<Complex>& a) {
    fft(a, true);
}
int main() {
    // Example usage:
    std::vector<Complex> input = {1, 2, 3, 4};
    // Compute FFT
    cooleyTukeyFFT(input);
    // Output FFT result
    for (const auto& val : input) {
        std::cout << val << " ";
    }
    std::cout << std::endl;
    // Compute inverse FFT
    inverseCooleyTukeyFFT(input);
    // Output inverse FFT result
    for (const auto& val : input) {
        std::cout << val << " ";
    }
    std::cout << std::endl;
    return 0;
}

输出

(10,0) (-2,2) (-2,0) (-2,-2) 
(1,0) (2,5.72119e-18) (3,0) (4,-5.72119e-18)

说明

包含必要的头文件

代码首先包含必要的 C++ 标准库头文件。这些头文件提供了处理复数 ()、数学运算 () 和输入/输出 () 的功能。

定义复数类型

使用 typedef 语句基于头文件中的 std::complex<double> 类型定义一个复数类型 (Complex)。此类型表示具有实部和虚部的数字。

FFT 函数

定义 fft 函数以递归方式执行 Cooley-Tukey 快速傅里叶变换。

它接受一个复数向量 (a) 和一个布尔标志 (invert)，该标志指示是计算正向 FFT 还是逆向 FFT。

该函数将输入分成偶数和奇数部分，对这些部分递归地计算 FFT，并使用周期性旋转因子组合结果。

包装器函数

两个包装器函数 cooleyTukeyFFT 和 inverseCooleyTukeyFFT 分别提供了一个方便的接口来计算 FFT 及其逆。它们使用适当的参数调用 fft 函数。

主函数

main 函数演示了在示例向量 (input) 上使用 FFT 函数。

它计算 FFT，打印结果，计算逆 FFT，并打印逆 FFT 结果。

输出到控制台

使用 std::cout 语句显示输出。更正后的代码包含了必要的 #include 指令，以解决与 std::cout 使用相关的编译错误。

执行流程

程序首先包含必要的库，定义复数类型，并声明函数。

在 main 函数中，使用示例向量来演示 FFT 计算。计算 FFT 及其逆，并将结果打印到控制台。

复杂度分析

时间复杂度

Cooley-Tukey FFT 算法由于其递归分治方法，时间复杂度为 O(n log n)。

在递归的每个级别，都有 O(n) 次操作，并且递归深度为 log(n)，因为输入被重复地减半。

空间复杂度

由于递归栈，该算法的空间复杂度为 O(n log n)。

在每次递归调用中，算法都需要额外的空间来存储临时数组（例如 a0 和 a1），并且最大递归深度为 log(n)。

总空间复杂度为 O(n log n)。

方法 2：使用迭代实现

通过采用基于循环的方法来迭代实现 Cooley-Tukey FFT 算法。将计算分解为多个阶段，每个阶段通过成对计算来管理多组蝶形运算。利用循环遍历这些阶段，优化算法以实现高效计算。这种迭代实现提供了 FFT 阶段的顺序进展，改善了缓存局部性并减少了内存访问模式。

程序

#include <iostream>
#include <complex>
#include <cmath>
#include <vector>
typedef std::complex<double> Complex;
void iterativeCooleyTukeyFFT(std::vector<Complex>& a) {
    int n = a.size();
    int levels = log2(n);
    for (int s = 1; s <= levels; ++s) {
        int m = 1 << s; // 2^s
        int halfM = m / 2;
        Complex wm = exp(Complex(0, -2.0 * M_PI / m));
        for (int k = 0; k < n; k += m) {
            Complex w = 1.0;
            for (int j = 0; j < halfM; ++j) {
                Complex t = w * a[k + j + halfM];
                Complex u = a[k + j];
                a[k + j] = u + t;
                a[k + j + halfM] = u - t;
                w *= wm;
            }
        }
    }
}
int main() {
    // Example usage:
    std::vector<Complex> input = {1, 2, 3, 4};
    // Compute FFT iteratively
    iterativeCooleyTukeyFFT(input);
    // Output FFT result
    for (const auto& val : input) {
        std::cout << val << " ";
    }
    std::cout << std::endl;
    return 0;
}

输出

(10,0) (-1,1) (-4,0) (-1,-1)

说明

包含必要的头文件

代码首先包含必要的 C++ 标准库头文件，例如用于输入/输出的和用于数学函数的。

定义复数类型

使用 typedef 语句基于头文件中的 std::complex 类型定义一个复数类型 (Complex)。此类型表示具有实部和虚部的数字。

迭代 Cooley-Tukey FFT 函数

iterativeCooleyTukeyFFT 函数迭代地执行 Cooley-Tukey FFT。

它使用循环遍历 FFT 算法的不同阶段。每个阶段处理蝶形运算组，并且循环变量控制这些组的大小和位置。

循环遍历阶段

由变量 s 控制的外部循环遍历 FFT 算法的各个阶段。阶段数由输入大小的 2 为底的对数确定。

循环遍历蝶形运算组

在每个阶段内，由变量 k 控制的内部循环遍历蝶形运算组。每个组代表一对复数，它们将进行计算。

旋转因子计算

在内循环中，算法计算旋转因子 (wm) 并使用它对组内的复数对执行蝶形运算。

输出到控制台

main 函数演示了在示例输入向量 (input) 上使用迭代 FFT 函数。

然后使用 std::cout 将 FFT 结果打印到控制台。

执行流程

程序首先包含必要的库，定义复数类型，并声明函数。

在 main 函数中，使用示例向量来演示迭代 FFT 计算。计算 FFT 结果并将其打印到控制台。

复杂度分析

时间复杂度

迭代 Cooley-Tukey FFT 算法的时间复杂度为 O(n log n)。

该算法遍历 log(n) 个阶段，在每个阶段内，它执行 O(n) 次操作（蝶形运算）。

空间复杂度

空间复杂度为 O(1) 辅助空间。

迭代方法使用的额外空间是恒定的，与输入大小无关。主要的修改数据结构是输入向量本身。

方法 3：使用位反转置换

将 位反转置换 整合到 Cooley-Tukey FFT 算法中，以重新排列输入数据，从而提高缓存局部性并最小化内存访问模式。该技术涉及以与内存结构对齐的方式重新排序元素，这可能会优化 FFT 计算期间的数据检索。位反转置换为改进内存访问模式提供了有效的策略，使得算法更能适应现代计算机体系结构。

程序

#include <iostream>
#include <complex>
#include <cmath>
#include <vector>
typedef std::complex<double> Complex;
// Function to perform bit-reversal permutation
void bitReversePermutation(std::vector<Complex>& a) {
    int n = a.size();
    int bits = log2(n);
    for (int i = 0; i < n; ++i) {
        int reversed = 0;
        for (int j = 0; j < bits; ++j) {
            reversed |= ((i >> j) & 1) << (bits - 1 - j);
        }
        if (reversed > i) {
            std::swap(a[i], a[reversed]);
        }
    }
}
// Iterative Cooley-Tukey FFT with bit-reversal permutation
void iterativeCooleyTukeyFFT(std::vector<Complex>& a) {
    bitReversePermutation(a);
    int n = a.size();
    int levels = log2(n);
    for (int s = 1; s <= levels; ++s) {
        int m = 1 << s; // 2^s
        int halfM = m / 2;
        Complex wm = exp(Complex(0, -2.0 * M_PI / m));
        for (int k = 0; k < n; k += m) {
            Complex w = 1.0;
            for (int j = 0; j < halfM; ++j) {
                Complex t = w * a[k + j + halfM];
                Complex u = a[k + j];
                a[k + j] = u + t;
                a[k + j + halfM] = u - t;
                w *= wm;
            }
        }
    }
}
int main() {
    // Example usage:
    std::vector<Complex> input = {1, 2, 3, 4, 5, 6, 7, 8};
    // Perform bit-reversal permutation
    bitReversePermutation(input);
    // Compute FFT iteratively with bit-reversed input
    iterativeCooleyTukeyFFT(input);
    // Output FFT result
    for (const auto& val : input) {
        std::cout << val << " ";
    }
    std::cout << std::endl;
    return 0;
}

输出

(36,0) (-1,2.41421) (-4,4) (-1,0.414214) (-16,0) (-1,-0.414214) (-4,-4) (-1,-2.41421)

说明

位反转置换

bitReversePermutation 函数使用位反转置换来重新排列输入数据。

它遍历输入向量，根据其位反转索引交换元素，并确保只执行一次交换以避免冗余。

带有位反转的迭代 Cooley-Tukey FFT

iterativeCooleyTukeyFFT 函数结合了位反转置换和 Cooley-Tukey FFT 算法。

它首先使用 bitReversePermutation 函数将位反转置换应用于输入向量。

然后，它迭代地对位反转的输入执行 Cooley-Tukey FFT 算法。

主函数

在 main 函数中，提供了示例输入，即样本向量 (input)。

调用 bitReversePermutation 以位反转顺序重新排列输入数据。

然后使用 iterativeCooleyTukeyFFT 函数对位反转的数据计算 FFT。

最后，将结果打印到控制台。

执行流程

程序首先对输入数据执行位反转置换。

之后，它将迭代 Cooley-Tukey FFT 算法应用于位反转的输入。

将结果打印到控制台。

位反转的目的

采用位反转置换是为了提高 FFT 计算期间的缓存局部性并减少内存访问模式。

它以与内存结构对齐的方式重新排列输入数据，从而可能提高性能。

提高缓存局部性

结合位反转和 FFT 计算旨在优化内存访问模式，以更好地利用缓存，这有助于更高效的处理。

复杂度分析

时间复杂度

位反转置换 (bitReversePermutation)

位反转置换的时间复杂度为 O(n log n)，其中 n 是输入向量的大小。

该函数遍历向量中的每个元素，并对每个元素执行位反转操作，该操作需要 O(log n) 的时间。

迭代 Cooley-Tukey FFT (iterativeCooleyTukeyFFT)

Cooley-Tukey FFT 算法的时间复杂度为 O(n log n)。

它涉及 log(n) 个阶段，每个阶段有 n/2 次操作（蝶形运算）。

主函数

main 函数涉及调用位反转置换和迭代 Cooley-Tukey FFT。

由于 Cooley-Tukey FFT 是主要的计算，所以 main 函数的总时间复杂度由其时间复杂度决定，即 O(n log n)。

空间复杂度

位反转置换 (bitReversePermutation)

位反转置换是就地操作的，直接修改输入向量。

空间复杂度为 O(1)，因为它使用恒定的额外空间。

迭代 Cooley-Tukey FFT (iterativeCooleyTukeyFFT)

迭代 Cooley-Tukey FFT 的空间复杂度为 O(1) 辅助空间。

它就地修改输入向量，而无需额外的与输入大小成比例的空间。

主函数

由于 main 函数不使用除输入向量之外的大量额外空间，因此其总空间复杂度为 O(1)。

方法 4：使用多线程进行并行化

通过使用多线程方法将计算分配给多个线程来并行化 Cooley-Tukey FFT 算法。它涉及将工作负载分解为并发任务，从而允许在不同处理单元上同时执行。可以使用 OpenMP 或 std::thread 等流行库来实现并行算法，从而提高整体性能。

程序

#include <iostream>
#include <complex>
#include <cmath>
#include <vector>
#include <omp.h>
typedef std::complex<double> Complex;
// Recursive FFT function (same as before)
void fft(std::vector<Complex>& a, bool invert) {
    int n = a.size();
    if (n <= 1) return;
    std::vector<Complex> a0(n / 2), a1(n / 2);
    for (int i = 0, j = 0; i < n; i += 2, ++j) {
        a0[j] = a[i];
        a1[j] = a[i + 1];
    }
    fft(a0, invert);
    fft(a1, invert);
    double angle = (invert ? 2 : -2) * M_PI / n;
    Complex w(1), wn(cos(angle), sin(angle));
    for (int i = 0; i < n / 2; ++i) {
        Complex t = w * a1[i];
        a[i] = a0[i] + t;
        a[i + n / 2] = a0[i] - t;
        if (invert) {
            a[i] /= 2;
            a[i + n / 2] /= 2;
        }
        w *= wn;
    }
}
// Parallel FFT function using OpenMP
void parallelCooleyTukeyFFT(std::vector<Complex>& a) {
    #pragma omp parallel
    {
        #pragma omp single nowait
        fft(a, false);
    }
}
int main() {
    // Example usage:
    std::vector<Complex> input = {1, 2, 3, 4};
    // Compute FFT in parallel
    parallelCooleyTukeyFFT(input);
    // Output FFT result
    for (const auto& val : input) {
        std::cout << val << " ";
    }
    std::cout << std::endl;
    return 0;
}

输出

(10,0) (-2,2) (-2,0) (-2,-2)

说明

递归 FFT 函数 (fft)

fft 函数是 Cooley-Tukey FFT 算法的递归实现。

它接受一个复数向量 (a) 和一个布尔标志 (invert) 来控制 FFT 的方向（正向或逆向）。

该函数递归地将 FFT 计算分解为更小的子问题，直到达到基本情况。

并行 FFT 函数 (parallelCooleyTukeyFFT)

parallelCooleyTukeyFFT 函数是使用 OpenMP 指令的 Cooley-Tukey FFT 算法的并行化版本。

它使用 #pragma omp parallel 创建一个线程池以进行并行执行。

#pragma omp single nowait 指令确保递归 FFT 仅由一个线程执行一次，而无需等待。其他线程继续工作。

nowait 子句用于允许线程独立继续，而无需等待单个线程完成。

主函数

在 main 函数中，提供了示例输入，即样本向量 (input)。

调用 parallelCooleyTukeyFFT 函数以并行方式使用 OpenMP 计算 FFT。

执行流程

程序首先包含必要的库并定义复数类型。

在 main 函数中，使用示例向量，并调用并行 FFT 函数。

然后，并行 FFT 函数利用多个线程并发执行 FFT 计算。

复杂度分析

时间复杂度

递归 FFT 函数 (fft)

递归 FFT 函数的时间复杂度为 O(n log n)，其中 n 是输入向量的大小。这与非并行版本相同。

FFT 算法涉及 log(n) 个阶段，每个阶段有 n/2 次操作（蝶形运算）。

并行 FFT 函数 (parallelCooleyTukeyFFT)

并行 FFT 函数的时间复杂度仍为 O(n log n)。

尽管并行化可以将工作负载分配给多个线程，但 FFT 算法的基本时间复杂度保持不变。

主函数

由于主要的计算是并行 FFT，因此程序的总时间复杂度为 O(n log n)。

空间复杂度

递归 FFT 函数 (fft)

递归 FFT 函数的空间复杂度为 O(log n)，这是由于递归栈。

在递归的每个级别，都会创建临时数组 (a0 和 a1)，但它们会被重用，因此递归栈决定了空间复杂度。

并行 FFT 函数 (parallelCooleyTukeyFFT)

并行 FFT 函数使用与非并行版本相同的空间复杂度。

每个线程操作输入向量的不同部分，但它们共享相同的内存空间，并且由于递归，空间复杂度仍为 O(log n)。

主函数

由于 FFT 函数中的递归栈，程序的总空间复杂度为 O(log n)。

下一主题std-to-underlying-function-in-cpp

C++ Cooley-Tukey FFT 算法

引言

方法 1：使用递归实现

程序

方法 2：使用迭代实现

程序

方法 3：使用位反转置换

程序

方法 4：使用多线程进行并行化

程序

联系信息

关注我们

教程

面试题

在线编译器

Python

Java

.Net Framework

AI, ML and Data Science

Cloud Technology

B.Tech and MCA

Web Technology

PHP

Software Testing

Technical Interview

Java Interview

Python

Web Interview

Database Interview

B.Tech / MCA

Important Interview

Software Testing Interview

Company Interviews

Online Compilers

Multiple Choice Questions

C++ 教程

C++ 控制语句

C++ 函数

C++ 数组

C++ 类和对象

C++ 构造函数

C++ 继承

C++ 多态

C++ 抽象

C++ 命名空间

C++ 模板

C++ 字符串

C++ 指针

信号处理

C++ 异常

C++ 文件与流

C++ STL 教程

面试题

选择题

C++ 程序

C++ STL Stack

C++ STL Bitset

C++ STL Deque

C++ STL List

C++ STL Map

C++ STL Math

C++ STL priority_queue

C++ STL Queue

C++ STL Multiset

C++ STL Multimap

C++ STL Set

C++ STD Strings

C++ STL Vector

C++ 操纵符

C++ STL Algorithms

C++ Algorithm

C++ Iterators

C++ 杂项

C++ Cooley-Tukey FFT 算法

引言

方法 1：使用递归实现

程序

方法 2：使用迭代实现

程序

方法 3：使用位反转置换

程序

方法 4：使用多线程进行并行化

程序

相关帖子

C++ 中的完美总计数

C++ 中的松鼠模拟

C++ 与 Eiffel 的区别

C++ 中 Tarjan 算法查找强连通分量

C++ std::common_type T1, T2 ::type 函数

C++ 与 Julia 的区别

C++ std::pmr::synchronized_pool_resource

C++ 中的 Enneacontahexagon 数

C++ 中分段筛（打印范围内的素数）

C++ 贝塞尔插值

订阅 Tpoint Tech