/***************************************************/
/* gccCALFSR.c                                     */
/***************************************************/
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#define Nmax 32768
void main_part(double band,char fnameR1[50],char fnameR2[50],char fnameW[50]);
void CALFSP(char fnameR[50],int *nn,double *dt,double fs0[Nmax],double fsp[Nmax],double band);
void FFT(int nn,double xr[Nmax],double xi[Nmax]);
void IACC(int nn,double dt,double ddy[Nmax],double dy[Nmax],double y[Nmax],
    double *ddymax,double *dymax,double *ymax);
void CRAC(int nn,double dt,double ddy[Nmax],double dy[Nmax],double y[Nmax],
    double ddymax,double dymax,double ymax);
void SWIN(int nn,double fs[(int)(Nmax/2)],double df,double band);
/*---------------------------------------------------------------------------*/

int main(int argc,char *argv[])
{
    char fnameR1[50]=""; /*分子入力ファイル名*/
    char fnameR2[50]=""; /*分母入力ファイル名*/
    char fnameW[50]=""; /*出力ファイル名*/
    double band;

    band=atof(argv[1]);      /*平滑化バンド幅*/
    strcpy(fnameR1,argv[2]); /*分子波形*/
    strcpy(fnameR2,argv[3]); /*分母波形*/
    strcpy(fnameW,argv[4]);  /*スペクトル比出力*/

    main_part(band,fnameR1,fnameR2,fnameW);

    return 0;
}
/*---------------------------------------------------------------------------*/

void main_part(double band,char fnameR1[50],char fnameR2[50],char fnameW[50])
{
    int i;
    double dt;
    char dat[256]="",sv[20]="";
    int nn,nd;
    FILE *fout;
    double frq[(int)(Nmax/2)+1],fsr[(int)(Nmax/2)+1];
    double fsp1[(int)(Nmax/2)+1],fs01[(int)(Nmax/2)+1];
    double fsp2[(int)(Nmax/2)+1],fs02[(int)(Nmax/2)+1];

    CALFSP(fnameR1,&nn,&dt,fs01,fsp1,band);
    CALFSP(fnameR2,&nn,&dt,fs02,fsp2,band);

    nd=(int)(nn/2);
    for(i=0;i<=nd;i++){
        frq[i]=(double)i/(double)nn/dt; /*振動数*/
        fsr[i]=fsp1[i]/fsp2[i]; /*スペクトル比*/
    }

    /*結果出力*/
    fout=fopen(fnameW,"w");
        strcpy(dat,"# in1=,");strcat(dat,fnameR1);strcat(dat,"\n");fputs(dat,fout);
        strcpy(dat,"# in2=,");strcat(dat,fnameR2);strcat(dat,"\n");fputs(dat,fout);
        strcpy(dat,"# out=,");strcat(dat,fnameW);strcat(dat,"\n");fputs(dat,fout);
        strcpy(dat,"# ");sprintf(sv,"%g",band);strcat(dat,"band=,");strcat(dat,sv);strcat(dat,"\n");fputs(dat,fout);
        strcpy(dat,"# ");sprintf(sv,"%d",nd+1);strcat(dat,"ndata=,");strcat(dat,sv);strcat(dat,"\n");fputs(dat,fout);
        strcpy(dat,"# ");strcat(dat,"frq,fs01,fs02,fsp1,fsp2,fsr\n");fputs(dat,fout);
    for(i=0;i<=nd;i++){
        sprintf(sv,"%g",frq[i]);strcpy(dat,sv);
        sprintf(sv,"%g",fs01[i]);strcat(dat,",");strcat(dat,sv);
        sprintf(sv,"%g",fs02[i]);strcat(dat,",");strcat(dat,sv);
        sprintf(sv,"%g",fsp1[i]);strcat(dat,",");strcat(dat,sv);
        sprintf(sv,"%g",fsp2[i]);strcat(dat,",");strcat(dat,sv);
        sprintf(sv,"%g",fsr[i]);strcat(dat,",");strcat(dat,sv);
        strcat(dat,"\n");
        fputs(dat,fout);
    }
    fclose(fout);


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

void CALFSP(char fnameR[50],int *nn,double *dt,double fs0[Nmax],double fsp[Nmax],double band)
{
    int i;
    double df;
    int ndata,nd;
    FILE *fin;
    char dat[256]="",sv[20]="";
    static double ddy[Nmax],dy[Nmax],y[Nmax];
    static double xr[Nmax],xi[Nmax];
    static double ddymax,dymax,ymax;

    /*データ入力*/
    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);

printf("nn=%d\n",*nn);

    for(i=1;i<=*nn;i++){ 
        ddy[i-1]=0.0;dy[i-1]=0.0;y[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);
            ddy[i-1]=atof(strtok(dat,","));
    }
    fclose(fin);

    /*平滑化フーリエスペクトル作成*/
    IACC(*nn,*dt,ddy,dy,y,&ddymax,&dymax,&ymax);     /*加速度時刻歴数値積分*/
    CRAC(*nn,*dt,ddy,dy,y,ddymax,dymax,ymax);        /*加速度時刻歴基線補正*/
    for(i=1;i<=*nn;i++){
        xr[i-1]=ddy[i-1];
        xi[i-1]=0.0;
    }
    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;
    }
    nd=(int)(*nn/ 2);
    for(i=1;i<=nd+1;i++){
        fsp[i-1]=*dt*(double)*nn*sqrt(xr[i-1]*xr[i-1]+xi[i-1]*xi[i-1]);
        fs0[i-1]=fsp[i-1];
    }
    df=1.0/ *dt/(double)*nn;
    SWIN(*nn,fsp,df,band); /*スペクトル平滑化*/

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

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;
    }
}
/*---------------------------------------------------------------------------*/

