C++实现二维快速傅里叶变换(FFT)

10 篇文章 1 订阅
3 篇文章 4 订阅

       上一篇文章里,我根据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)); //等于WN
	Cpx 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)); //等于WN
	Cpx 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)
{
	//第一遍fft
	for (int i = 0; i < height; i++)
	{
		vector<Cpx> tempData(width);
		//获取每行数据
		for (int j = 0; j < width; j++)
		{
			tempData[j] = src[i][j];
		}
		//一维FFT
		fft2(tempData, width, opt);
		//写入每行数据
		for (int j = 0; j < width; j++)
		{
			dst[i][j] = tempData[j];
		}
	}

	//第二遍fft
	for (int i = 0; i < width; i++)
	{
		vector<Cpx> tempData(height);
		//获取每列数据
		for (int j = 0; j < height; j++)
		{
			tempData[j] = dst[j][i];
		}
		//一维FFT
		fft2(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#实现 - 死猫 - 博客园

  • 25
    点赞
  • 105
    收藏
    觉得还不错? 一键收藏
  • 6
    评论

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

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

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

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

抵扣说明:

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

余额充值