#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#define Nmax 32768
void main_part(char fnameR[],char fnameW[]);
void FFT(int nn,double xr[Nmax],double xi[Nmax]);
/*---------------------------------------------------------------------------------*/

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];
    double xr[Nmax];
    double xi[Nmax];
    double cr[Nmax];
    double ci[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;xr[i-1]=0.0;xi[i-1]=0.0;
    }
    for(i=1;i<=ndata;i++){
        fgets(dat,sizeof dat,fin);data[i-1]=atof(strtok(dat,","));
        xr[i-1]=data[i-1];
    }
    fclose(fin);

    /*フーリエ変換・逆変換*/
    FFT(nn,xr,xi); /*フーリエ変換*/
    for(i=1;i<=nn;i++){ /*戻り値をデータ数で除す*/
        xr[i-1]=xr[i-1]/(double)nn;
        xi[i-1]=xi[i-1]/(double)nn;
        cr[i-1]=xr[i-1]; /*変換値実数部の記憶*/
        ci[i-1]=xi[i-1]; /*変換値虚数部の記憶*/
    }

    for(i=1;i<=nn;i++){ /*フーリエ逆変換用に虚数部の符号を反転*/
        xr[i-1]=xr[i-1];
        xi[i-1]=-xi[i-1];
    }
    FFT(nn,xr,xi); /*フーリエ逆変換*/

    /*データ書き出し*/
    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++){
        ck=sqrt(cr[i-1]*cr[i-1]+ci[i-1]*ci[i-1]);
        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",cr[i-1])   ;strcat(dat,",");strcat(dat,sv);
        sprintf(sv,"%e",ci[i-1])   ;strcat(dat,",");strcat(dat,sv);
        sprintf(sv,"%e",ck)        ;strcat(dat,",");strcat(dat,sv);
        sprintf(sv,"%e",xr[i-1])   ;strcat(dat,",");strcat(dat,sv);
        sprintf(sv,"%e",xi[i-1])   ;strcat(dat,",");strcat(dat,sv);
        strcat(dat,"\n");
        fputs(dat,fout);
    }
    fclose(fout);

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

void FFT(int nn,double xr[Nmax],double xi[Nmax])
{
    /**************************************************/
    /* 高速フーリエ変換・逆変換 */
    /**************************************************/
    /* nn  ：データ数（2の累乗）*/
    /* xr[]：入出力データ実数部 */
    /* xi[]：入出力データ虚数部 */

    int g,h,i,j,k,l,m,n,p,q;
    double a,b,xd;
    double s[(int)(nn/2)+1];
    double c[(int)(nn/2)+1];

    n=nn;
    i=0;j=0;k=0;l=0;p=0;h=0;g=0;q=0;
    m=(int)(log((double)n)/log(2.0)+1.0);
    a=0.0;
    b=M_PI*2.0/(double)n;
    for(i=0;i<=n/2;i++){
        s[i]=sin(a);
        c[i]=cos(a);
        a=a+b;
    }
    l=n;h=1;
    for(g=1;g<=m;g++){
        l=l/2;k=0;
        for(q=1;q<= h;q++){
            p=0;
            for(i=k;i<=l+k-1;i++){
                j=i+l;
                a=xr[i]-xr[j];
                b=xi[i]-xi[j];
                xr[i]=xr[i]+xr[j]; 
                xi[i]=xi[i]+xi[j];
                if(p==0){
                    xr[j]=a; 
                    xi[j]=b;
                }else{
                    xr[j]=a*c[p]+b*s[p];
                    xi[j]=b*c[p]-a*s[p];
                }
                p=p+h;
            }
            k=k+l+l;
        }
        h=h+h;
    }
    j=n/2;
    for(i=1;i<=n-1;i++){
        k=n;
        if(j<i){
            xd=xr[i];xr[i]=xr[j];xr[j]=xd;
            xd=xi[i];xi[i]=xi[j];xi[j]=xd;
        }
        k=k/2;
        while(j>= k){
            j=j-k;k=k/2;
            if(k==0)break;
        }
        j=j+k;
    }
}
/*---------------------------------------------------------------------------------*/