/*	
 *	clipdct.c
 *
 *	離散コサイン変換
 *
 *	CLiP - Common Library for P/ECE
 *	Copyright (C) 2001-2003 Naoyuki Sawa
 *
 *	* Mon Oct 20 04:14:00 JST 2003 Naoyuki Sawa
 *	- 新規作成。
 */
#include "clip.h"

/****************************************************************************
 *	DCT_REAL
 ****************************************************************************/

static void _dct_real(int N, const double* T, const double* src, double* dst);

DCT_REAL*
dct_real_new(int N)
{
	DCT_REAL* dct;
	int opt, k, n;
	double pi, ck;

	/* 次数が負値で指定されていたら、「扱いやすいDCT」を選択します。 */
	opt = N < 0;
	N = abs(N);

	/* DCT_REAL構造体のメモリを確保します。 */
	dct = (DCT_REAL*)calloc(1, sizeof(DCT_REAL));
	if(!dct) DIE();

	/* 次数を格納します。 */
	dct->N = N;

	/* エリアスセーフ用バッファを確保します。 */
	dct->v = (double*)calloc(N, sizeof(double));
	if(!dct->v) DIE();

	/* 変換行列のメモリを確保します。 */
	dct->T[0] = (double*)calloc(N * N, sizeof(double));
	if(!dct->T[0]) DIE();
	dct->T[1] = (double*)calloc(N * N, sizeof(double));
	if(!dct->T[1]) DIE();

	/* 変換行列を初期化します。 */
	pi = acos(-1);
	for(k = 0; k < N; k++) {
		if(!opt) {
			/* 正式なDCT。変換行列は正規直行。FDCT後の係数は各cos成分の振幅の√N倍(k=0)または√(N/2)倍(k<>0)です。 */
			ck = k == 0 ? sqrt(1.0 / N)	/* k = 0        */
			            : sqrt(2.0 / N);	/* k = 1〜(N-1) */
			for(n = 0; n < N; n++) {
				dct->T[0][k * N + n] = ck * cos(((2 * n + 1) * k * pi) / (2 * N)); /* FDCT */
				dct->T[1][n * N + k] = ck * cos(((2 * n + 1) * k * pi) / (2 * N)); /* IDCT */
			}
		} else {
			/* 扱いやすいDCT。変換行列は正規化されていない。FDCT後の係数は各cos成分の振幅そのもの(1倍)となります。 */
			ck = k == 0 ? 1.0 / N		/* k = 0        */
			            : 2.0 / N;		/* k = 1〜(N-1) */
			for(n = 0; n < N; n++) {
				dct->T[0][k * N + n] = ck * cos(((2 * n + 1) * k * pi) / (2 * N)); /* FDCT */
				dct->T[1][n * N + k] =      cos(((2 * n + 1) * k * pi) / (2 * N)); /* IDCT */
			}
		}
	}

	return dct;
}

void
dct_real_free(DCT_REAL* dct)
{
	/* エリアスセーフ用バッファを開放します。 */
	free(dct->v);

	/* 変換行列のメモリを開放します。 */
	free(dct->T[0]);
	free(dct->T[1]);

	/* DCT構造体のメモリを開放します。 */
	free(dct);
}

/* Forward Discrete Cosine Transform */
void
fdct_real(DCT_REAL* dct, const double* src, double* dst)
{
	if(src == dst) src = (const double*)memcpy(dct->v, src, sizeof(double) * dct->N);
	_dct_real(dct->N, dct->T[0/*FDCT*/], src, dst);
}

/* Inverse Discrete Cosine Transform */
void
idct_real(DCT_REAL* dct, const double* src, double* dst)
{
	if(src == dst) src = (const double*)memcpy(dct->v, src, sizeof(double) * dct->N);
	_dct_real(dct->N, dct->T[1/*IDCT*/], src, dst);
}

static void
_dct_real(int N, const double* T, const double* src, double* dst)
{
	int k, n;
	double tmp;
	const double* v;

	for(k = 0; k < N; k++) {
		v = src;
		tmp = 0;
		for(n = 0; n < N; n++) {
			tmp += *T++ * *v++;
		}
		*dst++ = tmp;
	}
}

/****************************************************************************
 *	DCT_HALF
 ****************************************************************************/

#define DCT_HALF_SHIFT	8			/* 変換行列の小数部ビット数(調整可) */
#define DCT_HALF_SCALE	(1 << DCT_HALF_SHIFT)	/* 変換行列の固定小数スケーリング値 */

static void _dct_half(int N, const short* T, const short* src, short* dst) __attribute__((unused));

