/***************************************************/ /* gccFSP.c */ /***************************************************/ #include #include #include #include #define Nmax 32768 void main_part(double band,char fnameR[50],char fnameW[50]); void FFT(int nn,double xr[Nmax],double xi[Nmax]); void SWIN(int nn,double fs[(int)(Nmax/2)],double df,double band); /*---------------------------------------------------------------------------*/ int main(int argc,char *argv[]) { char fnameR[50]=""; /*入力ファイル名*/ char fnameW[50]=""; /*出力ファイル名*/ double band; band=atof(argv[1]); /*平滑化バンド幅*/ strcpy(fnameR,argv[2]); /*入力波形*/ strcpy(fnameW,argv[3]); /*出力スペクトル*/ main_part(band,fnameR,fnameW); printf("リターンキーを押してください"); getchar(); return 0; } /*---------------------------------------------------------------------------*/ void main_part(double band,char fnameR[50],char fnameW[50]) { FILE *fin; FILE *fout; int i,nn,ndata,nd; double dt,df; char scom[256],dat[256]="",sv[20]=""; double ddy[Nmax],xr[Nmax],xi[Nmax]; double frq[(int)(Nmax/2)+1],fsp0[(int)(Nmax/2)+1],fsp1[(int)(Nmax/2)+1]; /*データ入力*/ fin=fopen(fnameR,"r"); fgets(dat,sizeof dat,fin);strcpy(scom,dat); 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= k){ j=j-k;k=k/2; if(k==0)break; } j=j+k; } } /*---------------------------------------------------------------------------*/ 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); } } } /*---------------------------------------------------------------------------*/