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 vcsFSP { public partial class Form1 : Form { public Form1() { InitializeComponent(); } private void Form1_Load(object sender, EventArgs e) { toolStripLabel1.Text = "バンド幅"; toolStripTextBox1.Text = "0.4"; } private void toolStripButton1_Click(object sender, EventArgs e) { //フーリエスペクトル比算定 //分子・分母にするデータのサンプリングピッチ・データ数は等しくとる string fnameR = ""; string fnameW = ""; System.IO.StreamReader sr; System.IO.StreamWriter sw; string dat = ""; string[] sbuf; char delim = ','; double dt, band; int i,nn,ndata,nd; double df; string scom=""; int ndmax = (int)Math.Pow(2, 16); double[] fsp0 = new double[ndmax + 1]; double[] fsp1 = new double[ndmax + 1]; double[] fq = new double[ndmax + 1]; for (i = 0; i <= ndmax; i++) { fsp0[i] = 0.0; fsp1[i] = 0.0; fq[i] = 0.0; } //平滑化バンド幅設定 band = double.Parse(toolStripTextBox1.Text); //入出力ファイル指定 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(); sbuf = dat.Split(delim); scom = sbuf[0]; 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); double[] ddy = new double[nn]; double[] xr = new double[nn]; double[] xi = new double[nn]; for (i = 1; i <= nn; i++) { //データ数を2の階乗にセット ddy[i - 1] = 0.0; } for (i = 1; i <= ndata; i++) { dat = sr.ReadLine(); sbuf = dat.Split(delim); ddy[i - 1] = double.Parse(sbuf[0]); } sr.Close(); //平滑化フーリエスペクトル作成 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++) { fsp0[i - 1] = dt * (double)nn * Math.Sqrt(xr[i - 1] * xr[i - 1] + xi[i - 1] * xi[i - 1]); fsp1[i - 1] = fsp0[i - 1]; } df = 1.0 / dt / (double)nn; SWIN(nn, fsp1, df, band); //スペクトル平滑化 for (i = 0; i <= nd; i++) { fq[i] = (double)i / (double)nn / dt; //振動数 } //結果書き出し sw = new System.IO.StreamWriter(fnameW, false, System.Text.Encoding.Default); dat = scom; sw.WriteLine(dat); dat = "バンド幅" + "," + band.ToString(); sw.WriteLine(dat); dat = "データ数" + "," + (nd + 1).ToString(); sw.WriteLine(dat); dat = "振動数,原波形SP,平滑化SP"; sw.WriteLine(dat); for (i = 0; i <= nd; i++) { dat = fq[i].ToString("E"); dat = dat + "," + fsp0[i].ToString("E"); dat = dat + "," + fsp1[i].ToString("E"); sw.WriteLine(dat); } sw.Close(); MessageBox.Show("完了"); } //-------------------------------------------------------------------- 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 < i) { xd = xr[i]; xr[i] = xr[j]; xr[j] = xd; xd = xi[i]; xi[i] = xi[j]; xi[j] = xd; } k = (int)(k / 2); while (j >= k) { j = j - k; k = (int)(k / 2); if (k == 0) break; } j = j + k; } } //--------------------------------------------------------------------------------- private void SWIN(int nn, double[] fs, double df, double band) { //************************************** //スペクトル平滑化 //************************************** //nn ：時刻歴全数 //fs() ：フーリエ・スペクトル //df ：スペクトルの振動数間隔（Hz） //band ：バンド幅（Hz） //nfold：スペクトル値総数 double pi = Math.PI; double tt, udf, dif, s; int i, j; int nfold, lmax, ll, ln, lt, le; double[] w = new double[nn]; double[] g = new double[nn]; double[] g1 = new double[nn]; double[] g2 = new double[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] * Math.Pow(Math.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] = Math.Sqrt(g[0] * tt); for (i = 2; i <= nfold - 1; i++) { fs[i - 1] = Math.Sqrt(g[i - 1] * tt / 2.0); } fs[nfold - 1] = Math.Sqrt(g[nfold - 1] * tt); } } } //-------------------------------------------------------------------- } }