FFT演算VBAマクロ：最小版

FFT演算VBAマクロ

FFT演算を実行するExcelVBAのマクロについて紹介します

今回紹介するマクロは，FFT演算の中身を理解しやすい事を意識し，最小限の実装を心がけています

注意事項

実用的なFFTソフトウェアで必要となりそうな処理を実装していません

• サンプリング周波数は定義していません
• 振幅の単位の物理量換算はしていません
• ウインドウ処理はしていません
• ２のべき乗のキリの良いデータ点数以上のデータは捨てています

これらの対応については今後別記事で紹介できればと思っています

コードの紹介

実装コンセプト

• 複素数の型を定義し複素数演算を関数化することで，バタフライ演算の部分のコードを簡潔にし，理解しやすくしています
• 演算に関するパラメータを１つの構造体としてまとめています
• FFT演算用のコードはクラスライブラリとして実装しています

複素数

複素数は実部と虚部を持った構造体として定義しています

'複素数型
Private Type Complex
    re As Double    '実数部
    Im As Double    '虚数部
End Type

複素数の演算

複素数の和，差，積を利用しますので，それぞれを関数化して用意しておきます

'複素数の掛け算
Private Function fComplexMulti(ByRef A As Complex, ByRef B As Complex) As Complex
    Dim rtn As Complex
    rtn.re = (A.re * B.re) - (A.Im * B.Im)
    rtn.Im = (A.re * B.Im) + (A.Im * B.re)
    fComplexMulti = rtn
End Function
'複素数の足し算
Private Function fComplexAdd(ByRef A As Complex, ByRef B As Complex) As Complex
    Dim rtn As Complex
    rtn.re = A.re + B.re
    rtn.Im = A.Im + B.Im
    fComplexAdd = rtn
End Function
'複素数の引き算
Private Function fComplexSub(ByRef A As Complex, ByRef B As Complex) As Complex
    Dim rtn As Complex
    rtn.re = A.re - B.re
    rtn.Im = A.Im - B.Im
    fComplexSub = rtn
End Function

パラメータ

FFT演算パラメータをまとめた構造体部分です，InputDataSizeをユーザー決定してやると，他のパラメータは導出できます．

'FFT解析のパラメータをまとめた構造体
Private Type Parameter
    InputDataSize  As Long      '受け取ったデータのサイズ（２のべき乗でなくても良い）
    AnalysisPoints As Long      '解析対象データ数（２のべき乗に制限したデータ点数）
    Radix As Long               '基数（２の何乗なのかを示す数）
    TwiddleFactor() As Complex  '回転因子の配列（基数に応じて生成する）
End Type

有効データ点数

注意事項に記述したとおり２のべき乗のキリの良い部分以降のデータの下図にのように捨てています

基数の算出

基数，即ち「２の何乗にしたらキリの良いデータ点数なのか？」を求めるには，入力されたデータ点数を２で割って何回割れたかをカウントすることで求めています．

'-------------------------------------------------------------------------------------
'基数を求める：データ数を繰り返し２で割って何回割れたかを求める，即ち２の何乗になるか
'-------------------------------------------------------------------------------------
Private Function RadixCalc(ByVal dataLen As Long)
    Dim rdx As Integer
    Dim temp As Double
    temp = dataLen
    rdx = 0
    Do Until temp <= 1
        temp = temp / 2
        rdx = rdx + 1
    Loop
    RadixCalc = rdx
End Function

回転因子の算出

基数の値を決定できると回転因子を算出することが可能です
半円πを有効データ点数で分割してSINとCOSを求め，複素数に当てはめたものが回転因子です

'-------------------------------------------------------------------------------------
'回転因子配列の生成
' 引数には基数を入れる
'　戻り値は回転因子の入った複素数型の配列
'-------------------------------------------------------------------------------------
Private Function fTwiddleFactor(ByVal iRadix As Integer) As Complex()
    Dim pointNum As Long    'データポイント数
    Dim twiddle() As Complex '回転因子配列
    Dim stepVal As Double
    Dim radian As Double
    
    '回転因子のポイント数は２の基数乗
    pointNum = (2 ^ (iRadix - 1))
    
    '半周期＝πをポイント数で割ると位相の分解能が求まる
    stepVal = WorksheetFunction.Pi() / (pointNum)
        
    '回転因子配列の生成
    ReDim twiddle(pointNum - 1)
    For i = 0 To pointNum - 1
        radian = stepVal * i          '角度をずらしてゆく
        twiddle(i).re = Cos(radian)   '実数
        twiddle(i).Im = Sin(radian)   '虚数
    Next
    
    '戻り値
    fTwiddleFactor = twiddle
End Function

前処理（ビットリバース）

FFT演算のコアとなるバタフライ演算を実行する前の前処理としてビットリバースを行います．
データのインデックスを２進数に換算し，MSB・LSBの順番をリバースしたインデックスを付与し，並べ替えを行うとビットリバースした配列を生成することが可能です．

バタフライ演算

先にビットリバース処理を行っているので，バタフライ演算の前側で回転因子をかけています．

