快速傅立叶变换（FFT）的C++实现

最新推荐文章于 2024-05-13 02:48:17 发布
最新推荐文章于 2024-05-13 02:48:17 发布
阅读量6.3k 收藏 12
点赞数
分类专栏： 乐程 文章标签： fft c++ float output permutation 算法

标准的离散傅立叶 DFT 变换形式如：

yk=Σj=0n-1 ajωn-kj  = A (ωn-k).

（ωnk 为复数 1 的第 k 个 n 次方根，且定义多项式 A (x) = Σj=0n-1 ajxj ）

而离散傅立叶逆变换 IDFT (Inverse DFT)形式如： aj=(Σk=0n-1 ykωnkj)/n .

以下不同颜色内容为引用并加以修正：

快速傅立叶变换（Fast Fourier Transform，FFT）是离散傅立叶变换（Discrete Fourier transform，DFT）的快速算法，它是根据离散傅立叶变换的奇、偶、虚、实等特性，对离散傅立叶变换的算法进行改进获得的。它对傅立叶变换的理论并没有新的发现，但是对于在计算机系统或者说数字系统中应用离散傅立叶变换，可以说是进了一大步。
设 Xn 为 N 项的复数序列，由 DFT 变换，任一 Xi 的计算都需要 N 次复数乘法和 N -1 次复数加法，而一次复数乘法等于四次实数乘法和两次实数加法，一次复数加法等于两次实数加法，即使把一次复数乘法和一次复数加法定义成一次“运算”（四次实数乘法和四次实数加法），那么求出 N 项复数序列的 Xi ，即 N 点 DFT 变换大约就需要 N2 次运算。当 N =1024 点甚至更多的时候，需要 N2 = 1048576 次运算，在 FFT 中，利用 ωn 的周期性和对称性，把一个 N 项序列（设 N 为偶数），分为两个 N / 2 项的子序列，每个 N / 2点 DFT 变换需要 (N / 2)2 次运算，再用 N 次运算把两个 N / 2点的 DFT 变换组合成一个 N 点的 DFT 变换。这样变换以后，总的运算次数就变成 N + 2 * (N / 2)2 = N + N2 / 2。继续上面的例子，N =1024 时，总的运算次数就变成了525312 次，节省了大约   50% 的运算量。而如果我们将这种“一分为二”的思想不断进行下去，直到分成两两一组的 DFT 运算单元，那么N 点的 DFT 变换就只需要 N * log2N 次的运算，N = 1024 点时，运算量仅有 10240 次，是先前的直接算法的1% ，点数越多，运算量的节约就越大，这就是 FFT 的优越性。

FFT 的实现可以自顶而下，采用递归，但是对于硬件实现成本高，对于软件实现都不够高效，改用迭代较好，自底而上地解决问题。感觉和归并排序的迭代版很类似，不过先要采用“位反转置换”的方法把 Xi 放到合适的位置，设 i 和 j 互为s = log2N 位二进制的回文数，假设 s = 3,   i = (110)2   = 6,   j = (011)2 = 3,   如果 i ≠ j , 那么 Xi 和 Xj 应该互换位置。（关于这个回文数的生成，是很有趣而且是很基本的操作，想当初偶初学 C++ 的时候就有这样的习题。）当“位反转置换”完成后，先将每一个Xi 看作是独立的多项式，然后两个两个地将它们合并成一个多项式（每个多项式有 2 项），合并实际上是“蝶形运算”（Butterfly Operation, 参考《算法导论》吧^_^），继续合并（第二次的每个多项式有 4 项），直到只剩下一个多项式（有 N 项），这样，合并的层数就是 log2N ，每层都有 N 次操作，所以总共有 N * log2N 次操作。迭代过程如下图所示，自底而上。

 

C++ 源程序，如下：

//
// 快速傅立叶变换 Fast Fourier Transform
// By rappizit@yahoo.com.cn
// 2007-07-20
// 版本 2.0
// 改进了《算法导论》的算法，旋转因子取 ωn-kj   (ωnkj 的共轭复数)
// 且只计算 n / 2 次，而未改进前需要计算 (n * lg n) / 2 次。
//

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

