Excel で FFTW を使えるようにしよう・その2

1、概要

前回は FFTW(高速フーリエ変換) のライブラリ libfftw3-3.dll と Excel を仲介する DLL の作成方法を解説しました。
今回はこの仲介 DLL を使って、Excel の複素関数が使える形式(x+yi)で出力するサンプルプログラムを作成したので、そのプログラムの解説と使い方を簡単に説明します。

Excel の分析ツールに付属するフーリエ解析(FFT) は、データ数が 4096 までの制限があったり２の整数乗でないと解析できませんが、FFTW はその制限が無く２の整数乗でなくても解析できて、 しかも Cooley, Tukey の FFT アルゴリズムより高速に動作します。
FFTW がどんなものかは　http://www.fftw.org　を参照ください。　

1-1、開発環境

・ OS : Windows7 Professional SP1、64bit版
・ Office : Office 2010、32bit 版又は 64bit 版

2、サンプルプログラムの入手

サンプルプログラム一式をここからダウンロード(file size:1,774kByte)ください。

・ダウンロードしたファイルを解凍すると、下記5つのファイルがあります。

1、vb_fftw.xlsm　→ Excel のサンプルプログラム
2、x86\fftw-vbif.dll　→　FFTW へ仲介する Sub と、x+yi からそれぞれのパートへ分割する c_split Function が入った　32bit 版の DLL
3、x86\libfftw3-3.dll　→ fftw.org から入手した 32bit 版の DLL、ここでは double 版を使用。
4、x64\fftw-vbif.dll　→　FFTW へ仲介する Sub と、x+yi からそれぞれのパートへ分割する c_split Function が入った 64bit 版の DLL
5、x64\libfftw3-3.dll　→ fftw.org から入手した 64bit 版の DLL、ここでは double 版を使用。

注) 32bit 版の DLL は 'x86'、64bit 版の DLL は 'x64' という名前のフォルダに入ってます。

3、Excel 側で使う準備

3-1、ランタイムライブラリの入手

注) Visual Studio 2017 をインストールしている場合、本項目の確認や設定は不要です。3-2項へ進んでください。

・コントロールパネルから 'プログラムと機能' を立ち上げ、
32bit 版の Excel を使われている方は 'Microsoft Visual C++2017 Redistributable (x86) - *' がインストールされているかどうか確認、
64bit 版の Excel を使われている方は 'Microsoft Visual C++2017 Redistributable (x64) - *' がインストールされているかどうか確認してください。

・インストールされてない場合は、下記 URL から、Visual Studio 2017 の Microsoft Visual C++ 再頒布可能パッケージを入手してインストールしてください。
https://support.microsoft.com/ja-jp/help/2977003/the-latest-supported-visual-c-downloads

3-2、DLL の環境設定

・Excel 側から DLL を呼び出せるようにパソコンの環境を整える。

・Windows OS が 32bit 版の場合：
x86\fftw-vbif.dll, x86\libfftw3-3.dll のファイルを、'C:\Windows\System32' へ入れる。

・Windows OS が 64bit 版の場合：
Excel が 64bit 版の場合は、x64\fftw-vbif.dll, x64\libfftw3-3.dll のファイルを 'C:\Windows\System32' へ入れる。
Excel が 32bit 版の場合は、x86\fftw-vbif.dll, x86\libfftw3-3.dll のファイルを 'C:\Windows\SysWOW64' へ入れる。

【上記と違うやり方】

・注) 上記の方法で DLL を登録した方は、ここの項目を飛ばして次へ行ってください。

・Windows のシステムフォルダへ入れるのが嫌な人は、環境変数の path に上記2つの DLL を入れたフォルダパスを追加する方法もあります。
Windows7 の場合は、マイコンピュータのプロパティ → システムの詳細設定 → 詳細設定 → 環境変数(N) から path の環境設定を編集してください。
例えば c:\work に入れた場合は、変数値(V) の最後に ';c:\work' を足してＯＫする。

・Excel 側から 'vb_fftw_dft_1d' の DLL をコールしたとき、下記のエラーが出る場合は環境設定が正しく出来てないので、ここの内容を良く見直してください。 Microsoft Visual C++ 再頒布可能パッケージをインストールしていない場合も同じエラーが出ます。

4、使い方の説明

プログラムの起動方法

