// MAKETSP.cpp 
//
#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#include <math.h>
/*---------------------------------------------------------
ＦＦＴのライブラリはこちらのサイト様のものを借用しました。
http://www.kurims.kyoto-u.ac.jp/~ooura/index-j.html      
---------------------------------------------------------  */
#include "fft4g.c"

/************** 定数定義 *************************************/
//TSPの全体のサンプル数
#define TSP_DEFAULT_TSP_SIZE		1 << 17
//正ＴＳＰ波形出力ファイル
#define TSP_DEFAULT_FILENAME		"TSP.pcm"
//逆ＴＳＰ波形出力ファイル
#define TSP_DEFAULT_FILENAME_INV	"TSP-INV.pcm"
//ＴＳＰ長
#define TSP_LENGTH					50000
//円周率
#define PI							3.14159265358979
//許容誤差
#define GOOD_ERROR					((1.0 / (1 << 25)) / (1 << 25))
//二分法打ち切り処理回数
#define MAX_LOOP					4000

/************** 構造体定義 *************************************/
//複素数
struct Complex {
	double real;
	double image;
};

/************** 関数プロトタイプ *************************************/
int makeTSP(Complex* pFreqArea, int iTSPSize, bool bInv);
int makeFile(Complex* pTsp, int iTSPSize,char* pchFileName);



/************** メイン関数 *************************************/
int main(int argc, char* argv[])
{
	char pchFileName[1024];
	int  iTSPSize;

	iTSPSize = TSP_DEFAULT_TSP_SIZE; 
	Complex* pTSP = new Complex[iTSPSize];

	strcpy(pchFileName,TSP_DEFAULT_FILENAME);
	makeTSP(pTSP , iTSPSize, false);
	makeFile(pTSP, iTSPSize, pchFileName);

	strcpy(pchFileName,TSP_DEFAULT_FILENAME_INV);
	makeTSP(pTSP, iTSPSize, true);
	makeFile(pTSP, iTSPSize, pchFileName);

	delete [] pTSP;


	return 0;

}

/**************************************************************
F
掃出線グラフ定義関数
x 　　... 時刻 (0<=x<=1)
戻り値...周波数
**************************************************************/
double F(double x) {

/*４点を通る線分の組み合わせ*/
	double t[4],hz[4];

	t[0]=0;   hz[0] = 0;
	t[1]=0.3; hz[1] = 220;
	t[2]=0.8; hz[2] = 3500;
	t[3]=1.0; hz[3] = 22050;
	
	for (int i=0;i<3 && !(t[i] <= x && t[i+1] > x) ;i++)
		;
	
	x= ((i==0) ? 
		(x * hz[1] / t[1]) : 
		(hz[i+1] - hz[i]) 
			/ (t[i+1] - t[i]) 
			* (x  -  t[i]) + hz[i])/22050;

	return x;        

/*指数関数グラフ
	return (pow(120,x*2)-1)/5000+ x*.2; */

/*従来型
	return x;*/
}