const int N = 1024;
const float PI = 3.1416;

inline void swap (float &a, float &b)
{
     float t;
     t = a;
     a = b;
     b = t;
}

void bitrp (float xreal [], float ximag [], int n)
{
     // 位反转置换 Bit-reversal Permutation
     int i, j, a, b, p;

     for (i = 1, p = 0; i < n; i *= 2)
{
         p ++;
}
     for (i = 0; i < n; i ++)
{
         a = i;
         b = 0;
         for (j = 0; j < p; j ++)
   {
             b = (b << 1) + (a & 1);     // b = b * 2 + a % 2;
             a >>= 1;         // a = a / 2;
   }
         if ( b > i)
   {
             swap (xreal [i], xreal [b]);
             swap (ximag [i], ximag [b]);
   }
}
}

void FFT(float xreal [], float ximag [], int n)
{
     // 快速傅立叶变换，将复数 x 变换后仍保存在 x 中，xreal, ximag 分别是 x 的实部和虚部
     float wreal [N / 2], wimag [N / 2], treal, timag, ureal, uimag, arg;
     int m, k, j, t, index1, index2;
    
     bitrp (xreal, ximag, n);

     // 计算 1 的前 n / 2 个 n 次方根的共轭复数 W'j = wreal [j] + i * wimag [j] , j = 0, 1,   , n / 2 - 1
     arg = - 2 * PI / n;
     treal = cos (arg);
     timag = sin (arg);
     wreal [0] = 1.0;
     wimag [0] = 0.0;
     for (j = 1; j < n / 2; j ++)
{
         wreal [j] = wreal [j - 1] * treal - wimag [j - 1] * timag;
         wimag [j] = wreal [j - 1] * timag + wimag [j - 1] * treal;
}

     for (m = 2; m <= n; m *= 2)
{
         for (k = 0; k < n; k += m)
   {
             for (j = 0; j < m / 2; j ++)
    {
                 index1 = k + j;
                 index2 = index1 + m / 2;
                 t = n * j / m;     // 旋转因子 w 的实部在 wreal [] 中的下标为 t
                 treal = wreal [t] * xreal [index2] - wimag [t] * ximag [index2];
                 timag = wreal [t] * ximag [index2] + wimag [t] * xreal [index2];
                 ureal = xreal [index1];
                 uimag = ximag [index1];
                 xreal [index1] = ureal + treal;
                 ximag [index1] = uimag + timag;
                 xreal [index2] = ureal - treal;
                 ximag [index2] = uimag - timag;
    }
   }
}
}

void   IFFT (float xreal [], float ximag [], int n)
{
     // 快速傅立叶逆变换
     float wreal [N / 2], wimag [N / 2], treal, timag, ureal, uimag, arg;
     int m, k, j, t, index1, index2;
    
     bitrp (xreal, ximag, n);

     // 计算 1 的前 n / 2 个 n 次方根 Wj = wreal [j] + i * wimag [j] , j = 0, 1,   , n / 2 - 1
     arg = 2 * PI / n;
     treal = cos (arg);
     timag = sin (arg);
     wreal [0] = 1.0;
     wimag [0] = 0.0;
     for (j = 1; j < n / 2; j ++)
{
         wreal [j] = wreal [j - 1] * treal - wimag [j - 1] * timag;
         wimag [j] = wreal [j - 1] * timag + wimag [j - 1] * treal;
}

     for (m = 2; m <= n; m *= 2)
{
         for (k = 0; k < n; k += m)
   {
             for (j = 0; j < m / 2; j ++)
    {
                 index1 = k + j;
                 index2 = index1 + m / 2;
                 t = n * j / m;     // 旋转因子 w 的实部在 wreal [] 中的下标为 t
                 treal = wreal [t] * xreal [index2] - wimag [t] * ximag [index2];
                 timag = wreal [t] * ximag [index2] + wimag [t] * xreal [index2];
                 ureal = xreal [index1];
                 uimag = ximag [index1];
                 xreal [index1] = ureal + treal;
                 ximag [index1] = uimag + timag;
                 xreal [index2] = ureal - treal;
                 ximag [index2] = uimag - timag;
    }
   }
}

     for (j=0; j < n; j ++)
{
         xreal [j] /= n;
         ximag [j] /= n;
}
}

