devbox
我家小花儿
这个屌丝很懒，什么也没留下！
关注作者
热门标签
热门文章
当前位置:   article > 正文

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

作者：我家小花儿 | 2024-02-18 05:01:09

赞

踩

c++ fft

       上一篇文章里，我根据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++代码如下

  1. void fft(vector<Cpx>& a, int lim, int opt)
  2. {
  3. 	if (lim == 1) return;
  4. 	vector<Cpx> a0(lim >> 1), a1(lim >> 1); // 初始化一半大小，存放偶数和奇数部分
  5. 	for (int i = 0; i < lim; i += 2)
  6. 		a0[i >> 1] = a[i], a1[i >> 1] = a[i + 1]; // 分成偶数部分和奇数部分
  7.  
  8. 	fft(a0, lim >> 1, opt); // 递归计算偶数部分
  9. 	fft(a1, lim >> 1, opt); // 递归计算偶数部分
  10.  
  11. 	Cpx wn(cos(2 * PI / lim), opt * -sin(2 * PI / lim)); //等于WN
  12. 	Cpx w(1, 0);
  13. 	for (int k = 0; k < (lim >> 1); k++) // 见蝶形图1运算过程
  14. 	{
  15. 		a[k] = a0[k] + w * a1[k];
  16. 		a[k + (lim >> 1)] = a0[k] - w * a1[k];
  17. 		w = w * wn;
  18. 	}
  19.  
  20. 	//for (int k = 0; k < (lim >> 1); k++) // 见蝶形图2，小优化一下，少一次乘法
  21. 	//{
  22. 	//	Cpx t = w * a1[k];
  23. 	//	a[k] = a0[k] + t;
  24. 	//	a[k + (lim >> 1)] = a0[k] - t;
  25. 	//	w = w * wn;
  26. 	//}
  27.  
  28. }

FFT方法二：

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

 

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

  1. //二进制逆序排列
  2. int ReverseBin(int a, int n)
  3. {
  4. 	int ret = 0;
  5. 	for (int i = 0; i < n; i++)
  6. 	{
  7. 		if (a&(1 << i)) ret |= (1 << (n - 1 - i));
  8. 	}
  9. 	return ret;
  10. }
  11.  
  12. void fft2(vector<Cpx>& a, int lim, int opt)
  13. {
  14. 	int index;
  15. 	vector<Cpx> tempA(lim);
  16. 	for (int i = 0; i < lim; i++)
  17. 	{
  18. 		index = ReverseBin(i, log2(lim));
  19. 		tempA[i] = a[index];
  20. 	}
  21.  
  22. 	vector<Cpx> WN(lim / 2);
  23. 	//生成WN表,避免重复计算
  24. 	for (int i = 0; i < lim / 2; i++)
  25. 	{
  26. 		WN[i] = Cpx(cos(2 * PI *i / lim), opt * -sin(2 * PI*i / lim));
  27. 	}
  28.  
  29. 	//蝶形运算
  30. 	int Index0, Index1;
  31. 	Cpx temp;
  32. 	for (int steplenght = 2; steplenght <= lim; steplenght *= 2)
  33. 	{
  34. 		for (int step = 0; step < lim / steplenght; step++)
  35. 		{
  36. 			for (int i = 0; i < steplenght / 2; i++)
  37. 			{
  38. 				Index0 = steplenght * step + i;
  39. 				Index1 = steplenght * step + i + steplenght / 2;
  40.  
  41. 				temp = tempA[Index1] * WN[lim / steplenght * i];
  42. 				tempA[Index1] = tempA[Index0] - temp;
  43. 				tempA[Index0] = tempA[Index0] + temp;
  44. 			}
  45. 		}
  46. 	}
  47. 	for (int i = 0; i < lim; i++)
  48. 	{
  49. 		if (opt == -1)
  50. 		{
  51. 			a[i] = tempA[i] / lim;
  52. 		}
  53. 		else
  54. 		{
  55. 			a[i] = tempA[i];
  56. 		}
  57. 	}
  58. }

二维FFT原理

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

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

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

例子

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

原图

二维FFT

