首页 > 编程语言 >Cooley-Tukey FFT算法的非递归实现

Cooley-Tukey FFT算法的非递归实现

时间：2024-09-08 23:27:43浏览次数：7  
标签：std Cooley point int FFT complex ++ Tukey

一维情况

#include <iostream>

#include <complex>
#include <cmath>

const double PI = 3.14159265358979323846;

// 交换位置 
void swap(std::complex<double>& a, std::complex<double>& b) {
    std::complex<double> temp = a;
    a = b;
    b = temp;
}

// 频率倒位 
void bitReverse(std::complex<double>* point, int N) {
    int j = 0;
    for (int i = 1; i < N - 1; ++i) {
        int k = N / 2;
        while (j >= k) {
            j -= k;
            k /= 2;
        }
        j += k;
        if (i < j) {
            swap(point[i], point[j]);
        }
    }
}

// Cooley-Tukey算法的非递归实现
void FFT(std::complex<double>* point, int N) {
    bitReverse(point, N);

    for (int s = 2; s <= N; s *= 2) { // 蝶形运算的Stage循环
        double angle = -2 * PI / s;
        std::complex<double> w(1, 0);
        std::complex<double> wn(cos(angle), sin(angle));

        for (int k = 0; k < N; k += s) {  // 一组蝶形运算group，注：每个group的蝶形因子乘数不同
            w = std::complex<double>(1, 0);
            for (int j = 0; j < s / 2; ++j) { // 一个蝶形运算  注：group内的蝶形运算
                std::complex<double> u = point[k + j];
                std::complex<double> v = point[k + j + s / 2] * w;
                point[k + j] = u + v;
                point[k + j + s / 2] = u - v;

                w *= wn;
            }
        }
    }
}

// FFT逆变换
void IFFT(std::complex<double>* point, int N) {
    // 共轭变换(即：实部不变，虚部取相反数)
    for (int i = 0; i < N; ++i) {
        point[i] = std::conj(point[i]);
    }

    // 应用FFT
    FFT(point, N);

    // 共轭变换(即：实部不变，虚部取相反数)并归一化
    for (int i = 0; i < N; ++i) {
        point[i] = std::conj(point[i]);
        point[i] /= N;
    }
}

// FFT变换测试
void FFTTest()
{
    const int N = 8;  // 输入序列的长度
    std::complex<double> input[N] = { 0, 1, 2, 3, 4, 5, 6, 7 };  // 输入序列

    // 将输入序列扩展为2的幂次方长度
    int paddedLength = (int)pow(2, ceil(log2(N)));
    std::complex<double>* paddedInput = new std::complex<double>[paddedLength];
    for (int i = 0; i < N; ++i) {
        paddedInput[i] = input[i];
    }
    // 补零到2的n次幂
    for (int i = N; i < paddedLength; ++i) {
        paddedInput[i] = 0;
    }

    // 执行FFT变换
    FFT(paddedInput, paddedLength);

    std::cout << "========== FFT变换 ========== " << std::endl;

    // 输出原始数据
    std::cout << "原始数据:----------" << std::endl;
    for (int k = 0; k < N; ++k) {
        std::cout << "f[" << k << "] = " << input[k] << std::endl;
    }

    // 输出变换结果
    std::cout << "变换结果:----------" << std::endl;
    for (int k = 0; k < paddedLength; ++k) {
        std::cout << "F[" << k << "] = " << paddedInput[k] << std::endl;
    }

    // 输出幅值结果
    std::cout << "幅值结果:----------" << std::endl;
    for (int k = 0; k < paddedLength; ++k) {
        std::cout << "Length[" << k << "] = " << std::abs(paddedInput[k]) << std::endl; // 计算复数的模
    }

    delete[] paddedInput;
}

// FFT逆变换测试
void IFFTTest()
{
    const int N = 8;  // 输入序列的长度
    std::complex<double> input[N] = {
        {28,0}, {-4,9.65685}, {-4,4}, {-4,1.65685},
        {-4,0}, {-4,-1.65685}, {-4,-4}, {-4,-9.65685}
    };  // 输入序列

    // 将输入序列扩展为2的幂次方长度
    int paddedLength = (int)pow(2, ceil(log2(N)));
    std::complex<double>* paddedInput = new std::complex<double>[paddedLength];
    for (int i = 0; i < N; ++i) {
        paddedInput[i] = input[i];
    }
    // 补零到2的n次幂
    for (int i = N; i < paddedLength; ++i) {
        paddedInput[i] = 0;
    }

    // 执行FFT逆变换
    IFFT(paddedInput, paddedLength);

    std::cout << "========== FFT逆变换 ========== " << std::endl;

    // 输出原始数据
    std::cout << "原始数据:----------" << std::endl;
    for (int k = 0; k < N; ++k) {
        std::cout << "F[" << k << "] = " << input[k] << std::endl;
    }

    // 输出变换结果
    std::cout << "变换结果:----------" << std::endl;
    for (int k = 0; k < paddedLength; ++k) {
        std::cout << "f[" << k << "] = " << paddedInput[k] << std::endl;
    }

    delete[] paddedInput;
}

