fft&ifft.cpp

使用fft实现的快速卷积

/**************************************************************************
 *  文件名：FFT&IFFT.cpp
 *
 *  函数：
 *
 *  FFT()				- 快速傅立叶变换
 *  IFFT()				- 快速傅立叶反变换
 *  FFT2()              - 图像的快速傅立叶变换
 *  IFFT2()             - 图像的快速傅立叶反变换
 *
 *************************************************************************/

#include "stdafx.h"
#include "FFT&IFFT.h"

#include <math.h>
#include <complex>
using namespace std;

// 常数π
#define PI 3.1415926535

/*************************************************************************
 *
 * 函数名称：
 *   FFT()
 *
 * 参数:
 *   complex<double> * TD	- 指向时域数组的指针
 *   complex<double> * FD	- 指向频域数组的指针
 *   r						- 2的幂数，即迭代次数，2的r次方=N
 *
 * 说明:
 *   该函数用来实现快速付立叶变换。
 *
 ************************************************************************/
void FFT(complex<double> * TD, complex<double> * FD, int r)
{
	// 付立叶变换点数
	long	count;
	
	// 循环变量
	int		i,j,k;
	
	// 中间变量
	int		bfsize,p;
	
	// 角度
	double	angle;
	
	complex<double> *W,*X1,*X2,*X;
	
	// 计算付立叶变换点数
	count = 1 << r;     //1左移r位，即2的r次方
	
	// 分配运算所需存储器
	W  = new complex<double>[count / 2];
	X1 = new complex<double>[count];
	X2 = new complex<double>[count];
	
	// 计算加权系数
	for(i = 0; i < count / 2; i++)
	{
		angle = -i * PI * 2 / count;
		W[i] = complex<double> (cos(angle), sin(angle));
	}
	
	// 将时域点写入X1
	memcpy(X1, TD, sizeof(complex<double>) * count);
	
	// 采用蝶形算法进行快速付立叶变换
	for(k = 0; k < r; k++)  //FFT的级数
	{
		for(j = 0; j < 1 << k; j++)   //组数
		{
			bfsize = 1 << (r-k);  //本级的跨度*2
			for(i = 0; i < bfsize / 2; i++)
			{
				p = j * bfsize;
				X2[i + p] = X1[i + p] + X1[i + p + bfsize / 2];
				X2[i + p + bfsize / 2] = (X1[i + p] - X1[i + p + bfsize / 2]) * W[i * (1<<k)];
			}
		}
		X  = X1;
		X1 = X2;
		X2 = X;
	}
	
	// 重新排序(输出是码位倒序的)
	for(j = 0; j < count; j++)
	{
		p = 0;
		for(i = 0; i < r; i++)
		{
			if (j&(1<<i))
			{
				p+=1<<(r-i-1);
			}
		}
		FD[j]=X1[p];
	}
	
	// 释放内存
	delete W;
	delete X1;
	delete X2;
}

/*************************************************************************
 *
 * 函数名称：
 *   IFFT()
 *
 * 参数:
 *   complex<double> * FD	- 指向频域值的指针
 *   complex<double> * TD	- 指向时域值的指针
 *   r						- 2的幂数，即迭代次数，2的r次方=N
 *
 * 说明:
 *   该函数用来实现快速付立叶反变换。
 *
 ************************************************************************/
void IFFT(complex<double> * FD, complex<double> * TD, int r)
{
	// 付立叶变换点数
	long	count;
	
	// 循环变量
	int		i;
	
	complex<double> *X;
	
	// 计算付立叶变换点数
	count = 1 << r;   //1左移r位，即2的r次方
	
	// 分配运算所需存储器
	X = new complex<double>[count];
	
	// 将频域点写入X
	memcpy(X, FD, sizeof(complex<double>) * count);

	// 求共轭
	for(i = 0; i < count; i++)
	{
		X[i] = complex<double> (X[i].real(), -X[i].imag());
	}
	
	// 调用快速付立叶变换
	FFT(X, TD, r);
	
	// 求时域点的共轭
	for(i = 0; i < count; i++)
	{
		TD[i] = complex<double> (TD[i].real() / count, -TD[i].imag() / count);
	}
	
	// 释放内存
	delete X;
}

