#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#define Nmax 32768

struct Complex{
    /* 構造体宣言 */
    double Real;
    double Imag;
};

struct Complex C_DOUBLE(double a,double b);
struct Complex C_PLUS(struct Complex z1,struct Complex z2);
struct Complex C_MINUS(struct Complex z1,struct Complex z2);
struct Complex C_TIMES(struct Complex z1,struct Complex z2);
struct Complex C_DIVIDED(struct Complex z1,struct Complex z2);
double D_ABS(struct Complex z);
struct Complex C_CONJ(struct Complex z);
struct Complex C_SQRT(struct Complex z);
struct Complex C_EXP(struct Complex z);
struct Complex C_COS(struct Complex z);
struct Complex C_SIN(struct Complex z);
struct Complex C_COSH(struct Complex z);
struct Complex C_SINH(struct Complex z);
void main_part(char fnameR[],char fnameW[]);
void CFFT(int n,struct Complex x[Nmax],int ind);
/*-------------------------------------------------------------------*/

int main(int argc,char *argv[])
{
    char fnameR[50]; /*入力ファイル名*/
    char fnameW[50]; /*出力ファイル名*/

    strcpy(fnameR,argv[1]);
    strcpy(fnameW,argv[2]);

    main_part(fnameR,fnameW);

    printf("リターンキーを押してください");
    getchar();
    return 0;
}
/*---------------------------------------------------------------------------*/

void main_part(char fnameR[],char fnameW[])
{
    int i;
    double dt;
    int ndata;
    int nn;
    FILE *fin,*fout;
    char dat[256];
    char sv[20];
    double ck;
    double data[Nmax];
    struct Complex x[Nmax];
    struct Complex y[Nmax];

    /*データ読み込み*/
    fin=fopen(fnameR,"r");
    fgets(dat,sizeof dat,fin);
    fgets(dat,sizeof dat,fin);strtok(dat,",");dt=atof(strtok(NULL,","));
    fgets(dat,sizeof dat,fin);strtok(dat,",");ndata=atoi(strtok(NULL,","));
    nn = 2;
    do{
        nn = nn * 2;
    }while( nn < ndata * 1);

    for(i=1;i<=nn;i++){ 
        data[i-1]=0.0;
        x[i-1]=C_DOUBLE(0.0,0.0);
    }
    for(i=1;i<=ndata;i++){
        fgets(dat,sizeof dat,fin);data[i-1]=atof(strtok(dat,","));
        x[i-1]=C_DOUBLE(data[i-1],0.0);
    }
    fclose(fin);

    /*フーリエ変換・逆変換*/
    CFFT(nn,x,-1); /* フーリエ変換 */
    for(i=1;i<=nn;i++){
        x[i-1]=C_DIVIDED(x[i-1],C_DOUBLE((double)nn,0)); /* 出力値をデータ数で除す */
        y[i-1]=x[i-1];
    }
    CFFT(nn,y,+1); /* フーリエ逆変換 */

    /*データ書き出し*/
    fout=fopen(fnameW,"w");
    fputs("m：番号\n",fout);
    fputs("data：波形時刻歴数値\n",fout);
    fputs("k：番号\n",fout);
    fputs("Re(Ck)：複素フーリエ係数実数部\n",fout);
    fputs("Im(Ck)：複素フーリエ係数虚数部\n",fout);
    fputs("abs(Ck)：複素フーリエ係数絶対値\n",fout);
    fputs("Re(xm)：フーリエ逆変換実数部\n",fout);
    fputs("Im(xm)：フーリエ逆変換虚数部\n",fout);
    fputs("m,data,k,Re(Ck),Im(Ck),abs(Ck),Re(xm),Im(xm)\n",fout);
    for(i=1;i<=nn;i++){
        sprintf(sv,"%d",i-1)          ;strcpy(dat,sv);
        sprintf(sv,"%e",data[i-1])    ;strcat(dat,",");strcat(dat,sv);
        sprintf(sv,"%d",i-1)          ;strcat(dat,",");strcat(dat,sv);
        sprintf(sv,"%e",x[i-1].Real)  ;strcat(dat,",");strcat(dat,sv);
        sprintf(sv,"%e",x[i-1].Imag)  ;strcat(dat,",");strcat(dat,sv);
        sprintf(sv,"%e",D_ABS(x[i-1]));strcat(dat,",");strcat(dat,sv);
        sprintf(sv,"%e",y[i-1].Real)  ;strcat(dat,",");strcat(dat,sv);
        sprintf(sv,"%e",y[i-1].Imag)  ;strcat(dat,",");strcat(dat,sv);
        strcat(dat,"\n");
        fputs(dat,fout);
    }
    fclose(fout);

}
/*-------------------------------------------------------------------*/

