/*************************************************************
Filename:	Fourier.h.
Writer:		Dai Zhiheng.
Time:		2005-3-27.
**************************************************************/
#include <math.h>
#include "myComplex.h"
#define PI	3.141593
/*************************************************************
函數(shù)功能:	產(chǎn)生倒序序列
			case inout r = 3; this function will format *a as
			{0、4、2、6、1、5、3、7}
**************************************************************/
void reverseOrder(int *a, int r)
{
	int count = 1 << r ;
	for(int j = 0; j < count; j++)
	{
		int p=0;
		for(int i = 0; i < r; i++)
		{
			if (j & (1<<i))
			{
					p += 1 << (r-i-1);
			}
		}
		*(a+j) = p;
	}
}
/*************************************************************
函數(shù)功能:	求矩陣的轉(zhuǎn)置
parameter:
			*f		數(shù)據(jù)指針
			m,n		矩陣維數(shù)
**************************************************************/
void complexArrayT( double *f, long m, long n)
{
	double *X = new  double [m*n*2];
	for(int i=0; i<m; i++)
		for(int j=0; j<n; j++)
		{
			*(X + i*2 + m*j*2) = *(f + i*n*2 +j*2);
			*(X + i*2 + m*j*2 +1) = *(f + i*n*2 +j*2 +1);
		}
	for(long t=0; t<m*n*2; t++)
		*(f+t) = *(X+t);
	delete [] X;
}
/*************************************************************
函數(shù)功能：	1-D,count行,r階Fourier變換
parameter:
			*f		時域數(shù)據(jù)的指針
			*F		頻域數(shù)據(jù)的指針
			 r		Fourier變換的階數(shù)基數(shù)
**************************************************************/
void DLineFourier(double *f, double *F, long r)
{
	int count = 1 << r ;
	int bfSize;
	int p;
	int *order = new int [count];
	double *W  = new double [count];
	double *X  = new double [count*count*2];
	double *X1 = new double [count*count*2];
	double *X2 = new double [count*count*2];
	//計算變換基W
	for(int i=0; i<count/2; i++)
	{
		double angle =-i*2*PI/count;
		*(W + i*2) = cos( angle ) ;
		*(W + i*2 + 1) = sin( angle ) ;
	}

	//產(chǎn)生倒序序列
	reverseOrder(order, r);
	//數(shù)據(jù)從 f （按每行）倒序到 X1
	for(int m=0; m< count; m++)
		for(int n=0; n< count; n++)
		{
			*(X1 + count*m*2 + n*2) = *(f + count*m*2 + (*(order+n))*2 );
			*(X1 + count*m*2 + n*2 + 1) = *(f + count*m*2 + (*(order+n))*2 +1 );
		}
		
	//Butterfly algorithm
	for(int l=0; l< count; l++)
	{
		for(int k=0; k<r; k++ )
		{	
			for(int j=0; j< 1<<(r-k-1); j++)
			{
				bfSize = 1<<k;
				for(int i=0; i<bfSize; i++)
				{
					p = j*bfSize*2;
					double *A  = new double[2];
					complexMUL(A, X1 + count*l*2 + i*2 + p*2 + bfSize*2, W + i*(1<<(r-k)) );
					complexADD(X2 + count*l*2 + i*2 + p*2 , X1 + count*l*2 + i*2 + p*2 , A);
					
					complexDEC(X2 + count*l*2 + i*2 + p*2 + bfSize*2, X1 + count*l*2 + i*2 + p*2, A);
					delete [] A;
				}
			}
			//交換數(shù)據(jù)
			for( int m=0; m< count*2; m++)
			{
				*(X + count*l*2 + m)	= *(X1 + count*l*2 + m);
			}
			for( m=0; m< count*2; m++)
			{
				*(X1 + count*l*2 + m)	= *(X2 + count*l*2 + m);
			}
			for( m=0; m< count*2; m++)
			{
				*(X2 + count*l*2 + m)	= *(X  + count*l*2 + m);
			}
		}
	}

	for(long t=0; t<count*count*2; t++)
	{
		*(F+t) = *(X1 + t);
	}
	delete [] W;
	delete [] order;
	delete [] X;
	delete [] X1;
	delete [] X2;
}
/*************************************************************
函數(shù)功能：	2-D,count行,count列，r階Fourier變換
parameter:
			*f		時域數(shù)據(jù)的指針
			*F		頻域數(shù)據(jù)的指針
			 r		Fourier變換的階數(shù)基數(shù)
**************************************************************/
void D2Fourier( double *f,  double *F, long r)
{
	int count = 1 << r ;
	double *X1 = new  double [count*count*2];
	double *X2 = new  double [count*count*2];
	//行的Fourier變換
	DLineFourier(f, X1, r);
	//轉(zhuǎn)置
	complexArrayT(X1, count, count);
	//列的Fourier變換
	DLineFourier(X1, X2, r);
	//轉(zhuǎn)置
	complexArrayT(X2, count, count);
	for(long m=0; m< count*count*2; m++)
		{
			*(F + m) = *(X2 + m);
		}
	delete [] X1;
	delete [] X2;
}
/*************************************************************
函數(shù)功能：	2-D圖象的r階Fourier變換
parameter:
			*f		時域數(shù)據(jù)的指針
			*F		頻域數(shù)據(jù)的指針
			 height	圖片的長度
			 width	圖片的寬度
			 r		Fourier變換的階數(shù)基數(shù)
**************************************************************/
void D2BMPFourier( double *f,  double *F, long height, long width, long r)
{
	int count = 1 << r ;
	double *X1 = new  double [count*count*2];
	double *X2 = new  double [count*count*2];
	for(int i=0; i<(height/count); i++)
		for(int j=0; j<(width/count); j++)
		{
			//對于每個count * count 塊
			for(int m=0; m< count; m++)
				for(int n=0; n< count*2; n++)
				{
					*(X1 + m*count*2 +n) = *(f + width*i*count*2 + width*m*2 + count*j*2 + n);
				}
			//塊的Fourier變換
			D2Fourier(X1, X2, r);
			//將數(shù)據(jù)輸出到頻域
			for(m=0; m< count; m++)
				for(int n=0; n< count*2; n++)
				{
					*(F + width*i*count*2 + width*m*2 + count*j*2 + n) = *(X2 + m*count*2 +n);
				}
		}
	delete [] X1;
	delete [] X2;
}
/*************************************************************
函數(shù)功能：	2-D,count行,count列，r階Fourier變換
parameter:
			*f		時域數(shù)據(jù)的指針
			*F		頻域數(shù)據(jù)的指針
			 r		Fourier變換的階數(shù)基數(shù)
**************************************************************/
void ID2Fourier( double *f,  double *F, long r)
{
	int count = 1 << r ;
	double *X1 = new  double [count*count*2];
	double *X2 = new  double [count*count*2];

	//行的Fourier變換
	//求復(fù)共軛
	for(int m=0; m< count; m++)
		for(int n=0; n< count; n++)
			{
				*(X1 + m*count*2 + n*2 ) = *(f + m*count*2 + n*2 );
				*(X1 + m*count*2 + n*2 +1) = -*(f + m*count*2 + n*2 +1);
			}
	DLineFourier(X1, X2, r);
	//求復(fù)共軛
	for(m=0; m< count; m++)
		for(int n=0; n< count; n++)
			{
				*(X1 + m*count*2 + n*2 ) = *(X2 + m*count*2 + n*2 )/count;
				*(X1 + m*count*2 + n*2 +1) = -*(X2 + m*count*2 + n*2 +1)/count;
			}

	//轉(zhuǎn)置
	complexArrayT(X1, count, count);

	//列的Fourier變換
	//求復(fù)共軛
	for(m=0; m< count; m++)
		for(int n=0; n< count; n++)
			{
				*(X2 + m*count*2 + n*2 ) = *(X1 + m*count*2 + n*2 );
				*(X2 + m*count*2 + n*2 +1) = -*(X1 + m*count*2 + n*2 +1);
			}
	DLineFourier(X2, X1, r);
	//求復(fù)共軛
	for(m=0; m< count; m++)
		for(int n=0; n< count; n++)
			{
				*(X2 + m*count*2 + n*2 ) = *(X1 + m*count*2 + n*2 )/count;
				*(X2 + m*count*2 + n*2 +1) = -*(X1 + m*count*2 + n*2 +1)/count;
			}	
	//轉(zhuǎn)置
	complexArrayT(X2, count, count);
	for(long t=0; t< count*count*2; t++)
		{
			*(F + t) = *(X2 + t);
		}
	delete [] X1;
	delete [] X2;
}
/*************************************************************
函數(shù)功能：	2-D圖象的r階Fourier逆變換
parameter:
			*F		頻域數(shù)據(jù)的指針
			*f		時域數(shù)據(jù)的指針
			 height	圖片的長度
			 width	圖片的寬度
			 r		Fourier變換的階數(shù)基數(shù)
**************************************************************/
void ID2BMPFourier( double *F,  double *f, long height, long width, long r)
{
	int count = 1 << r ;
	double *X1 = new double [count*count*2];
	double *X2 = new double [count*count*2];
	for(int i=0; i<(height/count); i++)
		for(int j=0; j<(width/count); j++)
		{
			//對于每個count * count 塊
			for(int m=0; m< count; m++)
				for(int n=0; n< count*2; n++)
				{
					*(X1 + m*count*2 + n)
						= *(F + width*i*count*2 + width*m*2 + count*j*2 + n);
				}

			//塊的Fourier變換
			ID2Fourier(X1, X2, r);

			//將數(shù)據(jù)輸出至?xí)r域
			for(m=0; m< count; m++)
				for(int n=0; n< count*2; n++)
				{
					*(f + width*i*count*2 + width*m*2 + count*j*2 + n)
						= (*(X2 + m*count*2 + n));
				}
		}
	delete [] X1;
	delete [] X2;
}