ニューメリカル･レシピのFFTプログラムをExcelのVBAに移植

ExcelのVBAで、サブルーチンに配列をそのまま渡せることを（いまさら）初めて知った。じゃあ、普通に配列渡して高速フーリエ変換（FFT)とか計算できるな、ってことで手元にあったニューメリカル・レシピ（in C)のCプログラムを移植してみた。

プログラムはこれです。↓

「FFT.bas」をダウンロード

使い方は、上のを標準モジュールとして読み込んだ後、配列data()を使って

call FFT(data nn, isign)

と呼ぶんだけど、ちょっと使い方が変わっていて

data(1) = データ1のreal part,

data(2) = データ1のimag part

data(3) = データ2のreal part

data(4) = データ2のimag part

・・・

とか並べます。データnnまで（つまり配列としては2*nnは必要）。当然nnは2^Mの形です。

isignは+1ならフーリエ変換、-1なら逆フーリエ変換ですよ。

まあ、

Option Explicit

Private Sub CommandButton1_Click()
    Dim data(16) As Double
    Dim i As Integer
       
    data(1) = 1#
    data(2) = 0#
    data(3) = 1#
    data(4) = 0#
    data(5) = 1#
    data(6) = 0#
    For i = 7 To 16
        data(i) = 0
    Next i

    Call FFT(data, 8, 1)
   
    For i = 1 To 16
        Worksheets("Sheet1").Cells(i, 5) = data(i)
    Next i
End Sub

Sub FFT(data() As Double, nn As Integer, isign As Integer)
    Dim n As Integer, mmax As Integer, m As Integer, j As Integer, istep As Integer, i As Integer
    Dim wtemp As Double, wr As Double, wpr As Double, wpi As Double, wi As Double, theta As Double
    Dim tempr As Double, tempi As Double
    Dim dummy As Double
   
    n = 2 * nn
    j = 1
    For i = 1 To n - 1 Step 2
        If j > i Then
            dummy = data(i)
            data(i) = data(j)
            data(j) = dummy
            
            dummy = data(i + 1)
            data(i + 1) = data(j + 1)
            data(j + 1) = dummy
        End If
        m = n / 2
        While (m >= 2 And j > m)
            j = j - m
            m = m / 2
        Wend
        j = j + m
    Next i
   
    mmax = 2
    While (n > mmax)
        istep = mmax * 2#
        theta = isign * (2# * 4# * Atn(1) / mmax)
        wtemp = Sin(0.5 * theta)
        wpr = -2# * wtemp * wtemp
        wpi = Sin(theta)
        wr = 1#
        wi = 0#
        For m = 1 To mmax - 1 Step 2
            For i = m To n Step istep
                j = i + mmax
                tempr = wr * data(j) - wi * data(j + 1)
                tempi = wr * data(j + 1) + wi * data(j)
                data(j) = data(i) - tempr
                data(j + 1) = data(i + 1) - tempi
                data(i) = data(i) + tempr
                data(i + 1) = data(i + 1) + tempi
            Next i
            wtemp = wr
            wr = wtemp * wpr - wi * wpi + wr
            wi = wi * wpr + wtemp * wpi + wi
               
        Next m
            
        mmax = istep
    Wend

End Sub

こんな感じで使えばいいかな。間違っていても責任はもてまへん。。。