void CFFT(int n,struct Complex x[Nmax],int ind)
{
    /************************************************************/
    /* 複素数計算による高速フーリエ変換 */
    /************************************************************/
    /* n：データ数（2の累乗） */
    /* x[]：入出力データ（複素数） */
    /* ind：-1=フーリエ変換，+1=フーリエ逆変換 */
    
    struct Complex temp,theta;
    int i,j,k,m,kmax,istep;
    double pi=M_PI;

    j=1;
    for(i=1;i<=n;i++){
        if(i<j){
            temp=x[j-1];
            x[j-1]=x[i-1];
            x[i-1]=temp;
        }
        m=(int)(n/2);
        do{
            if(j<=m)break;
            j=j-m;
            m=(int)(m/2);
        }while(m>=2);
        j=j+m;
    }
    kmax=1;
    do{
        istep=kmax*2;
        for(k=1;k<=kmax;k++){
            theta=C_DOUBLE(0.0,pi*(double)(ind*(k-1))/(double)kmax);
                for(i=k;i<=n;i=i+istep){
                    j=i+kmax;
                    temp=C_TIMES(x[j-1],C_EXP(theta));
                    x[j-1]=C_MINUS(x[i-1],temp);
                    x[i-1]=C_PLUS(x[i-1],temp);
                }
        }
        kmax=istep;
    }while(kmax<n);

}
/*-------------------------------------------------------------------*/

/**********************************************/
/* 複素数計算のための構造体宣言と関数定義開始 */
/**********************************************/
/* 定義済み関数 */
/* za=C_DOUBLE(a,b)   ：複素数実数部と虚数部への数値の代入（Re(z)=a，Im(z)=b）*/
/* za=C_PLUS(z1,z2)   ：加算（za=z1+z2）                  */
/* za=C_MINUS(z1,z2)  ：減算（za=z1-z2）                  */
/* za=C_TIMES(z1,z2)  ：乗算（za=z1*z2）                  */
/* za=C_DIVIDED(z1,z2)：除算（za=z1/z2）                  */
/* da=D_ABS(z)        ：絶対値（da=sqrt(Re(z)^2+Im(z)^2)）*/
/* za=C_CONJ(z)       ：共役複素数（z=a+ib→za=a-ib）     */
/* za=C_SQRT(z)       ：平方根（za=sqrt(z)）              */
/* za=C_EXP(z)        ：exp関数（za=exp(z)）              */
/* za=C_COS(z)        ：cos関数（za=cos(z)）              */
/* za=C_SIN(z)        ：sin関数（za=sin(z)）              */
/* za=C_COSH(z)       ：cosh関数（za=cosh(z)）            */
/* za=C_SINH(z)       ：sinh関数（za=sinh(z)）            */
/*----------------------------------------------------------------------*/

struct Complex C_DOUBLE(double a,double b)
{
    /* 複素数実数部と虚数部への数値の代入 */
    struct Complex zr;
    zr.Real=a;
    zr.Imag=b;
    return zr;
}
/*----------------------------------------------------------------------*/

struct Complex C_PLUS(struct Complex z1,struct Complex z2)
{
    /* 加算 */
    struct Complex zr;
    zr.Real=z1.Real+z2.Real;
    zr.Imag=z1.Imag+z2.Imag;
    return zr;
}
/*----------------------------------------------------------------------*/

struct Complex C_MINUS(struct Complex z1,struct Complex z2)
{
    /* 減算 */
    struct Complex zr;
    zr.Real=z1.Real-z2.Real;
    zr.Imag=z1.Imag-z2.Imag;
    return zr;
}

