www.pudn.com > spline_c++.rar > fft.cpp

#include "fft.h" 
#include "stdafx.h" 
#include "DIBAPI.h" 
 
#include  
#include  
#include  
using namespace std; 
 
// 常数π 
#define PI 3.1415926535 
 
VOID WINAPI FFT(complex * TD, complex * FD, int r) 
{ 
	// 付立叶变换点数 
	LONG	count; 
	 
	// 循环变量 
	int		i,j,k; 
	 
	// 中间变量 
	int		bfsize,p; 
	 
	// 角度 
	double	angle; 
	 
	complex *W,*X1,*X2,*X; 
	 
	// 计算付立叶变换点数 
	count = 1 << r; 
	 
	// 分配运算所需存储器 
	W  = new complex[count / 2]; 
	X1 = new complex[count]; 
	X2 = new complex[count]; 
	 
	// 计算加权系数 
	for(i = 0; i < count / 2; i++) 
	{ 
		angle = -i * PI * 2 / count; 
		W[i] = complex (cos(angle), sin(angle)); 
	} 
	 
	// 将时域点写入X1 
	memcpy(X1, TD, sizeof(complex) * count); 
	 
	// 采用蝶形算法进行快速付立叶变换 
	for(k = 0; k < r; k++) 
	{ 
		for(j = 0; j < 1 << k; j++) 
		{ 
			bfsize = 1 << (r-k); 
			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< * FD, complex * TD, int r) 
{ 
	// 付立叶变换点数 
	LONG	count; 
	 
	// 循环变量 
	int		i; 
	 
	complex *X; 
	 
	// 计算付立叶变换点数 
	count = 1 << r; 
	 
	// 分配运算所需存储器 
	X = new complex[count]; 
	 
	// 将频域点写入X 
	memcpy(X, FD, sizeof(complex) * count); 
	 
	// 求共轭 
	for(i = 0; i < count; i++) 
	{ 
		X[i] = complex (X[i].real(), -X[i].imag()); 
	} 
	 
	// 调用快速付立叶变换 
	FFT(X, TD, r); 
	 
	// 求时域点的共轭 
	for(i = 0; i < count; i++) 
	{ 
		TD[i] = complex (TD[i].real() / count, -TD[i].imag() / count); 
	} 
	 
	// 释放内存 
	delete X; 
} 
 
BOOL WINAPI Fourier(LPSTR lpDIBBits, LONG lWidth, LONG lHeight) 
{ 
	 
	// 指向源图像的指针 
	unsigned char*	lpSrc; 
	 
	// 中间变量 
	double	dTemp; 
	 
	// 循环变量 
	LONG	i; 
	LONG	j; 
	 
	// 进行付立叶变换的宽度和高度（2的整数次方） 
	LONG	w; 
	LONG	h; 
	 
	int		wp; 
	int		hp; 
	 
	// 图像每行的字节数 
	LONG	lLineBytes; 
	 
	// 计算图像每行的字节数 
	lLineBytes = WIDTHBYTES(lWidth * 8); 
	 
	// 赋初值 
	w = 1; 
	h = 1; 
	wp = 0; 
	hp = 0; 
	 
	// 计算进行付立叶变换的宽度和高度（2的整数次方） 
	while(w * 2 <= lWidth) 
	{ 
		w *= 2; 
		wp++; 
	} 
	 
	while(h * 2 <= lHeight) 
	{ 
		h *= 2; 
		hp++; 
	} 
	 
	// 分配内存 
	complex *TD = new complex[w * h]; 
	complex *FD = new complex[w * h]; 
	 
	// 行 
	for(i = 0; i < h; i++) 
	{ 
		// 列 
		for(j = 0; j < w; j++) 
		{ 
			// 指向DIB第i行，第j个象素的指针 
			lpSrc = (unsigned char*)lpDIBBits + lLineBytes * (lHeight - 1 - i) + j; 
			 
			// 给时域赋值 
			TD[j + w * i] = complex(*(lpSrc), 0); 
		} 
	} 
	 
	for(i = 0; i < h; i++) 
	{ 
		// 对y方向进行快速付立叶变换 
		FFT(&TD[w * i], &FD[w * i], wp); 
	} 
	 
	// 保存变换结果 
	for(i = 0; i < h; i++) 
	{ 
		for(j = 0; j < w; j++) 
		{ 
			TD[i + h * j] = FD[j + w * i]; 
		} 
	} 
	 
	for(i = 0; i < w; i++) 
	{ 
		// 对x方向进行快速付立叶变换 
		FFT(&TD[i * h], &FD[i * h], hp); 
	} 
	 
	// 行 
	for(i = 0; i < h; i++) 
	{ 
		// 列 
		for(j = 0; j < w; j++) 
		{ 
			// 计算频谱 
			dTemp = sqrt(FD[j * h + i].real() * FD[j * h + i].real() +  
				         FD[j * h + i].imag() * FD[j * h + i].imag()) / 100; 
			 
			// 判断是否超过255 
			if (dTemp > 255) 
			{ 
				// 对于超过的，直接设置为255 
				dTemp = 255; 
			} 
			 
			// 指向DIB第(i *TD = new complex[w * h]; 
	complex *FD = new complex[w * h]; 
	 
	// 行 
	for(i = 0; i < h; i++) 
	{ 
		// 列 
		for(j = 0; j < w; j++) 
		{ 
			// 指向DIB第i行，第j个象素的指针 
			lpSrc = (unsigned char*)lpDIBBits + lLineBytes * (lHeight - 1 - i) + j; 
			 
			// 给时域赋值 
			TD[j + w * i] = complex(*(lpSrc), 0); 
		} 
	} 
	 
	for(i = 0; i < h; i++) 
	{ 
		// 对y方向进行快速付立叶变换 
		IFFT(&TD[w * i], &FD[w * i], wp); 
	} 
	 
	// 保存变换结果 
	for(i = 0; i < h; i++) 
	{ 
		for(j = 0; j < w; j++) 
		{ 
			TD[i + h * j] = FD[j + w * i]; 
		} 
	} 
	 
	for(i = 0; i < w; i++) 
	{ 
		// 对x方向进行快速反付立叶变换 
		IFFT(&TD[i * h], &FD[i * h], hp); 
	} 
	 
	// 行 
	for(i = 0; i < h; i++) 
	{ 
		// 列 
		for(j = 0; j < w; j++) 
		{ 
			// 计算频谱 
			dTemp = sqrt(FD[j * h + i].real() * FD[j * h + i].real() +  
				         FD[j * h + i].imag() * FD[j * h + i].imag()) / 100; 
			 
			// 判断是否超过255 
			if (dTemp > 255) 
			{ 
				// 对于超过的，直接设置为255 
				dTemp = 255; 
			} 
			 
			// 指向DIB第(i0&&i0&&j