void IACC(int nn,double dt,double ddy[Nmax],double dy[Nmax],double y[Nmax],
    double *ddymax,double *dymax,double *ymax)
{
    /*加速度時刻歴数値積分               */
    /* nn     :加速度時刻歴データ総数    */
    /* dt     :時間間隔（入力値）        */
    /* ddy[]  :加速度時刻歴（入力値）    */
    /* dy[]   :速度時刻歴（計算出力）    */
    /* y[]    :変位時刻歴（計算出力値）  */
    /* ddymax :加速度最大値（計算出力値）*/
    /* dymax  :速度最大値（計算出力値）  */
    /* ymax   :変位最大値（計算出力値）  */

    int i;

    *ddymax=0.0;
    *dymax=0.0;
    *ymax=0.0;
    dy[0]=0.0;
    y[0]=0.0;
    for(i=1;i<=nn-1;i++){
        dy[i]=dy[i-1]+(ddy[i-1]+ddy[i])*dt/2.0;
        y[i]=y[i-1]+dy[i-1]*dt+(ddy[i-1]/3.0+ddy[i]/6.0)*dt*dt;
        if(*ddymax<fabs(ddy[i]))*ddymax=fabs(ddy[i]);
        if(*dymax<fabs(dy[i]))*dymax=fabs(dy[i]);
        if(*ymax<fabs(y[i]))*ymax=fabs(y[i]);
    }
}
/*---------------------------------------------------------------------------*/