コード全体

'***************************************************************************
'FFT演算クラス
'
'***************************************************************************
'////////////////////////////////////////////////
'ユーザー定義型
'////////////////////////////////////////////////
'複素数型
Private Type Complex
    re As Double    '実数部
    Im As Double    '虚数部
End Type

'FFT解析のパラメータをまとめた構造体
Private Type Parameter
    InputDataSize  As Long      '受け取ったデータのサイズ（２のべき乗でなくても良い）
    AnalysisPoints As Long      '解析対象データ数（２のべき乗に制限したデータ点数）
    Radix As Long               '基数（２の何乗なのかを示す数）
    TwiddleFactor() As Complex  '回転因子の配列（基数に応じて生成する）
End Type

'////////////////////////////////////////////////
'クラス変数
'////////////////////////////////////////////////
Private SignalData() As Double      '入力信号データ
Private Param As Parameter          'FFT解析に必要なパラメータ

'============================================================================================
'入力信号データを内部変数にセットし，データ数から決定可能なFFT解析パラメータを決定する
'============================================================================================
Public Sub SetSignalData(ByVal inputData)
    '入力データをクラス変数に保持
    SignalData = inputData
    
    '入力データ長
    Param.InputDataSize = UBound(SignalData) + 1
    '基数を求める
    Param.Radix = RadixCalc(Param.InputDataSize)
    '解析に有効な２のべき乗のデータポイント数を求める
    Param.AnalysisPoints = 2 ^ Param.Radix
    '基数から回転因子を作成
    Param.TwiddleFactor = fTwiddleFactor(Param.Radix)
End Sub

'============================================================================================
'FFT演算実行
'============================================================================================
Public Function Run()
    Dim rvsCmp() As Complex      '入力データをビットリバースした複素数の配列
    Dim bufData() As Complex
    Dim outData() As Complex
    Dim absData() As Double
    Dim stage As Integer
    '----------------------------------------------
    'ビットリバース処理と複素数化
    '----------------------------------------------
    ReDim rvsCmp(Param.AnalysisPoints - 1)
    For i = 0 To Param.AnalysisPoints - 1
        With rvsCmp(BitReverse(i, Param.Radix))
            .re = SignalData(i)
            .Im = 0             '入力信号に虚数は無いのでゼロ固定
        End With
    Next
    
    '----------------------------------------------
    'バタフライ演算
    '----------------------------------------------
    ReDim outData(Param.AnalysisPoints - 1)
    bufData = rvsCmp
    '基数分のステージ繰り返し
    For stage = 1 To Param.Radix
        '回転因子の乗算（バタフライ演算の前処理）
        outData = TwiddleFactorMulti(bufData, Param.Radix, stage, Param.TwiddleFactor)
        'バタフライ演算
        bufData = ButterflyOperation(outData, Param.Radix, stage)
    Next
        
    '----------------------------------------------
    '虚数を絶対値に変換（二乗平均）
    '----------------------------------------------
    ReDim absData(Param.AnalysisPoints - 1)
    For i = 0 To Param.AnalysisPoints - 1
        absData(i) = Sqr((bufData(i).re) ^ 2 + (bufData(i).Im) ^ 2)
    Next
    
    '並び替えた結果を出力する
    Run = absData
End Function

'-------------------------------------------------------------------------------------
'回転因子配列の生成
' 引数には基数を入れる
'　戻り値は回転因子の入った複素数型の配列
'-------------------------------------------------------------------------------------
Private Function fTwiddleFactor(ByVal iRadix As Integer) As Complex()
    Dim pointNum As Long    'データポイント数
    Dim twiddle() As Complex '回転因子配列
    Dim stepVal As Double
    Dim radian As Double
    
    '回転因子のポイント数は２の基数乗
    pointNum = (2 ^ (iRadix - 1))
    
    '半周期＝πをポイント数で割ると位相の分解能が求まる
    stepVal = WorksheetFunction.Pi() / (pointNum)
        
    '回転因子配列の生成
    ReDim twiddle(pointNum - 1)
    For i = 0 To pointNum - 1
        radian = stepVal * i          '角度をずらしてゆく
        twiddle(i).re = Cos(radian)   '実数
        twiddle(i).Im = Sin(radian)   '虚数
    Next
    
    '戻り値
    fTwiddleFactor = twiddle
End Function

'-------------------------------------------------------------------------------------
'基数を求める：データ数を繰り返し２で割って何回割れたかを求める，即ち２の何乗になるか
'-------------------------------------------------------------------------------------
Private Function RadixCalc(ByVal dataLen As Long)
    Dim rdx As Integer
    Dim temp As Double
    temp = dataLen
    rdx = 0
    Do Until temp <= 1
        temp = temp / 2
        rdx = rdx + 1
    Loop
    RadixCalc = rdx
End Function