int main() {
    
    // 执行FFT变换测试
    FFTTest();

    std::cout << std::endl;

    // 执行FFT逆变换测试
    IFFTTest();
    return 0;
}

 

输出结果如下：

========== FFT变换 ==========
原始数据:----------
f[0] = (0,0)
f[1] = (1,0)
f[2] = (2,0)
f[3] = (3,0)
f[4] = (4,0)
f[5] = (5,0)
f[6] = (6,0)
f[7] = (7,0)
变换结果:----------
F[0] = (28,0)
F[1] = (-4,9.65685)
F[2] = (-4,4)
F[3] = (-4,1.65685)
F[4] = (-4,0)
F[5] = (-4,-1.65685)
F[6] = (-4,-4)
F[7] = (-4,-9.65685)
幅值结果:----------
Length[0] = 28
Length[1] = 10.4525
Length[2] = 5.65685
Length[3] = 4.32957
Length[4] = 4
Length[5] = 4.32957
Length[6] = 5.65685
Length[7] = 10.4525

========== FFT逆变换 ==========
原始数据:----------
F[0] = (28,0)
F[1] = (-4,9.65685)
F[2] = (-4,4)
F[3] = (-4,1.65685)
F[4] = (-4,0)
F[5] = (-4,-1.65685)
F[6] = (-4,-4)
F[7] = (-4,-9.65685)
变换结果:----------
f[0] = (0,-0)
f[1] = (1,1.72255e-16)
f[2] = (2,4.44089e-16)
f[3] = (3,3.82857e-16)
f[4] = (4,-0)
f[5] = (5,-4.979e-17)
f[6] = (6,-4.44089e-16)
f[7] = (7,-5.05322e-16)

 

二维情况

二维FFT是在一维FFT基础上实现，实现过程为：

1.对二维输入数据的每一行进行FFT，变换结果仍然按行存入二维数组中。

2.在1的结果基础上，对每一列进行FFT，再存入原来数组，即得到二维FFT结果。

#include <iostream>
#include <complex>
#include <cmath>

const double PI = 3.14159265358979323846;

// 交换位置
void swap(std::complex<double>& a, std::complex<double>& b) {
    std::complex<double> temp = a;
    a = b;
    b = temp;
}

// 频率倒位
void bitReverse(std::complex<double>* point, int N) {
    int j = 0;
    for (int i = 1; i < N - 1; ++i) {
        int k = N / 2;
        while (j >= k) {
            j -= k;
            k /= 2;
        }
        j += k;
        if (i < j) {
            swap(point[i], point[j]);
        }
    }
}

// 一维Cooley-Tukey算法的非递归实现
void FFT(std::complex<double>* point, int N) {
    bitReverse(point, N);
    for (int s = 2; s <= N; s *= 2) {
        double angle = -2 * PI / s;
        std::complex<double> w(1, 0);
        std::complex<double> wn(cos(angle), sin(angle));
        for (int k = 0; k < N; k += s) {
            w = std::complex<double>(1, 0);
            for (int j = 0; j < s / 2; ++j) {
                std::complex<double> u = point[k + j];
                std::complex<double> v = point[k + j + s / 2] * w;
                point[k + j] = u + v;
                point[k + j + s / 2] = u - v;
                w *= wn;
            }
        }
    }
}

// 一维FFT逆变换
void IFFT(std::complex<double>* point, int N) {
    // 共轭变换(即：实部不变，虚部取相反数)
    for (int i = 0; i < N; ++i) {
        point[i] = std::conj(point[i]);
    }

    // 应用FFT
    FFT(point, N);

    // 共轭变换(即：实部不变，虚部取相反数)并归一化
    for (int i = 0; i < N; ++i) {
        point[i] = std::conj(point[i]);
        point[i] /= N;
    }
}

