FFTW を使用する C コードの自動生成と処理速度比較

# はじめに


`fft` を含む MATLAB コードも [MATLAB Coder](https://jp.mathworks.com/products/matlab-coder.html?s_eid=PSM_29435) で C コードに変換できるんですが、特に何も指定しないと FFTW ライブラリを使わない C コードが生成されます。

そうなると当然ながら処理速度はあまり期待できません。公式ページ [^0] にある通り、FFTW ライブラリを使ったコードを生成することは可能ですので、この記事では実際に試して処理速度の比較もしてみます。

[^0]: [生成されたスタンドアロン コードでの FFTW ライブラリ呼び出しを使用した高速フーリエ変換の高速化](https://jp.mathworks.com/help/coder/ug/speed-up-fast-fourier-transforms-in-generated-standalone-code-by-using-fftw-library-calls.html?s_eid=PSM_29435)


## やったこと

   -  `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 の再利用が効いている[^1]というのもありますが MATLAB 速いですね・・。[^2]

[^1]: 参考：[FFT アルゴリズムの探索メソッド比較](https://qiita.com/eigs/items/b2ba207c923efb490f36)）
[^2]: C の場合は dll  MATLAB から呼び出していることによる追加オーバーヘッドはあるかも。公式ページ：[共有ライブラリの C 関数の呼び出し](https://jp.mathworks.com/help/matlab/matlab_external/calling-functions-in-shared-libraries.html?s_eid=PSM_29435)）

## 動機


MATLAB で CFD を実装（[関連記事へのリンク](https://qiita.com/eigs/items/05ff8224a2301f8d9326#%E6%95%B0%E5%80%A4%E6%B5%81%E4%BD%93%E5%8A%9B%E5%AD%A6)）していて、C コード生成すれば多少は高速化できるのではと期待。`fft` を活用していたので、FFTW ライブラリを使ったコード生成が必要でした。
  
## 環境

   -  Windows 10 
   -  MATLAB R2020a (+MATLAB Coder) 
   -  Microsoft Visual C++ 2017 
   -  FFTW 3.3.5 

  
# コード生成対象


ここはシンプルにこんな関数。



```matlab:Code(Display)
function y = myfft(x)
%#codegen
tmp = fft(x);
y = abs(tmp);     
end
```



入力 `x` は double の1次元配列とします。`abs` を入れている理由は double 型の出力にしたかった為[^3]。

[^3]: 出力を複素数とした場合、生成した dll を loadlibrary で呼び出しがうまくできませんでした・・・（要勉強）




`abs` の処理時間も入ってしまうのが気になります。MATLAB で実行した場合下図の通りプロファイラーでみたところ影響は小さくない。（for ループで 50 回実行した際の時間）



<img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/149511/c48354f1-144c-f31f-9b04-7c59cc1816bc.png" alt="attach:cat" title="attach:cat" width=700px>





ですので、`fft` 部分の処理時間だけをちゃんと比較できるように、`abs` 部分だけ分けて（下記、複素数入力）コード生成したものの処理速度も計測しておきました。


```matlab:Code(Display)
function y = check_abs(tmp)
%#codegen
y = abs(tmp);     
end
```

  
# ステップ１：**FFTW ライブラリのインストール**


これは FFTW のページから。




[http://www.fftw.org/install/windows.html](http://www.fftw.org/install/windows.html)




Windows 向けにコンパイル済みのものが用意されていますのでこれを使います。 fftw-3.3.5-dll64.zip (3.1MB) をダウンロード、解凍して適当なところに置くだけ。以下に置いてみました。



```matlab:Code(Display)
C:\Users\minoue\fftw-3.3.5-dll64
```



以上！




・・・と思いきや、


> These DLLs were created by us, cross-compiled from GNU/Linux using [MinGW](http://www.mingw.org/); the 64-bit version is possible thanks to the [mingw-w64](http://mingw-w64.sourceforge.net/) 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](https://msdn.microsoft.com/en-us/library/7ykb2k5f.aspx) program included with VC++.




と書かれている通り、用意されている `.dll` ファイルに加えて `.lib` インポート ライブラリが必要だそうな。



```matlab:Code(Display)
     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` ってどこにあるんだ？




自分の環境ではこちらにありました。



```matlab:Code(Display)
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](https://stackoverflow.com/questions/5343932/lib-exe-visual-studio-generating-lib-files-from-dlls-and-def-files)


  
# ステップ２：FFT コールバック クラスの記述


これは公式ページ（[生成されたスタンドアロン コードでの FFTW ライブラリ呼び出しを使用した高速フーリエ変換の高速化](https://jp.mathworks.com/help/coder/ug/speed-up-fast-fourier-transforms-in-generated-standalone-code-by-using-fftw-library-calls.html?s_eid=PSM_29435)）に載っている `useMyFFTW` を転用。




変更点は以下の部分、FFTW ライブラリ等へのパスだけ。今回のケースだと



```matlab:Code(Display)
fftwLocation = '/usr/lib/fftw';
includePath = fullfile(fftwLocation, 'include');
buildInfo.addIncludePaths(includePath);
libPath = fullfile(fftwLocation, 'lib');
```



を



```matlab:Code(Display)
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` は検索パスのどこかにあれば大丈夫です。


<img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/149511/3e52458a-37e1-43d8-3159-093ed83722ff.png" alt="attach:cat" title="attach:cat" width=700px>




このアプリからスクリプト作成すると以下のような内容になっていました。2回目以降はこちら実行すれば C コード生成完了です。



```matlab:Code(Display)
%% クラス '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` の生成用）の場合は、下記に変更します。



```matlab:Code(Display)
ARGS{1}{1} = coder.typeof(1i,[Inf  1],[1 0]);
```



今回は FFTW を使ったバージョンを `myfft_wFFTW`、使わないバージョンを `myfft_noFFTW` として別々の名前で生成して、次の速度比較に使います。


参考まで出てきたコード（FFTW を使ったバージョン）の一部がこちら。

```c:生成コードfft.cの一部

    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`](https://jp.mathworks.com/help/matlab/ref/loadlibrary.html?s_eid=PSM_29435) を使用して、生成されたダイナミック ライブラリを MATLAB に読み込むことは推奨されず、結果として正しくない動作やクラッシュに至る可能性があります。" ([公式ページ](https://jp.mathworks.com/help/coder/ug/deploy-generated-code.html?s_eid=PSM_29435))




とのこと。




今回のケースは特に問題なさそうに見えますが非推奨である点は要注意です。通常生成されたコードの動作を MATLAB 上で検証する場合は MEX 関数として出力すればよいのですが、今回注目している `fft` については[公式ページ](https://jp.mathworks.com/help/matlab/ref/fft.html?s_eid=PSM_29435)にある通り


> MEX 出力の場合、MATLAB® Coder™ は MATLAB が FFT アルゴリズムに使用するライブラリを使用します。




MEX 出力してしまうと FFTW を使わないことができないので今回の比較用途には適しません。ということで**公式には推奨されていない方法**を使って無理やり検証してみます。


  
## loadlibrary で読み込み


まずは `loadlibrary` で自動生成した dll を読み込みます。不細工ですが、なぜか dll があるディレクトリに移動しないとロードできない（？）ので、cd しています。


<details><summary>loadlibrary で読み込み（折りたたみ）</summary><div>

```matlab:Code
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
```

</div></details>


## ポインタオブジェクトの作り方


ライブラリの関数シグネチャは `libfunctionsview` で確認してみます。



```matlab:Code
libfunctionsview myfft_wFFTW
```


<img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/149511/a7c31555-fef0-a1d3-b463-7a79662b83cc.png" alt="attach:cat" title="attach:cat" width=500px>




```matlab:Code
libfunctionsview check_abs
```

<img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/149511/cd763fa2-02a5-2939-5910-9157160bc7c4.png" alt="attach:cat" title="attach:cat" width=500px>




引数が `emxArray_real_TPtr` やら `emxArray_creal_TPtr `らしい。これ、`emxArray_real_T/emxArray_creal_T` のポインタ。生成されたコードをあさってみると `myfft_wFFTW_types.h` に定義がありました。それぞれ



```c:Code(Display)
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` 関数](https://jp.mathworks.com/help/matlab/ref/libpointer.html?s_eid=PSM_29435)で作ります。いろいろ試してみたところ以下で動きそうでした。この辺は要調査とは思いつつ、とりあえず前に進みます。（注：そもそも非推奨な行為）

参考：
   -  [公式ページ：ポインター引数](https://jp.mathworks.com/help/matlab/matlab_external/working-with-pointers.html?s_eid=PSM_29435) 
   -  [公式ページ：ポインターの受け渡し例](https://jp.mathworks.com/help/matlab/matlab_external/pass-pointers.html?s_eid=PSM_29435) 




### emxArray_real_TPtr


まずは double 型のやり取りに使う emxArray_real_TPtr から。



```matlab:Code
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
```


```text:Output
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 型のマッピング（複素数型）](https://jp.mathworks.com/help/coder/ug/mapping-matlab-types-to-cc-types.html#mw_d6bb3dc9-7621-4cbd-b360-2b4015f35fb1?s_eid=PSM_29435)）




いったん table 型でデータ作成してから `table2struct` 関数を使うことにします。



```matlab:Code
% 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` に与えます。



```matlab:Code
myEmxArray.data = table2struct(tmp);
myEmxArray.allocatedSize = N;

% ポインターオブジェクト作成
scc = libpointer('emxArray_creal_TPtr',myEmxArray);
scc.Value
```


```text:Output
ans = 
             data: [1x1 struct]
             size: [5 1]
    allocatedSize: 5
    numDimensions: 2
      canFreeData: 0

```

### 問題点


5x1 の構造体を与えたはずなのに、なぜか data が 1x1 のスカラー構造体になっています。




これは要追加調査ですが、以下での計算結果を比較してみるとデータの入力には問題はなさそう。出力として使った場合には計算結果を最初の1要素分しか確認できない。この事情から今回の調査では出力が複素数ではなく実数になる関数で速度比較しています。（もし解決方法分かる方いらっしゃれば教えてください。）


  
## 計算結果検証


ポインターオブジェクトを作る方法は上で説明しました。以下では関数化して使います。

<details><summary>ポインターオブジェクト作る関数（折りたたみ）</summary><div>


```matlab:Code(Display)
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
```



こんな感じ。必要なデータは入力引数として指定できるようにしておきます。

</div></details>



早速計算してみましょう。まずは同じ計算結果が返ってくるかどうかの検証。



```matlab:Code
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
```


```text:Output
ans = 3.1643e-16
```


```matlab:Code
norm(y2-y3)/N 
```


```text:Output
ans = 6.8226e-17
```



誤差は許容内ですね。`abs` にかかる処理時間計測用に用意した `check_abs` の計算結果も確認しておきます。




入力は複素数なので、`createEmxArray_creal_TPtr` を使います。



```matlab:Code
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
```


```text:Output
ans = 7.7579e-20
```



問題ないですね。


## 計算速度比較


それぞれ `timeit` で計測します。これは複数回実行して中央値を返す関数です。



```matlab:Code
f = @() calllib('myfft_noFFTW', 'myfft_noFFTW', dinput, doutput);
t_noFFTW = timeit(f)
```


```text:Output
t_noFFTW = 0.8093
```


```matlab:Code
f = @() calllib('myfft_wFFTW', 'myfft_wFFTW', dinput, doutput);
t_wFFTW = timeit(f)
```


```text:Output
t_wFFTW = 0.0876
```


```matlab:Code
f = @() calllib('check_abs', 'check_abs', cinput, doutput);
t_abs = timeit(f)
```


```text:Output
t_abs = 0.0107
```



`abs` にかかった時間を差し引いて `fft` にかかる処理時間は大体



   -  **C (FFTW 無し）: 0.8 秒**
   -  **C (FFTW 有り）: 0.08 秒**



であることが分かります。FFTW の効果は大きいですね。次に MATLAB でも計算してみると、MATLAB の `fft` にかかっている時間は 0.01 秒程度。



```matlab:Code
x = rand(N,1);
f = @() myfft(x);
t_MATLAB = timeit(f)
```


```text:Output
t_MATLAB = 0.0211
```


```matlab:Code
y = fft(x);
f = @() check_abs(y);
t_abs = timeit(f)
```


```text:Output
t_abs = 0.0118
```



`abs` にかかる時間は C コードと同程度ですね。


## Wisdom の影響


以前 [FFT アルゴリズムの探索メソッド比較](https://qiita.com/eigs/items/b2ba207c923efb490f36) でも触れた通り、`fft` 関数を実行するとまず FFT に使用するアルゴリズムを決定する処理が実行されます。MATLAB だと一度探索したアルゴリズム（Wisdom）は再利用されますのであまり心配はないのですが、生成された C コードを見ると fft 実行するたびに実施されています。




この辺も処理時間の差に影響しているだろうと思い、参考まで MATLAB で毎回アルゴリズム探索をした場合の処理時間を調べておきます。



```matlab:Code
fftw('dwisdom',[]);
```



というコマンドで wisdom が初期化されるので、次回 `fft` 実行時にアルゴリズム探索が走ります。



```matlab:Code
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 を初期化しない、すなわち初回だけアルゴリズム探索が走ります。



```matlab:Code
t2 = zeros(M,1);
for ii=1:M
    % fftw('dwisdom',[]);
    tic
    myfft(x);
    t2(ii,1) = toc;
end
```



処理時間を比較してみると・・ 毎回アルゴリズム探索したとしても 0.045 秒程度（`abs` 処理込み）。C (FFTW 有り）の倍くらいの速さがあります。



```matlab:Code
figure
histogram(t1)
hold on
histogram(t2)
hold off
xlabel('seconds')
legend('plan 毎回作成','plan 再利用')
```

<img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/149511/17fd3d3d-53ca-9ef3-fc9b-facc7de1e0f2.png" alt="attach:cat" title="attach:cat" width=500px>

# まとめ


`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 の再利用ができるように書き直した方がよいかもしれませんね。



```matlab:Code
% 後片付けして終わり。
clear
unloadlibrary myfft_noFFTW
unloadlibrary myfft_wFFTW
unloadlibrary check_abs
```

# Appendix: ヘルパー関数

```matlab:Code
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
```