Public Class Form1 Private Sub Form1_Load(ByVal sender As System.Object, ByVal e As System.EventArgs) Handles MyBase.Load ToolStripLabel1.Text = "バンド幅" ToolStripTextBox1.Text = "0.4" End Sub Private Sub ToolStripButton1_Click(ByVal sender As System.Object, ByVal e As System.EventArgs) Handles ToolStripButton1.Click Call main() MessageBox.Show("処理完了", "完了通知") End Sub Private Sub main() 'フーリエスペクトル比算定 '分子・分母にするデータのサンプリングピッチ・データ数は等しくとる Dim i As Integer Dim fnameR As String = "" Dim fnameW As String = "" Dim sr As System.IO.StreamReader Dim sw As System.IO.StreamWriter Dim dat As String Dim sbuf() As String Dim delim() As Char = {","c} Dim ndata As Integer Dim ddy() As Double Dim xr() As Double Dim xi() As Double Dim df As Double Dim scom As String Dim dt As Double Dim nn As Integer Dim nd As Integer Dim fsp0() As Double '原スペクトル Dim fsp1() As Double '平滑化スペクトル Dim fq() As Double '振動数 Dim band As Double Dim ndmax As Integer = CInt(2 ^ 16) ReDim fsp0(ndmax) ReDim fsp1(ndmax) ReDim fq(ndmax) For i = 0 To ndmax fsp0(i) = 0.0 fsp1(i) = 0.0 fq(i) = 0.0 Next i '平滑化バンド幅設定 band = CDbl(ToolStripTextBox1.Text) '入出力ファイル読み込み OpenFileDialog1.InitialDirectory = System.IO.Directory.GetCurrentDirectory() If OpenFileDialog1.ShowDialog() = DialogResult.OK Then fnameR = OpenFileDialog1.FileName If SaveFileDialog1.ShowDialog() = DialogResult.OK Then 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 = CDbl(sbuf(1)) dat = sr.ReadLine() : sbuf = dat.Split(delim) : ndata = CInt(sbuf(1)) nn = 2 Do nn = nn * 2 Loop While nn < ndata ReDim ddy(nn - 1) ReDim xr(nn - 1) ReDim xi(nn - 1) For i = 1 To nn 'データ数を2の階乗にセット ddy(i - 1) = 0.0 Next i For i = 1 To ndata dat = sr.ReadLine() : sbuf = dat.Split(delim) : ddy(i - 1) = CDbl(sbuf(0)) Next i sr.Close() '平滑化フーリエスペクトル作成 For i = 1 To nn xr(i - 1) = ddy(i - 1) xi(i - 1) = 0.0 Next i Call FFT(nn, xr, xi) 'フーリエ変換 For i = 1 To nn '戻り値をデータ数で除す xr(i - 1) = xr(i - 1) / CDbl(nn) xi(i - 1) = xi(i - 1) / CDbl(nn) Next i nd = CInt(nn / 2) For i = 1 To nd + 1 fsp0(i - 1) = dt * CDbl(nn) * Math.Sqrt(xr(i - 1) * xr(i - 1) + xi(i - 1) * xi(i - 1)) fsp1(i - 1) = fsp0(i - 1) Next i df = 1.0 / dt / CDbl(nn) Call SWIN(nn, fsp1, df, band) 'スペクトル平滑化 For i = 0 To nd fq(i) = CDbl(i) / CDbl(nn) / dt '振動数 Next i '結果書き出し 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 To nd dat = fq(i).ToString("E") dat = dat & "," & fsp0(i).ToString("E") dat = dat & "," & fsp1(i).ToString("E") sw.WriteLine(dat) Next i sw.Close() End Sub Private Sub FFT(ByVal nn As Integer, ByRef xr() As Double, ByRef xi() As Double) '実数計算高速フーリエ変換・逆変換 Dim g As Integer : Dim h As Integer : Dim i As Integer : Dim j As Integer : Dim k As Integer Dim l As Integer : Dim m As Integer : Dim n As Integer : Dim p As Integer : Dim q As Integer Dim a As Double : Dim b As Double : Dim xd As Double Dim s() As Double Dim c() As Double n = nn ReDim s(CInt(n / 2)) ReDim c(CInt(n / 2)) i = 0 : j = 0 : k = 0 : l = 0 : p = 0 : h = 0 : g = 0 : q = 0 m = CInt(System.Math.Log(CDbl(n)) / System.Math.Log(2.0)) a = 0.0 : b = Math.PI * 2.0 / CDbl(n) For i = 0 To CInt(n / 2) s(i) = System.Math.Sin(a) : c(i) = System.Math.Cos(a) : a = a + b Next i l = n : h = 1 For g = 1 To m l = CInt(l / 2) : k = 0 For q = 1 To h p = 0 For i = k To l + k - 1 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 Then xr(j) = a : xi(j) = b Else xr(j) = a * c(p) + b * s(p) : xi(j) = b * c(p) - a * s(p) End If p = p + h Next i k = k + l + l Next q h = h + h Next g j = CInt(n / 2) For i = 1 To n - 1 k = n If j < i Then xd = xr(i) : xr(i) = xr(j) : xr(j) = xd xd = xi(i) : xi(i) = xi(j) : xi(j) = xd End If k = CInt(k / 2) Do While j >= k j = j - k : k = CInt(k / 2) If k = 0 Then Exit Do Loop j = j + k Next i End Sub Private Sub SWIN(ByVal nn As Integer, ByRef fs() As Double, ByVal df As Double, ByVal band As Double) 'スペクトル平滑化 'nn ：時刻歴全数 'fs() ：フーリエ・スペクトル 'df ：スペクトルの振動数間隔（Hz） 'band ：バンド幅（Hz） 'nfold：スペクトル値総数 Dim w() As Double Dim g() As Double Dim g1() As Double Dim g2() As Double Dim tt As Double Dim udf As Double Dim dif As Double Dim s As Double Dim i As Integer Dim j As Integer Dim nfold As Integer Dim lmax As Integer Dim ll As Integer Dim ln As Integer Dim lt As Integer Dim le As Integer Dim pi As Double ReDim g(nn - 1) ReDim g1(nn - 1) ReDim g2(nn - 1) pi = Math.PI nfold = CInt(nn / 2) + 1 If band > 0.0 Then tt = 1.0 / df udf = 280.0 / 151.0 / band * df lmax = CInt(2.0 / udf) + 1 ReDim w(lmax - 1) If udf <= 0.5 Then w(0) = 0.75 * udf For i = 2 To lmax dif = 0.5 * pi * CDbl(i - 1) * udf w(i - 1) = w(0) * System.Math.Pow(System.Math.Sin(dif) / dif, 4.0) Next i g(0) = fs(0) * fs(0) / tt For i = 2 To nfold - 1 g(i - 1) = 2.0 * fs(i - 1) * fs(i - 1) / tt Next i 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 To lt g1(i - 1) = 0.0 Next i For i = 1 To nfold g1(ll - 2 + i) = g(i - 1) Next i For i = lmax To le s = w(0) * g1(i - 1) For j = 2 To lmax s = s + w(j - 1) * (g1(i - j) + g1(i + j - 2)) Next j g2(i - 1) = s Next i For j = 2 To lmax g2(ll + j - 2) = g2(ll + j - 2) + g2(ll - j) g2(ln - j) = g2(ln - j) + g2(ln + j - 2) Next j For i = 1 To nfold g(i - 1) = g2(ll - 2 + i) Next i fs(0) = System.Math.Sqrt(g(0) * tt) For i = 2 To nfold - 1 fs(i - 1) = System.Math.Sqrt(g(i - 1) * tt / 2.0) Next i fs(nfold - 1) = System.Math.Sqrt(g(nfold - 1) * tt) End If End If End Sub End Class