'-----------------------------------------------------------------------------------
'ビット反転:入力された数を２進数で表現したときにMSBからLSBの並び順を逆にする
' 例：入力 0010b →出力0100b
'-----------------------------------------------------------------------------------
Private Function BitReverse(ByVal InputNum As Integer, ByVal BitLen As Long) As Integer
    Dim outputNum As Integer
    'ﾋﾞｯﾄ反転動作
    For i = 0 To BitLen - 1
        If (InputNum And (2 ^ i)) <> 0 Then
            outputNum = outputNum + (2 ^ (BitLen - 1 - i))
        End If
    Next

    BitReverse = outputNum
End Function

'-----------------------------------------------------------------------------------
'バタフライ演算
'-----------------------------------------------------------------------------------
Private Function ButterflyOperation(ByRef iData() As Complex, ByVal iRdx As Integer, ByVal Stg) As Complex()
    Dim pear As Long
    Dim blok As Long
    
    Dim rtn() As Complex         '戻り値は複素数型配列
    ReDim rtn(UBound(iData))    '入力データと同じサイズの戻り値配列サイズを確保
    
    pear = 2 ^ (Stg - 1)    '1,2,4,8, ...
    blok = 2 ^ (iRdx - Stg)
    
    'バタフライ演算
    pt = 0
    For i = 1 To blok
        For j = 1 To pear '複素数で加算
            rtn(pt) = fComplexAdd(iData(pt), iData(pt + pear))
            pt = pt + 1
        Next
        For j = 1 To pear '複素数で減算
            rtn(pt) = fComplexSub(iData(pt), iData(pt - pear))
            pt = pt + 1
        Next
    Next
    ButterflyOperation = rtn
End Function

'----------------------------------------------------------------------------------------------
'配列に回転因子を掛ける処理，バタフライ演算の前処理
'----------------------------------------------------------------------------------------------
Private Function TwiddleFactorMulti(ByRef iData() As Complex, ByVal rdx As Integer, ByVal StageNum, ByRef twiddle() As Complex) As Complex()
    Dim BufTw() As Complex       '回転因子を掛けた結果を入れる複素数型のバッファ配列
    Dim blok    As Integer      'ブロックサイズ：バタフライ演算の塊が何個あるかを示す
    Dim pear    As Integer      'タスキ数：ブロックの中にたすき掛けが何個あるか
    Dim pt As Long
    
    blok = 2 ^ (rdx - StageNum)
    pear = 2 ^ (StageNum - 1)
    ReDim BufTw(UBound(iData))
    
    pt = 0
    For i = 1 To blok
        '前半：回転因子乗算なし
        For j = 1 To pear
            BufTw(pt) = iData(pt)
            pt = pt + 1
        Next
            
        '後半：回転因子乗算をする
        Ptw = 0
        For j = 1 To pear
            BufTw(pt) = fComplexMulti(iData(pt), twiddle(Ptw)) '複素数の掛け算
            Ptw = Ptw + blok    '回転因子のインデックスはブロックサイズ分飛ばす
            pt = pt + 1
        Next
    Next
    TwiddleFactorMulti = BufTw
End Function

'--------------------------------------------------------------------------------
'複素数の演算を実行する関数
'--------------------------------------------------------------------------------
'複素数の掛け算
Private Function fComplexMulti(ByRef A As Complex, ByRef B As Complex) As Complex
    Dim rtn As Complex
    rtn.re = (A.re * B.re) - (A.Im * B.Im)
    rtn.Im = (A.re * B.Im) + (A.Im * B.re)
    fComplexMulti = rtn
End Function
'複素数の足し算
Private Function fComplexAdd(ByRef A As Complex, ByRef B As Complex) As Complex
    Dim rtn As Complex
    rtn.re = A.re + B.re
    rtn.Im = A.Im + B.Im
    fComplexAdd = rtn
End Function
'複素数の引き算
Private Function fComplexSub(ByRef A As Complex, ByRef B As Complex) As Complex
    Dim rtn As Complex
    rtn.re = A.re - B.re
    rtn.Im = A.Im - B.Im
    fComplexSub = rtn
End Function

実行方法

以下３つのスッテプを通してFFT演算を実行できます．

1．クラスをインスタンスする
2. SetSignalDataメソッドの引数に入力信号となるDouble型の配列を与える
3. Runメソッドを実行するとFFT演算の結果が返ってきます

テスト用のコードはこんな感じです，ExcelのD列に信号を用意してE列にFFTの結果を書き出しています．

Sub FFT実行テスト()
    Dim FFT As clsFftBasic
    Set FFT = New clsFftBasic
    Const L = 4096
    
    Dim Data(L - 1) As Double
    
    'FFTに入力するデータをシートから取得
    For i = 0 To L - 1
        Data(i) = Range("D2").Offset(i, 0)
    Next
    
    'FFT実行
    FFT.SetSignalData (Data)
    f = FFT.Run()
    
    'FFT結果をセルに出力
    For i = 0 To (L / 2) - 1
         Range("E2").Offset(i, 0) = f(i)
    Next
    
End Sub

事項結果

今回は５つの正弦波とランダムノイズをミックスした4096ポイントの入力波形を用意し，Excelの分析ツールFFTと今回作成したコードの実行結果をならべて比較しました
小数の誤差はありますが，同様の結果を得ています