// 二维Cooley-Tukey算法的非递归实现
void FFT2D(std::complex<double>* point, int rows, int cols) {
    // 对每一行进行FFT变换
    for (int i = 0; i < rows; ++i) {
        FFT(&point[i * cols], cols);
    }

    // 对每一列进行FFT变换
    std::complex<double>* column = new std::complex<double>[rows];
    for (int j = 0; j < cols; ++j) {

        for (int i = 0; i < rows; ++i) {
            column[i] = point[i * cols + j];
        }
        FFT(column, rows);
        for (int i = 0; i < rows; ++i) {
            point[i * cols + j] = column[i];
        }
    }
    delete[] column;
}

// 二维FFT逆变换
void IFFT2D(std::complex<double>* point, int rows, int cols) {
    // 对每一行进行逆变换
    for (int i = 0; i < rows; ++i) {
        IFFT(&point[i * cols], cols);
    }

    // 对每一列进行逆变换
    std::complex<double>* column = new std::complex<double>[rows];
    for (int j = 0; j < cols; ++j) {

        for (int i = 0; i < rows; ++i) {
            column[i] = point[i * cols + j];
        }
        IFFT(column, rows);
        for (int i = 0; i < rows; ++i) {
            point[i * cols + j] = column[i];
        }
    }
    delete[] column;
}

// FFT变换测试
void FFTTest()
{
    const int rows = 4;  // 输入矩阵的行数
    const int cols = 4;  // 输入矩阵的列数
    std::complex<double> input[rows * cols] = {
        1, 2, 3, 4,
        5, 6, 7, 8,
        9, 10, 11, 12,
        13, 14, 15, 16
    };  // 输入矩阵

    // 将输入矩阵的行数和列数扩展为2的幂次方长度
    int paddedRows = (int)pow(2, ceil(log2(rows)));
    int paddedCols = (int)pow(2, ceil(log2(cols)));

    // 将长宽补零到2的n次幂
    std::complex<double>* paddedInput = new std::complex<double>[paddedRows * paddedCols];
    for (int i = 0; i < rows; ++i) {
        for (int j = 0; j < cols; ++j) {
            paddedInput[i * paddedCols + j] = input[i * cols + j];
        }
        for (int j = cols; j < paddedCols; ++j) {
            paddedInput[i * paddedCols + j] = 0;
        }
    }
    for (int i = rows; i < paddedRows; ++i) {
        for (int j = 0; j < paddedCols; ++j) {
            paddedInput[i * paddedCols + j] = 0;
        }
    }

    // 执行2D FFT变换
    FFT2D(paddedInput, paddedRows, paddedCols);

    std::cout << "========== FFT变换 ========== " << std::endl;

    // 输出原始数据
    std::cout << "原始数据:----------" << std::endl;
    for (int i = 0; i < rows; ++i) {
        for (int j = 0; j < cols; ++j) {
            std::cout << "f[" << i << "][" << j << "] = " << input[i * cols + j] << std::endl;
        }
    }

    // 输出变换结果
    std::cout << "变换结果:----------" << std::endl;
    for (int i = 0; i < paddedRows; ++i) {
        for (int j = 0; j < paddedCols; ++j) {
            std::cout << "F[" << i << "][" << j << "] = " << paddedInput[i * paddedCols + j] << std::endl;
        }
    }

    // 输出幅值结果
    std::cout << "幅值结果:----------" << std::endl;
    for (int i = 0; i < paddedRows; ++i) {
        for (int j = 0; j < paddedCols; ++j) {
            std::cout << "Length[" << i << "][" << j << "] = " << std::abs(paddedInput[i * paddedCols + j]) << std::endl; // 计算复数的模
        }
    }

    delete[] paddedInput;
}

