#include <stdio.h>
#include <math.h>
#include <stdlib.h>

#define STDIN  "<stdin>"
#define STDOUT "<stdout>"
#define TMPIN  ".in"
#define TMPOUT ".out"
#define N 2048 // 最大データ数
#ifdef M_PI
#define PI M_PI
#else
#define PI 3.14159265358979323846
#endif
#define Deg2Rad(x) (PI*x/180.)

//-----------------------------------------------------------------
void usage(char* cmd);
int file_data(char* in_file,double *x);
void out_data(char* out_file,double* psd,int data_num);
char* setgname(char* out_file);
void gnu_output(char* in_file,char* out_file);
void dft_spectrum(double *y, double *x, int data_num);
//-----------------------------------------------------------------

int main(int argc,char* argv[])
{
    int i,data_num;
    int use_simdata = 0,out_graph = 0;
    char *in_file = STDIN,*out_file = STDOUT;
    double x[N],psd[N/2];

    // read argument
    for(i=1;i<argc;i++){
        if(argv[i][0] == '-'){
            switch(argv[i][1]){
            case 'h': usage(argv[0]); return 0;
            case 'o':
                i++;
                if(argv[i] == NULL){
                    fprintf(stderr,"Output file not specified.\n");
                    usage(argv[0]);
                    return 0;
                }
                out_file = argv[i];
                break;
            case 'd':
                use_simdata = 1;
                break;
            case 'g':
                out_graph = 1;
                break;
            default:
                fprintf(stderr,"Unknown option '%c'\n",argv[i][1]);
                usage(argv[0]);
                return 0;
            }
        } else {
            in_file = argv[i];
        }
    }

    if(use_simdata) data_num = sim_data(x);
    else            data_num = file_data(in_file,x);
    if(data_num == 0) return;
    if(data_num)    dft_spectrum(x,psd,data_num);

    if(out_graph){
        if(!strcmp(in_file,STDIN)){
            in_file = TMPIN;
            out_data(in_file,x,data_num);
        }
        if(!strcmp(out_file,STDOUT))
            out_file = TMPOUT;
        out_data(out_file,psd,data_num/2);
        gnu_output(in_file,out_file);
        if(!strcmp(in_file,TMPIN))   remove(TMPIN);
        if(!strcmp(out_file,TMPOUT)) remove(TMPOUT);
    } else {
        out_data(out_file,psd,data_num/2);
    }
    
    return 0;
}

void usage(char* cmd)
{
    printf("usage: %s in_file [options]\n",cmd);
    printf("option:\n");
    printf("\t-h: display this information\n");
    printf("\t-o out_file: output to file.\n");
    printf("\t-d: use simulation data.\n");
    printf("\t-g: output to graph\n");
}

// ファイルのデータを格納する
int file_data(char* in_file,double *x)
{
    FILE* fp = stdin;
    int data_num = 0;
    char line[BUFSIZ];

    if(strcmp(in_file,STDIN))
        fp = fopen(in_file,"r");
    else
        printf("** input data (Ctrl-D end.)\n");
    if(fp == NULL) return 0;
    
    // ファイルからデータの読み取り
    while(fgets(line,BUFSIZ,fp)){
        sscanf(line,"%lf",&x[data_num]);
        data_num++;
        if(data_num > N){
            fprintf(stderr,"!!input data exceeded max read size\n");
            data_num = 0;
            break;
        }
    }
    fclose(fp);
    
    // データ数を返す
    return data_num;
}

// データを出力する
void out_data(char* out_file,double* psd,int data_num)
{
    FILE* fp = stdout;
    int k;
    if(strcmp(out_file,STDOUT))
        fp = fopen(out_file,"w+");
    else
        printf("** power spectrum\n");
    if(fp == NULL) return;
    
    for(k=0;k<data_num;k++)
        fprintf(fp,"%lf %lf\n",k/(double)data_num,psd[k]);
    fclose(fp);
}

// シミュレーションデータの作成
int sim_data(double* x)
{
    int k,data_num = N;
    for(k=0;k<N;k++){
        x[k] = sin(Deg2Rad(k)) + sin(Deg2Rad(5*k)) + sin(Deg2Rad(50*k));
        printf("%lf\n",x[k]);
    }
    return data_num;
}

char* setgname(char* out_file)
{
    int i;
    char* ofile = out_file;
    static char gname[256];
    if(!strcmp(out_file,TMPOUT))
        ofile = "output.png";
    for(i=strlen(ofile)-1;i>=0;i--){
        if(ofile[i]=='.') break;
    }
    if(i==-1) i = strlen(ofile);
    strncpy(gname,ofile,i);
    strcat(gname,".png");
    return gname;
}

// グラフ出力
void gnu_output(char* in_file,char* out_file)
{
    FILE* gfp = popen("gnuplot","w");
    fprintf(gfp,"set term png\n");
    fprintf(gfp,"set output \"%s\"\n",setgname(out_file));
    fprintf(gfp,"set multiplot\n");
    fprintf(gfp,"set size 1.0, 0.5\n");
    fprintf(gfp,"set origin 0, 0.5\n");
    fprintf(gfp,"plot \"%s\" w l 1\n",in_file);
    fprintf(gfp,"set origin 0, 0\n");
    fprintf(gfp,"plot \"%s\" t \"Power Spectrum by DFT\" w l 2\n",out_file);
    fflush(gfp);
    pclose(gfp);
}

// 離散フーリエ変換(Discrete Fourier Transform)
void dft_spectrum(double *x,double *psd,int data_num)
{
    int k,n;
    double Xr = 0.0,Xi = 0.0;
    // DFT
    for(k=0;k<data_num/2;k++){
        for(n=1;n<data_num;n++){
            Xr += x[n] * cos(2*PI*k*n/data_num);
            Xi += x[n] * sin(2*PI*k*n/data_num);
        }
        psd[k] = pow(Xr,2) + pow(Xi,2);
        Xr = 0.0; Xi = 0.0;
    }
}