#include "stdio.h"
#include "math.h"
typedef struct {//　■複素数構造体
    double Re;  //　　実数部(Real Part)
    double Im;  //　　虚数部(Image Part)
} Complex;
Complex setCom(double Re, double Im){	//■複素数設定
	Complex C; C.Re = Re; C.Im = Im;	return C;
}
Complex expJCom(double TH){				//■オイラーの公式
    return setCom(cos(TH),sin(TH));		//   exp(θ)=cosθ+j･sinθ
}										//   (j:虚数単位)
//----------------文字列出力---------------------------------------------
int strCom(Complex A, char str[]){		// ■文字列設定 
	  if(A.Im == 0) return(sprintf(str, "%lf", A.Re));
	  if(A.Im <  0) return(sprintf(str, "%lf - %lf j", A.Re, fabs(A.Im)));
	  else          return(sprintf(str, "%lf + %lf j", A.Re, A.Im));
}
void printCom(Complex A){				// ■文字列標準出力 
	  if(A.Im <  0) printf("%8.5lf - %8.5lf j", A.Re, fabs(A.Im));
	  else          printf("%8.5lf + %8.5lf j", A.Re, A.Im);
}
void fprintCom(FILE *fout, Complex A){	// ■文字列ファイル出力 
	  if(A.Im <  0) fprintf(fout,"%8.5lf - %8.5lf j", A.Re, fabs(A.Im));
	  else          fprintf(fout,"%8.5lf + %8.5lf j", A.Re, A.Im);
}
//----------------共通関数---------------------------------------------
double absCom(Complex A){				// ■複素数の絶対値
	if(A.Re == 0) return fabs(A.Im);	// sqrt(pow(A.Re, 2) + pow(A.Im, 2))の
    if(A.Im == 0) return fabs(A.Re);	// 指数部あふれを少なくするよう考慮
	double Re, Im; Re = fabs(A.Re); Im = fabs(A.Im);
	if(Im > Re) return Im * sqrt(1.0 + pow(Re / Im, 2));
	else        return Re * sqrt(1.0 + pow(Im / Re, 2));
}
double xDivAbs2(double X, double Y){	// ■ X/(pow(X,2)+pow(Y,2))計算
	if(X == 0) return 0;				// (pow(X,2)+pow(Y,2)の指数部あふれを
    if(Y == 0) return 1 / X;			// 少なくするよう考慮
	if(fabs(X) > fabs(Y)) return (1 / (X + Y / X * Y)); 
	double A = X / Y;     return (A /(1 + A * A) / Y);
}
//----------------複素数の加減算-----------------------------------------
Complex addCom(Complex A, Complex B){		/* 複素数 + 複素数 */
	return setCom(A.Re + B.Re, A.Im + B.Im);
}
Complex addCom(Complex A, double D){		/* 複素数 + 実数 　*/
	return setCom(A.Re + D, A.Im);
}
Complex addCom(double D, Complex A){		/* 実数 + 複素数　*/
	return setCom(A.Re + D, A.Im);
}
Complex subCom(Complex A, Complex B){		/* 複素数−複素数 　*/
	return setCom(A.Re - B.Re, A.Im - B.Im);
}
Complex subCom(Complex A, double D){		/* 複素数−実数　 　*/
	return setCom(A.Re - D, A.Im);
}
Complex subCom(double D, Complex A){		/* 実数−複素数　 　*/
	return setCom(D - A.Re, -A.Im);
}
Complex minusCom(Complex A){				/* −複素数　　　 　*/
	return setCom( -A.Re, -A.Im);
}
//----------------複素数の乗除算------------------------------------------
Complex multCom(Complex A, Complex B){/* 複素数 * 複素数 */
	return setCom(A.Re * B.Re - A.Im * B.Im, A.Im * B.Re + A.Re * B.Im);
}
Complex multCom(Complex A, double D){		/* 複素数　 * 実数 */
	return setCom(A.Re * D, A.Im * D);
}
Complex multCom(double D, Complex A){		/* 実数 * 複素数 */
	return setCom(A.Re * D, A.Im * D);
}
Complex revCom(Complex A){					/*   1 / 複素数    */
	return setCom(xDivAbs2(A.Re, A.Im), xDivAbs2(-A.Im, A.Re));
}
Complex divCom(Complex A, Complex B){		/* 複素数 / 複素数 */
	return multCom(A, revCom(B));
}
Complex divCom(Complex A, double D){		/* 複素数 / 実数 */
	return setCom(A.Re / D, A.Im / D);
}
Complex divCom(double D, Complex A){		/* 実数 / 複素数 */
	return multCom(D, revCom(A));
}
//----------------複素数関数------------------------------------------
Complex powCom(Complex A){					/* 複素数の2乗 */
	return setCom(A.Re * A.Re - A.Im * A.Im, A.Im * A.Re + A.Re * A.Im);
}
Complex powCom(Complex A, int N){			/* 複素数のN乗 */
	Complex P = setCom(1,0), X = A;
	for(int abN =abs(N); abN>0; abN >>= 1, X = powCom(X)) if(abN & 1) P=multCom(P,X);
	if(N<0) return revCom(P); else return P;
}
Complex sqrtCom(Complex A){					/* 平方根 */
	double S = absCom(A);
	return setCom(sqrt((S + A.Re)/2), sqrt((S - A.Re)/2));
}
Complex expCom(Complex A){					/* 自然指数 */
	double Ep=exp(A.Re), TH=A.Im;
	return setCom(cos(TH)*Ep, sin(TH)*Ep);
}
Complex logCom(Complex A){					/* 自然対数 */
	double Ap=absCom(A), D = A.Im / Ap;
	return setCom(log(Ap), (D>=-1 && D<=1) ? asin(D): 0);
}
Complex cosCom(Complex A){					/* cos */
	return setCom(cos(A.Re)*cosh(A.Im), -sin(A.Re)*sinh(A.Im));
}
Complex sinCom(Complex A){					/* sin */
	return setCom(sin(A.Re)*cosh(A.Im), cos(A.Re)*sinh(A.Im));
}
Complex tanCom(Complex A){					/* tan */
	double tanA  = tan(A.Re), tanhB  = tanh(A.Im) ;
	double tan2A = tanA*tanA, tanh2B = tanhB*tanhB;
	double D = 1 + tan2A * tanh2B;
	return setCom(tanA*(1-tanh2B)/D, tanhB*(tan2A+1)/D);
}
void DFT(int Max, Complex X[], Complex Y[], int Flag){//■フーリエ変換
	double mPi2 = -3.14159265358979 * 2, dMax = Max;
	if(Flag) mPi2 = -mPi2;
	for(int k = 0; k < Max; k++){
	  Y[k] = setCom(0,0);
	  for(int N = 0; N< Max; N++)
	     Y[k]= addCom(Y[k],multCom(X[N], expJCom(mPi2*(double)(N*k)/dMax)));
	}
	if(Flag) for(int k = 0; k < Max; k++) Y[k] = divCom(Y[k], dMax);
	Y[Max]=Y[0];
}
int main(){ //■フーリエ変換テストプログラム
	double Pi = 3.14159265358979 ; Complex X[181], Y1[181], Y2[181];
    for(int i=0;i<=181;i++){
		double TH=Pi*(double)i/90;
		X[i].Re = cos(TH)+ 0.5*cos(6*TH)+ 0.8 * sin(2*TH)+ 0.3 * sin(4*TH);
		X[i].Im = 0;
	}
	DFT(180, X, Y1, false);	DFT(180, Y1, Y2, true);
	for(int i=0;i<=180;i++) {
	  printf("\ni = %3d, ",i);printCom( X[i]); printf(", ");
	  printCom(Y1[i]); printf(", "); printCom(Y2[i]);
	}
	FILE *fout;
	if((fout=fopen("d:\\FurieTest.csv","wt"))==NULL) printf("Open Error");
	else{
		fprintf(fout,"番号, X(実数部), X(虚数部), Y1(実数部), Y1(虚数部), ");
		fprintf(fout,"Y2(実数部), Y2(虚数部)");
		for(int i=0;i<=180;i++) {
	       fprintf(fout, "\n%3d, %23.19lf, %23.19lf, ", i, X[i].Re, X[i].Im);
		   fprintf(fout, "%23.19lf, %23.19lf, %23.19lf, %23.19lf",
			              Y1[i].Re, Y1[i].Im, Y2[i].Re, Y2[i].Im);
		}
		fprintf(fout,"\n");	fclose(fout);
	}
	getchar();
	return 0;
}