// FFT逆变换测试
void IFFTTest()
{
    const int rows = 4;  // 输入矩阵的行数
    const int cols = 4;  // 输入矩阵的列数
    std::complex<double> input[rows * cols] = {
        {136,0}, {-8,8}, {-8,0}, {-8,-8},
        {-32,32}, {0,0}, {0,0}, {0,0},
        {-32,0}, {0,0}, {0,0}, {0,0},
        {-32,-32}, {0,0}, {0,0}, {0,0}
    };  // 输入矩阵

    // 将输入矩阵的行数和列数扩展为2的幂次方长度
    int paddedRows = (int)pow(2, ceil(log2(rows)));
    int paddedCols = (int)pow(2, ceil(log2(cols)));

    // 将长宽补零到2的n次幂
    std::complex<double>* paddedInput = new std::complex<double>[paddedRows * paddedCols];
    for (int i = 0; i < rows; ++i) {
        for (int j = 0; j < cols; ++j) {
            paddedInput[i * paddedCols + j] = input[i * cols + j];
        }
        for (int j = cols; j < paddedCols; ++j) {
            paddedInput[i * paddedCols + j] = 0;
        }
    }
    for (int i = rows; i < paddedRows; ++i) {
        for (int j = 0; j < paddedCols; ++j) {
            paddedInput[i * paddedCols + j] = 0;
        }
    }

    // 执行2D FFT变换
    IFFT2D(paddedInput, paddedRows, paddedCols);

    std::cout << "========== 执行FFT逆变换 ========== " << std::endl;

    // 输出原始数据
    std::cout << "原始数据:----------" << std::endl;
    for (int i = 0; i < rows; ++i) {
        for (int j = 0; j < cols; ++j) {
            std::cout << "F[" << i << "][" << j << "] = " << input[i * cols + j] << std::endl;
        }
    }

    // 输出变换结果
    std::cout << "变换结果:----------" << std::endl;
    for (int i = 0; i < paddedRows; ++i) {
        for (int j = 0; j < paddedCols; ++j) {
            std::cout << "f[" << i << "][" << j << "] = " << paddedInput[i * paddedCols + j] << std::endl;
        }
    }

    delete[] paddedInput;
}

int main() {

    // 执行FFT变换测试
    FFTTest();

    std::cout << std::endl;

    // 执行FFT逆变换测试
    IFFTTest();

    return 0;
}

 

输出结果如下：

========== FFT变换 ==========
原始数据:----------
f[0][0] = (1,0)
f[0][1] = (2,0)
f[0][2] = (3,0)
f[0][3] = (4,0)
f[1][0] = (5,0)
f[1][1] = (6,0)
f[1][2] = (7,0)
f[1][3] = (8,0)
f[2][0] = (9,0)
f[2][1] = (10,0)
f[2][2] = (11,0)
f[2][3] = (12,0)
f[3][0] = (13,0)
f[3][1] = (14,0)
f[3][2] = (15,0)
f[3][3] = (16,0)
变换结果:----------
F[0][0] = (136,0)
F[0][1] = (-8,8)
F[0][2] = (-8,0)
F[0][3] = (-8,-8)
F[1][0] = (-32,32)
F[1][1] = (0,0)
F[1][2] = (0,0)
F[1][3] = (0,0)
F[2][0] = (-32,0)
F[2][1] = (0,0)
F[2][2] = (0,0)
F[2][3] = (0,0)
F[3][0] = (-32,-32)
F[3][1] = (0,0)
F[3][2] = (0,0)
F[3][3] = (0,0)
幅值结果:----------
Length[0][0] = 136
Length[0][1] = 11.3137
Length[0][2] = 8
Length[0][3] = 11.3137
Length[1][0] = 45.2548
Length[1][1] = 0
Length[1][2] = 0
Length[1][3] = 0
Length[2][0] = 32
Length[2][1] = 0
Length[2][2] = 0
Length[2][3] = 0
Length[3][0] = 45.2548
Length[3][1] = 0
Length[3][2] = 0
Length[3][3] = 0

========== 执行FFT逆变换 ==========
原始数据:----------
F[0][0] = (136,0)
F[0][1] = (-8,8)
F[0][2] = (-8,0)
F[0][3] = (-8,-8)
F[1][0] = (-32,32)
F[1][1] = (0,0)
F[1][2] = (0,0)
F[1][3] = (0,0)
F[2][0] = (-32,0)
F[2][1] = (0,0)
F[2][2] = (0,0)
F[2][3] = (0,0)
F[3][0] = (-32,-32)
F[3][1] = (0,0)
F[3][2] = (0,0)
F[3][3] = (0,0)
变换结果:----------
f[0][0] = (1,-0)
f[0][1] = (2,6.12323e-17)
f[0][2] = (3,-0)
f[0][3] = (4,-6.12323e-17)
f[1][0] = (5,2.44929e-16)
f[1][1] = (6,3.06162e-16)
f[1][2] = (7,2.44929e-16)
f[1][3] = (8,1.83697e-16)
f[2][0] = (9,-0)
f[2][1] = (10,6.12323e-17)
f[2][2] = (11,-0)
f[2][3] = (12,-6.12323e-17)
f[3][0] = (13,-2.44929e-16)
f[3][1] = (14,-1.83697e-16)
f[3][2] = (15,-2.44929e-16)
f[3][3] = (16,-3.06162e-16)

 

扩展阅读

C++实现二维快速傅里叶变换（FFT）

 

标签：std,Cooley,point,int,FFT,complex,++,Tukey
From： https://www.cnblogs.com/kekec/p/18401750

