/************************************************************************** *  文件名：FreTrans.cpp * *  正交變換API函數庫： * *  FFT()				- 快速付立葉變換 *  IFFT()				- 快速付立葉反變換 *  DCT()				- 離散余弦變換 *  WALSH()				- 沃爾什－哈達瑪變換 * *  Fourier()			- 圖像的付立葉變換 *  DIBDct()			- 圖像的離散余弦變換 *  DIBWalsh()			- 圖像的沃爾什－哈達瑪變換 * *************************************************************************/#include "FreTrans.h"#include "stdafx.h"#include "DIBAPI.h"#include <math.h>#include <direct.h>#include <complex>using namespace std;// 常數π#define PI 3.1415926535/************************************************************************* * * 函數名稱： *   FFT() * * 參數: *   complex<double> * TD	- 指向時域數組的指針 *   complex<double> * FD	- 指向頻域數組的指針 *   r						－2的冪數，即迭代次數 * * 返回值: *   無。 * * 說明: *   該函數用來實現快速付立葉變換。 * ************************************************************************/VOID WINAPI 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;		// 分配運算所需存儲器	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++)	{		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<<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的冪數 * * 返回值: *   無。 * * 說明: *   該函數用來實現快速付立葉反變換。 * ************************************************************************/VOID WINAPI IFFT(complex<double> * FD, complex<double> * TD, int r){	// 付立葉變換點數	LONG	count;		// 循環變量	int		i;		complex<double> *X;		// 計算付立葉變換點數	count = 1 << 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;}/************************************************************************* * * 函數名稱： *   Fourier() * * 參數: *   LPSTR lpDIBBits    - 指向源DIB圖像指針 *   LONG  lWidth       - 源圖像寬度（象素數） *   LONG  lHeight      - 源圖像高度（象素數） * * 返回值: *   BOOL               - 成功返回TRUE，否則返回FALSE。 * * 說明: *   該函數用來對圖像進行付立葉變換。 * ************************************************************************/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<double> *TD = new complex<double>[w * h];	complex<double> *FD = new complex<double>[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<double>(*(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<h/2 ? i+h/2 : i-h/2)行，第(j<w/2 ? j+w/2 : j-w/2)個象素的指針			// 此處不直接取i和j，是為了將變換后的原點移到中心			//lpSrc = (unsigned char*)lpDIBBits + lLineBytes * (lHeight - 1 - i) + j;			lpSrc = (unsigned char*)lpDIBBits + lLineBytes * 				(lHeight - 1 - (i<h/2 ? i+h/2 : i-h/2)) + (j<w/2 ? j+w/2 : j-w/2);						// 更新源圖像			* (lpSrc) = (BYTE)(dTemp);		}	}		// 刪除臨時變量	delete TD;	delete FD;		// 返回	return TRUE;}/************************************************************************* * * 函數名稱： *   DCT() * * 參數: *   double * f				- 指向時域值的指針 *   double * F				- 指向頻域值的指針 *   r						－2的冪數 * * 返回值: *   無。 * * 說明: *   該函數用來實現快速離散余弦變換。該函數利用2N點的快速付立葉變換 * 來實現離散余弦變換。 * ************************************************************************/VOID WINAPI DCT(double *f, double *F, int r){	// 離散余弦變換點數	LONG	count;		// 循環變量	int		i;		// 中間變量	double	dTemp;		complex<double> *X;		// 計算離散余弦變換點數	count = 1<<r;		// 分配內存	X = new complex<double>[count*2];		// 賦初值為0	memset(X, 0, sizeof(complex<double>) * count * 2);		// 將時域點寫入數組X	for(i=0;i<count;i++)	{		X[i] = complex<double> (f[i], 0);	}		// 調用快速付立葉變換	FFT(X,X,r+1);		// 調整系數	dTemp = 1/sqrt(count);		// 求F[0]	F[0] = X[0].real() * dTemp;		dTemp *= sqrt(2);		// 求F[u]		for(i = 1; i < count; i++)	{		F[i]=(X[i].real() * cos(i*PI/(count*2)) + X[i].imag() * sin(i*PI/(count*2))) * dTemp;	}		// 釋放內存	delete X;}/************************************************************************* * * 函數名稱： *   IDCT() * * 參數: *   double * F				- 指向頻域值的指針 *   double * f				- 指向時域值的指針 *   r						－2的冪數 * * 返回值: *   無。 * * 說明: *   該函數用來實現快速離散余弦反變換。該函數也利用2N點的快速付立葉變換 * 來實現離散余弦反變換。 * ************************************************************************/VOID WINAPI IDCT(double *F, double *f, int r){	// 離散余弦反變換點數	LONG	count;		// 循環變量	int		i;		// 中間變量	double	dTemp, d0;		complex<double> *X;		// 計算離散余弦變換點數	count = 1<<r;		// 分配內存	X = new complex<double>[count*2];		// 賦初值為0	memset(X, 0, sizeof(complex<double>) * count * 2);		// 將頻域變換后點寫入數組X	for(i=0;i<count;i++)	{		X[i] = complex<double> (F[i] * cos(i*PI/(count*2)), F[i] * sin(i*PI/(count*2)));	}		// 調用快速付立葉反變換	IFFT(X,X,r+1);		// 調整系數	dTemp = sqrt(2.0/count);	d0 = (sqrt(1.0/count) - dTemp) * F[0];		// 計算f(x)	for(i = 0; i < count; i++)	{		f[i] = d0 + X[i].real()* dTemp * 2 * count;	}		// 釋放內存	delete X;}/************************************************************************* * * 函數名稱： *   DIBDct() * * 參數: *   LPSTR lpDIBBits    - 指向源DIB圖像指針 *   LONG  lWidth       - 源圖像寬度（象素數） *   LONG  lHeight      - 源圖像高度（象素數） * * 返回值: *   BOOL               - 成功返回TRUE，否則返回FALSE。 * * 說明: *   該函數用來對圖像進行離散余弦變換。 * ************************************************************************/BOOL WINAPI DIBDct(LPSTR lpDIBBits, LONG lWidth, LONG lHeight){		// 指向源圖像的指針	unsigned char*	lpSrc;		// 循環變量	LONG	i;	LONG	j;		// 進行付立葉變換的寬度和高度（2的整數次方）	LONG	w;	LONG	h;		// 中間變量	double	dTemp;