Option Explicit On
Option Strict On
Public Class Form1
    '***********************************************************
    '複素数計算のための構造体宣言とFunctionプロシージャ開始
    '***********************************************************
    '定義済み関数
    'za=C_DOUBLE(a,b)   ：複素数実数部と虚数部への数値の代入（Re(z)=a，Im(z)=b）
    'za=C_PLUS(z1,z2)   ：加算（za=z1+z2）
    'za=C_MINUS(z1,z2)  ：減算（za=z1-z2）
    'za=C_TIMES(z1,z2)  ：乗算（za=z1*z2）
    'za=C_DIVIDED(z1,z2)：除算（za=z1/z2）
    'da=D_ABS(z)        ：絶対値（da=sqrt(Re(z)^2+Im(z)^2)）
    'za=C_CONJ(z)       ：共役複素数（z=a+ib→za=a-ib）
    'za=C_SQRT(z)       ：平方根（za=sqrt(z)）
    'za=C_EXP(z)        ：exp関数（za=exp(z)）

    Private Structure Complex
        '構造体宣言
        Dim Real As Double
        Dim Imag As Double
    End Structure

    Private Function C_DOUBLE(ByVal a As Double, ByVal b As Double) As Complex
        '複素数実数部と虚数部への数値の代入
        C_DOUBLE.Real = a
        C_DOUBLE.Imag = b
    End Function

    Private Function C_PLUS(ByVal z1 As Complex, ByVal z2 As Complex) As Complex
        '加算
        C_PLUS.Real = z1.Real + z2.Real
        C_PLUS.Imag = z1.Imag + z2.Imag
    End Function

    Private Function C_MINUS(ByVal z1 As Complex, ByVal z2 As Complex) As Complex
        '減算
        C_MINUS.Real = z1.Real - z2.Real
        C_MINUS.Imag = z1.Imag - z2.Imag
    End Function

    Private Function C_TIMES(ByVal z1 As Complex, ByVal z2 As Complex) As Complex
        '乗算
        C_TIMES.Real = z1.Real * z2.Real - z1.Imag * z2.Imag
        C_TIMES.Imag = z1.Real * z2.Imag + z2.Real * z1.Imag
    End Function

    Private Function C_DIVIDED(ByVal z1 As Complex, ByVal z2 As Complex) As Complex
        '除算
        If z2.Real = 0.0 And z2.Imag = 0.0 Then
            MessageBox.Show("分母が0です！")
            Me.Close()
        Else
            C_DIVIDED.Real = (z1.Real * z2.Real + z1.Imag * z2.Imag) / (z2.Real ^ 2 + z2.Imag ^ 2)
            C_DIVIDED.Imag = (z2.Real * z1.Imag - z1.Real * z2.Imag) / (z2.Real ^ 2 + z2.Imag ^ 2)
        End If
    End Function

    Private Function D_ABS(ByVal z As Complex) As Double
        '絶対値
        D_ABS = Math.Sqrt(z.Real ^ 2 + z.Imag ^ 2)
    End Function

    Private Function C_CONJ(ByVal z As Complex) As Complex
        '共役複素数
        C_CONJ.Real = z.Real
        C_CONJ.Imag = -z.Imag
    End Function

    Private Function C_SQRT(ByVal z As Complex) As Complex
        '平方根
        If z.Real = 0.0 And z.Imag = 0.0 Then
            C_SQRT.Real = 0.0
            C_SQRT.Imag = 0.0
        Else
            Dim rr As Double
            Dim phi As Double
            rr = Math.Sqrt(z.Real ^ 2 + z.Imag ^ 2)
            phi = Math.Acos(z.Real / rr)
            If z.Imag < 0.0 Then phi = 2.0 * Math.PI - phi 'phiの範囲を0〜2πとする（acosの戻り値は0〜π）
            C_SQRT.Real = rr ^ (0.5) * Math.Cos(0.5 * phi)
            C_SQRT.Imag = rr ^ (0.5) * Math.Sin(0.5 * phi)
        End If
    End Function

    Private Function C_EXP(ByVal z As Complex) As Complex
        'EXP関数
        C_EXP.Real = Math.Exp(z.Real) * Math.Cos(z.Imag)
        C_EXP.Imag = Math.Exp(z.Real) * Math.Sin(z.Imag)
    End Function
    '***********************************************************
    '複素数計算のための構造体宣言とFunctionプロシージャ終了
    '***********************************************************


    Private Sub CFFT(ByVal n As Integer, ByRef x() As Complex, ByVal ind As Integer)
        '***********************************************************
        '複素数計算による高速フーリエ変換
        '***********************************************************
        'n  ：データ数（2の累乗） 
        'x()：入出力データ（複素数） 
        'ind：-1=フーリエ変換，+1=フーリエ逆変換 
        Dim temp As Complex
        Dim theta As Complex
        Dim j As Integer
        Dim k As Integer
        Dim m As Integer
        Dim kmax As Integer
        Dim istep As Integer
        Dim pi As Double = Math.PI

        j = 1
        For i = 1 To n
            If i < j Then
                temp = x(j - 1)
                x(j - 1) = x(i - 1)
                x(i - 1) = temp
            End If
            m = CInt(n / 2)
            Do
                If j <= m Then Exit Do
                j = j - m
                m = CInt(m / 2)
            Loop While m >= 2
            j = j + m
        Next i
        kmax = 1
        Do
            istep = kmax * 2
            For k = 1 To kmax
                theta = C_DOUBLE(0.0, pi * CDbl(ind * (k - 1)) / CDbl(kmax))
                For i = k To n Step istep
                    j = i + kmax
                    temp = C_TIMES(x(j - 1), C_EXP(theta))
                    x(j - 1) = C_MINUS(x(i - 1), temp)
                    x(i - 1) = C_PLUS(x(i - 1), temp)
                Next i
            Next k
            kmax = istep
        Loop While kmax < n

    End Sub


    Private Sub Button1_Click(ByVal sender As System.Object, ByVal e As System.EventArgs) Handles Button1.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 x() As Complex
        Dim y() As Complex

        Call LVSET() 'ListView設定

        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 x(nn - 1)
        ReDim y(nn - 1)

        For i = 1 To nn 'データ数を2の累乗にセット
            data(i - 1) = 0.0
            x(i - 1) = C_DOUBLE(0.0, 0.0)
        Next i
        For i = 1 To ndata
            dat = sr.ReadLine : sbuf = dat.Split(delim) : data(i - 1) = CDbl(sbuf(0))
            data(i - 1) = data(i - 1) + 0.0
            x(i - 1) = C_DOUBLE(data(i - 1), 0.0)
        Next i
        sr.Close()
        Call CFFT(nn, x, -1) 'フーリエ変換
        For i = 1 To nn
            x(i - 1) = C_DIVIDED(x(i - 1), C_DOUBLE(CDbl(nn), 0)) '出力値をデータ数で除す
            y(i - 1) = x(i - 1)
        Next i
        Call CFFT(nn, y, +1) 'フーリエ逆変換

        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(x(i - 1).Real.ToString("0.000"))
            lvitem.SubItems.Add(x(i - 1).Imag.ToString("0.000"))
            lvitem.SubItems.Add(D_ABS(x(i - 1)).ToString("0.000"))
            lvitem.SubItems.Add(y(i - 1).Real.ToString("0.000"))
            lvitem.SubItems.Add(y(i - 1).Imag.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 & "," & x(i - 1).Real.ToString("0.000")
            dat = dat & "," & x(i - 1).Imag.ToString("0.000")
            dat = dat & "," & D_ABS(x(i - 1)).ToString("0.000")
            dat = dat & "," & y(i - 1).Real.ToString("0.000")
            dat = dat & "," & y(i - 1).Imag.ToString("0.000")
            sw.WriteLine(dat)
        Next i
        sw.Close()

    End Sub

    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

End Class