cuFFTパラメータ編

## 概要
* cuFFTで主に使用するパラメータの紹介

## はじめに
最初に言います。
**「cuFFTまじでむずい！！」**
少し扱う機会があったので、勉強をしてみたのですが最初使い方が本当にわかりませんでした。
今もわからない部分はありますが、忘れたくないので勉強したことを記事にしていこうと思います。
間違っている部分があれば、遠慮なく教えていただけると幸いです。

## 関連記事

>* [cufftComplexに画素値を代入してみる](https://qiita.com/toroi-ex/items/66d3985c844c71f9cd79)

## cuFFTとは
まずcuFFTとは何かと言いますと
「GPU上で離散フーリエ変換を行うCUDAのライブラリ」
です。
CUDAのドキュメントにもIntroductionのところに説明が書いてあります。

>**英語：**
The cuFFT library provides a simple interface for computing FFTs on an NVIDIA GPU, which allows users to quickly leverage the floating-point power and parallelism of the GPU in a highly optimized and tested FFT library.

>**日本語(DeepL)：**
cuFFTライブラリは、NVIDIA GPU上でFFTを計算するためのシンプルなインターフェースを提供し、高度に最適化されテストされたFFTライブラリでGPUの浮動小数点演算能力と並列性を迅速に活用することを可能にします。

* [cuFFTドキュメント](https://docs.nvidia.com/cuda/cufft/index.html#introduction)

## cuFFTで主に使うであろうもの
cuFFTの大まかな流れとしては

1. cuFFT planの保存とアクセスに使用するハンドルを定義
1. planの詳細なパラメータを定義
1. FFTを実行

みたいな感じです。途中のホストからデバイスにデータをコピーするといった処理は省略していますが、ざっくり書くと上記のようになります。
少し補足をすると、「plan」とは「CUFFTプランの保存とアクセスに使用されるハンドル型」です。わかりやすく言い換えると、フーリエ変換をするときにこのplanを介して行うみたいな感じです。

### plan関数

#### cufftPlan1D
1次元信号をフーリエ変換するときに定義する「plan」で、パラメータとしてはこんな感じ。

```.cpp
cufftPlan1d(cufftHandle *plan, int nx, cufftType type, int batch);
```
|Input||
|:--                |:--|
|plan              |ハンドル|
|nx              |データのサイズ（信号の長さ）|
|type              |フーリエ変換の方法を選択[（下記参照）](#typeについて)|
|batch       |Inputのフーリエ変換計算回数：基本はnxと同じ値|

#### cufftPlan2D
2次元のデータをフーリエ変換するときに定義する「plan」で、パラメータとしてはこんな感じ。
画像のフーリエ変換などに使用するはず。

```.cpp
cufftPlan2d(cufftHandle *plan, int nx, int ny, cufftType type);
```

|Input||
|:--                |:--|
|plan              |ハンドル|
|nx              |データのサイズ（データの大きさ：画像だとwidth）|
|ny|データのサイズ（データの大きさ：画像だとheight）|
|type              |フーリエ変換の方法を選択[（下記参照）](#typeについて)|

※`cufftPlan2d`にはbatchのオプションはないので、もしbatchを使いたければ`cufftPlanMany`を使用すること

#### cufftPlan3d
3次元データをフーリエ変換するときに定義する「plan」で、パラメータとしてはこんな感じ。

```.cpp
cufftPlan3d(cufftHandle *plan, int nx, int ny, int nz, cufftType type);
```

|Input||
|:--                |:--|
|plan              |ハンドル|
|nx              |xのデータのサイズ|
|ny|yのデータのサイズ|
|nz|zのデータのサイズ|
|type              |フーリエ変換の方法を選択[（下記参照）](#typeについて)|

#### cufftPlanMany
これはわかりにくいので詳しく説明します。まずパラメータがこんな感じです。

```.cpp
cufftPlanMany(cufftHandle *plan, int rank, int *n, int *inembed,
        int istride, int idist, int *onembed, int ostride,
        int odist, cufftType type, int batch);
```

|Input||
|:--                |:--|
|plan              |ハンドル|
|rank|入力データの次元（1次元～3次元）|
|n              |入力データの大きさ・フーリエ変換のカーネルサイズ<font color="gray">（おそらくこの解釈で合ってるはず）</font>|
|inembed|入力するpixelの幅（画像であればwidth） or NULL|
|istride|入力するブロックの間隔 or NULL|
|idist|入力するブロックの大きさ or NULL|
|onembed|出力するpixelの幅（画像であればwidth） or NULL|
|ostride|出力するブロックの間隔 or NULL|
|odist|出力するブロックの大きさ or NULL|
|type              |フーリエ変換の方法を選択[（下記参照）](#typeについて)|
|batch|Inputのフーリエ変換計算回数|

とパラメータがたくさんあります。`cufftPlanMany`は**Advanced Data Layout**というもので、ユーザーが自由にデータの入出力を調整することができます。例えば、640×300の画像データを7×7に分割して入力したり、640×300の画像データを7×7の分割データで出力したりすることができます。
簡潔に言えば、`cufftPlan1d`とかの設定をより細かくできるみたいな感じです。
それでは、それぞれのパラメータの説明に入ります。

##### rank
これは単純に、入力データが何次元かを指定するものです。1次元～3次元で選択します。

##### n
このnはわかりにくいです。まずは公式ドキュメントの内容を見てみます。

> 英語：
Array of size rank, describing the size of each dimension, n[0] being the size of the outermost and n[rank-1] innermost (contiguous) dimension of a transform.

> 日本語（DeepL）：
大きさ rank の配列で，各次元の大きさを表す．n[0] は変換の一番外側の次元， n[rank-1] は一番内側の（連続した）次元の大きさである．

「うーん。。。つまりどういうことだ？」
最初見たときに↑みたいな感じになりました。
この`n`が意味するのはフーリエ変換に使用する**kernel サイズ**を意味しているみたいです。
例えば、2次元データである画像に対してフーリエ変換を行うときに、

```.cpp
int rank = 2;
int n[rank] = {7, 7}; // nの配列サイズはrank
```

と設定してあげれば、7×7のカーネルでフーリエ変換ができるということです。

##### inembed
これはフーリエ変換を行うブロックのPixel幅を設定します。
公式ドキュメントでは以下のように説明されています。

> **英語：**
The inembed and onembed parameters define the number of elements in each dimension in the input array and the output array respectively. The inembed[rank-1] contains the number of elements in the least significant (innermost) dimension of the input data excluding the istride elements; the number of total elements in the least significant dimension of the input array is then istride*inembed[rank-1]. The inembed[0] or onembed[0] corresponds to the most significant (that is, the outermost) dimension and is effectively ignored since the idist or odist parameter provides this information instead. Note that the size of each dimension of the transform should be less than or equal to the inembed and onembed values for the corresponding dimension, that is n[i] ≤ inembed[i], n[i] ≤ onembed[i]，where i∈{0,…,rank−1}

> **日本語（DeepL）：**
inembedとonembedパラメータは，それぞれ入力配列と出力配列の各次元の要素数を定義します．inembed[rank-1] には，入力データの最下位（最内周）次元の要素数（istride 要素を除く）が入ります．したがって，入力配列の最下位次元の総要素数は istride*inembed[rank-1] となります．inembed[0] または onembed[0] は，最も大きい（つまり一番外側の）次元に対応し， idist または odist パラメータが代わりにこの情報を提供するので，実質的には無視されます．変換の各次元のサイズは，対応する次元の inembed と onembed の値以下でなければならないことに注意してください．

ドキュメントでは`inembed`は入力配列の要素数を定義すると記載されています。
またまた画像を例えに用いますと、`inembed`は入力する画像の幅を設定するということになります。
ここで、公式ドキュメントに記載してある式を以下に示します。

> **1D：**
input[ b * idist + x * istride]
output[ b * odist + x * ostride]

> **2D：**
input[b * idist + (x * inembed[1] + y) * istride]
output[b * odist + (x * onembed[1] + y) * ostride]

> **3D：**
input[b * idist + ((x * inembed[1] + y) * inembed[2] + z) * istride]
output[b * odist + ((x * onembed[1] + y) * onembed[2] + z) * ostride]

上記の式はInput、Outputされる値を表したものになります。
この式を見るとお気づきの方もいるかもしれませんが、`inembed[0]`の値は無視されています。なので、`inembed`で意味のある値は**inembed[1]**以上の値になります。
ただ、`that is n[i] ≤ inembed[i], n[i] ≤ onembed[i]`という記載があるのが少し気にかかりはします。inembed[0]は意味を持たないのに、なぜ条件に含まれているのかが不思議です。

`inembed`は下記のように定義します。

上記の式の意味合いなどは別の記事を記載いたしますので少々お待ちを。。。<font color="gray">（長くなってしまいそうなので）</font>

```.cpp
int inembed[rank] = {1, width};
```

`n`と`inembed`を2Dで使う際のイメージ図を作ってみました。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/654472/3af620c5-d3c9-54fb-d578-65a7fc23ba3b.png)

##### idist
これは、入力する要素の数を定義します。1個ずつでよければ1を設定し、7×7で入力したければ49を設定します。ただ、公式ドキュメントでは

> idist または odist パラメータが代わりにこの情報を提供するので，実質的には無視されます

と書かれているため、重要なのは**odist**の定義の方になります。

##### istride
この`istride`は深層学習とかで畳込み層にて使用するstrideとは意味合いが少し違います。
これは図を見てもらえばわかりやすいかなと思うので先に図を載せます。1の場合と2の場合を載せています。
図の補足として、四角内の数字は入力順を表しています。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/654472/4c321068-177a-b984-6863-c2328bb117b6.png)

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/654472/c710b794-7455-ce4a-a13d-3b5fbbafe8a4.png)

言葉で表すと、**入力する形を一つのブロック(kernel×kernel)とし、そのブロックの入力間隔をistrideで調整する**ということになるのかなと思います。
図では

* `istride=1`であれば、入力はブロック順の通り
* `istride=2`であれば、入力は一個飛ばし

のようになります。
それでは次に出力に使用するパラメータの説明に入っていきましょう。

##### onembed
**onembed**は**inembed**とそんなに違いはありません。出力のwidthを定義する感じです。**inembed**と同様に0番目の値は無視されるので、2Dで使用すると

```.cpp
int onembed[rank] = {1, width}
```

のようになります。

##### odist
`odist`は出力する要素の数を定義します。

```.cpp
int odist = 1; // 要素1個を1ブロックとする場合

int odist = 49; // 要素49個を1ブロックとする場合（7×7）
```

##### ostride
これも`istride`とほぼ同じ役割です。出力するブロックの間隔を定義します。

##### batch
フーリエ変換の計算回数を定義します。正確には並列処理で実行する個数の定義を行います。batchはGPUのスペックが許す限りは大きくしていいのかなと思います。

> 英語：
where batch denotes the number of transforms that will be executed in parallel, rank is the number of dimensions of the input data

> 日本語（DeepL）：
ここで、batch は並列に実行される変換の数を表し、rank は入力データの次元数を表す

#### typeについて
planを作成する際にはフーリエ変換の方法を選択する必要があります。

|type||
|:--                |:--|
|CUFFT_R2C|Float型の実数からFloat型の複素数に変換|
|CUFFT_C2R|Float型の複素数からFloat型の実数に変換|
|CUFFT_C2C           |Float型の複素数からFloat型の複素数に変換|
|CUFFT_D2Z|Double型の実数からDouble型の複素数に変換|
|CUFFT_Z2D|Double型の複素数からDouble型の実数に変換|
|CUFFT_Z2Z|Double型の複素数からDouble型の複素数に変換|

上記の中から実装環境に合わせて選択します。

#### Plan作成のエラーチェック
上記で紹介してきたplanには戻り値があり、それでplanの作成に成功したか失敗したか知ることができます。
まず戻り値の内容としては公式ドキュメントより以下のようになってます。

```
typedef enum cufftResult_t {
    CUFFT_SUCCESS        = 0,  //  The cuFFT operation was successful
    CUFFT_INVALID_PLAN   = 1,  //  cuFFT was passed an invalid plan handle
    CUFFT_ALLOC_FAILED   = 2,  //  cuFFT failed to allocate GPU or CPU memory
    CUFFT_INVALID_TYPE   = 3,  //  No longer used
    CUFFT_INVALID_VALUE  = 4,  //  User specified an invalid pointer or parameter
    CUFFT_INTERNAL_ERROR = 5,  //  Driver or internal cuFFT library error
    CUFFT_EXEC_FAILED    = 6,  //  Failed to execute an FFT on the GPU
    CUFFT_SETUP_FAILED   = 7,  //  The cuFFT library failed to initialize
    CUFFT_INVALID_SIZE   = 8,  //  User specified an invalid transform size
    CUFFT_UNALIGNED_DATA = 9,  //  No longer used
    CUFFT_INCOMPLETE_PARAMETER_LIST = 10, //  Missing parameters in call
    CUFFT_INVALID_DEVICE = 11, //  Execution of a plan was on different GPU than plan creation
    CUFFT_PARSE_ERROR    = 12, //  Internal plan database error 
    CUFFT_NO_WORKSPACE   = 13  //  No workspace has been provided prior to plan execution
    CUFFT_NOT_IMPLEMENTED = 14, // Function does not implement functionality for parameters given.
    CUFFT_LICENSE_ERROR  = 15, // Used in previous versions.
    CUFFT_NOT_SUPPORTED  = 16  // Operation is not supported for parameters given.
} cufftResult;
```

肝心のどうやって値が返ってくるかわかりにくい。。笑
取得できる値としては、0,1,2,...の数字が取得できます。この数字なんですが、`cufftResult`という型なので

```.cpp
cufftResult cufftR = cufftPlanMany(&plan, rank, n, inembed, istride, idist, onembed, ostride, odist, CUFFT_C2C, batch);
```

のようする必要があります。あとは、Swich文やif文でリザルトの内容を表示するものを作ってあげればわかりやすくなるかなと思います。
もちろん`CUFFT_SUCCESS `であればplanの作成は成功してます。

### フーリエ変換関数
最後にフーリエ変換を実際に行う関数の説明を行っていきます。この関数もplanと同じ方法で戻り値の確認をすることができます。

#### cufftExecR2C

実数から虚数に変換するフーリエ変換です。

```.cpp
cufftExecR2C(cufftHandle plan, cufftReal *idata, cufftComplex *odata);
```

|Input||
|:--                |:--|
|plan              |ハンドル|
|idata|入力データ|
|odata             |出力データ|

出力データのところにidataを指定するとインプレイス変換され、idataに出力データが保存されます。

#### cufftExecC2R

虚数から実数に変換するフーリエ変換です。R2Cでフーリエ変換した後に、C2Rで逆フーリエ変換するパターンが多いのかなと思います。

```.cpp
cufftExecC2R(cufftHandle plan, cufftComplex *idata, cufftReal *odata);
```

|Input||
|:--                |:--|
|plan              |ハンドル|
|idata|入力データ|
|odata             |出力データ|

出力データのところにidataを指定するとインプレイス変換され、idataに出力データが保存されます。

#### cufftExecC2C

複素数から複素数に変換するフーリエ変換です。

```.cpp
cufftExecC2C(cufftHandle plan, cufftComplex *idata, cufftComplex *odata, int direction);
```

|Input||
|:--                |:--|
|plan              |ハンドル|
|idata|入力データ|
|odata             |出力データ|
|direction|CUFFT_FORWARD or CUFFT_INVERSEのどちらかを指定する（FORWARDがフーリエ変換、INVERSEが逆フーリエ変換）|

出力データのところにidataを指定するとインプレイス変換され、idataに出力データが保存されます。
`cufftComplex`の型は、xに実数部、yに虚数部が設定されています。↓非常に簡略化してるけどこんな感じで値は参照可能。

```.cpp
cufftComplex *idata;
x = idata[10].x;
y = idata[10].y;
```

あとはダブル型の関数もありますが、全く一緒なのでリンクだけ貼っておきます。


* [cuFFT Execution](https://docs.nvidia.com/cuda/cufft/index.html#plan-execution)


## 終わりに
コードの実装も書くと記事が長くなるので一旦ここでこの記事は締めます。
これからもこれ必要だなと感じたら説明を加えて行きます。