using System; using System.Collections.Generic; using System.ComponentModel; using System.Data; using System.Drawing; using System.Linq; using System.Text; using System.Windows.Forms; namespace vcsFFT { public partial class Form1 : Form { public Form1() { InitializeComponent(); } //--------------------------------------------------------------------------------- private void toolStripButton1_Click(object sender, EventArgs e) { int i; double dt; int ndata; int nn; System.IO.StreamReader sr; System.IO.StreamWriter sw; char delim = ','; string dat; string[] sbuf; string fname1=""; string fname2=""; openFileDialog1.InitialDirectory = System.IO.Directory.GetCurrentDirectory(); if (openFileDialog1.ShowDialog() == DialogResult.OK) { fname1 = openFileDialog1.FileName; } sr = new System.IO.StreamReader(fname1, System.Text.Encoding.Default); dat = sr.ReadLine(); dat = sr.ReadLine(); sbuf = dat.Split(delim) ; dt = double.Parse(sbuf[1]); dat = sr.ReadLine(); sbuf = dat.Split(delim) ; ndata = int.Parse(sbuf[1]); nn = 2; do{ nn = nn * 2; }while( nn < ndata * 1); //必要に応じて後続の0を追加 double[] data=new double[nn]; double[] xr=new double[nn]; double[] xi=new double[nn]; double[] cr=new double[nn]; double[] ci=new double[nn]; 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++){ dat=sr.ReadLine();sbuf=dat.Split(delim); data[i-1]=double.Parse(sbuf[0]); xr[i-1]=data[i-1]; } sr.Close(); 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); //フーリエ逆変換 if(saveFileDialog1.ShowDialog()==DialogResult.OK){fname2=saveFileDialog1.FileName;} sw=new System.IO.StreamWriter(fname2,false,System.Text.Encoding.Default); sw.WriteLine("m:番号"); sw.WriteLine("data:波形時刻歴数値"); sw.WriteLine("k:番号"); sw.WriteLine("Re(Ck):複素フーリエ係数実数部"); sw.WriteLine("Im(Ck):複素フーリエ係数虚数部"); sw.WriteLine("abs(Ck):複素フーリエ係数絶対値"); sw.WriteLine("Re(xm):フーリエ逆変換実数部"); sw.WriteLine("Im(xm):フーリエ逆変換虚数部"); sw.WriteLine("m,data,k,Re(Ck),Im(Ck),abs(Ck),Re(xm),Im(xm)"); for(i=1;i<= nn;i++){ dat=(i-1).ToString("0"); dat=dat+","+data[i-1].ToString("0.000"); dat=dat+","+(i-1).ToString("0"); dat=dat+","+cr[i-1].ToString("0.000"); dat=dat+","+ci[i-1].ToString("0.000"); dat=dat+","+(Math.Sqrt(cr[i-1]*cr[i-1]+ci[i-1]*ci[i-1])).ToString("0.000"); dat=dat+","+xr[i-1].ToString("0.000"); dat=dat+","+xi[i-1].ToString("0.000"); sw.WriteLine(dat); } sw.Close(); } //--------------------------------------------------------------------------------- private void FFT(int nn, double[] xr, double[] xi) { //************************************************** //高速フーリエ変換・逆変換 //************************************************** //nn :データ数(2の累乗) //xr[]:入出力データ実数部 //xi[]:入出力データ虚数部 int g,h,i,j,k,l,m,n,p,q; double a,b,xd; n=nn; double[] s = new double[(int)(n / 2) + 1]; double[] c = new double[(int)(n / 2) + 1]; i=0;j=0;k=0;l=0;p=0;h=0;g=0;q=0; m=(int)(Math.Log((double)n)/Math.Log(2.0)+1.0); a=0.0; b=Math.PI*2.0/(double)n; for(i=0;i<=(int)(n/2);i++){ s[i]=Math.Sin(a); c[i]=Math.Cos(a); a=a+b; } l=n;h=1; for(g=1;g<=m;g++){ l=(int)(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=(int)(n/2); for(i=1;i<=n-1;i++){ k=n; if(j= k){ j=j-k;k=(int)(k/2); if(k==0)break; } j=j+k; } } //--------------------------------------------------------------------------------- } }