逆变换

 C++实现二维FFT代码

  1. #include <iostream>
  2. #include<opencv2/opencv.hpp>
  3.  
  4. using namespace cv;
  5. using namespace std;
  6. const int height = 512, width = 512;
  7.  
  8. const double PI = acos(-1); // pi值
  9.  
  10.  
  11. struct Cpx // 定义一个复数结构体和复数运算法则
  12. {
  13. 	double r, i;
  14. 	Cpx() : r(0), i(0) {}
  15. 	Cpx(double _r, double _i) : r(_r), i(_i) {}
  16. };
  17. Cpx operator + (Cpx a, Cpx b) { return Cpx(a.r + b.r, a.i + b.i); }
  18. Cpx operator - (Cpx a, Cpx b) { return Cpx(a.r - b.r, a.i - b.i); }
  19. Cpx operator * (Cpx a, Cpx b) { return Cpx(a.r * b.r - a.i * b.i, a.r * b.i + a.i * b.r); }
  20. Cpx operator / (Cpx a, int b) { return Cpx(a.r*1.0 / b, a.i*1.0 / b); }
  21.  
  22. Mat BGR2GRAY(Mat img)
  23. {
  24. 	int w = img.cols;
  25. 	int h = img.rows;
  26. 	Mat grayImg(h, w, CV_8UC1);
  27.  
  28. 	uchar *p = grayImg.ptr<uchar>(0);
  29. 	Vec3b *pImg = img.ptr<Vec3b>(0);
  30.  
  31. 	for (int i = 0; i < w*h; ++i)
  32. 	{
  33. 		p[i] = 0.2126*pImg[i][2] + 0.7152*pImg[i][1] + 0.0722*pImg[i][0];
  34. 	}
  35. 	return grayImg;
  36. }
  37.  
  38. Mat Resize(Mat img)
  39. {
  40. 	int w = img.cols;
  41. 	int h = img.rows;
  42. 	Mat out(height, width, CV_8UC1);
  43. 	uchar *p = out.ptr<uchar>(0);
  44. 	uchar *p2 = img.ptr<uchar>(0);
  45. 	int x_before, y_before;
  46. 	for (int y = 0; y < height; y++) {
  47. 		for (int x = 0; x < width; x++)
  48. 		{
  49. 			x_before = (int)x*w*1.0 / width;
  50. 			y_before = (int)y*h*1.0 / height;
  51. 			p[y * width + x] = p2[y_before * w + x_before];
  52. 		}
  53. 	}
  54. 	return out;
  55.  
  56. }
  57.  
  58. //自动缩放到0-255范围内,并变换象限，将低频移至中间
  59. int AutoScale(Mat src, Mat out)
  60. {
  61. 	int w = src.cols;
  62. 	int h = src.rows;
  63. 	double *p = src.ptr<double>(0);
  64. 	uchar *pOut = out.ptr<uchar>(0);
  65. 	double max = p[0];
  66. 	double min = p[0];
  67. 	for (int i = 0; i < w*h; i++)
  68. 	{
  69. 		if (p[i] > max) max = p[i];
  70. 		if (p[i] < min) min = p[i];
  71. 	}
  72.  
  73. 	double scale = 255.0 / (max - min);
  74.  
  75. 	for (int i = 0; i < w*h; i++)
  76. 	{
  77. 		int j = i + w * h / 2 + w / 2;
  78. 		if (j > w*h) j = j - w * h;   //低频移至中间
  79. 		pOut[i] = (uchar)((p[j] - min)*scale);
  80. 	}
  81. 	return 0;
  82. }
  83.  
  84. void fft(vector<Cpx>& a, int lim, int opt)
  85. {
  86. 	if (lim == 1) return;
  87. 	vector<Cpx> a0(lim >> 1), a1(lim >> 1); // 初始化一半大小，存放偶数和奇数部分
  88. 	for (int i = 0; i < lim; i += 2)
  89. 		a0[i >> 1] = a[i], a1[i >> 1] = a[i + 1]; // 分成偶数部分和奇数部分
  90.  
  91. 	fft(a0, lim >> 1, opt); // 递归计算偶数部分
  92. 	fft(a1, lim >> 1, opt); // 递归计算偶数部分
  93.  
  94. 	Cpx wn(cos(2 * PI / lim), opt * -sin(2 * PI / lim)); //等于WN
  95. 	Cpx w(1, 0);
  96. 	for (int k = 0; k < (lim >> 1); k++) // 见蝶形图1运算过程
  97. 	{
  98. 		a[k] = a0[k] + w * a1[k];
  99. 		a[k + (lim >> 1)] = a0[k] - w * a1[k];
  100. 		w = w * wn;
  101. 	}
  102.  
  103. 	//for (int k = 0; k < (lim >> 1); k++) // 见蝶形图2，小优化一下，少一次乘法
  104. 	//{
  105. 	//	Cpx t = w * a1[k];
  106. 	//	a[k] = a0[k] + t;
  107. 	//	a[k + (lim >> 1)] = a0[k] - t;
  108. 	//	w = w * wn;
  109. 	//}
  110.  
  111. }
  112.  
  113. //二进制逆序排列
  114. int ReverseBin(int a, int n)
  115. {
  116. 	int ret = 0;
  117. 	for (int i = 0; i < n; i++)
  118. 	{
  119. 		if (a&(1 << i)) ret |= (1 << (n - 1 - i));
  120. 	}
  121. 	return ret;
  122. }
  123.  
  124. void fft2(vector<Cpx>& a, int lim, int opt)
  125. {
  126. 	int index;
  127. 	vector<Cpx> tempA(lim);
  128. 	for (int i = 0; i < lim; i++)
  129. 	{
  130. 		index = ReverseBin(i, log2(lim));
  131. 		tempA[i] = a[index];
  132. 	}
  133.  
  134. 	vector<Cpx> WN(lim / 2);
  135. 	//生成WN表,避免重复计算
  136. 	for (int i = 0; i < lim / 2; i++)
  137. 	{
  138. 		WN[i] = Cpx(cos(2 * PI *i / lim), opt * -sin(2 * PI*i / lim));
  139. 	}
  140.  
  141. 	//蝶形运算
  142. 	int Index0, Index1;
  143. 	Cpx temp;
  144. 	for (int steplenght = 2; steplenght <= lim; steplenght *= 2)
  145. 	{
  146. 		for (int step = 0; step < lim / steplenght; step++)
  147. 		{
  148. 			for (int i = 0; i < steplenght / 2; i++)
  149. 			{
  150. 				Index0 = steplenght * step + i;
  151. 				Index1 = steplenght * step + i + steplenght / 2;
  152.  
  153. 				temp = tempA[Index1] * WN[lim / steplenght * i];
  154. 				tempA[Index1] = tempA[Index0] - temp;
  155. 				tempA[Index0] = tempA[Index0] + temp;
  156. 			}
  157. 		}
  158. 	}
  159. 	for (int i = 0; i < lim; i++)
  160. 	{
  161. 		if (opt == -1)
  162. 		{
  163. 			a[i] = tempA[i] / lim;
  164. 		}
  165. 		else
  166. 		{
  167. 			a[i] = tempA[i];
  168. 		}
  169. 	}
  170. }
  171.  
  172. void FFT2D(Cpx(*src)[width], Cpx(*dst)[width], int opt)
  173. {
  174. 	//第一遍fft
  175. 	for (int i = 0; i < height; i++)
  176. 	{
  177. 		vector<Cpx> tempData(width);
  178. 		//获取每行数据
  179. 		for (int j = 0; j < width; j++)
  180. 		{
  181. 			tempData[j] = src[i][j];
  182. 		}
  183. 		//一维FFT
  184. 		fft2(tempData, width, opt);
  185. 		//写入每行数据
  186. 		for (int j = 0; j < width; j++)
  187. 		{
  188. 			dst[i][j] = tempData[j];
  189. 		}
  190. 	}
  191.  
  192. 	//第二遍fft
  193. 	for (int i = 0; i < width; i++)
  194. 	{
  195. 		vector<Cpx> tempData(height);
  196. 		//获取每列数据
  197. 		for (int j = 0; j < height; j++)
  198. 		{
  199. 			tempData[j] = dst[j][i];
  200. 		}
  201. 		//一维FFT
  202. 		fft2(tempData, height, opt);
  203. 		//写入每列数据
  204. 		for (int j = 0; j < height; j++)
  205. 		{
  206. 			dst[j][i] = tempData[j];
  207. 		}
  208. 	}
  209. }
  210.  
  211. void Mat2Cpx(Mat src, Cpx(*dst)[width])
  212. {
  213. 	//这里Mat里的数据得是unchar类型
  214. 	uchar *p = src.ptr<uchar>(0);
  215. 	for (int i = 0; i < height; i++)
  216. 	{
  217. 		for (int j = 0; j < width; j++)
  218. 		{
  219. 			dst[i][j] = Cpx(p[i*width + j],0);
  220. 		}
  221. 	}
  222. }
  223.  
  224. void Cpx2Mat(Cpx(*src)[width], Mat dst)
  225. {
  226. 	//这里Mat里的数据得是unchar类型
  227. 	uchar *p = dst.ptr<uchar>(0);
  228. 	for (int i = 0; i < height; i++)
  229. 	{
  230. 		for (int j = 0; j < width; j++)
  231. 		{
  232. 			double g = sqrt(src[i][j].r*src[i][j].r);
  233. 			p[j + i * width] = (uchar)g;
  234. 		}
  235. 	}
  236. }
  237.  
  238. void Cpx2MatDouble(Cpx(*src)[width], Mat dst)
  239. {
  240. 	//这里Mat里的数据得是double类型
  241. 	double *p = dst.ptr<double>(0);
  242. 	for (int i = 0; i < height; i++)
  243. 	{
  244. 		for (int j = 0; j < width; j++)
  245. 		{
  246. 			double g = sqrt(src[i][j].r*src[i][j].r+ src[i][j].i*src[i][j].i);
  247. 			g = log(g + 1);  //转换为对数尺度
  248. 			p[j + i * width] = (double)g;
  249. 		}
  250. 	}
  251. }
  252.  
  253.  
  254. int main()
  255. {
  256. 	Mat img = imread("d:\\1.jpg");
  257. 	imshow("img", img);
  258.  
  259. 	Mat gray = BGR2GRAY(img);
  260. 	imshow("gray", gray);
  261. 	imwrite("gray.jpg", gray);
  262.  
  263. 	Mat imgRez = Resize(gray);
  264. 	imshow("imgRez", imgRez);
  265. 	imwrite("imgRez.jpg", imgRez);
  266.  
  267.  
  268. 	Cpx(*src)[width] = new Cpx[height][width];
  269.  
  270. 	Mat2Cpx(imgRez, src);
  271.  
  272. 	Cpx(*dst)[width] = new Cpx[height][width];
  273.  
  274.  
  275.  
  276. 	double t1 = getTickCount();
  277. 	FFT2D(src, dst, 1);
  278. 	double t2 = getTickCount();
  279.  
  280. 	double t1t2 = (t2 - t1) / getTickFrequency();
  281. 	std::cout << "DFT耗时: " << t1t2 << "秒" << std::endl;
  282.  
  283.  
  284. 	Cpx(*dst2)[width] = new Cpx[height][width];
  285. 	double t3 = getTickCount();
  286. 	FFT2D(dst, dst2, -1);
  287. 	double t4 = getTickCount();
  288.  
  289. 	double t3t4 = (t4 - t3) / getTickFrequency();
  290. 	std::cout << "DFT耗时: " << t3t4 << "秒" << std::endl;
  291.  
  292. 	Mat out2 = Mat::zeros(height, width, CV_8UC1);
  293. 	Cpx2Mat(dst2, out2);
  294. 	imshow("out2", out2);
  295. 	imwrite("out2.jpg", out2);
  296.  
  297. 	Mat out = Mat::zeros(height, width, CV_64F);
  298. 	Cpx2MatDouble(dst, out);
  299. 	Mat out3 = Mat::zeros(height, width, CV_8UC1);
  300. 	AutoScale(out, out3);
  301. 	imshow("out3", out3);
  302. 	imwrite("out3.jpg", out3);
  303.  
  304.  
  305. 	waitKey(0);
  306. }

OK！

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

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

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

声明：本文内容由网友自发贡献，不代表【wpsshop博客】立场，版权归原作者所有，本站不承担相应法律责任。如您发现有侵权的内容，请联系我们。转载请注明出处：https://www.wpsshop.cn/w/我家小花儿/article/detail/104973
推荐阅读
相关标签

Copyright © 2003-2013 www.wpsshop.cn 版权所有，并保留所有权利。

  

闽ICP备14008679号