・ vb_fftw.xlsm を立ち上げると、TOOL タブが表示されますので、その中の FFTW をクリックするとプログラムが立ち上がります。
このリボンは vb_fftw.xlsm のシートを作業対象としている時には出て、他のワークブックへ行くと表示されなくなります。
他のワークブックでもこのリボンを表示したい場合は、ファイルを vb_fftw.xlam のアドイン形式に保存しなおして、アドイン登録して使ってください。
下記は vb_fftw.xlsm を立ち上げ、TOOL タブをクリックしたときの画面です。

パラメータの与え方

・ FFTW をクリックすると、下記のダイアログが立ち上がります。

・ IN には入力元の先頭セルを指定する。 → データの終わりは (xlDown) で検出。
・ OUT には出力先の先頭セルを指定する。
→ 分析ツールのように囲うのが面倒なのでこのようにした、ワークシート判別しているので、IN,OUT は同じワークブック内の違うシートタブを選択しても正常に動作します。
・WINDOW には窓関数の種類を選択する、rectangular を選択すると窓関数は何も処理しない。
・INV. には逆変換したい場合チェックを入れる。
・上記を入力し終わったら、 GO ボタンを押す。

実際の使用例１

・ 前回のシートで作成した、サンプリング周波数 50kHz、生成周波数 1kHz の Sin 5波、ポイント数 250 ポイントのデータで FFT かけてみる。

・ Excel の B 列に元となる Sin 5 波の波形を入れる。
・ FFTW 変換 → IN:B3、OUT:C3、WINDOW:rectangular、Inv.:チェック入れないで GO ボタンを押す。
・ Excel の H3 に '=imabs(C3)' と入れて行コピーする。
・ 下記はグラフ範囲 H3:H128 を指定して表示したグラフです。

実際の使用例2

・ 実際の使用例１の FFT 結果に逆 FFT かけ、波形が元に戻るか確認。

・ 先ほどのデータ、グラフと H 列を削除。
・ FFTW 逆変換 → IN:C3、OUT:H3、WINDOW:rectangular、Inv.:チェック入れて GO ボタンを押す。
・ Excel の I3 に '=imreal(H3)' と入れて行コピーする。
・ 下記はグラフ範囲 B3:B252 と I3:I252 を重ねて表示したグラフです、ぴったり重なった。

実際の使用例3

・ 折り返しが繋がってない波形で FFT かけてみる。

・ 先ほどの Sin 5波、ポイント数 250 ポイントのデータを 240 ポイントに削って窓関数を使って FFT かけてみる。
・ 違いが解りやすいように FFT 後は、 '=20*Log10(imabs())' に式を変更し 20 ポイント分グラフ表示 
・ グラフ青：250 ポイント,Rectangular窓、グラフ赤：240 ポイント,Rectangular窓、グラフ緑：240 ポイント,Hanning窓。
・ 折り返しが繋がってない波形では、窓関数入れるとスペクトルは訛ったが S/N は上がるようだ。

5、プログラムの説明

5-1、DLL 側の変更点

・ Excel VBA で x+yi 形式の文字列分割ルーチンを作成しようとしたが、VBA 標準では正規表現が使えないのでこれまた難儀、実行スピードのことも考えて DLL 側に実装することにした。
fftw-vbif.cpp に足したコードは下記の通り。

#include "stdlib.h"
#include "fftw-vbif.h"

void __stdcall c_split(char **moji, complex *data)
{

	int		n;
	char	*stopstring;

	n = strlen(*moji);
	n--;

	// データの最後が 'i' の場合
	if (*(*moji + n) == 'i') {

		data[0].real = strtod(*moji, &stopstring);

		if (*stopstring == 'i') {
			data[0].image = data[0].real;
			data[0].real = 0.0;
		}
		data[0].image = strtod(stopstring, &stopstring);
	}

	else {
		// Real part のみの場合
		data[0].real = strtod(*moji, &stopstring);
		data[0].image = 0.0;
	}

}
	

・ またヘッダーファイル fftw_vbif.h も新たにプロジェクトへ追加し、下記コードを追加。

// public データ格納用構造体
typedef struct {
	double	real;		// Real part
	double	image;		// Image part
} complex;
	

・ その他エクスポート定義ファイル fftw-vbif.def も編集しましたが、自分で DLL ビルドしてみたい方は、前回の記事を参考に変更ください。

5-2、VBA 側のコード説明

