はじめに
fft を含む MATLAB コードも MATLAB Coder で C コードに変換できるんですが、特に何も指定しないと FFTW ライブラリを使わない C コードが生成されます。
そうなると当然ながら処理速度はあまり期待できません。公式ページ 1 にある通り、FFTW ライブラリを使ったコードを生成することは可能ですので、この記事では実際に試して処理速度の比較もしてみます。
やったこと
-
fftを使った関数の FFTW ライブラリを使った C コード自動生成 - ざっくり速度比較：MATLAB vs 自動生成 C コード
結果
サイズ 1,000,000 x 1 の double 型配列を fft で処理するのにかかる時間を計測。
| 処理時間 | |
|---|---|
| MATLAB | 0.01 秒 |
| C（FFTW あり） | 0.08 秒 |
| C（FFTW なし） | 0.8 秒 |
Wisdom の再利用が効いている2というのもありますが MATLAB 速いですね・・。3
動機
MATLAB で CFD を実装（関連記事へのリンク）していて、C コード生成すれば多少は高速化できるのではと期待。fft を活用していたので、FFTW ライブラリを使ったコード生成が必要でした。
環境
- Windows 10
- MATLAB R2020a (+MATLAB Coder)
- Microsoft Visual C++ 2017
- FFTW 3.3.5
コード生成対象
ここはシンプルにこんな関数。
function y = myfft(x)
%#codegen
tmp = fft(x);
y = abs(tmp);
end
入力 x は double の1次元配列とします。abs を入れている理由は double 型の出力にしたかった為4。
abs の処理時間も入ってしまうのが気になります。MATLAB で実行した場合下図の通りプロファイラーでみたところ影響は小さくない。（for ループで 50 回実行した際の時間）
ですので、fft 部分の処理時間だけをちゃんと比較できるように、abs 部分だけ分けて（下記、複素数入力）コード生成したものの処理速度も計測しておきました。
function y = check_abs(tmp)
%#codegen
y = abs(tmp);
end
ステップ１：FFTW ライブラリのインストール
これは FFTW のページから。
http://www.fftw.org/install/windows.html
Windows 向けにコンパイル済みのものが用意されていますのでこれを使います。 fftw-3.3.5-dll64.zip (3.1MB) をダウンロード、解凍して適当なところに置くだけ。以下に置いてみました。
C:\Users\minoue\fftw-3.3.5-dll64
以上！
・・・と思いきや、
These DLLs were created by us, cross-compiled from GNU/Linux using MinGW; the 64-bit version is possible thanks to the mingw-w64 project. You should be able to call them from any compiler. In order to link to them from Visual C++, you will need to create .lib "import libraries" using the lib.exe program included with VC++.
と書かれている通り、用意されている .dll ファイルに加えて .lib インポート ライブラリが必要だそうな。
lib /def:libfftw3-3.def
lib /def:libfftw3f-3.def
lib /def:libfftw3l-3.def
を fftw のルートディレクトリ（C:\Users\minoue\fftw-3.3.5-dll64）で実行すれば OKです。ただ lib ってどこにあるんだ？
自分の環境ではこちらにありました。
C:\Program Files (x86)\Microsoft Visual Studio\2017\Professional\VC\Tools\MSVC\14.11.25503\bin\HostX64\x64
PATH に追加したうえで上のコードをコマンドプロンプトから実行して完了。（MATLAB 上から ! を使ってやってもいい。）
参考：Stack overflow: lib.exe, Visual Studio, generating .lib files from dll's and def files
ステップ２：FFT コールバック クラスの記述
これは公式ページ（生成されたスタンドアロン コードでの FFTW ライブラリ呼び出しを使用した高速フーリエ変換の高速化）に載っている useMyFFTW を転用。
変更点は以下の部分、FFTW ライブラリ等へのパスだけ。今回のケースだと
fftwLocation = '/usr/lib/fftw';
includePath = fullfile(fftwLocation, 'include');
buildInfo.addIncludePaths(includePath);
libPath = fullfile(fftwLocation, 'lib');
を
fftwLocation = 'C:\Users\minoue\fftw-3.3.5-dll64';
includePath = fullfile(fftwLocation);
buildInfo.addIncludePaths(includePath);
libPath = fullfile(fftwLocation);
に変えて保存します。
ステップ３：コード生成
MATLAB Coder アプリで実施すると楽。ここでは MATLAB で検証するために dll を作成します。
【詳細設定】の【カスタム FFT ライブラリのコールバック】に先程こしらえた useMyFFTW と記載して、あとは通常通り。useMyFFTW.m は検索パスのどこかにあれば大丈夫です。
このアプリからスクリプト作成すると以下のような内容になっていました。2回目以降はこちら実行すれば C コード生成完了です。
%% クラス 'coder.CodeConfig' の構成オブジェクトを作成します。
cfg = coder.config('dll','ecoder',false);
cfg.GenerateReport = true;
cfg.ReportPotentialDifferences = false;
cfg.CustomFFTCallback = 'useMyFFTW'; % FFTW を使わない場合はコメントアウト
%% エントリ ポイント 'myfft' の引数のタイプを定義します。
RGS = cell(1,1);
ARGS{1} = cell(1,1);
ARGS{1}{1} = coder.typeof(0,[Inf 1],[1 0]);
%% MATLAB Coder を起動します。
codegen -config cfg myfft -args ARGS{1}
入力引数は double 型のスカラー値を想定していますが、入力引数が複素数（check_abs.m の生成用）の場合は、下記に変更します。
ARGS{1}{1} = coder.typeof(1i,[Inf 1],[1 0]);
今回は FFTW を使ったバージョンを myfft_wFFTW、使わないバージョンを myfft_noFFTW として別々の名前で生成して、次の速度比較に使います。
参考まで出てきたコード（FFTW を使ったバージョン）の一部がこちら。
plan = fftw_plan_guru_dft_r2c(1, &iodims, 1, &howmany_iodims, (double *)
&x->data[0], (fftw_complex *)&y->data[0], (unsigned int)
FFTW_PRESERVE_INPUT | (unsigned int)FFTW_ESTIMATE);
fftw_execute(plan);
fftw_destroy_plan(plan);
FFTW のお作法通り、プラン作成、実行、プラン破棄やってますね。
速度検証 MATLAB vs C（FFTW有り）vs C（FFTW無し）
先程作った dll を MATLAB 上から実行して処理速度を確認してみます。ただ、公式には
"関数
loadlibraryを使用して、生成されたダイナミック ライブラリを MATLAB に読み込むことは推奨されず、結果として正しくない動作やクラッシュに至る可能性があります。" (公式ページ)
とのこと。
今回のケースは特に問題なさそうに見えますが非推奨である点は要注意です。通常生成されたコードの動作を MATLAB 上で検証する場合は MEX 関数として出力すればよいのですが、今回注目している fft については公式ページにある通り
MEX 出力の場合、MATLAB® Coder™ は MATLAB が FFT アルゴリズムに使用するライブラリを使用します。
MEX 出力してしまうと FFTW を使わないことができないので今回の比較用途には適しません。ということで公式には推奨されていない方法を使って無理やり検証してみます。
loadlibrary で読み込み
まずは loadlibrary で自動生成した dll を読み込みます。不細工ですが、なぜか dll があるディレクトリに移動しないとロードできない（？）ので、cd しています。
loadlibrary で読み込み（折りたたみ）
if ~libisloaded('myfft_noFFTW')
currdir = pwd;
cd('codegen\dll\myfft_noFFTW');
[notfound,warnings] = loadlibrary('myfft_noFFTW','myfft_noFFTW.h');
cd(currdir);
end
if ~libisloaded('myfft_wFFTW')
currdir = pwd;
cd('codegen\dll\myfft_wFFTW');
[notfound,warnings] = loadlibrary('myfft_wFFTW','myfft_wFFTW.h');
cd(currdir);
end
if ~libisloaded('check_abs')
currdir = pwd;
cd('codegen\dll\check_abs');
[notfound,warnings] = loadlibrary('check_abs','check_abs.h');
cd(currdir);
end
ポインタオブジェクトの作り方
ライブラリの関数シグネチャは libfunctionsview で確認してみます。
libfunctionsview myfft_wFFTW
libfunctionsview check_abs
引数が emxArray_real_TPtr やら emxArray_creal_TPtr らしい。これ、emxArray_real_T/emxArray_creal_T のポインタ。生成されたコードをあさってみると myfft_wFFTW_types.h に定義がありました。それぞれ
typedef struct {
creal_T *data;
int *size;
int allocatedSize;
int numDimensions;
boolean_T canFreeData;
} emxArray_creal_T;
struct emxArray_real_T
{
double *data;
int *size;
int allocatedSize;
int numDimensions;
boolean_T canFreeData;
};
とのこと。
今回は共有 C ライブラリを実行する際に入力・出力はポインター オブジェクトで与えます。ポインタオブジェクトは libpointer 関数で作ります。いろいろ試してみたところ以下で動きそうでした。この辺は要調査とは思いつつ、とりあえず前に進みます。（注：そもそも非推奨な行為）
参考：
emxArray_real_TPtr
まずは double 型のやり取りに使う emxArray_real_TPtr から。
N = 5; % fft にかけるデータ長
% emxArray_real_TPtr を定義するための情報を構造体に纏めます。
myEmxArray.size = [N 1]; % 'size' specifies the dimensions of the input array
myEmxArray.numDimensions = 2;
myEmxArray.data = ones(N,1); %'data' holds the input array
myEmxArray.allocatedSize = N; % 'allocatedSize' is the number of elements in the array
% ポインターオブジェクト作成
sc = libpointer('emxArray_real_TPtr',myEmxArray);
sc.Value
ans =
data: [5x1 double]
size: [5 1]
allocatedSize: 5
numDimensions: 2
canFreeData: 0
入力として使う場合は入力データを data で与えて、出力として使う場合は結果は data に入って帰ってくることになります。
emxArray_creal_TPtr
次に複素数型のやり取りに使う emxArray_creal_TPtr。
上の emxArray_real_TPtr では data は double 型の配列でよかったのですが、複素数は少しややこしい。1つの複素数毎に 1 つの構造体（メンバーに re と im を持つ）を用意する必要があるそうな（公式ページ：生成コードの型への MATLAB 型のマッピング（複素数型））
いったん table 型でデータ作成してから table2struct 関数を使うことにします。
% emxArray_creal_TPt を定義するための情報を構造体に纏めます。
clear myEmxArray
myEmxArray.size = [N 1];
myEmxArray.numDimensions = 2;
tmp = table(ones(N,1),ones(N,1),'VariableNames',["re","im"])
| re | im | |
|---|---|---|
| 1 | 1 | 1 |
| 2 | 1 | 1 |
| 3 | 1 | 1 |
| 4 | 1 | 1 |
| 5 | 1 | 1 |
こんな感じ。5x1 の構造体に変換して data に与えます。
myEmxArray.data = table2struct(tmp);
myEmxArray.allocatedSize = N;
% ポインターオブジェクト作成
scc = libpointer('emxArray_creal_TPtr',myEmxArray);
scc.Value
ans =
data: [1x1 struct]
size: [5 1]
allocatedSize: 5
numDimensions: 2
canFreeData: 0
問題点
5x1 の構造体を与えたはずなのに、なぜか data が 1x1 のスカラー構造体になっています。
これは要追加調査ですが、以下での計算結果を比較してみるとデータの入力には問題はなさそう。出力として使った場合には計算結果を最初の1要素分しか確認できない。この事情から今回の調査では出力が複素数ではなく実数になる関数で速度比較しています。（もし解決方法分かる方いらっしゃれば教えてください。）
計算結果検証
ポインターオブジェクトを作る方法は上で説明しました。以下では関数化して使います。
ポインターオブジェクト作る関数（折りたたみ）
function pointer = createEmxArray_real_TPtr(data)
% emxArray_real_TPtr を定義するための情報を構造体に纏めます。
myEmxArray.size = size(data); % 'size' specifies the dimensions of the input array
myEmxArray.numDimensions = length(size(data));
myEmxArray.data = data; %'data' holds the input array
myEmxArray.allocatedSize = numel(data); % 'allocatedSize' is the number of elements in the array
% ポインターオブジェクト作成
pointer = libpointer('emxArray_real_TPtr',myEmxArray);
end
function pointer = createEmxArray_creal_TPtr(data)
% emxArray_creal_TPtr を定義するための情報を構造体に纏めます。
myEmxArray.size = size(data);
myEmxArray.numDimensions = length(size(data));
tmp = table(real(data),imag(data),'VariableNames',["re","im"]);
myEmxArray.data = table2struct(tmp);
myEmxArray.allocatedSize = numel(data);
% ポインターオブジェクト作成
pointer = libpointer('emxArray_creal_TPtr',myEmxArray);
end
こんな感じ。必要なデータは入力引数として指定できるようにしておきます。
早速計算してみましょう。まずは同じ計算結果が返ってくるかどうかの検証。
N = 1e6; % 配列長さ。
x = rand(N,1); % 乱数で検証
dinput = createEmxArray_real_TPtr(x); % 入力（double）
doutput = createEmxArray_real_TPtr(zeros(N,1)); % 出力（double）
% FFTW 使わないバージョン
calllib('myfft_noFFTW', 'myfft_noFFTW', dinput, doutput);
y1 = doutput.Value.data; % 結果
% FFTW 使うバージョン
calllib('myfft_wFFTW', 'myfft_wFFTW', dinput, doutput);
y2 = doutput.Value.data; % 結果
% MATLAB 関数による計算
y3 = myfft(x);
% 誤差
norm(y1-y3)/N
ans = 3.1643e-16
norm(y2-y3)/N
ans = 6.8226e-17
誤差は許容内ですね。abs にかかる処理時間計測用に用意した check_abs の計算結果も確認しておきます。
入力は複素数なので、createEmxArray_creal_TPtr を使います。
x = complex(rand(N,1),rand(N,1)); % 乱数
cinput = createEmxArray_creal_TPtr(x); % 入力（複素数）
calllib('check_abs', 'check_abs', cinput, doutput);
abs1 = doutput.Value.data;
% MATLAB 関数
abs2 = check_abs(x);
% 誤差
norm(abs1-abs2)/N
ans = 7.7579e-20
問題ないですね。
計算速度比較
それぞれ timeit で計測します。これは複数回実行して中央値を返す関数です。
f = @() calllib('myfft_noFFTW', 'myfft_noFFTW', dinput, doutput);
t_noFFTW = timeit(f)
t_noFFTW = 0.8093
f = @() calllib('myfft_wFFTW', 'myfft_wFFTW', dinput, doutput);
t_wFFTW = timeit(f)
t_wFFTW = 0.0876
f = @() calllib('check_abs', 'check_abs', cinput, doutput);
t_abs = timeit(f)
t_abs = 0.0107
abs にかかった時間を差し引いて fft にかかる処理時間は大体
- C (FFTW 無し）: 0.8 秒
- C (FFTW 有り）: 0.08 秒
であることが分かります。FFTW の効果は大きいですね。次に MATLAB でも計算してみると、MATLAB の fft にかかっている時間は 0.01 秒程度。
x = rand(N,1);
f = @() myfft(x);
t_MATLAB = timeit(f)
t_MATLAB = 0.0211
y = fft(x);
f = @() check_abs(y);
t_abs = timeit(f)
t_abs = 0.0118
abs にかかる時間は C コードと同程度ですね。
Wisdom の影響
以前 FFT アルゴリズムの探索メソッド比較 でも触れた通り、fft 関数を実行するとまず FFT に使用するアルゴリズムを決定する処理が実行されます。MATLAB だと一度探索したアルゴリズム（Wisdom）は再利用されますのであまり心配はないのですが、生成された C コードを見ると fft 実行するたびに実施されています。
この辺も処理時間の差に影響しているだろうと思い、参考まで MATLAB で毎回アルゴリズム探索をした場合の処理時間を調べておきます。
fftw('dwisdom',[]);
というコマンドで wisdom が初期化されるので、次回 fft 実行時にアルゴリズム探索が走ります。
M = 100;
x = rand(N,1);
t1 = zeros(M,1);
for ii=1:M
fftw('dwisdom',[]);
tic
myfft(x);
t1(ii,1) = toc;
end
今度は wisdom を初期化しない、すなわち初回だけアルゴリズム探索が走ります。
t2 = zeros(M,1);
for ii=1:M
% fftw('dwisdom',[]);
tic
myfft(x);
t2(ii,1) = toc;
end
処理時間を比較してみると・・ 毎回アルゴリズム探索したとしても 0.045 秒程度（abs 処理込み）。C (FFTW 有り）の倍くらいの速さがあります。
figure
histogram(t1)
hold on
histogram(t2)
hold off
xlabel('seconds')
legend('plan 毎回作成','plan 再利用')
まとめ
fft を含むコードを何も指定せず MATLAB Coder で C コード生成すると FFTW ライブラリを使わない C コードが生成されます。そうなると当然ながら処理速度はあまり期待できません。この記事では FFTW ライブラリを使ったコード生成を試して、MATLAB 上で検証を行ってみました。
-
fftを使った関数の FFTW ライブラリを使った C コード自動生成 - ざっくり速度比較：MATLAB vs 自動生成 C コード
関数呼び出しのオーバーヘッド、アルゴリズム探索にかかる時間など考慮すべき点は他にもたくさんあるとは思いますが、FFTW ライブラリの力は明白、そして少なくとも MATLAB 上で使うなら MATLAB で計算するのが一番速いという結果でした。
| 処理時間 | |
|---|---|
| MATLAB | 0.01 秒 |
| C（FFTW あり） | 0.08 秒 |
| C（FFTW なし） | 0.8 秒 |
C コード生成した場合には、もし同じサイズの FFT を繰り返す処理があるのであれば wisdom の再利用ができるように書き直した方がよいかもしれませんね。
% 後片付けして終わり。
clear
unloadlibrary myfft_noFFTW
unloadlibrary myfft_wFFTW
unloadlibrary check_abs
Appendix: ヘルパー関数
function pointer = createEmxArray_real_TPtr(data)
% emxArray_real_TPtr を定義するための情報を構造体に纏めます。
myEmxArray.size = size(data); % 'size' specifies the dimensions of the input array
myEmxArray.numDimensions = length(size(data));
myEmxArray.data = data; %'data' holds the input array
myEmxArray.allocatedSize = numel(data); % 'allocatedSize' is the number of elements in the array
% ポインターオブジェクト作成
pointer = libpointer('emxArray_real_TPtr',myEmxArray);
end
function pointer = createEmxArray_creal_TPtr(data)
% emxArray_creal_TPtr を定義するための情報を構造体に纏めます。
myEmxArray.size = size(data);
myEmxArray.numDimensions = length(size(data));
tmp = table(real(data),imag(data),'VariableNames',["re","im"]);
myEmxArray.data = table2struct(tmp);
myEmxArray.allocatedSize = numel(data);
% ポインターオブジェクト作成
pointer = libpointer('emxArray_creal_TPtr',myEmxArray);
end
-
参考：FFT アルゴリズムの探索メソッド比較） ↩
-
C の場合は dll MATLAB から呼び出していることによる追加オーバーヘッドはあるかも。公式ページ：共有ライブラリの C 関数の呼び出し） ↩
-
出力を複素数とした場合、生成した dll を loadlibrary で呼び出しがうまくできませんでした・・・（要勉強） ↩