#include #include #include #include #define Nmax 8192 /*gnuplot用構造体宣言*/ /*scom,sxjiku,syjiku :タイトル(凡例),x軸ラベル,y軸ラベル*/ /*slogx,xmin,xmax,dx,sxfmt :x軸種,x軸最小値,x軸最大値,x軸増分,x軸書式*/ /*slogy,ymin,ymax,dy,syfmt :y軸種,y軸最小値,y軸最大値,y軸増分,y軸書式*/ struct GPdata{ char scom[100];char sxjiku[100];char syjiku[100]; char slogx[3];double xmin;double xmax;double dx;char sxfmt[20]; char slogy[3];double ymin;double ymax;double dy;char syfmt[20]; }; /*関数プロトタイプ*/ void main_part(char fnameR[50],char fnameW[50],char fnameP[50]); void FFT(int nn,double xr[Nmax],double xi[Nmax]); void GPLOT(char fnameP[50],struct GPdata gp); /*---------------------------------------------------------------------------------*/ int main(int argc,char *argv[]) { char fnameR[50]; /*入力ファイル名*/ char fnameW[50]; /*出力ファイル名*/ char fnameP[50]; /*出力画像ファイル名*/ strcpy(fnameR,argv[1]); strcpy(fnameW,argv[2]); strcpy(fnameP,argv[3]); main_part(fnameR,fnameW,fnameP); printf("リターンキーを押してください"); getchar(); return 0; } /*---------------------------------------------------------------------------*/ void main_part(char fnameR[],char fnameW[],char fnameP[]) { int i,k; double dt; int ndata,nn,mm; FILE *fin,*fout; char dat[256]; char sv[20]; double datax[Nmax]; double datay[Nmax]; double func[Nmax]; double xr[Nmax]; double xi[Nmax]; double ak[(int)(Nmax/2)+1]; double bk[(int)(Nmax/2)+1]; struct GPdata gp={ "","","", "",0.0,0.0,0.0,"", "",0.0,0.0,0.0,"" }; /*データ入力*/ fin=fopen(fnameR,"r"); fgets(dat,sizeof dat,fin); strcpy(sv,strtok(dat,",")); strncat(gp.scom,sv,strlen(sv)-1); fgets(dat,sizeof dat,fin); strcpy(sv,strtok(dat,",")); strncat(gp.sxjiku,sv,strlen(sv)); strcpy(sv,strtok(NULL,",")); strncat(gp.syjiku,sv,strlen(sv)-1); fgets(dat,sizeof dat,fin); strcpy(gp.slogx,strtok(dat,",")); gp.xmin=atof(strtok(NULL,",")); gp.xmax=atof(strtok(NULL,",")); if(strncmp(gp.slogx,"N",1)==0){ gp.dx=atof(strtok(NULL,",")); strcpy(sv,strtok(NULL,",")); strncat(gp.sxfmt,sv,strlen(sv)-1); } fgets(dat,sizeof dat,fin); strcpy(gp.slogy,strtok(dat,",")); gp.ymin=atof(strtok(NULL,",")); gp.ymax=atof(strtok(NULL,",")); if(strncmp(gp.slogy,"N",1)==0){ gp.dy=atof(strtok(NULL,",")); strcpy(sv,strtok(NULL,",")); strncat(gp.syfmt,sv,strlen(sv)-1); } fgets(dat,sizeof dat,fin); ndata=atof(strtok(dat,",")); mm=atoi(strtok(NULL,",")); printf("ndata=%d\n",ndata+1); for(i=0;i<=ndata-1;i++){ fgets(dat,sizeof dat,fin); datax[i]=atof(strtok(dat,",")); datay[i]=atof(strtok(NULL,",")); } fclose(fin); /*フーリエ変換*/ nn = 2; do{ nn = nn * 2; }while( nn < ndata * 1); for(i=1;i<=nn;i++){ xr[i-1]=0.0;xi[i-1]=0.0; } for(i=1;i<=ndata;i++){ xr[i-1]=datay[i-1]-datay[ndata-1]; } 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; } /*回帰推定値*/ dt=datax[1]-datax[0]; for(i=0;i<=(int)(nn/2);i++){ ak[i]= 2.0*xr[i]; bk[i]=-2.0*xi[i]; } for(i=0;i<=ndata-1;i++){ func[i]=0.5*ak[0]; for(k=1;k<=mm;k++){ func[i]=func[i]+ak[k]*cos(2.0*M_PI*(double)(k*i)/(double)nn/dt); func[i]=func[i]+bk[k]*sin(2.0*M_PI*(double)(k*i)/(double)nn/dt); } func[i]=func[i]+datay[ndata-1]; } /*結果出力*/ fout=fopen(fnameW,"w"); strcat(gp.scom,"\n"); fputs(gp.scom,fout); sprintf(sv,"%d",ndata);strcpy(dat,sv);strcat(dat,"\n"); fputs(dat,fout); fputs("x,y,func\n",fout); for(i=0;i<=ndata-1;i++){ sprintf(sv,"%e",datax[i]) ;strcpy(dat,sv); sprintf(sv,"%e",datay[i]) ;strcat(dat,",");strcat(dat,sv); sprintf(sv,"%e",func[i]) ;strcat(dat,",");strcat(dat,sv); strcat(dat,"\n"); fputs(dat,fout); } fclose(fout); /*入力データプロット用一時ファイル書き出し*/ fout=fopen("temp0.prn","w"); for(i=0;i<=ndata-1;i++){ fprintf(fout,"%f %f\n",datax[i],datay[i]); } fclose(fout); /*回帰曲線プロット用一時ファイル書き出し*/ fout=fopen("temp1.prn","w"); for(i=0;i<=ndata-1;i++){ fprintf(fout,"%f %f\n",datax[i],func[i]); } fclose(fout); /*残差プロット用一時ファイル書き出し*/ fout=fopen("temp2.prn","w"); for(i=0;i<=ndata-1;i++){ fprintf(fout,"%f %f\n",datax[i],datay[i]-func[i]); } fclose(fout); /*gnuplotによる作図*/ GPLOT(fnameP,gp); } /*---------------------------------------------------------------------------------*/ 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= k){ j=j-k;k=k/2; if(k==0)break; } j=j+k; } } /*---------------------------------------------------------------------------------*/ void GPLOT(char fnameP[50],struct GPdata gp) { FILE *fp; /****************************/ /*gnuplot制御用ファイル作成*/ /****************************/ fp=fopen("gpl_temp.txt","w"); /*出力ファイル定義*/ fprintf(fp,"set terminal png\n"); fprintf(fp,"set output \"%s\"\n",fnameP); /*軸ラベル定義*/ fprintf(fp,"set xlabel \"%s\"\n",gp.sxjiku); fprintf(fp,"set ylabel \"%s\"\n",gp.syjiku); /*x軸定義*/ if(strncmp(gp.slogx,"L",1)==0){ fprintf(fp,"set logscale x\n"); fprintf(fp,"set mxtics 10\n"); } if(strncmp(gp.slogx,"N",1)==0){ fprintf(fp,"set xtics %f\n",gp.dx); fprintf(fp,"set format x %s\n",gp.sxfmt); } /*y軸定義*/ if(strncmp(gp.slogy,"L",1)==0){ fprintf(fp,"set logscale y\n"); fprintf(fp,"set mytics 10\n"); } if(strncmp(gp.slogy,"N",1)==0){ fprintf(fp,"set ytics %f\n",gp.dy); fprintf(fp,"set format y %s\n",gp.syfmt); } /*軸範囲定義*/ fprintf(fp,"set xrange [%f:%f]\n",gp.xmin,gp.xmax); fprintf(fp,"set yrange [%f:%f]\n",gp.ymin,gp.ymax); /*軸グリッド描画設定*/ fprintf(fp,"set grid xtics ytics mxtics mytics\n"); /*プロット実行(プロットは各種定義が終わって最後に実行)*/ fprintf(fp,"plot \\\n"); fprintf(fp,"\"temp0.prn\" with points pointtype 1 title \"Input data\", \\\n"); fprintf(fp,"\"temp1.prn\" with lines linetype 3 title \"Estimated\", \\\n"); fprintf(fp,"\"temp2.prn\" with points pointtype 2 title \"Residual\"\n"); fclose(fp); /****************************/ /*gnuplot制御*/ /****************************/ fp=popen("gnuplot.exe","w"); fprintf(fp,"load \"gpl_temp.txt\"\n"); fflush(fp); pclose(fp); } /*---------------------------------------------------------------------------*/