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 vcsCFFT { struct Complex { //構造体宣言 public double Real; public double Imag; } public partial class Form1:Form { public Form1() { InitializeComponent(); } //------------------------------------------------------------------- private void toolStripButton1_Click(object sender, EventArgs e) { int i, ndata, nn; double dt; System.IO.StreamReader sr; System.IO.StreamWriter sw; string dat; string[] sbuf; char delim = ','; string fnameR = ""; string fnameW = ""; //入出力ファイル指定 openFileDialog1.InitialDirectory = System.IO.Directory.GetCurrentDirectory(); if (openFileDialog1.ShowDialog() == DialogResult.OK) { fnameR = openFileDialog1.FileName; } if (saveFileDialog1.ShowDialog() == DialogResult.OK) { fnameW = saveFileDialog1.FileName; } sr = new System.IO.StreamReader(fnameR, 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]; Complex[] x = new Complex[nn]; Complex[] y = new Complex[nn]; for (i = 1; i <= nn; i++){ //データ数を2の累乗にセット data[i - 1] = 0.0; x[i - 1] = C_DOUBLE(0.0, 0.0); } for (i = 1; i <= ndata; i++){ dat = sr.ReadLine(); sbuf = dat.Split(delim); data[i - 1] = double.Parse(sbuf[0]); x[i - 1] = C_DOUBLE(data[i - 1], 0.0); } sr.Close(); CFFT(nn, x, -1);//フーリエ変換 for (i = 1; i <= nn; i++){ x[i - 1] = C_DIVIDED(x[i - 1], C_DOUBLE((double)nn, 0));//出力値をデータ数で除す y[i - 1] = x[i - 1]; } CFFT(nn, y, +1);//フーリエ逆変換 sw = new System.IO.StreamWriter(fnameW, 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 + "," + x[i - 1].Real.ToString("0.000"); dat = dat + "," + x[i - 1].Imag.ToString("0.000"); dat = dat + "," + D_ABS(x[i - 1]).ToString("0.000"); dat = dat + "," + y[i - 1].Real.ToString("0.000"); dat = dat + "," + y[i - 1].Imag.ToString("0.000"); sw.WriteLine(dat); } sw.Close(); } //------------------------------------------------------------------- //*********************************************************** //複素数計算のための構造体宣言とFunctionプロシージャ開始 //*********************************************************** //定義済み関数 //za=C_DOUBLE(a,b):複素数実数部と虚数部への数値の代入(Re(z)=a,Im(z)=b) //za=C_PLUS(z1,z2):加算(za=z1+z2) //za=C_MINUS(z1,z2):減算(za=z1-z2) //za=C_TIMES(z1,z2):乗算(za=z1*z2) //za=C_DIVIDED(z1,z2):除算(za=z1/z2) //da=D_ABS(z):絶対値(da=sqrt(Re(z)^2+Im(z)^2)) //za=C_CONJ(z):共役複素数(z=a+ib→za=a-ib) //za=C_SQRT(z):平方根(za=sqrt(z)) //za=C_EXP(z):exp関数(za=exp(z)) //------------------------------------------------------------------- private Complex C_DOUBLE(double a,double b) { //複素数実数部と虚数部への数値の代入 Complex zr; zr.Real=a; zr.Imag=b; return zr; } //------------------------------------------------------------------- private Complex C_PLUS(Complex z1,Complex z2) { //加算 Complex zr; zr.Real=z1.Real+z2.Real; zr.Imag=z1.Imag+z2.Imag; return zr; } //------------------------------------------------------------------- private Complex C_MINUS(Complex z1,Complex z2) { //減算 Complex zr; zr.Real=z1.Real-z2.Real; zr.Imag=z1.Imag-z2.Imag; return zr; } //------------------------------------------------------------------- private Complex C_TIMES(Complex z1,Complex z2) { //乗算 Complex zr; zr.Real=z1.Real*z2.Real-z1.Imag*z2.Imag; zr.Imag=z1.Real*z2.Imag+z2.Real*z1.Imag; return zr; } //------------------------------------------------------------------- private Complex C_DIVIDED(Complex z1,Complex z2) { //除算 Complex zr; zr.Real=0.0; zr.Imag=0.0; if(z2.Real==0.0&&z2.Imag==0.0){ MessageBox.Show("分母が0です!"); this.Close(); }else{ zr.Real=(z1.Real*z2.Real+z1.Imag*z2.Imag)/(z2.Real*z2.Real+z2.Imag*z2.Imag); zr.Imag=(z2.Real*z1.Imag-z1.Real*z2.Imag)/(z2.Real*z2.Real+z2.Imag*z2.Imag); } return zr; } //------------------------------------------------------------------- private double D_ABS(Complex z) { //絶対値 return Math.Sqrt(z.Real*z.Real+z.Imag*z.Imag); } //------------------------------------------------------------------- private Complex C_CONJ(Complex z) { //共役複素数 Complex zr; zr.Real=z.Real; zr.Imag=-z.Imag; return zr; } //------------------------------------------------------------------- private Complex C_SQRT(Complex z) { //平方根 Complex zr; if(z.Real==0.0&&z.Imag==0.0){ zr.Real=0.0; zr.Imag=0.0; }else{ double rr; double phi; rr=Math.Sqrt(z.Real*z.Real+z.Imag*z.Imag); phi=Math.Acos(z.Real/rr); if (z.Imag < 0.0) phi = 2.0 * Math.PI - phi; /* phiの範囲を0〜2πとする(acosの戻り値は0〜π) */ zr.Real=Math.Sqrt(rr)*Math.Cos(0.5*phi); zr.Imag=Math.Sqrt(rr)*Math.Sin(0.5*phi); } return zr; } //------------------------------------------------------------------- private Complex C_EXP(Complex z) { //EXP関数 Complex zr; zr.Real=Math.Exp(z.Real)*Math.Cos(z.Imag); zr.Imag=Math.Exp(z.Real)*Math.Sin(z.Imag); return zr; } //*********************************************************** //複素数計算のための構造体宣言とFunctionプロシージャ終了 //*********************************************************** //------------------------------------------------------------------- private void CFFT(int n,Complex[]x,int ind) { //*********************************************************** //複素数計算による高速フーリエ変換 //*********************************************************** //n:データ数(2の累乗) //x[]:入出力データ(複素数) //ind:-1=フーリエ変換,+1=フーリエ逆変換 Complex temp,theta; int i,j,k,m,kmax,istep; double pi=Math.PI; j=1; for(i=1;i<=n;i++){ if(i=2); j=j+m; } kmax=1; do{ istep=kmax*2; for(k=1;k<=kmax;k++){ theta=C_DOUBLE(0.0,pi*(double)(ind*(k-1))/(double)kmax); for(i=k;i<=n;i=i+istep){ j=i+kmax; temp=C_TIMES(x[j-1],C_EXP(theta)); x[j-1]=C_MINUS(x[i-1],temp); x[i-1]=C_PLUS(x[i-1],temp); } } kmax=istep; }while(kmax