void FFT_test ()
{
     char inputfile [] = "input.txt";     // 从文件 input.txt 中读入原始数据
     char outputfile [] = "output.txt";     // 将结果输出到文件 output.txt 中
     float xreal [N] = {}, ximag [N] = {};
     int n, i;
     FILE *input, *output;

     if (!(input = fopen (inputfile, "r")))
{
         printf ("Cannot open file. ");
         exit (1);
}
     if (!(output = fopen (outputfile, "w")))
{
         printf ("Cannot open file. ");
         exit (1);
}

     i = 0;
     while ((fscanf (input, "%f%f", xreal + i, ximag + i)) != EOF)
{
         i ++;
}
     n = i;     // 要求 n 为 2 的整数幂
     while (i > 1)
{
         if (i % 2)
   {
             fprintf (output, "%d is not a power of 2! ", n);
             exit (1);
   }
         i /= 2;
}

     FFT (xreal, ximag, n);
     fprintf (output, "FFT:     i      real imag ");
     for (i = 0; i < n; i ++)
{
         fprintf (output, "%4d     %8.4f     %8.4f ", i, xreal [i], ximag [i]);
}
     fprintf (output, "================================= ");

     IFFT (xreal, ximag, n);
     fprintf (output, "IFFT:     i      real imag ");
     for (i = 0; i < n; i ++)
{
         fprintf (output, "%4d     %8.4f     %8.4f ", i, xreal [i], ximag [i]);
}

     if ( fclose (input))
{
         printf ("File close error. ");
         exit (1);
}
     if ( fclose (output))
{
         printf ("File close error. ");
         exit (1);
}
}

int main ()
{
     FFT_test ();

     return 0;
}


// sample: input.txt


     0.5751          0
     0.4514          0
     0.0439          0
     0.0272          0
     0.3127          0
     0.0129          0
     0.3840          0
     0.6831          0
     0.0928          0
     0.0353          0
     0.6124          0
     0.6085          0
     0.0158          0
     0.0164          0
     0.1901          0
     0.5869          0


// sample: output.txt

FFT:
    i            real          imag
    0       4.6485       0.0000
    1       0.0176       0.3122
    2       1.1114       0.0429
    3       1.6776      -0.1353
    4      -0.2340       1.3897
    5       0.3652      -1.2589
    6      -0.4325       0.2073
    7      -0.1312       0.3763
    8      -0.1949       0.0000
    9      -0.1312      -0.3763
   10      -0.4326      -0.2073
   11       0.3652       1.2589
   12      -0.2340      -1.3897
   13       1.6776       0.1353
   14       1.1113      -0.0429
   15       0.0176      -0.3122
=================================
IFFT:
    i            real          imag
    0       0.5751      -0.0000
    1       0.4514       0.0000
    2       0.0439      -0.0000
    3       0.0272      -0.0000
    4       0.3127      -0.0000
    5       0.0129      -0.0000
    6       0.3840      -0.0000
    7       0.6831       0.0000
    8       0.0928       0.0000
    9       0.0353      -0.0000
   10       0.6124       0.0000
   11       0.6085       0.0000
   12       0.0158       0.0000
   13       0.0164       0.0000
   14       0.1901       0.0000
   15       0.5869      -0.0000


**************************************************/

lishizelibin
关注
  • 0
    点赞
  • 踩
  • 12
    收藏
    觉得还不错? 一键收藏
  • 0
    评论