DCT_HALF*
dct_half_new(int N)
{
	DCT_HALF* dct;
	int opt, k, n;
	double pi, ck;

	/* 次数が負値で指定されていたら、「扱いやすいDCT」を選択します。 */
	opt = N < 0;
	N = abs(N);

	/* DCT_HALF構造体のメモリを確保します。 */
	dct = (DCT_HALF*)calloc(1, sizeof(DCT_HALF));
	if(!dct) DIE();

	/* 次数を格納します。 */
	dct->N = N;

	/* エリアスセーフ用バッファを確保します。 */
	dct->v = (short*)calloc(N, sizeof(short));
	if(!dct->v) DIE();

	/* 変換行列のメモリを確保します。 */
	dct->T[0] = (short*)calloc(N * N, sizeof(short));
	if(!dct->T[0]) DIE();
	dct->T[1] = (short*)calloc(N * N, sizeof(short));
	if(!dct->T[1]) DIE();

	/* 変換行列を初期化します。 */
	pi = acos(-1);
	for(k = 0; k < N; k++) {
		if(!opt) {
			/* 正式なDCT。変換行列は正規直行。FDCT後の係数は各cos成分の振幅の√N倍(k=0)または√(N/2)倍(k<>0)です。 */
			ck = k == 0 ? sqrt(1.0 / N)	/* k = 0        */
			            : sqrt(2.0 / N);	/* k = 1〜(N-1) */
			for(n = 0; n < N; n++) {
				dct->T[0][k * N + n] = (short)(ck * cos(((2 * n + 1) * k * pi) / (2 * N)) * DCT_HALF_SCALE); /* FDCT(DCT-II ) */
				dct->T[1][n * N + k] = (short)(ck * cos(((2 * n + 1) * k * pi) / (2 * N)) * DCT_HALF_SCALE); /* IDCT(DCT-III) */
			}
		} else {
			/* 扱いやすいDCT。変換行列は正規化されていない。FDCT後の係数は各cos成分の振幅そのもの(1倍)となります。 */
			ck = k == 0 ? 1.0 / N		/* k = 0        */
			            : 2.0 / N;		/* k = 1〜(N-1) */
			for(n = 0; n < N; n++) {
				dct->T[0][k * N + n] = (short)(ck * cos(((2 * n + 1) * k * pi) / (2 * N)) * DCT_HALF_SCALE); /* FDCT */
				dct->T[1][n * N + k] = (short)(     cos(((2 * n + 1) * k * pi) / (2 * N)) * DCT_HALF_SCALE); /* IDCT */
			}
		}
	}

	return dct;
}

void
dct_half_free(DCT_HALF* dct)
{
	/* エリアスセーフ用バッファを開放します。 */
	free(dct->v);

	/* 変換行列のメモリを開放します。 */
	free(dct->T[0]);
	free(dct->T[1]);

	/* DCT構造体のメモリを開放します。 */
	free(dct);
}

/* Forward Discrete Cosine Transform */
void
fdct_half(DCT_HALF* dct, const short* src, short* dst)
{
	if(src == dst) src = (const short*)memcpy(dct->v, src, sizeof(short) * dct->N);
	_dct_half(dct->N, dct->T[0/*FDCT*/], src, dst);
}

/* Inverse Discrete Cosine Transform */
void
idct_half(DCT_HALF* dct, const short* src, short* dst)
{
	if(src == dst) src = (const short*)memcpy(dct->v, src, sizeof(short) * dct->N);
	_dct_half(dct->N, dct->T[1/*IDCT*/], src, dst);
}

static void
_dct_half(int N, const short* T, const short* src, short* dst)
{
#if 0
//	/* 低速版 */
//	int k, n;
//	int tmp;
//	const short* v;
//	for(k = 0; k < N; k++) {
//		v = src;
//		tmp = 0;
//		for(n = 0; n < N; n++) {
//			tmp += *T++ * *v++;
//		}
//		*dst++ = tmp >> DCT_HALF_SHIFT;
//	}
#else
	/* 高速版 */
	asm("
	; 注意!!
	; mac命令は自動的に%ahr:%alrをリセット★しません!!★。
	; mac命令実行前に明示的に0をロードする必要があります。
	;
	ld.w %alr, %r8		; 初回のmac用の%alrリセット
	ld.w %r7, %r12		; %r7 <= N (ループカウンタ)
	ld.w %r6, %r13		; %r6 <= T
_dct_half_L10:
	;
	ld.w %r5, %r14		; %r5 <= src
	ld.w %r4, %r12		; %r4 <= N
	mac %r4			; 畳み込み演算
	ld.w %r4, %alr
	xsra %r4, 8		; 固定小数桁合わせ
	ld.h [%r15]+, %r4	; 結果格納
	;
	xsub %r7, %r7, 1	; ループ処理
	xjrne.d _dct_half_L10
	ld.w %alr, %r8		; (delay) 二回目以降のmac用の%alrリセット
	;
	ret
	");
#endif
}