/*************************************************************************
 *
 * 函数名称：
 *   FFT2()
 *
 * 参数:
 *   complex<double> * TD	- 指向时域数组的指针
 *   complex<double> * FD	- 指向频域数组的指针
 *   width, height			- 2维变换的宽度和高度
 *
 * 说明:
 *   该函数用来实现对图像的快速付立叶变换。
 *
 ************************************************************************/
void FFT2(complex<double> * TD, complex<double> * FD, long width, long height)
{
    // 求宽度、高度的2的幂(这里暂时要求输入图像的宽度和高度都是2的整数次方)
    int wr=0;  // 2的幂数，即FFT迭代次数，2的r次方=width
    wr = log(width)/log(2);

    int hr=0;  // 2的幂数，即FFT迭代次数，2的r次方=height
    hr = log(height)/log(2);
	
	// 循环变量
	long i, j;
	for(i = 0; i < height; i++)
	{
		// 对x方向进行快速付立叶变换
		FFT(&TD[width * i], &FD[width * i], wr);  //每一行的FFT
	}
	
	// 保存变换结果,行列互换便于y方向变换
	for(i = 0; i < height; i++)
	{
		for(j = 0; j < width; j++)
		{
			TD[i + height * j] = FD[j + width * i];
		}
	}

	//临时变量
	complex<double> *FDtemp;

	FDtemp = new complex<double>[width * height];  //保存频域变换结果
	
	for(i = 0; i < width; i++)
	{
		// 对y方向进行快速付立叶变换
		FFT(&TD[i * height], &FDtemp[i * height], hr);  //每一列的FFT
	}

	//将FDtemp再行列互换一次才是最终结果
    for(i = 0; i < height; i++)
	{
		for(j = 0; j < width; j++)
		{
			FD[i * width + j] = FDtemp[j * height + i];
		}
	}

	//释放内存
	delete FDtemp;
}

/*************************************************************************
 *
 * 函数名称：
 *   IFFT2()
 *
 * 参数:
 *   complex<double> * FD	- 指向频域值的指针
 *   complex<double> * TD	- 指向时域值的指针
 *   width, height			- 2维变换的宽度和高度
 *
 * 说明:
 *   该函数用来实现对图像的快速付立叶反变换。
 *
 ************************************************************************/
void IFFT2(complex<double> * FD, complex<double> * TD, long width, long height)
{
    // 求宽度、高度的2的幂(这里暂时要求输入图像的宽度和高度都是2的整数次方)
    int wr=0;  // 2的幂数，即FFT迭代次数，2的r次方=width
    wr = log(width)/log(2);

    int hr=0;  // 2的幂数，即FFT迭代次数，2的r次方=height
    hr = log(height)/log(2);
	
	// 循环变量
	long i, j;

	for(i = 0; i < height; i++)
	{
		// 对x方向进行快速付立叶反变换
		IFFT(&FD[i * width], &TD[i * width], wr);  //每一行的IFFT
	}

	// 保存变换结果,行列互换便于y方向变换
	for(i = 0; i < height; i++)
	{
		for(j = 0; j < width; j++)
		{
			FD[i + height * j] = TD[j + width * i];
		}
	}
	
    //临时变量
	complex<double> *TDtemp;

	TDtemp = new complex<double>[width * height];  //保存时域变换结果

	for(i = 0; i < width; i++)
	{
		// 对y方向进行快速付立叶反变换
		IFFT(&FD[height * i], &TDtemp[height * i], hr);  //每一列的IFFT
	}

	//将TDtemp再行列互换一次才是最终结果
    for(i = 0; i < height; i++)
	{
		for(j = 0; j < width; j++)
		{
			TD[i * width + j] = TDtemp[j * height + i];
		}
	}

	//释放内存
	delete TDtemp;
}