Black Engine

/* File		:	struct.h
 * Creator	:	BlackEngine
 * Date		:	2006/03/28
 * Version	:	0.0.1 ( created )
 * Descript	:	Definition of Complex Number , and etc
 */
 
typedef struct _COMPLEX {
	double	real;
	double	imag;
} COMPLEX, *pCOMPLEX;

/* File		:	fft.h
 * Creator	:	BlackEngine
 * Date		:	2006/03/28
 * Version	:	0.0.1 ( created )
 * Descript	:	Perform the FFT Algorithm
 */
#ifndef _BLACKENGINE_FFT_H_
#define _BLACKENGINE_FFT_H_
 
#include "struct.h"
/* TODO : at 03/29 , decision the return values of FFT Function. */
int fft_2d(COMPLEX ** ,int ,int ,int );
int fft(int ,int ,double *,double *);
#endif /* ifndef _BLACKENGINE_FFT_H_ */

/* File		:	fft.c
 * Creator	:	BlackEngine
 * Date		:	2006/03/28
 * Version	:	0.0.1 ( created )
 * Descript	:	Perform the FFT Algorithm ( 2-D Image FFT )
 */
#include 
#include 
#include
 
#include "fft.h"
 
#define TRUE 1
#define FALSE 0
/*#define NULL 0*/
 
/* Powerof2 function : compare the integer with 2^m
 * and find maximum m , that is 2^m <= nx
 * twopm = 2^m
 */
int Powerof2(int nx, int* m, int* twopm) {
	int pwr;
	*m = 0;
	for ( pwr = 1 ; pwr < nx ; pwr = pwr*2) {
		*m = *m + 1;
	}
	*twopm = pwr;
	return(TRUE);
}
 
int fft_2d(COMPLEX **c, int nx, int ny, int dir) {
	int i,j;
	int m,twopm;
	double *real,*imag;
 
	/* Transform the rows */
	real = (double *)malloc(nx * sizeof(double));
	imag = (double *)malloc(nx * sizeof(double));
	if (real == NULL || imag == NULL)
		return(FALSE);
	if (!Powerof2(nx,&m,&twopm) || twopm != nx)
		return(FALSE);
	for (j=0;j<<= m;
 
	/* N의 절반값을 half에 저장한다. */
	half = N >> 1;
 
	/* Bit Reverse를 하는 과정 */
	for ( i = 0 , j = 0 ; i < N - 1 ; i++ ) {
		if ( i < j )
		{
			swap_x = x[i];
			swap_y = y[i];
			x[i] = x[j];
			y[i] = y[j];
			x[j] = swap_x;
			y[j] = swap_y;
		}
		step = half;
 
		while ( step <= j )
		{
			j = j - k;
			k >>= 1;
		}
		j = j + step;
	}
 
	/* 초기의 값을 입력하는 부분,
	 * factor의 역할은 각 단계(Stage)에서 omega의 값을 변경하기 위한
	 * 기본 값이다. 즉 omega N= 4일 경우에는 각도가 pi/2 만큼씩 변하므로
	 * exp(-j pi/2) 의 값을 가지고 있는 것이 factor이다
	 * factor는 각 stage가 올라갈 때 마다 변경되어야 한다.
	 * 그 각도가 절반으로 줄어야 하므로
	 * cos(t) + j sin(t) 가 cos(t/2) + jsin(t/2) 로 되어야 한다.
	 * cos(t) = 2cos^2(t/2) -1 이므로 cos(t/2) = sqrt((1+cos(t))/2)이고
	 * sin(t/2) = sqrt(1-cos^2(t/2))이므로 sqrt((1-cos(t))/2) 이다.
	 */
	factor_x = -1.0;
	factor_y = 0;
	butterfly = 1;
 
	for ( stage = 0 ; stage < m ; stage++ )
	{
		/* step 은 butterfly가 두개의 샘플을 가지고 작업하므로
		 * 현재 위치에서 몇번째 샘플을 가지고 해야 하는 가를
		 * 가리키는 변수이다.
		 * butterfly는 같은 omega 값에의해 반복되는 butterfly가
		 * 존재하므로 그 위치를 가리키기 위한 값이다.
		 */
		step = butterfly;
		butterfly <<= 1;
		/* 각 단계에서 omega는 항상 1+j0으로 시작한다. */
		omega_x = 1.0;
		omega_y = 0;
		for ( numofbutter = 0 ; numofbutter < step ; numofbutter++)
		{
			for( i = numofbutter; i < N; i = i + butterfly)
			{
				/* omega_N(k) 일때 omega_N(k+N/2) = -omega_N(k)
				 * 의 성질을 이용한다. 그래서 두번째 샘플에 미리
				 * omega_N(k)를 곱한 후 결과의 처음 샘플에는
				 * 그냥 더하고 두번째 샘플에는 빼준다(-이므로)
				 */
				swap_x = x[i+step]* omega_x - y[i+step]*omega_y;
				swap_y = y[i+step]*omega_x + x[i+step]*omega_y;
 
				x[i+step] = x[i] - swap_x;
				y[i+step] = y[i] - swap_y;
 
				x[i] = x[i] + swap_x;
				y[i] = y[i] + swap_y;
 
			}
			/* 공통된 omega를 사용하는 샘플을 다 곱한 후
			 * omega를 factor의 각도만큼 증가시켜야 한다.
			 * omega(k+1) 의 값은 omega와 factor의 곱을 통해
			 * 구할 수 있다.
			 */
			swap_x = omega_x * factor_x - omega_y * factor_y;
			swap_y = omega_y * factor_x + omega_x * factor_y;
			omega_x = swap_x;
			omega_y = swap_y;
 
		}
		/* 한 단계가 끝났으므로 factor를 기존 값에서 각도가
		 * 절반인 값으로 구해야 한다.
		 */
		factor_y = sqrt((1.0+factor_x)/2.0);
		if( dir == 1)
			factor_y = -factor_y;
		factor_x = sqrt((1.0-factor_x)/2.0);
	}
	/* Inverse FFT 의 경우에는 sampling을 해줘야 한다.
	 * 다른 소스에서는 FFT의 값이 크게 나오므로 fft시에 sampling을
	 * 하고 inverse에서 안하는 경우도 있다 ( 상관없음 )
	 */
	if( dir == -1 )
	{
		for ( i = 0 ; i < N ; i++)
		{
			x[i] /= N;
			y[i] /= N;
		}
	}
 
	return TRUE;
}