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
        Call LVSET()
    End Sub

    Private Sub LVSET()
        ListView1.Clear()
        ListView1.View = View.Details
        ListView1.Columns.Add("m", 60, HorizontalAlignment.Right)
        ListView1.Columns.Add("data", 60, HorizontalAlignment.Right)
        ListView1.Columns.Add("k", 60, HorizontalAlignment.Right)
        ListView1.Columns.Add("Re(Ck)", 60, HorizontalAlignment.Right)
        ListView1.Columns.Add("Im(Ck)", 60, HorizontalAlignment.Right)
        ListView1.Columns.Add("abs(Ck)", 60, HorizontalAlignment.Right)
        ListView1.Columns.Add("Re(xm)", 60, HorizontalAlignment.Right)
        ListView1.Columns.Add("Im(xm)", 60, HorizontalAlignment.Right)
    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 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 fname1 As String = ""
        Dim fname2 As String = ""
        Dim data() As Double
        Dim xr() As Double
        Dim xi() As Double
        Dim cr() As Double
        Dim ci() As Double

        Dim ctl As Control
        For Each ctl In Me.Controls
            Console.WriteLine(ctl.Name)
        Next

        Call LVSET()

        If OpenFileDialog1.ShowDialog() = Windows.Forms.DialogResult.OK Then fname1 = OpenFileDialog1.FileName
        sr = New System.IO.StreamReader(fname1, System.Text.Encoding.Default)
        dat = sr.ReadLine
        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 * 1 '必要に応じて後続の0を追加
        ReDim data(nn - 1)
        ReDim xr(nn - 1)
        ReDim xi(nn - 1)
        ReDim cr(nn - 1)
        ReDim ci(nn - 1)

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

        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)
            cr(i - 1) = xr(i - 1) '変換値実数部の記憶
            ci(i - 1) = xi(i - 1) '変換値虚数部の記憶
        Next i

        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 nn
            Debug.WriteLine(CStr(i - 1) & " " & CStr(cr(i - 1)) & " " & CStr(ci(i - 1)))
        Next

        For i = 1 To nn
            Dim lvitem As ListViewItem = New ListViewItem((i - 1).ToString)
            lvitem.SubItems.Add(data(i - 1).ToString("0.000"))
            lvitem.SubItems.Add((i - 1).ToString)
            lvitem.SubItems.Add(cr(i - 1).ToString("0.000"))
            lvitem.SubItems.Add(ci(i - 1).ToString("0.000"))
            lvitem.SubItems.Add((Math.Sqrt(cr(i - 1) ^ 2 + ci(i - 1) ^ 2)).ToString("0.000"))
            lvitem.SubItems.Add(xr(i - 1).ToString("0.000"))
            lvitem.SubItems.Add(xi(i - 1).ToString("0.000"))
            ListView1.Items.Add(lvitem)
        Next i

        If SaveFileDialog1.ShowDialog() = Windows.Forms.DialogResult.OK Then fname2 = SaveFileDialog1.FileName
        sw = New System.IO.StreamWriter(fname2, 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 To nn
            dat = (i - 1).ToString
            dat = dat & "," & data(i - 1).ToString("0.000")
            dat = dat & "," & (i - 1).ToString
            dat = dat & "," & cr(i - 1).ToString("0.000")
            dat = dat & "," & ci(i - 1).ToString("0.000")
            dat = dat & "," & (Math.Sqrt(cr(i - 1) ^ 2 + ci(i - 1) ^ 2)).ToString("0.000")
            dat = dat & "," & xr(i - 1).ToString("0.000")
            dat = dat & "," & xi(i - 1).ToString("0.000")
            sw.WriteLine(dat)
        Next i
        sw.Close()

    End Sub


    Private Sub FFT(ByVal nn As Integer, ByRef xr() As Double, ByRef xi() As Double)
        '**************************************************
        '高速フーリエ変換・逆変換
        '**************************************************
        'nn  ：データ数（2の累乗）
        'xr()：入出力データ実数部
        'xi()：入出力データ虚数部
        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

End Class