Option Explicit On
Option Strict On
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 nnn As Integer
        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 strcom1() As String
        Dim strcom2() As String
        Dim rn1() As String
        Dim rn2() As String
        Dim wn0() As String
        Dim ncase As Integer
        Dim icase As Integer

        Dim dt As Double
        Dim nn As Integer
        Dim nd As Integer
        Dim fs01() As Double '分子原スペクトル
        Dim fs02() As Double '分子原スペクトル
        Dim fsp1() As Double '分子平滑化スペクトル
        Dim fsp2() As Double '分母平滑化スペクトル
        Dim fspr() As Double '平滑化スペクトル比
        Dim fq() As Double '振動数
        Dim band As Double

        Dim ndmax As Integer = CInt(2 ^ 16)
        ReDim fs01(ndmax)
        ReDim fs02(ndmax)
        ReDim fsp1(ndmax)
        ReDim fsp2(ndmax)
        ReDim fq(ndmax)

        For i = 0 To ndmax
            fs01(i) = 0.0
            fs02(i) = 0.0
            fsp1(i) = 0.0
            fsp2(i) = 0.0
            fq(i) = 0.0
        Next i

        '平滑化バンド幅設定
        band = CDbl(ToolStripTextBox1.Text)

        '入出力ファイル読み込み
        OpenFileDialog1.InitialDirectory = System.IO.Directory.GetCurrentDirectory()
        If OpenFileDialog1.ShowDialog() = Windows.Forms.DialogResult.OK Then fnameR = OpenFileDialog1.FileName
        sr = New System.IO.StreamReader(fnameR, System.Text.Encoding.Default)
        dat = sr.ReadLine() : sbuf = dat.Split(delim) : ncase = CInt(sbuf(0))
        ReDim strcom1(ncase - 1)
        ReDim strcom2(ncase - 1)
        ReDim rn1(ncase - 1)
        ReDim rn2(ncase - 1)
        ReDim wn0(ncase - 1)
        For icase = 1 To ncase
            dat = sr.ReadLine() : sbuf = dat.Split(delim)
            strcom1(icase - 1) = sbuf(0) '分子スペクトル名
            strcom2(icase - 1) = sbuf(1) '分母スペクトル名
            rn1(icase - 1) = sbuf(2) '分子加速度時刻歴ファイル名
            rn2(icase - 1) = sbuf(3) '分母加速度時刻歴ファイル名
            wn0(icase - 1) = sbuf(4) '結果出力ファイル名
        Next icase
        sr.Close()

        For icase = 1 To ncase
            For nnn = 1 To 2
                '原波形読み込みと平滑化フーリエスペクトル計算
                Select Case nnn
                    Case 1 '分子側加速度時刻歴入力・スペクトル計算
                        fnameR = System.IO.Path.GetDirectoryName(fnameR) & "\" & rn1(icase - 1)
                        Call CALFSP(fnameR, dt, nn, band, fs01, fsp1)
                    Case 2 '分母側加速度時刻歴入力・スペクトル計算
                        fnameR = System.IO.Path.GetDirectoryName(fnameR) & "\" & rn2(icase - 1)
                        Call CALFSP(fnameR, dt, nn, band, fs02, fsp2)
                End Select
            Next nnn
            'スペクトル比計算・振動数設定
            nd = CInt(nn / 2)
            ReDim fspr(nd)
            For i = 0 To nd
                fspr(i) = fsp1(i) / fsp2(i)
                fq(i) = CDbl(i) / CDbl(nn) / dt '振動数
            Next i

            '結果書き出し
            fnameW = System.IO.Path.GetDirectoryName(fnameR) & "\" & wn0(icase - 1)
            sw = New System.IO.StreamWriter(fnameW, False, System.Text.Encoding.Default)
            dat = "分子ファイル名" & "," & strcom1(icase - 1) & "," & rn1(icase - 1) : sw.WriteLine(dat)
            dat = "分母ファイル名" & "," & strcom2(icase - 1) & "," & rn2(icase - 1) : sw.WriteLine(dat)
            dat = "出力ファイル名" & "," & "sp比" & "," & wn0(icase - 1) : sw.WriteLine(dat)
            dat = "バンド幅" & "," & band.ToString : sw.WriteLine(dat)
            dat = "データ数" & "," & (nd + 1).ToString : sw.WriteLine(dat)
            dat = "振動数,分子原SP,分母原SP,平滑化分子SP,平滑化分母SP,平滑化SP比" : sw.WriteLine(dat)
            For i = 0 To nd
                dat = fq(i).ToString("E")
                dat = dat & "," & fs01(i).ToString("E")
                dat = dat & "," & fs02(i).ToString("E")
                dat = dat & "," & fsp1(i).ToString("E")
                dat = dat & "," & fsp2(i).ToString("E")
                dat = dat & "," & fspr(i).ToString("E")
                sw.WriteLine(dat)
            Next i
            sw.Close()
        Next icase

        '出力ファイル一覧出力
        fnameW = System.IO.Path.GetDirectoryName(fnameR) & "\fileoutlist.txt"
        sw = New System.IO.StreamWriter(fnameW, False, System.Text.Encoding.Default)
        dat = ncase.ToString("0") : sw.WriteLine(dat)
        For icase = 1 To ncase
            dat = wn0(icase - 1) : sw.WriteLine(dat)
        Next icase
        sw.Close()

    End Sub


    Private Sub CALFSP(ByVal fnameR As String, ByRef dt As Double, ByRef nn As Integer, _
                       ByVal band As Double, ByRef fs0() As Double, ByRef fsp() As Double)

        '加速度時刻歴波形読み込みとスペクトル計算制御
        Dim sr As System.IO.StreamReader
        Dim dat As String
        Dim sbuf() As String
        Dim delim() As Char = {","c}
        Dim ndata As Integer
        Dim nd As Integer
        Dim ddy() As Double
        Dim dy() As Double
        Dim y() As Double
        Dim ddymax As Double
        Dim dymax As Double
        Dim ymax As Double
        Dim xr() As Double
        Dim xi() As Double
        Dim df As Double

        '原波形読み込み
        sr = New System.IO.StreamReader(fnameR, System.Text.Encoding.Default)
        dat = sr.ReadLine() : sbuf = dat.Split(delim)
        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
        nd = CInt(nn / 2)
        ReDim ddy(nn - 1)
        ReDim dy(nn - 1)
        ReDim y(nn - 1)
        ReDim xr(nn - 1)
        ReDim xi(nn - 1)

        For i = 1 To nn 'データ数を2の階乗にセット
            ddy(i - 1) = 0.0 : dy(i - 1) = 0.0 : y(i - 1) = 0.0
        Next
        For i = 1 To ndata
            dat = sr.ReadLine() : sbuf = dat.Split(delim) : ddy(i - 1) = CDbl(sbuf(0))
        Next i
        sr.Close()

        '平滑化フーリエスペクトル作成
        Call IACC(ndata, dt, ddy, dy, y, ddymax, dymax, ymax) '数値積分
        Call CRAC(ndata, dt, ddy, dy, y, ddymax, dymax, ymax) '基線補正
        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
        For i = 1 To nd + 1
            fs0(i - 1) = dt * CDbl(nn) * Math.Sqrt(xr(i - 1) * xr(i - 1) + xi(i - 1) * xi(i - 1)) '分子フーリエスペクトル
            fsp(i - 1) = fs0(i - 1)
        Next i
        df = 1.0 / dt / CDbl(nn)
        Call SWIN(nn, fsp, df, band) 'スペクトル平滑化

    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 IACC(ByVal nn As Integer, ByVal dt As Double, ByRef ddy() As Double, ByRef dy() As Double, ByRef y() As Double, _
                 ByRef ddymax As Double, ByRef dymax As Double, ByRef ymax As Double)

        '加速度時刻歴数値積分
        'nn:加速度時刻歴データ総数
        'dt    :時間間隔（入力値）
        'ddy() :加速度時刻歴（入力値）
        'dy()  :速度時刻歴（計算出力）
        'y()   :変位時刻歴（計算出力値）
        'ddymax:加速度最大値（計算出力値）
        'dymax :速度最大値（計算出力値）
        'ymax  :変位最大値（計算出力値）

        Dim i As Integer

        ddymax = 0.0
        dymax = 0.0
        ymax = 0.0

        dy(0) = 0.0
        y(0) = 0.0
        For i = 1 To nn - 1
            dy(i) = dy(i - 1) + (ddy(i - 1) + ddy(i)) * dt / 2.0
            y(i) = y(i - 1) + dy(i - 1) * dt + (ddy(i - 1) / 3.0 + ddy(i) / 6.0) * dt * dt
            If (ddymax < System.Math.Abs(ddy(i))) Then ddymax = System.Math.Abs(ddy(i))
            If (dymax < System.Math.Abs(dy(i))) Then dymax = System.Math.Abs(dy(i))
            If (ymax < System.Math.Abs(y(i))) Then ymax = System.Math.Abs(y(i))
        Next i

    End Sub


    Private Sub CRAC(ByVal nn As Integer, ByVal dt As Double, ByRef ddy() As Double, ByRef dy() As Double, ByRef y() As Double, _
                     ByVal ddymax As Double, ByVal dymax As Double, ByVal ymax As Double)

        '加速度時刻歴基線補正
        'nn:加速度時刻歴データ総数
        'dt    :時間間隔（入力値）
        'ddy() :加速度時刻歴（入力値・書換出力値）
        'dy()  :速度時刻歴（入力値）
        'y()   :変位時刻歴（入力値）
        'ddymax:加速度最大値（入力値）
        'dymax :速度最大値（入力値）
        'ymax  :変位最大値（入力値）

        Dim i As Integer
        Dim tt As Double
        Dim t As Double
        Dim sum As Double
        Dim a1 As Double
        Dim a0 As Double
        Dim acmax As Double
        Dim coef As Double

        tt = CDbl(nn - 1) * dt
        t = 0.0
        For i = 0 To nn - 1
            y(i) = y(i) * (3.0 * tt - 2.0 * t) * t * t
            t = t + dt
        Next i
        sum = (y(0) + y(nn - 1)) / 2.0
        For i = 1 To nn - 2
            sum = sum + y(i)
        Next i
        sum = sum * dt
        a1 = 28.0 / 13.0 / tt / tt * (2.0 * dy(nn - 1) - 15.0 / System.Math.Pow(tt, 5.0) * sum)
        a0 = dy(nn - 1) / tt - a1 / 2.0 * tt
        t = 0.0
        acmax = 0.0
        For i = 0 To nn - 1
            ddy(i) = ddy(i) - a0 - a1 * t
            If (acmax < System.Math.Abs(ddy(i))) Then acmax = System.Math.Abs(ddy(i))
            t = t + dt
        Next i
        coef = ddymax / acmax
        For i = 0 To nn - 1
            ddy(i) = ddy(i) * coef
        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