/*--- fftImage.c --- ２次元ＦＦＴの結果を画像化する ---------------------------
		（X_SIZE，Y_SIZEが２のべき乗の場合に限る）
                 井上誠喜他，1999を元に，田中が翻案，2002．
-----------------------------------------------------------------------------*/
#include	<stdio.h>
#include	<stdlib.h>
#include	<math.h>

#define X_SIZE 256	/* 画像の横方向の大きさ		*/
#define Y_SIZE 256	/* 画像の縦方向の大きさ		*/
#define X_EXP	8	/* X_EXP = log2(X_SIZE)		*/
#define Y_EXP	8	/* Y_EXP = log2(Y_SIZE)		*/
#define HIGH   255	/* ２値画像のハイ・レベル	*/
#define LOW     0	/* ２値画像のロー・レベル	*/
#define	PI     3.1415926535

#define	OPT		1	/*　OPT = 1　光学的ＤＦＴ（直流分が中央）	*/
			    	/*　OPT = 0　通常のＤＦＴ（直流分が左端）	*/

int fft2 (double a_rl[Y_SIZE][X_SIZE], double a_im[Y_SIZE][X_SIZE], int inv);
void fft1core(double a_rl[], double a_im[], int length, 
	int ex, double sin_tbl[], double cos_tbl[], double buf[]);
void cstb(int length, int inv, double sin_tbl[], double cos_tbl[]);
void birv(double a[], int length, int ex, double b[]);
void rvmtx1(double a[Y_SIZE][X_SIZE], double b[X_SIZE][Y_SIZE], 
	int xsize, int ysize);
void rvmtx2(double a[X_SIZE][Y_SIZE], double b[Y_SIZE][X_SIZE], 
	int xsize, int ysize);

void main( void )
{
	FILE		*fpr,*fpw;
	FILE		*fopen();

	int		image_in[Y_SIZE][X_SIZE];
	double  	*ar;		/* データ実数部（入出力兼用）*/
	double  	*ai;		/* データ虚数部（入出力兼用）*/
	double 	norm, max;
	double	data;
	long	i, j;


	ar = (double *)malloc((size_t)Y_SIZE*X_SIZE*sizeof(double));
	ai = (double *)malloc((size_t)Y_SIZE*X_SIZE*sizeof(double));
	if ((ar == NULL) || (ai == NULL)) return -1;

	/* 原画像を読み込み，２次元ＦＦＴの入力データに変換する */
	fpw = fopen("frequency-myarea.gray","wb");
	fpr = fopen("myarea.gray","rb");

	for (i = 0; i < Y_SIZE; i++) {
		for (j = 0; j < X_SIZE; j++) {
			fseek(fpr,X_SIZE*i + j,0);
			image_in[i][j] = (int)fgetc(fpr);
			ar[X_SIZE*i + j] = image_in[i][j];
			/* printf ("%f\n", ar[X_SIZE*i + j]); */ 
			ai[X_SIZE*i + j] = 0.0;
		}
	}
	/* ２次元ＦＦＴを実行する */
	if (fft2((double (*)[X_SIZE])ar, (double (*)[X_SIZE])ai, 1) == -1) return -1;
	/* ＦＦＴ結果を画像化する */
	max = 0;
	for (i = 0; i < Y_SIZE; i++) {
		for (j = 0; j < X_SIZE; j++) {
			norm = ar[X_SIZE*i + j]*ar[X_SIZE*i + j] 
			     + ai[X_SIZE*i + j]*ai[X_SIZE*i + j];	/* 振幅成分計算 */
			if (norm != 0.0) norm = log(norm) / 2.0;
			else norm = 0.0;
			ar[X_SIZE*i + j] = (double)norm;
			if (norm > max) max = norm;
		}
	}
	for (i = 0; i < Y_SIZE; i++) {
		for (j = 0; j < X_SIZE; j++) {
			ar[X_SIZE*i + j] = (double)(ar[X_SIZE*i + j]*255 / max);
		}
	}
	/* ＦＦＴ結果を画像データに変換する */
	for (i = 0; i < Y_SIZE; i++) {
		for (j = 0; j < X_SIZE; j++) {
			data = ar[X_SIZE*i + j];
			if (data > 255) data = 255;
			if (data <   0) data = 0;
			fputc(data,fpw);
		}
	}
	free(ar);
	free(ai);
	return 0;
	fclose(fpr);
	fclose(fpw);
}

