Option Explicit

Sub wavelet()
'ガボールウェーブレット変換
'結果をGraphRで読めるようにcsvファイル出力
'coded by 伊津野和行　ver.1 2015/9/3  ver.5 2018/3/20
'参考：http://www.softist.com/　のC++コード
'たたみ込み積分をフーリエスペクトルどうしのかけ算で計算

Dim sampleFreq As Double      'サンプリング周波数 Hz
Dim sigma As Double   '減衰係数σ
Dim n As Long      '計算するデータ点数　2の累乗でなければならない．
Dim nn As Long     '実際のデータ点数
Dim m As Long      '周波数の個数
Dim fmin As Double, fmax As Double  '最低および最高周波数(Hz)
Dim df                '周波数の計算刻み幅
Dim f As Double       '計算対象の周波数
Dim i As Long, j As Long, k As Long  'カウンター
Dim rc As Integer     'ファイルを開いたときのリターンコード
Dim data_r() As Double, data_i() As Double  '入力波形の動的配列　実部と虚部
Dim gabor_r() As Double, gabor_i() As Double 'ガボール・フィルタの実部と虚部
Dim wavelet() As Double 'ウェーブレットを記憶する動的配列
Dim timeAxis() As Double     '時間 sec
Dim fileNumber As Long, csvFile As String  'csvファイルを出力するための変数
Dim amp As Double       '出力がHzなら1，kHzなら1000, MHzなら1000000
Dim sec As Double       '出力が秒なら1，ミリ秒なら1000
Dim temp_r As Double, temp_i As Double '一時記憶変数
Dim ind As Integer 'FFTのスイッチ　-1：フーリエ変換　+1：逆フーリエ変換
Dim unit(2) As String    '出力データの単位を記録するための文字変数
Dim nyquist As Double    'ナイキスト周波数
Dim temp As Double       '2の階乗を計算するための一時記憶変数

'データの入力
    Sheets("data").Select
    nn = [b3]
    If (nn And (nn - 1)) <> 0 Then
        temp = Log(nn) / Log(2)
        n = 2 ^ Application.RoundUp(temp, 0)
    Else
        n = nn
    End If
    
    sigma = [b1]
    sampleFreq = [b2]
    nyquist = sampleFreq / 2#
    m = [b4]
    fmin = [b5]
    fmax = [b6]
    If fmax > nyquist Then
        MsgBox "B6セルはサンプリング周波数の半分以下でなければなりません。", vbCritical, "Error"
        End
    End If
    ReDim data_r(n): ReDim data_i(n)
    ReDim gabor_r(n): ReDim gabor_i(n)
    ReDim wavelet(n)
    ReDim timeAxis(n)
    For i = 1 To n
     timeAxis(i) = CDbl(i) / sampleFreq
     data_r(i) = 0#
     data_i(i) = 0#
    Next i
    j = n / 2 + 1   'ガウス関数の時間軸に合わせるためデータの後半を先に読み込み，あとから前半部を読み込んでいる．
    For i = 1 To nn
        If (j > n) Then
            k = j - n
        Else
            k = j
        End If
        data_r(k) = Cells(i + 1, "D")
        j = j + 1
    Next i
    If [b7] = "Hz" Then
        amp = 1#
        unit(2) = " y-axis: Hz"
    ElseIf [b7] = "kHz" Then
        amp = 1000#
        unit(2) = " y-axis: kHz"
    Else
        amp = 1000000#
        unit(2) = " y-axis: MHz"
    End If
    If [b8] = "秒" Then
       sec = 1#
       unit(1) = " x-axis: s"
    ElseIf [b8] = "ミリ秒" Then
       sec = 1000#
       unit(1) = " x-axis: ms"
    Else
       sec = 1000000#
       unit(1) = " x-axis: micro-s"
    End If

'Graph-R用にCSVファイルを出力　まずは4行目までを出力
    csvFile = ActiveWorkbook.Path & "\wavelet.csv"
    fileNumber = FreeFile
    Open csvFile For Output As #fileNumber
    Write #fileNumber, "DataFormat", Val(1)
    Write #fileNumber, "Gabor Wavelet Transform   Units" + unit(1) + unit(2)
    Write #fileNumber, Now
    Write #fileNumber, " ";
    For j = 1 To n - 1
        Write #fileNumber, Val(timeAxis(j) * sec);
    Next j
    Write #fileNumber, Val(timeAxis(n) * sec)

'入力波のフーリエ変換
    ind = -1
    Call fast(data_r, data_i, n, ind)

