概要

 高速フーリエ変換(FFT)をVBAで実装し、動作することを確認した。2015/6/7 プログラムの一部間違いを修正。


背景と目的

 以前から、高速フーリエ変換(FFT)を行うVBAのプログラムが欲しかったのだがWeb上を探してもコメントが書かれていなかったり、データサイズが固定だったりとなかなか希望に合うものが見つからなかったので、自分でプログラムを実装してみる。


詳細
1.方法

 アルゴリズムは、最も有名なCooley-Tukey型アルゴリズムとする。実装にあたっては、手元にあった『C言語ではじめる音のプログラミング』という本(図1)の2.3節を参考にする。要するに、C言語で書かれたサンプルをVBAで書き直すようなものである。

IMG_0800
図1 『C言語ではじめる音のプログラミング』という本

2.実装

 『C言語ではじめる音のプログラミング』を参考に、以下のようなプログラムを実装した。基本的な構成としては、
・バタフライ演算をデータサイズに合わせて繰り返し実行するループ
・アルゴリズム上順序がばらばらになってしまう計算結果データを並び替える処理
を行っている。
 入力は、xという変数として出力は参照渡しの引数y_reとy_imにそれぞれ実部、虚部を戻している。また、引数Nはデータサイズで、2のべき乗でなければ処理をしないようにした。また、窓関数をかけるオプションも考えたが、FFT処理関数の中に仕込む必要はないと考え、実装しないことにした。

Public Function FFT(x() As Double, ByRef y_re() As Double, ByRef y_im() As Double, N As Long) As Boolean

    '引数
    'x          入力データ
    'y_re       戻り値実部
    'y_im       戻り値虚部
    'N          FFTサイズ
    
    '戻り値
    'TRUE   FFT成功
    'FALSE  データ数が2のべき乗でないので失敗
    
    '変数
    Dim PI As Double 'π
    Dim StageCount As Integer 'ステージ数
    Dim BlockCount As Long 'ブロック数
    Dim NodeCount As Long 'ノード数
    Dim stage As Integer 'ステージ番号
    Dim block As Long 'ブロック番号
    Dim node As Long 'ノード番号
    Dim n1 As Long '計算用
    Dim n2 As Long '計算用
    Dim i As Long '計算用
    Dim r As Double '計算用
    Dim index() As Long 'インデックス並べ替え用
    Dim re1 As Double '計算用
    Dim im1 As Double '計算用
    Dim re2 As Double '計算用
    Dim im2 As Double '計算用
    
    'FFTサイズ確認
    If CLng(N And (N - 1)) <> 0 Then
        FFT = False
        Exit Function '2のべき乗でなければ終わり
    End If
    
    '計算準備
    PI = Math.Atn(1) * 4 'π
    StageCount = CInt(Round(Log(N) / Log(2), 0)) 'ステージ数
    y_re = x '実部
    ReDim y_im(N) '虚部
    
    'バタフライ演算の繰り返しループ
    For stage = 0 To StageCount - 1

        BlockCount = 2 ^ stage 'ブロック数
        NodeCount = N / BlockCount 'ノード数
        r = 2 * PI / N * BlockCount '定数
        
        For block = 0 To BlockCount - 1
            For node = 0 To NodeCount / 2 - 1
                'バタフライ演算
                n1 = node + NodeCount * block
                n2 = n1 + NodeCount / 2
                re1 = y_re(n1)
                im1 = y_im(n1)
                re2 = y_re(n2)
                im2 = y_im(n2)
                y_re(n1) = re1 + re2
                y_im(n1) = im1 + im2
                
                '2015/06/07 符号を間違えていたので修正
                y_re(n2) = (re1 - re2) * Cos(r * node) + (im1 - im2) * Sin(r * node)
                y_im(n2) = -(re1 - re2) * Sin(r * node) + (im1 - im2) * Cos(r * node)
                'y_re(n2) = (re1 - re2) * Cos(r * node) - (im1 - im2) * Sin(r * node)
                'y_im(n2) = (re1 - re2) * Sin(r * node) + (im1 - im2) * Cos(r * node)
            Next node
        Next block
    Next stage
    
    '並び替え処理
    ReDim index(N)
    For stage = 0 To StageCount - 1
        For i = 0 To 2 ^ stage - 1
            index(2 ^ stage + i) = index(i) + 2 ^ (StageCount - stage - 1)
        Next i
    Next stage
    For i = 0 To N - 1
        If index(i) > i Then
            re1 = y_re(index(i))
            im1 = y_im(index(i))
            y_re(index(i)) = y_re(i)
            y_im(index(i)) = y_im(i)
            y_re(i) = re1
            y_im(i) = im1
        End If
    Next i
    
    FFT = True

End Function
3.動作確認と実行速度

 上記のプログラムをN=8のデータを使って動作させたところ正しく動作していることが確認できた。また、実行速度は私の環境ではN=2^16で体感的に1秒以内であったので十分実用になりそうであった。


まとめと今後の課題

 高速フーリエ変換(FFT)をVBAで実装し、動作することを確認できた。今後の信号処理に活用したい。