专栏目录
快速傅立叶变换：FFT算法原理及C++实现
白马负金羁
01-07 2万+
傅立叶变换以高等数学（微积分）中的傅立叶级数为基础发展而来，它是信号处理（特别是图像处理）中非常重要的一种时频变换手段，具有重要应用。在图像编码、压缩、降噪、数字水印方面都有重要意义。此外，快速傅立叶变换算法还位列20世纪十大算法之列，它是“动态规划”策略在算法设计中的杰出代表
C++实现fft2函数功能（附完整源码）
希望我的博客，能帮上你解决学习中工作中所遇到的问题
06-26 250
C++实现fft2函数功能（附完整源码）
C++实现快速傅里叶变换
庐州牧
05-07 4227
使用C++实现快速傅里叶变换代码，使用Qt展示测试信号的变换结果。
【FFT实战篇】C++实现：利用快速傅里叶变换快速计算（多项式）乘法
GoWilli的CSDN
03-22 2410
快速傅里叶变换FFT全流程解析，洛谷P3803 FFT多项式乘法
2024年C C++最全快速傅立叶变换（FFT）的C++实现与Matlab实验_c+，从零开始学数据结构和算法
最新发布
2401_84975937的博客
05-13 440
/ 且只计算 n / 2 次，而未改进前需要计算 (n * lg n) / 2 次。// 改进了《算法导论》的算法，旋转因子取 ωn-kj (ωnkj 的共轭复数)// 位反转置换 Bit-reversal Permutation。// 快速傅立叶变换 Fast Fourier Transform。
图像的傅里叶变换和逆变换C++版
Dezeming的博客
01-04 6482
首先把整体程序分段，每段都是单独的功能。最后再把完整代码给出： 头文件：包含OpenCV和标准库输出调试 #include<opencv2/core/core.hpp> #include<opencv2/highgui/highgui.hpp> #include<opencv2/imgproc/imgproc.hpp> using namespace cv; #include <iostream> using std::cout; using std::
FFT.rar_快速傅立叶变换c++实现
09-20
总之，“FFT.rar_快速傅立叶变换c++实现”提供了一种C++实现的快速傅立叶变换算法，对于理解和应用FFT具有很高的价值。开发者可以通过阅读和理解源代码，进一步学习FFT的工作原理，将其应用于自己的项目中，提高计算...
一维快速傅里叶变换FFT的C++实现
03-12
一维快速傅里叶变换FFT的C++实现，里面是FFT1.cpp函数，用于进行一维数组的FFT。有详细的注释和说明。
快速傅里叶变换(FFT)算法C++实现代码.docx
10-25
快速傅里叶变换（FFT）算法C++实现代码 快速傅里叶变换（FFT）算法是傅里叶变换的一种快速算法，其中傅里叶变换是将时域信号转换为频域信号的过程。快速傅里叶变换算法可以将信号从时域转换到频域，或者从频域转换...
快速傅里叶变换的C++实现
11-08
快速傅里叶变换 (fast Fourier transform), 即利用计算机计算离散傅里叶变换（DFT)的高效、快速计算方法的统称，简称FFT。快速傅里叶变换是1965年由J.W.库利和T.W.图基提出的。采用这种算法能使计算机计算离散傅里叶...
fft.zip_fft c源代码_快速傅立叶变换 c++
09-19
快速傅立叶变换（FFT）是一种在数字信号处理领域中广泛应用的算法，它极大地提高了计算离散傅立叶变换（DFT）的效率。在给定的“fft.zip”压缩包中，包含了一个名为“fft.txt”的文件，这很可能是用C语言编写的快速...
快速傅里叶变换FFT算法完整实现C++代码
03-10
快速傅里叶变换FFT算法完整实现C++代码 附带详细注释和数十组测试数据
FFT快速傅里叶变换的C++程序类
07-01
FFT快速傅里叶变换的C++程序,已封装成类,里面有详细的使用说明,特别适合C++Builder等标准编译器使用，还准备了一个小小的PPT(是别人做了ppt,我在网上下后做了修改,适合自学者编写FFT程序时参考)
自己编的FFT以及FFT2函数，比较慢……
04-18
自己编的FFT以及FFT2函数，从复数定义以及复数的运算开始写的，比较慢……
基于c++的快速傅里叶变换
03-26
C++实现快速傅里叶变换(FFT)算法，附带源码
快速傅里叶变换(FFT)算法C++实现代码
简单代码的专栏
06-25 1万+
#include #define DOUBLE_PI   6.283185307179586476925286766559// 快速傅里叶变换// data 长度为 (2 * 2^n), data 的偶位为实数部分, data 的奇位为虚数部分// isInverse表示是否为逆变换void FFT(double * data, int n, bool isInverse = false){ 
OpenCV之傅里叶变换（C++实现）
热门推荐
xddwz的博客
12-09 1万+
一、傅立叶变换 离散傅立叶变换(Discrete Fourier Transform, DFT)，是指傅立叶变换在时域和频域上都呈现离散的形式，将时域信号的采样变换为在离散时间傅立叶变换频域的采样。 简单来说，对一张图像进行离散傅立叶变换就是将它分解成正弦和余弦两部分，也就是将图像从空间域转换到频域。这一转换的理论基础是：任一函数都可以表示成无数个正弦和余弦函数的和的形式。 二维图像的傅立叶变化用数学公式可以表示为： 式中的f是空间域值，F是频域值。变换之后的频域值是复数，因此，显示傅立叶变换之
C++实现快速傅里叶变换（FFT）
baidu_41901258的博客
03-23 1万+
由于项目要从MATLAB搬到VS上，开始认真研究怎么在C++中实现FFT，更准确的来说是DFT和IDFT。 FFT的公式人人都知道，但要在计算机上实现需要一点小技巧：蝶形运算。 蝶形运算的原理就是通过递归缩小进行DFT的序列长度。根据Wn因子的特点，我们需要把序列分为奇数下标序列和偶数下标序列，经过了不停的拆分和重新组合，对得到的仅有的2点进行DFT，不断扩展序列长度、递归计算得到最后的结果...
FFT c++语言实现
08-19 2983
/*  * test.cpp  *  *  Created on: 2013-8-18  *      Author: zhijian  */ #include   #include   #define N 128 //N 点 #define RATE 128 //抽样频率 #define LEN 7 //1 #define PI 3.14159265
c++实现快速傅里叶变换（FFT）
05-29
快速傅里叶变换（FFT）是一种高效的计算离散傅里叶变换（DFT）的算法。下面是使用 C++ 实现 FFT 的示例代码： ```c++ #include <iostream> #include <cmath> #include <complex> using namespace std; const double PI = acos(-1.0); void FFT(complex<double>* a, int n, int inv) { if (n == 1) return; complex<double> a0[n / 2], a1[n / 2]; for (int i = 0; i < n / 2; i++) { a0[i] = a[i * 2]; a1[i] = a[i * 2 + 1]; } FFT(a0, n / 2, inv); FFT(a1, n / 2, inv); complex<double> w = 1, wn = exp(complex<double>(0, inv * 2 * PI / n)); for (int i = 0; i < n / 2; i++) { a[i] = a0[i] + w * a1[i]; a[i + n / 2] = a0[i] - w * a1[i]; w *= wn; } } int main() { int n = 8; complex<double> a[n] = {1, 2, 3, 4, 5, 6, 7, 8}; FFT(a, n, 1); for (int i = 0; i < n; i++) { cout << a[i] << endl; } FFT(a, n, -1); for (int i = 0; i < n; i++) { cout << a[i] / n << endl; } return 0; } ``` 上述代码中，FFT 函数实现了对复数数组 a 的 DFT 或 IDFT 计算。其中，inv 参数表示计算 DFT 还是 IDFT，1 表示 DFT，-1 表示 IDFT。在 FFT 函数内部，首先进行递归，将 a 数组划分为两个长度为 n/2 的数组，然后对这两个数组分别进行 FFT 计算。最后，按照蝴蝶操作的方式将两个数组合并为一个长度为 n 的数组，完成 DFT 或 IDFT 计算。在 main 函数中，我们构造了一个长度为 8 的复数数组并对其进行了 DFT 和 IDFT 计算，并输出了结果。

“相关推荐”对你有帮助么？

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

当前余额3.43元 前往充值 >
需支付：10.00元
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付元
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明：

1.余额是钱包充值的虚拟货币，按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载，可以购买VIP、付费专栏及课程。

余额充值