相关文章

  • RFFT：数据与代码已开源，京东推出广告图生成新方法 | ECCV 2024
    论文将多模态可靠反馈网络（RFNet）结合到一个循环生成图片过程中，可以增加可用的广告图片数量。为了进一步提高生产效率，利用RFNet反馈进行创新的一致条件正则化，对扩散模型进行微调（RFFT），显著增加生成图片的可用率，减少了循环生成中的尝试次数，并提供了高效的生产过程，而不牺牲视觉吸引力。......
  • 如何编译FFTW3库：静态库与动态库的编译指南
    目录1.下载并解压FFTW3库2.配置编译选项3.编译并安装库4.验证编译结果5.在项目中使用FFTW3库6.总结FFTW3（FastestFourierTransformintheWest）是一个广泛使用的高性能傅里叶变换库。它支持多种优化，适用于多线程计算和SIMD指令，是处理大规模数据傅里叶变换的理......
  • 在 Ubuntu 中查找库的位置：以 FFTW3 库为例
    目录一、为什么要查找库的位置？二、查找库位置的常用方法三、实践示例四、总结在Ubuntu或其他基于Linux的操作系统中，开发者常常需要查找已安装库的位置，以便进行编译、链接或配置环境变量。本文将详细介绍如何在Ubuntu中查找库的位置，并以常用的FFTW3库为例进行......
  • 音频去噪：使用Python和FFT增强音质
    根据定义,声音去噪是从音频信号中去除不需要的噪音或干扰,以提高其质量和清晰度的过程。这涉及识别和隔离噪音成分(通常以不规则或高频元素为特征),并将其过滤掉,同时保持原始声音的完整性。声音去噪目标是改善聆听体验以及音频分析和处理的准确性。过滤掉噪音对于高保真音频来说......
  • 小学生都能看懂的FFT
    小学生都能看懂的FFT!!作者(小学生)其实学这个学了两个月，但我相信你只要努力，就能成功好的，废话不多说，正片开始FFT章节1:了解FFT是干嘛的oiwiki：FFT支持在O(n......
  • python 实现FFT快速傅立叶变换算法
    FFT快速傅里叶变换介绍FFT（快速傅里叶变换）是计算离散傅里叶变换（DFT）及其逆变换的一种高效算法。DFT是一种将信号从时域转换到频域的数学工具，而FFT通过减少计算量来加速这一过程。FFT的基本思想FFT利用了DFT中的对称性和周期性，通过分而治之的策略将DFT分解为更小的DFT，从而显......
  • Contest5385 - FFT
    ContestA签到题BFFT/NTT快速傅立叶P3803【模板】多项式乘法（FFT）&P1919【模板】高精度乘法|A*BProblem升级版参考资料：炫酷反演魔术-博客-vfleaking的博客题解P3803【【模板】多项式乘法（FFT）】-洛谷专栏FFT总体思路FFT处理循环卷积问题，而卷积问题通用的......
  • STM32H7 HAL库CubeMX 双重ADC模式同步采样详细配置+FFT计算相位差
    前言在电赛备赛期间琢磨了一下ADC同步采样的实现方式，本来是打算直接用AD7606来着，但是搞了半天也没把驱动整出来...考虑到AD7606本身采样率也拉不到太高，于是就花了几天时间把片上ADC配出来了。查资料的时候我发现关于STM32双重ADC模式的资料是真的少，用FFT算两路信号相位差的实例代......
  • MATLAB仿真：数字信号处理用FFT对信号分析
    目录1．实验目的2实验原理3．实验仪器及设备4．实验步骤及内容（1）对以下序列进行谱分析。（2）对以下周期序列进行谱分析。 （3）对模拟周期信号进行谱分析1．实验目的学习用FFT 对连续信号和时域离散信号进行谱分析的方法，了解可能出现的分析误差及其原因，以便正确应用 FFT。2......
  • scipy.fft (Python) 结果和 FFTW (C) 结果之间的微小差异
    我正在尝试使用C中的FFTW从Python中的一些已知工作代码重新创建结果。我发现结果中有一些小错误。scipy.fft我的输入数据是真实的3d，尺寸=(294,294,294)。我的scipy.fft调用如下所示：我的fftw代码如下所示这个：complex_data_out=scipy.fft.fftn......

赞助商

阅读排行

网站主页关于我们联系我们网站地图

本网站内容转载自其他媒体，侵权联系[admin##ips99.com]。

Copyright © 2020-2023 IPS99 版权所有 IPS99