工作と競馬

電子工作、プログラミング、木工といった工作の記録記事、競馬に関する考察記事を掲載するブログ

カテゴリ:4.ソフトウェア > 4.1 ExcelVBA

概要

 高速フーリエ変換(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で実装し、動作することを確認できた。今後の信号処理に活用したい。

概要

 ExcelVBAにおける複素数クラスモジュールを自作し、複素数計算のコーディングができることを確認した。(2015/03/01 4.動作確認に説明を追加)


背景と目的

 私はExcelVBAをよく使用するが、電気回路やディジタル信号処理の計算を行うときに複素数を扱えると便利だろうということに気づいた。そこで、ExcelVBAで複素数の計算を出来るようにしてみる。


詳細
1.構想

 始めに思いついたのは、実部と虚部の2つの要素を持つ構造体だったが、複素数は極形式で表したりすることもあるし、四則演算やべき乗といった計算はつき物なので、これらをまとめたクラスを作成するのが、いいのではないかと考えた。そこで、今回は、複素数クラスを作成することにする。

2.設計

 複素数クラスは以下のように設計した。プロパティは、実部、虚部だけでなく絶対値と偏角を持たせた。

  • Real:実部
  • Imaginary:虚部
  • Absolute:絶対値
  • Argument:偏角

 また、メソッドは以下のとおり。よく使う演算を用意した。

  • imsum:加算
  • imsub:引き算
  • improduct:乗算
  • imdiv:除算
  • impower:べき乗
3.実装

 VBAのクラスモジュールを使用し、実装したソースコードは以下のとおり。

Option Explicit '変数の宣言を強制

'変数の宣言部----------------------------------------------------------------------------------------------------------
Private m_re As Double '実部
Private m_im As Double '虚部
Private m_abs As Double '絶対値
Private m_arg As Double '偏角[±πrad]
Private pi As Double

'プロパティの宣言部----------------------------------------------------------------------------------------------------
Public Property Get Real() As Double
    Real = m_re
End Property
Public Property Get Imaginary() As Double
    Imaginary = m_im
End Property
Public Property Get Absolute() As Double
    Absolute = m_abs
End Property
Public Property Get Argument() As Double
    Argument = m_arg
End Property

'コンストラクタ
Private Sub Class_Initialize()

    SetReIm 1, 0
    
    pi = Math.Atn(1) * 4

End Sub
'デストラクタ
Private Sub Class_Terminate()

End Sub
'初期化メソッド
Public Sub Init(Optional a_re As Double = 1, Optional a_im As Double = 0)

    SetReIm a_re, a_im

End Sub

'値をセットする(複素形式)
Public Sub SetReIm(a_re As Double, a_im As Double)

    '値の格納
    m_re = a_re
    m_im = a_im

    '絶対値の計算
    Call imabs
    
    '偏角の計算
    Call imarg
    
End Sub

'値をセットする(極形式)
Public Sub SetAbsArg(a_abs As Double, a_arg As Double)

    '値の格納
    m_abs = a_abs
    If a_arg > 0 Then
        m_arg = a_arg - Int((a_arg + pi) / (2 * pi)) * 2 * pi
    ElseIf a_arg < 0 Then
        m_arg = a_arg - Int((a_arg - pi) / (2 * pi) + 1) * 2 * pi
    Else
        m_arg = a_arg
    End If

    '絶対値の計算
    Call imreal
    
    '偏角の計算
    Call imimaginary
    
End Sub

'各種演算
'和
Public Function imsum(a_cn As ComplexNumber) As ComplexNumber

    Set imsum = New ComplexNumber

    imsum.SetReIm m_re + a_cn.Real, m_im + a_cn.Imaginary

End Function

'積
Public Function improduct(a_cn As ComplexNumber) As ComplexNumber

    Set improduct = New ComplexNumber

    improduct.SetAbsArg m_abs * a_cn.Absolute, m_arg + a_cn.Argument

End Function

'差
Public Function imsub(a_cn As ComplexNumber) As ComplexNumber

    Set imsub = New ComplexNumber

    imsub.SetReIm m_re - a_cn.Real, m_im - a_cn.Imaginary

End Function

'商
Public Function imdiv(a_cn As ComplexNumber) As ComplexNumber

    Set imdiv = New ComplexNumber

    imdiv.SetAbsArg m_abs / a_cn.Absolute, m_arg - a_cn.Argument

End Function

'べき乗
Public Function impower(a_n As Double) As ComplexNumber

    Set impower = New ComplexNumber

    impower.SetAbsArg m_abs ^ a_n, m_arg * a_n

End Function

'絶対値を計算する
Private Sub imabs()

    m_abs = Math.Sqr(m_re ^ 2 + m_im ^ 2)

End Sub

'偏角(±πで表示)を計算する
Private Sub imarg()

    '結果は0~±π/2の値を返す
    'atn関数は実部虚部の正負が判断できないから
    '±π/2~±πまでの間は実部、虚部の符号を考えて
    '場合わけする必要がある

    If m_re > 0 Then
        '実部が正のとき(0~±π/2)
        m_arg = Math.Atn(m_im / m_re)
    ElseIf m_re < 0 Then
        '実部が負のとき(±π/2~±π)
        If m_im > 0 Then
            '虚部が正のとき(+π/2<θ<+π)
            m_arg = Math.Atn(m_im / m_re) + Math.Atn(1) * 4
        Else
            '虚部が負または0のとき(-π/2>θ>=-π)
            m_arg = Math.Atn(m_im / m_re) - Math.Atn(1) * 4
        End If
    Else
        '実部が0のとき(0除算を避ける)
        If m_im > 0 Then
            '虚部が正のとき(+π/2)
            m_arg = Math.Atn(1) * 2
        ElseIf m_im < 0 Then
            '虚部が負のとき(-π/2)
            m_arg = -Math.Atn(1) * 2
        Else
            '虚部が0のとき(本当は不定形だが0とする)
            m_arg = 0
        End If
    End If

End Sub

'実部を計算する
Private Sub imreal()

    m_re = m_abs * Math.Cos(m_arg)

End Sub

'虚部を計算する
Private Sub imimaginary()

    m_im = m_abs * Math.Sin(m_arg)

End Sub

4.動作確認

 動作確認は、以下のプロシージャを作成して行った。(念のため、手順を書くと、同一プロジェクト内に、上記ComplexNumberクラスモジュールをインポートおよび標準モジュールを追加。標準モジュール内に以下のコードを書き、実行する。詳細は図1を参照)

Public Sub test_cn()

    'クラスのオブジェクトを作成
    Dim cn1 As New ComplexNumber
    Dim cn2 As New ComplexNumber
    Dim cn3 As New ComplexNumber
    
    'cn1,cn2にそれぞれ大きさ1、偏角π/4を与える
    cn1.SetAbsArg 1, Math.Atn(1) '1/√2 + j1/√2
    cn2.SetAbsArg 1, Math.Atn(1) '1/√2 + j1/√2
    
    '足し算のメソッドimsumを試す
    Set cn3 = cn1.imsum(cn2) 'cn3 = √2 + j√2となる
    
    '結果確認
    'イミディエイトウインドウに、cn3の各プロパティを出力して確かめる
    Debug.Print cn3.Real '実部(√2)
    Debug.Print cn3.Imaginary '虚部(√2)
    
    Debug.Print cn3.Absolute '絶対値(2)
    Debug.Print cn3.Argument * 180 / Math.Atn(1) / 4 '偏角(π/4)

End Sub
20150301224646
図1 動作確認した結果


まとめと今後の課題

 ExcelVBAで使用できる複素数クラスモジュールを作成し、動作確認できた。今後はこれを用いてディジタル信号処理、電気回路の計算を行いたい。

このページのトップヘ