Option Explicit On
Option Strict On
Public Class Form1

    Private Sub Form1_Load(ByVal sender As Object, ByVal e As System.EventArgs) Handles Me.Load
        ToolStripLabel1.Text = "Band of period (sec)"
        ToolStripTextBox1.Text = "-1"
        ToolStripTextBox2.Text = "-1"
    End Sub

    Private Sub ToolStripButton1_Click(ByVal sender As System.Object, ByVal e As System.EventArgs) Handles ToolStripButton1.Click

        Dim i As Integer
        Dim dt As Double
        Dim ndata As Integer
        Dim nn As Integer
        Dim nd As Integer
        Dim wave() As Double
        Dim tdraw() As Double
        Dim xr() As Double
        Dim xi() As Double
        Dim datax() As Double
        Dim datay() As Double
        Dim fk As Double
        Dim fq As Double
        Dim tt As Double
        Dim rand As New Random
        Dim pi As Double = Math.PI
        Dim icase As Integer
        Dim tcutdw As Double
        Dim tcutup As Double

        dt = 0.01 '時間刻み
        ndata = 50000 'データ数
        ReDim wave(ndata - 1)
        ReDim tdraw(ndata - 1)
        nn = 2
        Do
            nn = nn * 2
        Loop While nn < ndata
        ReDim xr(nn - 1)
        ReDim xi(nn - 1)
        ReDim datax(CInt(nn / 2))
        ReDim datay(CInt(nn / 2))
        For i = 1 To nn 'データ数を2の階乗にセット
            xr(i - 1) = 0.0 : xi(i - 1) = 0.0
        Next
        '試験波形作成
        For i = 1 To ndata
            tt = CDbl(i - 1) * dt
            wave(i - 1) = 0.0
            'wave(i - 1) = 1.0 * (rand.NextDouble() - 0.5) * 2.0 'rand.NextDouble()：0〜1の乱数
            wave(i - 1) = wave(i - 1) + 1.0 * Math.Cos(2.0 * pi / 0.1 * tt)
            wave(i - 1) = wave(i - 1) + 0.5 * Math.Cos(2.0 * pi / 1.0 * tt)
            wave(i - 1) = wave(i - 1) + 1.0 * Math.Cos(2.0 * pi / 10.0 * tt)
            wave(i - 1) = wave(i - 1) + 2.0 * Math.Cos(2.0 * pi / 100.0 * tt)
            xr(i - 1) = wave(i - 1)
            tdraw(i - 1) = tt
        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 CInt(nn / 2) + 1
            fq = CDbl(i - 1) / CDbl(nn) / dt '振動数
            If i = 1 Then fq = CDbl(1) / CDbl(2 * nn) / dt '振動数
            fk = dt * CDbl(nn) * Math.Sqrt(xr(i - 1) * xr(i - 1) + xi(i - 1) * xi(i - 1)) 'フーリエスペクトル
            datax(i - 1) = 1.0 / fq '周期
            datay(i - 1) = fk
        Next i
        '原波形出力
        nd = ndata - 1
        icase = 0
        SCRDRAW(icase, nd, tdraw, wave)
        '原波形フーリエスペクトル出力
        nd = CInt(nn / 2)
        icase = 1
        SCRDRAW(icase, nd, datax, datay)

        tcutup = CDbl(nn) * dt : tcutdw = 2.0 * dt '理論最大周期・最小周期
        ToolStripLabel1.Text = tcutdw.ToString("0.00") & "〜" & tcutup.ToString("0") & "(sec)"
        '*************************************************************
        '周波数処理
        tcutdw = CDbl(ToolStripTextBox1.Text) '周期下限
        tcutup = CDbl(ToolStripTextBox2.Text) '周期上限
        If tcutdw < 0.0 Then tcutdw = 2.0 * dt
        If tcutup < 0.0 Then tcutup = CDbl(nn) * dt
        For i = 2 To CInt(nn / 2)  '特定周波数領域のカット
            fq = CDbl(i - 1) / CDbl(nn) / dt '振動数
            If (fq < (1.0 / tcutup)) Or ((1.0 / tcutdw) < fq) Then
                xr(i - 1) = 0.0 : xr(nn - i + 1) = 0.0
                xi(i - 1) = 0.0 : xi(nn - i + 1) = 0.0
            End If
        Next i
        i = CInt(nn / 2) + 1 : fq = CDbl(i - 1) / CDbl(nn) / dt
        If (fq < (1.0 / tcutup)) Or ((1.0 / tcutdw) < fq) Then
            xr(i - 1) = 0.0 : xi(i - 1) = 0.0
        End If
        '*************************************************************
        For i = 1 To nn 'フーリエ逆変換用に虚数部の符号を反転
            xr(i - 1) = xr(i - 1)
            xi(i - 1) = -xi(i - 1)
        Next i
        Call FFT(nn, xr, xi) 'フーリエ逆変換
        '処理波形記憶
        For i = 1 To ndata
            wave(i - 1) = xr(i - 1)
        Next i
        For i = 1 To nn '配列初期化
            xr(i - 1) = 0.0 : xi(i - 1) = 0.0
        Next
        For i = 1 To ndata
            xr(i - 1) = wave(i - 1)
        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 CInt(nn / 2) + 1
            fq = CDbl(i - 1) / CDbl(nn) / dt '振動数
            If i = 1 Then fq = CDbl(1) / CDbl(2 * nn) / dt '振動数
            fk = dt * CDbl(nn) * Math.Sqrt(xr(i - 1) * xr(i - 1) + xi(i - 1) * xi(i - 1)) 'フーリエスペクトル
            datax(i - 1) = 1.0 / fq '周期
            datay(i - 1) = fk
        Next i
        '処理波形出力
        nd = ndata - 1
        icase = 2
        SCRDRAW(icase, nd, tdraw, wave)
        '処理波形フーリエスペクトル出力
        nd = CInt(nn / 2)
        icase = 3
        SCRDRAW(icase, nd, datax, datay)

    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) + 1.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 SCRDRAW(ByVal icase As Integer, ByVal nd As Integer, ByRef datax() As Double, ByRef datay() As Double)
        Dim kpt As Integer
        Dim HSIZE As Integer
        Dim DHL As Integer
        Dim DHR As Integer
        Dim VSIZE As Integer
        Dim DVL As Integer
        Dim DVU As Integer
        Dim sxjiku As String = ""
        Dim xmin As Double
        Dim xmax As Double
        Dim dx As Double
        Dim slogx As String = ""
        Dim syjiku As String = ""
        Dim ymin As Double
        Dim ymax As Double
        Dim dy As Double
        Dim slogy As String = ""
        Dim kxi0 As Integer
        Dim kxf0 As Integer
        Dim kyi0 As Integer
        Dim kyf0 As Integer
        Dim i As Integer
        Dim pi As Double = Math.PI
        Dim bmp As Bitmap
        Dim g As Graphics
        Dim fname0 As String
        Dim fname1 As String
        Dim fname2 As String
        Dim fname3 As String

        Select Case icase
            Case 0, 2
                HSIZE = 800 : DHL = 80 : DHR = 20
                VSIZE = 200 : DVL = 50 : DVU = 20
                sxjiku = "時刻 (sec)"
                syjiku = "応力 (MPa)"
                xmin = 0.0 : xmax = 500.0 : dx = 50.0 : slogx = "N"
                ymin = -5.0 : ymax = 5.0 : dy = 1.0 : slogy = "N"
            Case 1, 3
                HSIZE = 500 : DHL = 80 : DHR = 20
                VSIZE = 400 : DVL = 50 : DVU = 20
                sxjiku = "周期  T(sec)"
                syjiku = "フーリエ・スペクトル  (MPa*sec)"
                xmin = -2.0 : xmax = 4.0 : dx = 1.0 : slogx = "L"
                ymin = 0.0 : ymax = 500.0 : dy = 50.0 : slogy = "N"
        End Select
        kxi0 = DHL
        kxf0 = HSIZE - DHR
        kyi0 = DVU
        kyf0 = VSIZE - DVL
        If slogx = "L" Then
            For i = 0 To nd
                If datax(i) < System.Math.Pow(10.0, xmin) Then datax(i) = System.Math.Pow(10.0, xmin)
                datax(i) = System.Math.Log(datax(i)) / System.Math.Log(10.0)
            Next
        End If
        If slogy = "L" Then
            For i = 0 To nd
                If datay(i) < System.Math.Pow(10.0, ymin) Then datay(i) = System.Math.Pow(10.0, ymin)
                datay(i) = System.Math.Log(datay(i)) / System.Math.Log(10.0)
            Next
        End If
        '画面描画
        If icase = 3 Then
            kpt = 1
            PictureBox2.Size = New Size(kpt * HSIZE, kpt * VSIZE)
            PictureBox2.Visible = True
            PictureBox2.Left = 0
            PictureBox2.Top = ToolStrip1.Height
            Me.ClientSize = New Size(PictureBox2.Width, PictureBox2.Height + ToolStrip1.Height)
            bmp = New Bitmap(PictureBox2.Width, PictureBox2.Height)
            PictureBox2.Image = bmp
            g = Graphics.FromImage(PictureBox2.Image)
            g.FillRectangle(Brushes.White, 0, 0, PictureBox2.Width, PictureBox2.Height)
            Call PLOT(g, kpt, sxjiku, syjiku, xmin, xmax, dx, ymin, ymax, dy, kxi0, kxf0, kyi0, kyf0, nd, datax, datay, slogx, slogy)
            g.Dispose()
        End If
        '画像書き出し
        kpt = 2
        PictureBox1.Size = New Size(kpt * HSIZE, kpt * VSIZE)
        PictureBox1.Visible = False
        bmp = New Bitmap(PictureBox1.Width, PictureBox1.Height)
        PictureBox1.Image = bmp
        g = Graphics.FromImage(PictureBox1.Image)
        g.FillRectangle(Brushes.White, 0, 0, PictureBox1.Width, PictureBox1.Height)
        Call PLOT(g, kpt, sxjiku, syjiku, xmin, xmax, dx, ymin, ymax, dy, kxi0, kxf0, kyi0, kyf0, nd, datax, datay, slogx, slogy)
        g.Dispose()

        fname0 = "c:\vb2010pro\vb4FFTex0\data-file\out0.png"
        fname1 = "c:\vb2010pro\vb4FFTex0\data-file\out1.png"
        fname2 = "c:\vb2010pro\vb4FFTex0\data-file\out2.png"
        fname3 = "c:\vb2010pro\vb4FFTex0\data-file\out3.png"
        Select Case icase
            Case 0 : PictureBox1.Image.Save(fname0, System.Drawing.Imaging.ImageFormat.Png)
            Case 1 : PictureBox1.Image.Save(fname1, System.Drawing.Imaging.ImageFormat.Png)
            Case 2 : PictureBox1.Image.Save(fname2, System.Drawing.Imaging.ImageFormat.Png)
            Case 3 : PictureBox1.Image.Save(fname3, System.Drawing.Imaging.ImageFormat.Png)
        End Select

    End Sub

    Private Sub PLOT(ByVal g As Graphics, ByVal kpt As Integer, ByVal sxjiku As String, ByVal syjiku As String, _
                     ByVal xmin As Double, ByVal xmax As Double, ByVal dx As Double, _
                     ByVal ymin As Double, ByVal ymax As Double, ByVal dy As Double, _
                     ByVal kxi0 As Integer, ByVal kxf0 As Integer, ByVal kyi0 As Integer, ByVal kyf0 As Integer, _
                     ByVal nd As Integer, ByRef datax() As Double, ByRef datay() As Double, _
                     ByVal slogx As String, ByVal slogy As String)

        Dim i As Integer
        Dim kxi As Integer : Dim kxf As Integer
        Dim kyi As Integer : Dim kyf As Integer
        Dim xx As Double : Dim yy As Double
        Dim kxx1 As Integer : Dim kyy1 As Integer
        Dim kxx2 As Integer : Dim kyy2 As Integer
        Dim str As String
        Dim f As New Font("ＭＳ ゴシック", kpt * 10)
        Dim TextSize1 As New System.Drawing.SizeF
        Dim TextSize2 As New System.Drawing.SizeF

        kxi = kpt * kxi0 : kxf = kpt * kxf0 : kyi = kpt * kyi0 : kyf = kpt * kyf0

        '座標軸
        Dim LPen As New System.Drawing.Pen(System.Drawing.Color.Black)
        LPen.DashStyle = Drawing2D.DashStyle.Dot
        If slogx = "N" Then Call DRN_XJIKU(g, f, LPen, kpt, xmin, xmax, dx, kxi, kxf, kyi, kyf)
        If slogy = "N" Then Call DRN_YJIKU(g, f, LPen, kpt, ymin, ymax, dy, kxi, kxf, kyi, kyf)
        If slogx = "L" Then Call DRL_XJIKU(g, f, LPen, kpt, xmin, xmax, kxi, kxf, kyi, kyf)
        If slogy = "L" Then Call DRL_YJIKU(g, f, LPen, kpt, ymin, ymax, kxi, kxf, kyi, kyf)
        '枠線描画
        g.DrawRectangle(New Pen(Color.Black, 2), kxi, kyi, kxf - kxi, kyf - kyi)
        'ｙ軸名描画
        str = syjiku
        TextSize2 = g.MeasureString(str, f)
        TextSize1 = g.MeasureString("-1.0", f)
        kxx1 = CInt(kxi - kpt * 10 - TextSize1.Width - TextSize2.Height)
        kyy1 = CInt((kyi + kyf) / 2 + TextSize2.Width / 2)
        Call INC_STR(g, f, str, kxx1, kyy1, -90)
        'ｘ軸名描画
        str = sxjiku
        TextSize2 = g.MeasureString(str, f)
        kxx1 = CInt((kxi + kxf) / 2 - TextSize2.Width / 2)
        kyy1 = CInt(kyf + kpt * 10 + TextSize1.Height)
        Call INC_STR(g, f, str, kxx1, kyy1, 0)
        'データプロット
        'dc = kpt * 5
        'For i = 0 To nd
        'xx = datax(i)
        'yy = datay(i)
        'kxx = kxi + CInt((xx - xmin) * (kxf - kxi) / (xmax - xmin)) - CInt(dc / 2)
        'kyy = kyf - CInt((yy - ymin) * (kyf - kyi) / (ymax - ymin)) - CInt(dc / 2)
        'g.DrawEllipse(New Pen(Color.Red, 1), kxx, kyy, dc, dc)
        'Next
        'データを線で連結
        xx = datax(0)
        yy = datay(0)
        kxx1 = kxi + CInt((xx - xmin) * (kxf - kxi) / (xmax - xmin))
        kyy1 = kyf - CInt((yy - ymin) * (kyf - kyi) / (ymax - ymin))
        For i = 1 To nd
            xx = datax(i)
            yy = datay(i)
            kxx2 = kxi + CInt((xx - xmin) * (kxf - kxi) / (xmax - xmin))
            kyy2 = kyf - CInt((yy - ymin) * (kyf - kyi) / (ymax - ymin))
            g.DrawLine(New Pen(Color.Blue, 2), kxx1, kyy1, kxx2, kyy2)
            kxx1 = kxx2
            kyy1 = kyy2
        Next i
        f.Dispose()

    End Sub


    Private Sub INC_STR(ByVal g As Graphics, ByVal f As System.Drawing.Font, ByVal str As String, _
                        ByVal kxx As Integer, ByVal kyy As Integer, ByVal ang As Single)
        '軸名描画
        g.ScaleTransform(1.0, 1.0)                '横・縦の表示比率を設定
        g.TranslateTransform(kxx, kyy)            '表示位置の設定（表示位置を原点とする座標移動）
        g.RotateTransform(ang)                    '表示角度を指定
        g.DrawString(str, f, Brushes.Black, 0, 0) '描画実行
        g.ResetTransform()                        '単位行列にリセット
    End Sub


    Private Sub DRN_XJIKU(ByVal g As Graphics, ByVal f As System.Drawing.Font, ByVal LPen As System.Drawing.Pen, _
                          ByVal kpt As Integer, _
                          ByVal xmin As Double, ByVal xmax As Double, ByVal dx As Double, _
                          ByVal kxi As Integer, ByVal kxf As Integer, _
                          ByVal kyi As Integer, ByVal kyf As Integer)
        Dim i As Integer
        Dim ix As Integer
        Dim xx As Double
        Dim wv As Double
        Dim str As String
        Dim kxx As Integer
        Dim TextSize1 As New System.Drawing.SizeF

        '普通ｘ軸描画
        ix = CInt((xmax - xmin) / dx)
        For i = 0 To ix
            xx = xmin + CDbl(i) * dx
            kxx = kxi + CInt((xx - xmin) * (kxf - kxi) / (xmax - xmin))
            wv = xmin + CDbl(i) * dx
            str = wv.ToString("0")
            If CInt(dx * 1000.0) Mod 10 <> 0 Then str = wv.ToString("0.000")
            If CInt(dx * 1000.0) Mod 10 = 0 And CInt(dx * 100.0) Mod 10 <> 0 Then str = wv.ToString("0.00")
            If (CInt(dx * 1000.0) Mod 10 = 0 And CInt(dx * 100.0) Mod 10 = 0) And CInt(dx * 10.0) Mod 10 <> 0 Then str = wv.ToString("0.0")
            TextSize1 = g.MeasureString(str, f)
            g.DrawLine(LPen, kxx, kyi, kxx, kyf)
            g.DrawString(str, f, Brushes.Black, kxx - TextSize1.Width / 2, kyf + kpt * 3)
        Next i
    End Sub


    Private Sub DRN_YJIKU(ByVal g As Graphics, ByVal f As System.Drawing.Font, ByVal LPen As System.Drawing.Pen, _
                          ByVal kpt As Integer, _
                          ByVal ymin As Double, ByVal ymax As Double, ByVal dy As Double, _
                          ByVal kxi As Integer, ByVal kxf As Integer, _
                          ByVal kyi As Integer, ByVal kyf As Integer)
        Dim i As Integer
        Dim iy As Integer
        Dim yy As Double
        Dim wv As Double
        Dim str As String
        Dim kyy As Integer
        Dim TextSize1 As New System.Drawing.SizeF

        '普通ｙ軸描画
        iy = CInt((ymax - ymin) / dy)
        For i = 0 To iy
            yy = ymin + CDbl(i) * dy
            kyy = kyf - CInt((yy - ymin) * (kyf - kyi) / (ymax - ymin))
            wv = ymin + CDbl(i) * dy
            str = wv.ToString("0")
            If CInt(dy * 1000.0) Mod 10 <> 0 Then str = wv.ToString("0.000")
            If CInt(dy * 1000.0) Mod 10 = 0 And CInt(dy * 100.0) Mod 10 <> 0 Then str = wv.ToString("0.00")
            If (CInt(dy * 1000.0) Mod 10 = 0 And CInt(dy * 100.0) Mod 10 = 0) And CInt(dy * 10.0) Mod 10 <> 0 Then str = wv.ToString("0.0")
            TextSize1 = g.MeasureString(str, f)
            g.DrawLine(LPen, kxi, kyy, kxf, kyy)
            g.DrawString(str, f, Brushes.Black, kxi - TextSize1.Width - kpt * 2, kyy - TextSize1.Height / 2)
        Next i
    End Sub


    Private Sub DRL_XJIKU(ByVal g As Graphics, ByVal f As System.Drawing.Font, ByVal LPen As System.Drawing.Pen, _
                          ByVal kpt As Integer, _
                          ByVal xmin As Double, ByVal xmax As Double, _
                          ByVal kxi As Integer, ByVal kxf As Integer, _
                          ByVal kyi As Integer, ByVal kyf As Integer)
        Dim i As Integer
        Dim j As Integer
        Dim xx As Double
        Dim str As String
        Dim kxx As Integer
        Dim TextSize1 As New System.Drawing.SizeF

        '対数ｘ軸描画
        For i = CInt(xmin) To CInt(xmax) - 1
            For j = 1 To 9
                xx = System.Math.Log(CDbl(j) * System.Math.Pow(10.0, CDbl(i))) / System.Math.Log(10.0)
                kxx = kxi + CInt((xx - xmin) * (kxf - kxi) / (xmax - xmin))
                g.DrawLine(LPen, kxx, kyi, kxx, kyf)
            Next j
        Next i
        For i = CInt(xmin) To CInt(xmax)
            xx = CDbl(i)
            kxx = kxi + CInt((xx - xmin) * (kxf - kxi) / (xmax - xmin))
            str = CStr(System.Math.Pow(10.0, CDbl(i)))
            TextSize1 = g.MeasureString(str, f)
            g.DrawString(str, f, Brushes.Black, kxx - TextSize1.Width / 2, kyf + kpt * 3)
        Next i
    End Sub


    Private Sub DRL_YJIKU(ByVal g As Graphics, ByVal f As System.Drawing.Font, ByVal LPen As System.Drawing.Pen, _
                          ByVal kpt As Integer, _
                          ByVal ymin As Double, ByVal ymax As Double, _
                          ByVal kxi As Integer, ByVal kxf As Integer, _
                          ByVal kyi As Integer, ByVal kyf As Integer)
        Dim i As Integer
        Dim j As Integer
        Dim yy As Double
        Dim str As String
        Dim kyy As Integer
        Dim TextSize1 As New System.Drawing.SizeF

        '対数ｙ軸描画
        For i = CInt(ymin) To CInt(ymax) - 1
            For j = 1 To 9
                yy = System.Math.Log(CDbl(j) * System.Math.Pow(10.0, CDbl(i))) / System.Math.Log(10.0)
                kyy = kyf - CInt((yy - ymin) * (kyf - kyi) / (ymax - ymin))
                g.DrawLine(LPen, kxi, kyy, kxf, kyy)
            Next j
        Next i
        For i = CInt(ymin) To CInt(ymax)
            yy = CDbl(i)
            kyy = kyf - CInt((yy - ymin) * (kyf - kyi) / (ymax - ymin))
            str = CStr(System.Math.Pow(10.0, CDbl(i)))
            TextSize1 = g.MeasureString(str, f)
            g.DrawString(str, f, Brushes.Black, kxi - TextSize1.Width - kpt * 2, kyy - TextSize1.Height / 2)
        Next i
    End Sub


End Class