/**************************************************************
makeTSP  
ＴＳＰ生成関数
pFreqArea  .. [OUT] 実時間波形を返すバュファ
iTSPSize   .. [IN]  ＴＳＰのサイズ
bInv       .. [IN]  falseならば正ＴＳＰ  trueならば逆ＴＳＰ
**************************************************************/
int makeTSP(Complex* pFreqArea, int iTSPSize, bool bInv) {
	
	int		i;
	double	Q ;	//ＴＳＰサイズ
	double	A;//振幅重み係数
	int		iHalfTSPSize;	//ＴＳＰサイズの半分
	double	iHalfTSPSizeInv;//iHalfTSPSizeの逆数
	double	direction;		//位相関数内の変数の符号
	double	*pPhaseGraph;	//位相回転量格納配列
	double	*pAmpGraph;		//振幅格納配列


	iHalfTSPSize = iTSPSize / 2;
	iHalfTSPSizeInv = 1.0 / iHalfTSPSize; 


	pPhaseGraph = new double[iHalfTSPSize+1];
	pAmpGraph = new double[iHalfTSPSize+1];


	double	x_y,bef_ix,bef_ix2, a, b, ix, nx, //逆関数計算用
			lpf_ix, lpf_ix2, lpf_ix3,bef_lpf_ix,lpf_start_cutoff;//フィルタ用
	
	bef_ix = bef_ix2 = ix =lpf_ix= lpf_ix2=bef_lpf_ix=lpf_ix3=0;
	lpf_start_cutoff  = 1;

	/*********************************************************
	* 0 ≦ n ≦ 1におけるF()の逆関数の積分関数を求める
	* F()は原点を通る単調増加する連続関数とする
	*********************************************************/
	pPhaseGraph[0] = 0; 
	for (i=1; i <= iHalfTSPSize;i++) {

		x_y = i * iHalfTSPSizeInv;

		nx = x_y - F(x_y);

		//ixが定まっていない時は仮にx_yと同じとする
		if (ix == 0) ix = x_y;
		
		if (nx != 0) { 
			//初期値の推定（単調増加なので、かならず1つ前のループの値よりも大きくなると仮定）
			if (nx > 0) {
				for (a = b = ix; F(b) < x_y; b =b * 1.1 + .01)
					a = b;
			}
			else  {
				a = bef_ix;
				b = x_y;
			}
			int cnt = 0;
			//二分法によりixを収束させる
			while (fabs(nx) > GOOD_ERROR && cnt++ < MAX_LOOP) {
				ix = (a + b) * 0.5;
				nx = x_y - F(ix);

				if (nx == 0.0) break;
				if (nx > 0.0) {
					a = ix;
				} else {
					b = ix;
				}			
			}	
		} else {
			ix = x_y;
		}
		
		//ローパスフィルタ
		if (i > 1) {
			lpf_ix += ( ix - lpf_ix) * lpf_start_cutoff;
			lpf_ix2 += ( lpf_ix - lpf_ix2) * lpf_start_cutoff;
			lpf_ix3 += ( lpf_ix2 - lpf_ix3) * lpf_start_cutoff;
			//フィルタの強さの設定（０Ｈｚ付近はフィルタをかけない）
			if (lpf_start_cutoff > 0.005) lpf_start_cutoff *= .99;
		} else
		{
			lpf_ix = lpf_ix2 = lpf_ix3 = ix;
		}
		
		//位相シフト量グラフ
		pPhaseGraph[i] = pPhaseGraph[i-1] + lpf_ix3;
		//増幅量グラフ
		pAmpGraph[i] = lpf_ix3 - bef_lpf_ix;
		if (i==1) pAmpGraph[0] = lpf_ix3; //[0]の時は[1]の値と同じ値を設定して近似させる

		bef_lpf_ix = lpf_ix3;
		bef_ix = ix;
	}
	
	//F()は単調増加するのでpPhaseGraph[iHalfTSPSize]が最大値
	Q = 2 * PI *  TSP_LENGTH  / ( pPhaseGraph[iHalfTSPSize]);

	direction = bInv ? 1.0 : -1.0;

	//周波数領域の設定
	for (i = 0; i <= iHalfTSPSize; i++) {
		//0 ~ iHalfTSPSize
		
		double jexp; 

		jexp =  pPhaseGraph[i]  * Q;

		//振幅はf()の微分の平方根
		A = sqrt(pAmpGraph[i]);

		if (bInv) A = 1.0 / A; //逆ＴＳＰの場合
		
		//位相関数
		pFreqArea[i].real = cos( direction * jexp ) * A;
		pFreqArea[i].image = -cos( PI / 2.0 + direction * jexp) * A;

		//ナイキスト周波数より上の帯域の設定
		if (i > 0 ) {  
			pFreqArea[iTSPSize - i].real  =   pFreqArea[i].real;
			pFreqArea[iTSPSize - i].image = - pFreqArea[i].image;	
		}
	}
	
	delete [] pPhaseGraph;
	delete [] pAmpGraph;
	
	//時間領域に展開
	int *ip = new int[3+(int)sqrt(iTSPSize)]; //バタフライワーク
	double *table = new double[iTSPSize/2]; //サインワーク
	ip[0] = 0;
	//FFTライブラリ呼び出し
	cdft(iTSPSize * 2, 1, (double*)pFreqArea, ip, table);

	delete [] ip;
	delete [] table;

	return 0;
}


int makeFile(Complex* pTsp, int iTSPSize,char* pchFileName) {

	FILE *fp;
	float data;
	fp = fopen(pchFileName,"wb");

	int i;
	
	for (i =0 ; i < iTSPSize ; i++) {
		int shift;

		shift =(i +  (iTSPSize - TSP_LENGTH)  / 2) % iTSPSize;
		data = (float)pTsp[i].real;
		fwrite(&data, sizeof(float), 1, fp);
	}

	fclose(fp);


	return 0;
}