上一篇文章里，我根据DFT公式用C++实现了二维离散傅里叶变换。但跑一张300*300的图片都要好几分钟，速度实在太慢了。我研究了下快速傅里叶变换，在网上找了一些资料，然后用C++实现了二维快速傅里叶变换，速度超级快，512*512的图片在0.1秒内就能处理完。

FFT原理

我们知道一维离散傅里叶变换公式如下

其中。

用矩阵表示即为：

这样就有N*N次乘法计算，很耗时。可以通过一些方式减小计算量。

由于是个周期函数，有如下两个性质

这样

从上面两个式子可以看出，

我们只要分别计算奇数项的DFT值

和偶数项的DFT值，

然后再通过一些加减，得到整体DFT的值X(k),乘法计算量大致就是(N/2)*(N/2)+N/2)*(N/2)=N*N/2。是之前的一半。

依次类推，奇数项和偶数项计算DFT值，也可以通过上面的方法，这样乘法计算量就每次除以2的降低。只要N满足是2的n次方，减到最后，计算量就减少2的n次方倍。

例如以下8个点计算FFT

FFT方法一：

分析可知，可以通过递归实现这种算法。

C++代码如下

void fft(vector<Cpx>& a, int lim, int opt)
{if (lim == 1) return;vector<Cpx> a0(lim >> 1), a1(lim >> 1); // 初始化一半大小，存放偶数和奇数部分for (int i = 0; i < lim; i += 2)a0[i >> 1] = a[i], a1[i >> 1] = a[i + 1]; // 分成偶数部分和奇数部分fft(a0, lim >> 1, opt); // 递归计算偶数部分fft(a1, lim >> 1, opt); // 递归计算偶数部分Cpx wn(cos(2 * PI / lim), opt * -sin(2 * PI / lim)); //等于WNCpx w(1, 0);for (int k = 0; k < (lim >> 1); k++) // 见蝶形图1运算过程{a[k] = a0[k] + w * a1[k];a[k + (lim >> 1)] = a0[k] - w * a1[k];w = w * wn;}//for (int k = 0; k < (lim >> 1); k++) // 见蝶形图2，小优化一下，少一次乘法//{//   Cpx t = w * a1[k];//   a[k] = a0[k] + t;//   a[k + (lim >> 1)] = a0[k] - t;//    w = w * wn;//}}

FFT方法二：

当然，再次观察这个蝶形图，也可以有另外一种常规算法

从蝶形图可以看出，从x到X，可以写成3次循环实现。当然关于x的排序，可以推理出为二进制逆序。

//二进制逆序排列
int ReverseBin(int a, int n)
{int ret = 0;for (int i = 0; i < n; i++){if (a&(1 << i)) ret |= (1 << (n - 1 - i));}return ret;
}void fft2(vector<Cpx>& a, int lim, int opt)
{int index;vector<Cpx> tempA(lim);for (int i = 0; i < lim; i++){index = ReverseBin(i, log2(lim));tempA[i] = a[index];}vector<Cpx> WN(lim / 2);//生成WN表,避免重复计算for (int i = 0; i < lim / 2; i++){WN[i] = Cpx(cos(2 * PI *i / lim), opt * -sin(2 * PI*i / lim));}//蝶形运算int Index0, Index1;Cpx temp;for (int steplenght = 2; steplenght <= lim; steplenght *= 2){for (int step = 0; step < lim / steplenght; step++){for (int i = 0; i < steplenght / 2; i++){Index0 = steplenght * step + i;Index1 = steplenght * step + i + steplenght / 2;temp = tempA[Index1] * WN[lim / steplenght * i];tempA[Index1] = tempA[Index0] - temp;tempA[Index0] = tempA[Index0] + temp;}}}for (int i = 0; i < lim; i++){if (opt == -1){a[i] = tempA[i] / lim;}else{a[i] = tempA[i];}}
}

二维FFT原理

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

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

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

例子

我这边用一张图做了变换处理，

原图

二维FFT

逆变换

C++实现二维FFT代码

#include <iostream>
#include<opencv2/opencv.hpp>using namespace cv;
using namespace std;
const int height = 512, width = 512;const double PI = acos(-1); // pi值struct Cpx // 定义一个复数结构体和复数运算法则
{double r, i;Cpx() : r(0), i(0) {}Cpx(double _r, double _i) : r(_r), i(_i) {}
};
Cpx operator + (Cpx a, Cpx b) { return Cpx(a.r + b.r, a.i + b.i); }
Cpx operator - (Cpx a, Cpx b) { return Cpx(a.r - b.r, a.i - b.i); }
Cpx operator * (Cpx a, Cpx b) { return Cpx(a.r * b.r - a.i * b.i, a.r * b.i + a.i * b.r); }
Cpx operator / (Cpx a, int b) { return Cpx(a.r*1.0 / b, a.i*1.0 / b); }Mat BGR2GRAY(Mat img)
{int w = img.cols;int h = img.rows;Mat grayImg(h, w, CV_8UC1);uchar *p = grayImg.ptr<uchar>(0);Vec3b *pImg = img.ptr<Vec3b>(0);for (int i = 0; i < w*h; ++i){p[i] = 0.2126*pImg[i][2] + 0.7152*pImg[i][1] + 0.0722*pImg[i][0];}return grayImg;
}Mat Resize(Mat img)
{int w = img.cols;int h = img.rows;Mat out(height, width, CV_8UC1);uchar *p = out.ptr<uchar>(0);uchar *p2 = img.ptr<uchar>(0);int x_before, y_before;for (int y = 0; y < height; y++) {for (int x = 0; x < width; x++){x_before = (int)x*w*1.0 / width;y_before = (int)y*h*1.0 / height;p[y * width + x] = p2[y_before * w + x_before];}}return out;}//自动缩放到0-255范围内,并变换象限，将低频移至中间
int AutoScale(Mat src, Mat out)
{int w = src.cols;int h = src.rows;double *p = src.ptr<double>(0);uchar *pOut = out.ptr<uchar>(0);double max = p[0];double min = p[0];for (int i = 0; i < w*h; i++){if (p[i] > max) max = p[i];if (p[i] < min) min = p[i];}double scale = 255.0 / (max - min);for (int i = 0; i < w*h; i++){int j = i + w * h / 2 + w / 2;if (j > w*h) j = j - w * h;   //低频移至中间pOut[i] = (uchar)((p[j] - min)*scale);}return 0;
}void fft(vector<Cpx>& a, int lim, int opt)
{if (lim == 1) return;vector<Cpx> a0(lim >> 1), a1(lim >> 1); // 初始化一半大小，存放偶数和奇数部分for (int i = 0; i < lim; i += 2)a0[i >> 1] = a[i], a1[i >> 1] = a[i + 1]; // 分成偶数部分和奇数部分fft(a0, lim >> 1, opt); // 递归计算偶数部分fft(a1, lim >> 1, opt); // 递归计算偶数部分Cpx wn(cos(2 * PI / lim), opt * -sin(2 * PI / lim)); //等于WNCpx w(1, 0);for (int k = 0; k < (lim >> 1); k++) // 见蝶形图1运算过程{a[k] = a0[k] + w * a1[k];a[k + (lim >> 1)] = a0[k] - w * a1[k];w = w * wn;}//for (int k = 0; k < (lim >> 1); k++) // 见蝶形图2，小优化一下，少一次乘法//{//   Cpx t = w * a1[k];//   a[k] = a0[k] + t;//   a[k + (lim >> 1)] = a0[k] - t;//    w = w * wn;//}}//二进制逆序排列
int ReverseBin(int a, int n)
{int ret = 0;for (int i = 0; i < n; i++){if (a&(1 << i)) ret |= (1 << (n - 1 - i));}return ret;
}void fft2(vector<Cpx>& a, int lim, int opt)
{int index;vector<Cpx> tempA(lim);for (int i = 0; i < lim; i++){index = ReverseBin(i, log2(lim));tempA[i] = a[index];}vector<Cpx> WN(lim / 2);//生成WN表,避免重复计算for (int i = 0; i < lim / 2; i++){WN[i] = Cpx(cos(2 * PI *i / lim), opt * -sin(2 * PI*i / lim));}//蝶形运算int Index0, Index1;Cpx temp;for (int steplenght = 2; steplenght <= lim; steplenght *= 2){for (int step = 0; step < lim / steplenght; step++){for (int i = 0; i < steplenght / 2; i++){Index0 = steplenght * step + i;Index1 = steplenght * step + i + steplenght / 2;temp = tempA[Index1] * WN[lim / steplenght * i];tempA[Index1] = tempA[Index0] - temp;tempA[Index0] = tempA[Index0] + temp;}}}for (int i = 0; i < lim; i++){if (opt == -1){a[i] = tempA[i] / lim;}else{a[i] = tempA[i];}}
}void FFT2D(Cpx(*src)[width], Cpx(*dst)[width], int opt)
{//第一遍fftfor (int i = 0; i < height; i++){vector<Cpx> tempData(width);//获取每行数据for (int j = 0; j < width; j++){tempData[j] = src[i][j];}//一维FFTfft2(tempData, width, opt);//写入每行数据for (int j = 0; j < width; j++){dst[i][j] = tempData[j];}}//第二遍fftfor (int i = 0; i < width; i++){vector<Cpx> tempData(height);//获取每列数据for (int j = 0; j < height; j++){tempData[j] = dst[j][i];}//一维FFTfft2(tempData, height, opt);//写入每列数据for (int j = 0; j < height; j++){dst[j][i] = tempData[j];}}
}void Mat2Cpx(Mat src, Cpx(*dst)[width])
{//这里Mat里的数据得是unchar类型uchar *p = src.ptr<uchar>(0);for (int i = 0; i < height; i++){for (int j = 0; j < width; j++){dst[i][j] = Cpx(p[i*width + j],0);}}
}void Cpx2Mat(Cpx(*src)[width], Mat dst)
{//这里Mat里的数据得是unchar类型uchar *p = dst.ptr<uchar>(0);for (int i = 0; i < height; i++){for (int j = 0; j < width; j++){double g = sqrt(src[i][j].r*src[i][j].r);p[j + i * width] = (uchar)g;}}
}void Cpx2MatDouble(Cpx(*src)[width], Mat dst)
{//这里Mat里的数据得是double类型double *p = dst.ptr<double>(0);for (int i = 0; i < height; i++){for (int j = 0; j < width; j++){double g = sqrt(src[i][j].r*src[i][j].r+ src[i][j].i*src[i][j].i);g = log(g + 1);  //转换为对数尺度p[j + i * width] = (double)g;}}
}int main()
{Mat img = imread("d:\\1.jpg");imshow("img", img);Mat gray = BGR2GRAY(img);imshow("gray", gray);imwrite("gray.jpg", gray);Mat imgRez = Resize(gray);imshow("imgRez", imgRez);imwrite("imgRez.jpg", imgRez);Cpx(*src)[width] = new Cpx[height][width];Mat2Cpx(imgRez, src);Cpx(*dst)[width] = new Cpx[height][width];double t1 = getTickCount();FFT2D(src, dst, 1);double t2 = getTickCount();double t1t2 = (t2 - t1) / getTickFrequency();std::cout << "DFT耗时: " << t1t2 << "秒" << std::endl;Cpx(*dst2)[width] = new Cpx[height][width];double t3 = getTickCount();FFT2D(dst, dst2, -1);double t4 = getTickCount();double t3t4 = (t4 - t3) / getTickFrequency();std::cout << "DFT耗时: " << t3t4 << "秒" << std::endl;Mat out2 = Mat::zeros(height, width, CV_8UC1);Cpx2Mat(dst2, out2);imshow("out2", out2);imwrite("out2.jpg", out2);Mat out = Mat::zeros(height, width, CV_64F);Cpx2MatDouble(dst, out);Mat out3 = Mat::zeros(height, width, CV_8UC1);AutoScale(out, out3);imshow("out3", out3);imwrite("out3.jpg", out3);waitKey(0);
}

OK！

参考资料：FFT原理及C++与MATLAB混合编程详细介绍 - 羽扇纶巾o0 - 博客园

二维离散傅里叶（DFT）与快速傅里叶（FFT）的实现_asmartplum的博客-CSDN博客_二维fft

FFT与二维FFT的C#实现 - 死猫 - 博客园

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

快速傅里叶变换 java_二维快速傅里叶变换的java实现
图像处理与模式识别课作业,没学过信号与系统(哭晕). 恶补了两天冈萨里斯的书,写一下实现原理及过程看了网络上很多版本的概念讲解,道理我都懂,但是在将算法迁移到傅里叶变换的实现上时,出现了一些问题.接 ...
傅里叶变换一维和二维快速傅里叶变换（代码和性能的优化）
1.介绍. 在类FourierUtils的fftProgress方法中,有这个代码段,我们可以将Complext.euler(flag * i)提前计算好,设置大小为2次幂N,如果没有的话,也要调节到 ...
matlab fft2怎么移动频率对称,fft2 二维快速傅里叶变换（Matlab）
1.语法: Y = fft2(X) Y = fft2(X,m,n) 2.说明: Y = fft2(X) 使用快速傅里叶变换算法返回矩阵的二维傅里叶变换,这等同于计算 fft(fft(X).').'.如 ...
fft2 二维快速傅里叶变换（Matlab）
1.语法: Y = fft2(X) Y = fft2(X,m,n) 2.说明: Y = fft2(X) 使用快速傅里叶变换算法返回矩阵的二维傅里叶变换,这等同于计算 fft(fft(X ...
C语言二维数组范德蒙,浅谈范德蒙德(Vandermonde)方阵的逆矩阵的求法以及快速傅里叶变换(FFT)中IDFT的原理...
浅谈范德蒙德(Vandermonde)方阵的逆矩阵与拉格朗日(Lagrange)插值的关系以及快速傅里叶变换(FFT)中IDFT的原理标签: 行列式矩阵线性代数 FFT 拉格朗日插值只要稍微看 ...
基于python的快速傅里叶变换FFT（二）
基于python的快速傅里叶变换FFT(二) 本文在上一篇博客的基础上进一步探究正弦函数及其FFT变换. 知识点 FFT变换,其实就是快速离散傅里叶变换,傅立叶变换是数字信号处理领域一种很重要的算 ...
快速傅里叶变换python_基于python的快速傅里叶变换FFT（二）
基于python的快速傅里叶变换FFT(二) 本文在上一篇博客的基础上进一步探究正弦函数及其FFT变换. 知识点 FFT变换,其实就是快速离散傅里叶变换,傅立叶变换是数字信号处理领域一种很重要的算法. ...
【OpenCV 例程200篇】73. 二维连续傅里叶变换
[OpenCV 例程200篇]73. 二维连续傅里叶变换欢迎关注『OpenCV 例程200篇』系列,持续更新中欢迎关注『Python小白的OpenCV学习课』系列,持续更新中 2.1 二维 ...
离散傅里叶变换 (DFT)、快速傅里叶变换 (FFT)
目录离散傅里叶变换 (DFT) 离散傅里叶变换的基离散傅里叶变换快速傅里叶变换 (FFT) 卷积线性时不变系统傅里叶级数参考文献离散傅里叶变换 (DFT) 离散傅里叶变换的基对于周期为 ...

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

FFT原理

二维FFT原理

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

最新文章

热门文章