void CRAC(int nn,double dt,double ddy[Nmax],double dy[Nmax],double y[Nmax],
    double ddymax,double dymax,double ymax)
{
    /*加速度時刻歴基線補正                       */
    /* nn     :加速度時刻歴データ総数            */
    /* dt     :時間間隔（入力値）                */
    /* ddy[]  :加速度時刻歴（入力値・書換出力値）*/
    /* dy[]   :速度時刻歴（入力値）              */
    /* y[]    :変位時刻歴（入力値）              */
    /* ddymax :加速度最大値（入力値）            */
    /* dymax  :速度最大値（入力値）              */
    /* ymax   :変位最大値（入力値）              */

    int i;
    double tt,t,sum,a1,a0,acmax,coef;

    tt=(double)(nn-1)*dt;
    t=0.0;
    for(i=0;i<=nn-1;i++){
        y[i]=y[i]*(3.0*tt-2.0*t)*t*t;
        t=t+dt;
    }
    sum=(y[0]+y[nn-1])/2.0;
    for(i=1;i<=nn-2;i++){
        sum=sum+y[i];
    }
    sum=sum*dt;
    a1=28.0/13.0/tt/tt*(2.0*dy[nn-1]-15.0/pow(tt,5.0)*sum);
    a0=dy[nn-1]/tt-a1/2.0*tt;
    t=0.0;
    acmax=0.0;
    for(i=0;i<=nn-1;i++){
        ddy[i]=ddy[i]-a0-a1*t;
        if(acmax<fabs(ddy[i]))acmax=fabs(ddy[i]);
        t=t+dt;
    }
    coef=ddymax/acmax;
    for(i=0;i<=nn-1;i++){
        ddy[i]=ddy[i]*coef;
    }
}
/*---------------------------------------------------------------------------*/

void SWIN(int nn,double fs[(int)(Nmax/2)],double df,double band)
{
    /**********************************/
    /*スペクトル平滑化                */
    /**********************************/
    /*nn：時刻歴全数                  */
    /*fs[]：フーリエ・スペクトル      */
    /*df：スペクトルの振動数間隔（Hz）*/
    /*band：バンド幅（Hz）            */
    /*nfold：スペクトル値総数         */

    double pi=M_PI;
    double tt,udf,dif,s;
    int i,j;
    int nfold,lmax,ll,ln,lt,le;
    double w[nn];
    double g[nn];
    double g1[nn];
    double g2[nn];

    nfold=(int)(nn/2)+1;

    if(band>0.0){
        tt=1.0/df;
        udf=280.0/151.0/band*df;
        lmax=(int)(2.0/udf)+1;
        if(udf<=0.5){
            w[0]=0.75*udf;
            for(i=2;i<=lmax;i++){
                dif=0.5*pi*(double)(i-1)*udf;
                w[i-1]=w[0]*pow(sin(dif)/dif,4.0);
            }
            g[0]=fs[0]*fs[0]/tt;
            for(i=2;i<=nfold-1;i++){
                g[i-1]=2.0*fs[i-1]*fs[i-1]/tt;
            }
            g[nfold-1]=fs[nfold-1]*fs[nfold-1]/tt;
            ll=lmax*2-1;
            ln=ll-1+nfold;
            lt=(ll-1)*2+nfold;
            le=lt-lmax+1;
            for(i=1;i<=lt;i++){
                g1[i-1]=0.0;
            }
            for(i=1;i<=nfold;i++){
                g1[ll-2+i]=g[i-1];
            }
            for(i=lmax;i<=le;i++){
                s=w[0]*g1[i-1];
                for(j=2;j<=lmax;j++){
                    s=s+w[j-1]*(g1[i-j]+g1[i+j-2]);
                }
                g2[i-1]=s;
            }
            for(j=2;j<=lmax;j++){
                g2[ll+j-2]=g2[ll+j-2]+g2[ll-j];
                g2[ln-j]=g2[ln-j]+g2[ln+j-2];
            }
            for(i=1;i<=nfold;i++){
                g[i-1]=g2[ll-2+i];
            }
            fs[0]=sqrt(g[0]*tt);
            for(i=2;i<=nfold-1;i++){
                fs[i-1]=sqrt(g[i-1]*tt/2.0);
            }
            fs[nfold-1]=sqrt(g[nfold-1]*tt);
        }
    }
}
/*---------------------------------------------------------------------------*/