・ 今回 DLL 側に足した c_split sub の使い方を以下に説明します。

手続き宣言部分

・変数 str は VBA 側 String、 C 言語側 char 型なので普通に書けば良い。
変数 in は complex 型なので、Any 型と書いておく。

#If VBA7 And Win64 Then
    Declare PtrSafe Sub c_split Lib "fftw-vbif.dll" (str As String, data As Any)
#Else
    Declare Sub c_split Lib "fftw-vbif.dll" (str As String, data As Any)
#End If
	

・ 64bit 版の Excel で実行する場合は、Declare ステートメントの次に PtrSafe 属性を付けて宣言してください。
#If 文中の VBA7 は Office2010 以降か Office2007 以前か、Win64 は Excel が 64bit 版か 32bit 版か調べられる定数です。
上記の様に書いておくと、Excel のビット数やバージョンが何であれ、共通のコードで実行出来るようになります。

手続きコール部分

・VBA のサンプルコードを下記に示す。
手続きで宣言した str 部分は、 C言語側、VBA 側ともアドレス渡しになっているので普通に書けばよい。
手続きで宣言した in 部分は、C言語側 complex 型のアドレス渡しになっているので in(i) のように書けば、in(i).real, in(i).image に分割された数値が入る。

Dim in(1000) As complex
Dim str As String
Dim i As Long

i = 10

str = "1+2i"
Call c_split( str, in(i))
	

5-3、サンプルコードの説明

Excel のサンプルプログラム、vb_fftw.xlsm のコード説明

・DLL に実装した c_split sub 以外は至って普通のコードなので、直接 VBE を立ち上げて中身見てください。　... 全然説明になってない。
処理の流れは、Dialog 立ち上げて必要なパラメータを入力してもらった後、シートからデータ読込み → 窓関数処理 → FFTW or 逆 FFTW 実行 → シートへデータ書き込みの順に処理が進み、サブルーチンから抜けます。
FFT の使い方は他の文献に譲りますが、正変換はそのまま、逆変換は最終列を 1/n にするとレベルが合います、本サンプルプログラムに実装したコードを下記に示します。
また VBA コード部、シートへの書き込み方法とか、コードの書き方次第でプログラムの実行スピードがとても遅くなったりします、次回この辺り高速に動くプログラムの書き方など説明できればと思ってます。

' FFTW 実行
If dlg.inv = False Then

    ' 正変換
    Call vb_fftw_dft_1d(n, fftw_in(0), fftw_out(0), FFTW_FORWARD, FFTW_ESTIMATE)

Else
    
    ' 逆変換
    Call vb_fftw_dft_1d(n, fftw_in(0), fftw_out(0), FFTW_BACKWARD, FFTW_ESTIMATE)
    
    ' レベル調整(最終列 1/n)
    For i = 0 To n - 1
        fftw_out(i).real = fftw_out(i).real / n
        fftw_out(i).image = fftw_out(i).image / n
    Next i

End If	
	

6、プロジェクトファイル

今回作成した DLL のプロジェクトファイル一式はここにあります。(file size:40kByte)
Visual Studio 2017 を使って、x86 or x64,Release で Build すると fftw-vbif.dll　が作られ、どこのパソコンでも使えるようになります。
(Visual Studio の入っていないパソコン環境で、この DLL を使う場合は Visual C++2017 のランタイムライブラリを入れる必要があります)

7、最後に

このソフトの実行スピードは自分のパソコンで、書き込み先のセルが空の場合は 25万ポイントの計算で 3秒程度、書き込み先のセルのデータを上書きする場合は 6秒程度に伸びます。
FFTW の実行時間は 1秒もかからず、Excel の Range オブジェクトにアクセスする時間の方が長いようです。
このツールもちょっとお勉強するにはちょうど良いですが、仕事で使おうとすると使いづらいかもしれません。皆さんもどんどん自分に合わせて使いやすいように改良して頂ければと思います。

8、参考文献、参考 URL

1、中村尚五著：ビギナーズ デジタルフーリエ変換，東京電機大学出版局(1989)
2、坂巻佳壽美著：見てわかる デジタル信号処理，工業調査会(1998)
3、http://www.fftw.org(本家 FFTW のサイト)
4、http://skomo.o.oo7.jp/f46/hp46_1.htm(同じような DLL の作成方法について書かれてます)

戻る

このサイトについて