'ガボールウェーブレット変換
    df = (fmax - fmin) / (m - 1)
    f = fmin
    For i = 1 To m
        [b13] = i: Calculate
        If f < 0.0001 Then
            Write #fileNumber, Val(f / amp);
            For j = 1 To n - 1
                Write #fileNumber, Val(0);
            Next j
            Write #fileNumber, Val(0)
        Else
            'ガボール・フィルタのフーリエスペクトル
            Call GaborTransformFFT(n, f, sampleFreq, sigma, gabor_r, gabor_i)
            '入力波形のスペクトル×フィルターのスペクトル
            For j = 1 To n
                temp_r = data_r(j) * gabor_r(j) - data_i(j) * gabor_i(j)
                gabor_i(j) = data_i(j) * gabor_r(j) + data_r(j) * gabor_i(j)
                gabor_r(j) = temp_r
            Next j
            'フーリエ逆変換してウェーブレットに
            ind = 1
            Call fast(gabor_r, gabor_i, n, ind)
            'nで割る必要あり．絶対値を求めて記憶．
            For j = 1 To n
                temp_r = gabor_r(j) / n
                temp_i = gabor_i(j) / n
                wavelet(j) = Sqr((temp_r * temp_r + temp_i * temp_i) * f / sampleFreq)
            Next j
            'CSVファイルに保存
            Write #fileNumber, Val(f / amp);
            For j = 1 To n - 1
                Write #fileNumber, Val(wavelet(j));
            Next j
            Write #fileNumber, Val(wavelet(n))
        End If
        f = f + df
    Next i
    Close fileNumber
    [b14] = "計算終了"
    
End Sub

Sub fast(xr() As Double, xi() As Double, nn As Long, ind As Integer)

'Fast Fourier transform (FFT)
'大崎順彦：新・地震動のスペクトル解析入門，鹿島出版会
'のFortranコードをVBAに移植．
'xr(): 変数の実数部, xi(): 変数の虚数部
'nn: xr()とxi()の宣言された大きさ．2の累乗でなければならない．
'ind: フーリエ変換の時は-1，フーリエ逆変換の時は1とする．

Dim tr As Double, ti As Double
Dim theta As Double
Dim i As Long, j As Long, k As Long
Dim m As Long, kmax As Long, istep As Long
Dim PI As Double
PI = Atn(1#) * 4#

j = 1
For i = 1 To nn
  If i < j Then
    tr = xr(j): ti = xi(j)
    xr(j) = xr(i): xi(j) = xi(i)
    xr(i) = tr: xi(i) = ti
  End If
  m = nn / 2
  Do While (j > m)
    j = j - m
    m = m / 2
    If (m < 2) Then Exit Do
  Loop
  j = j + m
Next i
kmax = 1
Do While kmax < nn
  istep = kmax * 2
  For k = 1 To kmax
    theta = PI * CDbl(ind * (k - 1)) / CDbl(kmax)
    For i = k To nn Step istep
      j = i + kmax
      tr = xr(j) * Cos(theta) - xi(j) * Sin(theta)
      ti = xr(j) * Sin(theta) + xi(j) * Cos(theta)
      xr(j) = xr(i) - tr: xi(j) = xi(i) - ti
      xr(i) = xr(i) + tr: xi(i) = xi(i) + ti
    Next i
  Next k
  kmax = istep
Loop
End Sub

Sub GaborTransformFFT(n As Long, f As Double, sampleFreq As Double, _
                      sigma As Double, gabor_r() As Double, gabor_i() As Double)

'ガボール関数を作り，そのフーリエスペクトルを求める
Dim pi2 As Double                    '2π
Dim t As Double, omegaT As Double
Dim gauss As Double
Dim i As Long, j As Long, ind As Integer

pi2 = Atn(1#) * 8#

For i = 1 To n
  j = i - n / 2
  t = j / sampleFreq * f
  'ガウス関数
  gauss = Exp(-t * t / 2# / sigma / sigma) / (Sqr(pi2) * sigma)
  'ωt
  omegaT = pi2 * t
  gabor_r(i) = gauss * Cos(omegaT)
  gabor_i(i) = gauss * Sin(omegaT)
Next i

'フーリエ変換
ind = -1
Call fast(gabor_r, gabor_i, n, ind)

End Sub

Sub init()
'D列のデータを消去する
 Dim lastrow As Long
 lastrow = Cells(Rows.Count, 4).End(xlUp).Row
 Range(Cells(2, 4), Cells(lastrow, 4)).Clear
End Sub