// MAKETSP.cpp : コンソール アプリケーション用のエントリ ポイントの定義
//


#include "stdafx.h"
#include "string.h"
#include "stdlib.h"
#include "math.h"
#include "fft4g.c"

//定数定義
#define TSP_DEFAULT_TSP_SIZE	1 << 17
#define TSP_DEFAULT_FILENAME	"TSP1.pcm"
#define TSP_MAX_TSP_POW			17
#define TSP_MIN_TSP_POW			3
#define TSP_M_MUL				.45
#define PI						3.14159265358979

/*
enum FileType {
	type_16=1,
	type_32,
	type_f,
};
*/

//構造体定義
struct Complex {
	double real;
	double image;
};

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



//メイン関数
//呼出形式：プログラム名 [出力ファイル名] [ＴＳＰサイズ(2の指数)] [ファイル形式 (16 24 32 32F)]
int main(int argc, char* argv[])
{
	char pchFileName[1024];
	int  iTSPSize;

	//予定外の引数の数
	if (argc < 1 || argc > 4) {
		printf("引数の個数が違います。");
		return -1;
	}

	//ファイル名が指定されている場合
	if (argc >= 2) {
		strcpy(pchFileName,argv[1]);
	} else {
		strcpy(pchFileName,TSP_DEFAULT_FILENAME);
	}

	//TSPサイズが指定されている場合
	if (argc >= 3) {
		iTSPSize = atoi(argv[2]);

		if (iTSPSize > TSP_MAX_TSP_POW || iTSPSize < TSP_MIN_TSP_POW) {
			printf("ＴＳＰの長さは2の(%d〜%d)乗サンプルです。", TSP_MIN_TSP_POW, TSP_MAX_TSP_POW);
			return -1;
		}
		if ((iTSPSize % 2) == 1) {
			printf("ＴＳＰの長さは偶数です。");
			return -1;
		}
		iTSPSize = 2 << (iTSPSize-1);
	} else {
		iTSPSize = TSP_DEFAULT_TSP_SIZE; 
	}
/*	
	//ファイル形式が指定されている場合
	if (argc == 4) {
		if (strcmp(argv[3],"16") == 0) {
			file_type = type_16;
		} 
		else if (strcmp(argv[3],"32") == 0) {
			file_type = type_32;
		}
		else if (strcmp(argv[3],"F") == 0 || strcmp(argv[3],"f") == 0) {
			file_type = type_f;
		} 
	}
*/

	double* pTSP = new double[iTSPSize*2];
	
	makeTSP((Complex*)pTSP, iTSPSize);
	makeFile((Complex*)pTSP, iTSPSize, pchFileName);

	delete [] pTSP;

	return 0;

}



int makeTSP(Complex* pFreqArea, int iTSPSize) {
	
	int i;
	double dsize ;
	double sht;
	double powerr;
	dsize = iTSPSize;
	dsize =pow(iTSPSize,1.4);
//	dsize =iTSPSize * iTSPSize* iTSPSize;
	dsize = 1.0 / dsize;
	int m;
	m = (int)(iTSPSize * TSP_M_MUL);

	sht=0;
	powerr = 1;
	//初期化
	for(i=0;i<iTSPSize;i++)  pFreqArea[i].real = pFreqArea[i].image = 0;
	//周波数領域の設定
	for (i = 0; i <= iTSPSize / 2; i++) {
		//0 ~ iTSPSize / 2
		
		double jexp; 
		

		jexp = 4 * m * PI * pow((double)i, 1.4)  * dsize;


		powerr =  pow(i+1 , .31) / 10;
		pFreqArea[i].real = cos( jexp ) * powerr;
		pFreqArea[i].image = -cos( PI / 2.0 + jexp) * powerr;

		//iTSPSize / 2 + 1 ~ iTSPSize / 2
		if (i > 0 && i < iTSPSize / 2) {
			pFreqArea[iTSPSize - i].real  =   pFreqArea[i].real;
			pFreqArea[iTSPSize - i].image = - pFreqArea[i].image;	
		}
	}

		
	//時間領域に展開
	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 + (int)((0.5 - TSP_M_MUL ) * iTSPSize)) % iTSPSize;
		data = (float)pTsp[i].real;
		fwrite(&data, sizeof(float), 1, fp);
	}

	fclose(fp);


	return 0;
}