/*----------------------------------------------------------------------*/
struct Complex C_TIMES(struct Complex z1,struct Complex z2)
{
    /* 乗算 */
    struct Complex zr;
    zr.Real=z1.Real*z2.Real-z1.Imag*z2.Imag;
    zr.Imag=z1.Real*z2.Imag+z2.Real*z1.Imag;
    return zr;
}

/*----------------------------------------------------------------------*/
struct Complex C_DIVIDED(struct Complex z1,struct Complex z2)
{
    /* 除算 */
    struct Complex zr;
    zr.Real=0.0;
    zr.Imag=0.0;
    if(z2.Real==0.0&&z2.Imag==0.0){
        printf("分母が0です！");
        zr.Real=0.0;
        zr.Imag=0.0;
    }else{
        zr.Real=(z1.Real*z2.Real+z1.Imag*z2.Imag)/(z2.Real*z2.Real+z2.Imag*z2.Imag);
        zr.Imag=(z2.Real*z1.Imag-z1.Real*z2.Imag)/(z2.Real*z2.Real+z2.Imag*z2.Imag);
    }
    return zr;
}

/*----------------------------------------------------------------------*/
double D_ABS(struct Complex z)
{
    /* 絶対値 */
    return sqrt(z.Real*z.Real+z.Imag*z.Imag);
}

/*----------------------------------------------------------------------*/
struct Complex C_CONJ(struct Complex z)
{
    /* 共役複素数 */
    struct Complex zr;
    zr.Real=z.Real;
    zr.Imag=-z.Imag;
    return zr;
}

/*----------------------------------------------------------------------*/
struct Complex C_SQRT(struct Complex z)
{
    /* 平方根 */
    struct Complex zr;
    double rr;
    double phi;
    if(z.Real==0.0&&z.Imag==0.0){
        zr.Real=0.0;
        zr.Imag=0.0;
    }else{
        rr=sqrt(z.Real*z.Real+z.Imag*z.Imag);
        phi=acos(z.Real/rr);
        if(z.Imag<0.0)phi=2.0*M_PI-phi; /* phiの範囲を0〜2πとする（acosの戻り値は0〜π） */
        zr.Real=pow(rr,0.5)*cos(0.5*phi);
        zr.Imag=pow(rr,0.5)*sin(0.5*phi);
    }
    return zr;
}

/*----------------------------------------------------------------------*/
struct Complex C_EXP(struct Complex z)
{
    /* EXP関数 */
    struct Complex zr;
    zr.Real=exp(z.Real)*cos(z.Imag);
    zr.Imag=exp(z.Real)*sin(z.Imag);
    return zr;
}

/*----------------------------------------------------------------------*/
struct Complex C_COS(struct Complex z)
{
    /* cos関数 */
    struct Complex zr;
    zr.Real=0.5*(exp(-z.Imag)+exp(z.Imag))*cos(z.Real);
    zr.Imag=0.5*(exp(-z.Imag)-exp(z.Imag))*sin(z.Real);
    return zr;
}

/*----------------------------------------------------------------------*/
struct Complex C_SIN(struct Complex z)
{
    /* sin関数 */
    struct Complex zr;
    zr.Real=0.5*(exp(-z.Imag)+exp(z.Imag))*sin(z.Real);
    zr.Imag=-0.5*(exp(-z.Imag)-exp(z.Imag))*cos(z.Real);
    return zr;
}

/*----------------------------------------------------------------------*/
struct Complex C_COSH(struct Complex z)
{
    /* cosh関数 */
    struct Complex zr;
    zr.Real=0.5*(exp(z.Real)+exp(-z.Real))*cos(z.Imag);
    zr.Imag=0.5*(exp(z.Real)-exp(-z.Real))*sin(z.Imag);
    return zr;
}

/*----------------------------------------------------------------------*/
struct Complex C_SINH(struct Complex z)
{
    /* sinh関数 */
    struct Complex zr;
    zr.Real=0.5*(exp(z.Real)-exp(-z.Real))*cos(z.Imag);
    zr.Imag=0.5*(exp(z.Real)+exp(-z.Real))*sin(z.Imag);
    return zr;
}

/**********************************************/
/* 複素数計算のための構造体宣言と関数定義終了 */
/**********************************************/