/*--- fft2 --- ２次元ＦＦＴの実行 ---------------------------------------------
		（X_SIZE，Y_SIZEが２のべき乗の場合に限る）
	a_rl:	データ実数部（入出力兼用）
	a_im:	データ虚数部（入出力兼用）
	inv:	1: ＤＦＴ，-1: 逆ＤＦＴ
-----------------------------------------------------------------------------*/
int fft2 (double a_rl[Y_SIZE][X_SIZE], double a_im[Y_SIZE][X_SIZE], int inv)
{
	double	*b_rl;		/* データ転置作業用配列（実数部）*/
	double	*b_im;		/* データ転置作業用配列（虚数部）*/
	double	*hsin_tbl;	/* 水平用SIN計算用テーブル		*/
	double	*hcos_tbl;	/* 水平用COS計算用テーブル		*/
	double	*vsin_tbl;	/* 垂直用SIN計算用テーブル		*/
	double	*vcos_tbl;	/* 垂直用COS計算用テーブル		*/
	double	*buf_x;		/* 作業用バッファ（水平方向）	*/
	double	*buf_y;		/* 作業用バッファ（垂直方向）	*/
	int		i;

	b_rl = (double *)calloc((size_t)X_SIZE*Y_SIZE, sizeof(double));
	b_im = (double *)calloc((size_t)X_SIZE*Y_SIZE, sizeof(double));
	hsin_tbl = (double *)calloc((size_t)X_SIZE, sizeof(double));
	hcos_tbl = (double *)calloc((size_t)X_SIZE, sizeof(double));
	vsin_tbl = (double *)calloc((size_t)Y_SIZE, sizeof(double));
	vcos_tbl = (double *)calloc((size_t)Y_SIZE, sizeof(double));
	buf_x = (double *)malloc((size_t)X_SIZE*sizeof(double));
	buf_y = (double *)malloc((size_t)Y_SIZE*sizeof(double));
	if ((b_rl == NULL) || (b_im == NULL)
		|| (hsin_tbl == NULL) || (hcos_tbl == NULL)
		|| (vsin_tbl == NULL) || (vcos_tbl == NULL)
		|| (buf_x == NULL) || (buf_y == NULL)) {
		return -1;
	}
	cstb(X_SIZE, inv, hsin_tbl, hcos_tbl);	/* 水平用SIN,COSテーブル作成	*/
	cstb(Y_SIZE, inv, vsin_tbl, vcos_tbl);	/* 垂直用SIN,COSテーブル作成	*/
	/* 水平方向のＦＦＴ */
	for (i = 0; i < Y_SIZE; i++) {
		fft1core(&a_rl[(long)i][0], &a_im[(long)i][0],
					X_SIZE, X_EXP, hsin_tbl, hcos_tbl, buf_x);
	}
	/* ２次元データの転置 */
	rvmtx1((double (*)[X_SIZE])a_rl, (double (*)[X_SIZE])b_rl, X_SIZE, Y_SIZE);
	rvmtx1((double (*)[X_SIZE])a_im, (double (*)[X_SIZE])b_im, X_SIZE, Y_SIZE);
	/* 垂直方向のＦＦＴ */
	for (i = 0; i < X_SIZE; i++) {
		fft1core(&b_rl[(long)Y_SIZE*i], &b_im[(long)Y_SIZE*i], 
					Y_SIZE, Y_EXP, vsin_tbl, vcos_tbl, buf_y);
	}
	/* ２次元データの転置 */
	rvmtx2((double (*)[Y_SIZE])b_rl, (double (*)[Y_SIZE])a_rl, X_SIZE, Y_SIZE);
	rvmtx2((double (*)[Y_SIZE])b_im, (double (*)[Y_SIZE])a_im, X_SIZE, Y_SIZE);
	free(b_rl);
	free(b_im);
	free(hsin_tbl);
	free(hcos_tbl);
	free(vsin_tbl);
	free(vcos_tbl);
}

