fourier.h

实现jpeg的压缩编码功能 在demo2的工程文件中即可查看

/*************************************************************
Filename:	Fourier.h.
Writer:		Dai Zhiheng.
Time:		2005-3-27.
**************************************************************/
#include <math.h>
#include "myComplex.h"
#define PI	3.141593
/*************************************************************
函数功能:	产生倒序序列
			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;
	}
}
/*************************************************************
函数功能:	求矩阵的转置
parameter:
			*f		数据指针
			m,n		矩阵维数
**************************************************************/
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;
}
/*************************************************************
函数功能：	1-D,count行,r阶Fourier变换
parameter:
			*f		时域数据的指针
			*F		频域数据的指针
			 r		Fourier变换的阶数基数
**************************************************************/
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 ) ;
	}

	//产生倒序序列
	reverseOrder(order, r);
	//数据从 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;
				}
			}
			//交换数据
			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;
}
/*************************************************************
函数功能：	2-D,count行,count列，r阶Fourier变换
parameter:
			*f		时域数据的指针
			*F		频域数据的指针
			 r		Fourier变换的阶数基数
**************************************************************/
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);
	//转置
	complexArrayT(X1, count, count);
	//列的Fourier变换
	DLineFourier(X1, X2, r);
	//转置
	complexArrayT(X2, count, count);
	for(long m=0; m< count*count*2; m++)
		{
			*(F + m) = *(X2 + m);
		}
	delete [] X1;
	delete [] X2;
}
/*************************************************************
函数功能：	2-D图象的r阶Fourier变换
parameter:
			*f		时域数据的指针
			*F		频域数据的指针
			 height	图片的长度
			 width	图片的宽度
			 r		Fourier变换的阶数基数
**************************************************************/
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);
			//将数据输出到频域
			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;
}
/*************************************************************
函数功能：	2-D,count行,count列，r阶Fourier变换
parameter:
			*f		时域数据的指针
			*F		频域数据的指针
			 r		Fourier变换的阶数基数
**************************************************************/
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变换
	//求复共轭
	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);
	//求复共轭
	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;
			}

	//转置
	complexArrayT(X1, count, count);

	//列的Fourier变换
	//求复共轭
	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);
	//求复共轭
	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;
			}	
	//转置
	complexArrayT(X2, count, count);
	for(long t=0; t< count*count*2; t++)
		{
			*(F + t) = *(X2 + t);
		}
	delete [] X1;
	delete [] X2;
}
/*************************************************************
函数功能：	2-D图象的r阶Fourier逆变换
parameter:
			*F		频域数据的指针
			*f		时域数据的指针
			 height	图片的长度
			 width	图片的宽度
			 r		Fourier变换的阶数基数
**************************************************************/
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);

			//将数据输出至时域
			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;
}