/*--- fft1core --- １次元ＦＦＴの計算の核になる部分 ---------------------------
	a_rl:	データ実数部（入出力兼用）
	a_im:	データ虚数部（入出力兼用）
	ex:		データ個数を２のべき乗で与える(データ個数 = 2のex乗個）
	sin_tbl:	SIN計算用テーブル
	cos_tbl:	COS計算用テーブル
-----------------------------------------------------------------------------*/
void fft1core(double a_rl[], double a_im[], int length, 
	int ex, double sin_tbl[], double cos_tbl[], double buf[])
{
	int		i, j, k, w, j1, j2;
	int		numb, lenb, timb;
	double   xr, xi, yr, yi, nrml;

	if (OPT == 1) {
		for (i = 1; i < length; i+=2) {
			a_rl[i] = -a_rl[i];
			a_im[i] = -a_im[i];
		}
	}
	numb = 1;
	lenb = length;
	for (i = 0; i < ex; i++) {
		lenb /= 2;
		timb = 0;
		for (j = 0; j < numb; j++) {
			w = 0;
			for (k = 0; k < lenb; k++) {
				j1 = timb + k;
				j2 = j1 + lenb;
				xr = a_rl[j1];
				xi = a_im[j1];
				yr = a_rl[j2];
				yi = a_im[j2];
				a_rl[j1] = xr + yr;
				a_im[j1] = xi + yi;
				xr = xr - yr;
				xi = xi - yi;
				a_rl[j2] = xr*cos_tbl[w] - xi*sin_tbl[w];
				a_im[j2] = xr*sin_tbl[w] + xi*cos_tbl[w];
				w += numb;
			}
			timb += (2*lenb);
		}
		numb *= 2;
	}
	birv(a_rl, length, ex, buf);		/*　実数データの並べ換え　*/
	birv(a_im, length, ex, buf);		/*　虚数データの並べ換え　*/
	if (OPT == 1) {
		for (i = 1; i < length; i+=2) {
			a_rl[i] = -a_rl[i];
			a_im[i] = -a_im[i];
		}
	}
	nrml = (double)(1.0 / sqrt((double)length));
	for (i = 0; i < length; i++) {
		a_rl[i] *= nrml;
		a_im[i] *= nrml;
	}
}

/*--- cstb --- SIN,COSテーブル作成 --------------------------------------------
	length:		データ個数
	inv:		1: ＤＦＴ, -1: 逆ＤＦＴ
	sin_tbl:	SIN計算用テーブル
	cos_tbl:	COS計算用テーブル
-----------------------------------------------------------------------------*/
void cstb(int length, int inv, double sin_tbl[], double cos_tbl[])
{
	int		i;
	double  	xx, arg;

	xx = (double)(((-PI) * 2.0) / (double)length);
	if (inv < 0) xx = -xx;
	for (i = 0; i < length; i++) {
		arg = (double)i * xx;
		sin_tbl[i] = (double)sin(arg);
		cos_tbl[i] = (double)cos(arg);
	}
}

/*--- birv --- データの並べ換え -----------------------------------------------
	a:		データの配列
	length:	データ個数
	ex:		データ個数を２のべき乗で与える(length = 2のex乗個）
	b:		作業用バッファ
-----------------------------------------------------------------------------*/
void birv(double a[], int length, int ex, double b[])
{
	int	i, ii, k, bit;

	for (i = 0; i < length; i++) {
		for (k = 0, ii=i, bit=0; ; bit<<=1, ii>>=1) {
			bit = (ii & 1) | bit;
			if (++k == ex) break;
		}
		b[i] = a[bit];
	}
	for (i = 0; i < length; i++) 
		a[i] = b[i];
}

/*--- rvmtx1 --- ２次元データの転置 -------------------------------------------
	a:		２次元入力データ
	b:		２次元出力データ
	xsize:	水平データ個数
	ysize:	垂直データ個数
-----------------------------------------------------------------------------*/
void rvmtx1(double a[Y_SIZE][X_SIZE], double b[X_SIZE][Y_SIZE], 
	int xsize, int ysize)
{
	int  i, j;

	for (i = 0; i < ysize; i++) {
		for (j = 0; j < xsize; j++)
			b[j][i] = a[i][j];
	}
}

/*--- rvmtx2 --- ２次元データの転置 -------------------------------------------
	a:		２次元入力データ
	b:		２次元出力データ
	xsize:	水平データ個数
	ysize:	垂直データ個数
-----------------------------------------------------------------------------*/
void rvmtx2(double a[X_SIZE][Y_SIZE], double b[Y_SIZE][X_SIZE], 
	int xsize, int ysize)
{
	int i, j;

	for (i = 0; i < ysize; i++) {
		for (j = 0; j < xsize; j++)
			b[i][j] = a[j][i];
	}
}