二次元データの離散フーリエ変換(DFT)と、その高速版である二次元高速フーリエ変換(FFT)について、JavaScriptコードでの解説記事です。
この記事内の各種のDFT2D/FFT2DのJavaScript実装は、以下のgistのファイルfft2d.jsに置きました。

https://gist.github.com/bellbind/b7583ca138a6a1a14763756cd75ccc86

二次元版のDFTとFFTの実装では、一次元版のDFTとFFTの実装を使用します。
一次元版のDFT/FFTの解説は、「高速フーリエ変換FFTを理解する」にあります。
この中の複素数値演算、DFTの計算と、FFTの仕組みは二次元版でもほぼ同じものを使用します。

さらに1次元DFT/FFTから2次元DFT/FFTに拡張するのと同様の手法で、３次元のDFT/FFTも2次元のDFT/FFTを用いて実装できます。上記gistのfft3d.jsは、3次元版のDFT/FFTのJavaScript実装です。

1. ユーティリティ関数の用意

まず、2次元DFT/FFTの実装コードで使用するユーティリティ関数を用意します。

1.1. 複素数値演算の関数

ここではコードの理解のしやすさのため、浮動小数値をフラットな配列には並べません。
フーリエ変換する対象である複素数値は、浮動小数値２個のJavaScript配列で表現し、データではこの複素数値のJavaScript配列を用います。

関数として実装したこれらの複素数演算は、1次元版フーリエ変換で用いるものとまったく同じです。

const zero = () => [0, 0];
const expi = t => [Math.cos(t), Math.sin(t)];
const add = ([ax, ay], [bx, by]) => [ax + bx, ay + by];
const sub = ([ax, ay], [bx, by]) => [ax - bx, ay - by];
const mul = ([ax, ay], [bx, by]) => [ax * bx - ay * by, ax * by + ay * bx];
const divr = ([ax, ay], r) => [ax / r, ay / r];

expi(t)は実数$t$を引数に、複素数 $\exp(t i)$を返す関数です。

1.2. 複素数ベクトル演算の関数

複素数値の配列を、多段の複素数ベクトルとして、1階複素数ベクトルと、2階複素数ベクトル(複素数ベクトルを要素としたベクトル)のためのDFT/FFTで使用する演算を関数として用意します。
実装コード内でループではなく、ベクトル演算を用いることで、1次元版と2次元版のフーリエ変換実装の比較がしやすくなります。

const v1mul = (a1d, b) => a1d.map(a => mul(a, b));
const v1add = (a1d, b1d) => a1d.map((a, i) => add(a, b1d[i]));
const v1sub = (a1d, b1d) => a1d.map((a, i) => sub(a, b1d[i]));
const v2add = (a2d, b2d) => a2d.map((a1d, i) => v1add(a1d, b2d[i]));
const v2sub = (a2d, b2d) => a2d.map((a1d, i) => v1sub(a1d, b2d[i]));

const v1sum = c1d => c1d.reduce(add, zero());
const v2sum = c2d => c2d.reduce(v1add, c2d[0].map(zero));

全体に渡って、可読性のためベクトルの次元(1d、2d)を変数名につけています。

これらは、以下の演算を行います:

v1mul: 1階複素数ベクトルの各要素ごとに、同じ複素数値をかける関数(複素数スカラーの積)
v1add: 1階複素数ベクトルの各要素ごとの和
v1sub: 1階複素数ベクトルの各要素ごとの差
v2add: 2階複素数ベクトルの各要素ごとの和
v2sub: 2階複素数ベクトルの各要素ごとの差
v1sum: 1階複素数ベクトルの要素である複素数値の総和(複素数値)
v2sum: ２階複素数ベクトルの要素である1階複素数ベクトルの総和(1階複素数ベクトル)

そして、v1で始まる1階ベクトル演算のコードと、v2で始まる2階ベクトル演算のコードは、ほぼ同じ記述になります。

注: v2sumは、全複素数値の総和ではなく、多階ベクトルの階数を1つ減らす操作です。

1.3. 行列演算の関数

2次元データのDFT/FFTでは、二次元配列データ(2階ベクトル)を行列とみなして転置操作(transpose)を利用する実装もあるので、用意しておきます。

const transpose = c2d => c2d[0].map((_, j) => c2d.map((_, i) => c2d[i][j]));

転置は、N行M列の行列を、M行N列の行列に変換する操作です。

1.4. DFT/IDFT関数のファクトリー関数

フーリエ変換とフーリエ逆変換の実装処理はほぼ同じものであり、ただ掛け合わせるexpiに渡す角度の符号と、結果の複素数値を要素数で割る後処理だけが違います。

ここでは、この角度と二次元データを渡す、フーリエ変換コア実装の関数(dft2dc(t, c2d))を受け取り、フーリエ変換と逆フーリエ変換の関数を生成して返すmakeDft2dを用意しておき、以降、各種フーリエ変換実装のコア部分だけを扱うようにします。

const makeDft2d = dft2dc => [
  f2d => dft2dc(-2 * Math.PI, f2d),
  F2d => {
    const f2d = dft2dc(2 * Math.PI, F2d);
    return f2d.map(c1d => c1d.map(c => divr(c, f2d.length * c1d.length)));
  },
];

2. ベーシックな2次元DFT実装

1次元のDFT/IDFTの定義は、$f$を変換対象データ、$F$をフーリエ変換後データとし、それぞれの結果の各要素ごとに、

\begin{align}
F(u) =& \sum_{I=0}^{M} f(I) \exp(-2\pi i\frac{u I}{M}) \\
f(s) =& \frac{1}{M} \sum_{I=0}^{M} F(I) \exp(2\pi i\frac{s I}{M}) \\
\end{align}

でした。DFT全体では、$M$ 個の和で求める $F(u)$ が $M$ 個あるので、計算量のオーダーは $M^2$ となりました。

そして、２次元DFT/IDFの場合は、結果の各要素ごとに、

\begin{align}
F(u, v) =& \sum_{J=0}^{N} \sum_{I=0}^{M} f(I, J) \exp(-2\pi i(\frac{u I}{M} + \frac{v J}{N})) \\
f(s, t) =& \frac{1}{M N} \sum_{J=0}^{N} \sum_{I=0}^{M} F(I, J) \exp(2\pi i(\frac{s I}{M} + \frac{t J}{N})) \\
\end{align}

となります。２次元DFT全体では、$M \times N$ 個の和で求める $F(u,v)$ が $M \times N$ 個あるので、計算量のオーダーは $M^2 N^2$ となります。

DFTとIDFTの仕組みは、$\exp$の角度の符号とIDFTのときに割る部分を除けば、同じ形の式になります。
そこで、角度を $T$ 、入力データを $c(u, v)$ 、出力データを $R(u, v)$ とすると、２次元DFT/IDFTの共通部分は以下のようになります:

$$
R(u, v) = \sum_{J=0}^{N} \sum_{I=0}^{M} c(I, J) \exp(T i(\frac{u I}{M} + \frac{v J}{N}))
$$

この共通部分dft2d0c(t, c2d)として、そのままJavaScriptコードにして書くと、以下ようになります:

const dft2d0c = (t, c2d) => c2d.map((c1d, v) => c1d.map(
  (_, u) => v1sum(c2d.map((c1d, j) => v1sum(c1d.map(
    (c, i) => mul(c, expi(t * (u * i / c1d.length + v * j / c2d.length)))))))));

const [dft2d0, idft2d0] = makeDft2d(dft2d0c);

2段のv1sum(map())が、定義での２つの $\sum$ に対応します。その外側の2段のmapで全要素についての結果を求めます。

2.1. 実行結果

まずDFT/IDFTの実行例で使用する、二次元複素数データ(2階複素数ベクトル)の表示と比較の関数print2dとcmp2dを用意します。

const print2d = c2d => c2d.forEach(c1d => {
  console.log("  " + c1d.reduce((s, [cx, cy]) => `${s}[${cx} + ${cy}i], `, ""));
});

const cmp2d = (title, a2d, b2d, eps = 1e-12) => {
  a2d.forEach((a1d, j) => a1d.forEach(([ax, ay], i) => {
    const [bx, by] = b2d[j][i];
    console.assert(
      Math.abs(ax - bx) < eps && Math.abs(ay - by) < eps,
      `${title}: [${j}][${i}]`);
  }));
};

例として8x8の二次元複素数データを用意します。値には、添字の整数値をもとにsinを使って生成した数値を用いました。

const v2d0 = [...Array(8).keys()].map(j => [...Array(8).keys()].map(
  i => [Math.sin(5 * i + 3 * j), 0]));

console.log("[f0]");
print2d(v2d0);

[f0]
  [0 + 0i], [-0.9589242746631385 + 0i], [-0.5440211108893698 + 0i], [0.6502878401571168 + 0i], [0.9129452507276277 + 0i], [-0.13235175009777303 + 0i], [-0.9880316240928618 + 0i], [-0.428182669496151 + 0i], 
  [0.1411200080598672 + 0i], [0.9893582466233818 + 0i], [0.4201670368266409 + 0i], [-0.7509872467716762 + 0i], [-0.8462204041751706 + 0i], [0.27090578830786904 + 0i], [0.9999118601072672 + 0i], [0.2963685787093853 + 0i], 
  [-0.27941549819892586 + 0i], [-0.9999902065507035 + 0i], [-0.2879033166650653 + 0i], [0.8366556385360561 + 0i], [0.7625584504796028 + 0i], [-0.404037645323065 + 0i], [-0.9917788534431158 + 0i], [-0.158622668804709 + 0i], 
  [0.4121184852417566 + 0i], [0.9906073556948704 + 0i], [0.14987720966295234 + 0i], [-0.9055783620066238 + 0i], [-0.6636338842129675 + 0i], [0.5290826861200238 + 0i], [0.9637953862840878 + 0i], [0.017701925105413577 + 0i], 
  [-0.5365729180004349 + 0i], [-0.9613974918795568 + 0i], [-0.008851309290403876 + 0i], [0.956375928404503 + 0i], [0.5514266812416906 + 0i], [-0.6435381333569995 + 0i], [-0.9165215479156338 + 0i], [0.123573122745224 + 0i], 
  [0.6502878401571168 + 0i], [0.9129452507276277 + 0i], [-0.13235175009777303 + 0i], [-0.9880316240928618 + 0i], [-0.428182669496151 + 0i], [0.7451131604793488 + 0i], [0.8509035245341184 + 0i], [-0.26237485370392877 + 0i], 
  [-0.7509872467716762 + 0i], [-0.8462204041751706 + 0i], [0.27090578830786904 + 0i], [0.9999118601072672 + 0i], [0.2963685787093853 + 0i], [-0.8317747426285983 + 0i], [-0.7682546613236668 + 0i], [0.39592515018183416 + 0i], 
  [0.8366556385360561 + 0i], [0.7625584504796028 + 0i], [-0.404037645323065 + 0i], [-0.9917788534431158 + 0i], [-0.158622668804709 + 0i], [0.9017883476488092 + 0i], [0.6702291758433747 + 0i], [-0.5215510020869119 + 0i],

この値v2d0に、dft2d0c実装によるDFTとIDFTをかけ、そのフーリエ変換値の表示と、逆変換結果と元の値との比較を行います。

const F2d0 = dft2d0(v2d0), f2d0 = idft2d0(F2d0);

console.log("[F]");
print2d(F2d0);
cmp2d("idft2d0", v2d0, f2d0);

[F]
  [-0.22229879303422662 + 0i], [-0.5829478362350251 + 0.49946881512451435i], [1.6158074809677128 + -1.0544318738657057i], [0.6760817853439278 + -0.21346789979905423i], [0.5900664050710727 + 2.241523206033985e-16i], [0.676081785343928 + 0.2134678997990545i], [1.6158074809677119 + 1.0544318738657061i], [-0.5829478362350273 + -0.499468815124517i], 
  [-0.21287180658860183 + -0.11538090860711825i], [-0.7584032823126821 + 0.23478295922095072i], [2.00640184828377 + -0.3839209751643051i], [0.7582086658433942 + 0.06030647478566076i], [0.5925982026948169 + 0.23974247206021915i], [0.5975583383569237 + 0.49438127886705846i], [1.2128641500717299 + 1.7602061254268897i], [-0.3825161873697558 + -0.7808584007869516i], 
  [-0.16712422884347958 + -0.2797205014555071i], [-0.9654928015319779 + -0.11358054181791366i], [2.541564819192021 + 0.5106522342275537i], [0.8813610100676346 + 0.438011330141586i], [0.6048845902158891 + 0.581213004078648i], [0.49189285016154494 + 0.906730520794004i], [0.6177737384375472 + 2.8259064480072i], [-0.05422136193782978 + -1.21028228948894i], 
  [0.10732646598272466 + -0.6921977533614565i], [-1.273805659144174 + -0.84651756443816i], [3.7802212514289124 + 2.457388998368781i], [1.2209707817448012 + 1.3255728696286773i], [0.6785935767517792 + 1.4382726098168974i], [0.25719092331458027 + 2.002132047333954i], [-0.9804017887341514 + 5.799275821635207i], [0.9812312076828965 + -2.4295200260174825i], 
  [-11.138588774903457 + -4.482527616176913e-15i], [-15.043956868557988 + -9.17184867739379i], [8.765752386031306 + 19.362749574878002i], [-1.4102755770682913 + 3.9199549904830135i], [-2.3417121856333027 + -3.789391776621495e-16i], [-1.4102755770682853 + -3.919954990483025i], [8.765752386031343 + -19.36274957487798i], [-15.04395686855801 + 9.171848677393761i], 
  [0.10732646598272844 + 0.6921977533614586i], [0.9812312076828968 + 2.4295200260174923i], [-0.9804017887341425 + -5.799275821635225i], [0.25719092331458415 + -2.0021320473339532i], [0.6785935767517806 + -1.4382726098168899i], [1.2209707817448074 + -1.3255728696286715i], [3.7802212514289257 + -2.4573889983687702i], [-1.2738056591441682 + 0.8465175644381502i], 
  [-0.16712422884348274 + 0.2797205014555082i], [-0.05422136193784044 + 1.2102822894889438i], [0.6177737384375543 + -2.8259064480072045i], [0.49189285016154427 + -0.9067305207940063i], [0.6048845902158861 + -0.5812130040786413i], [0.8813610100676356 + -0.4380113301416003i], [2.5415648191920175 + -0.5106522342275472i], [-0.9654928015319677 + 0.11358054181792498i], 
  [-0.2128718065885955 + 0.11538090860712158i], [-0.3825161873697558 + 0.7808584007869519i], [1.2128641500717323 + -1.7602061254268893i], [0.5975583383569242 + -0.4943812788670471i], [0.5925982026948178 + -0.23974247206021976i], [0.7582086658433971 + -0.060306474785652944i], [2.0064018482837875 + 0.3839209751643107i], [-0.7584032823126798 + -0.23478295922095205i],

逆変換結果の比較では警告が出ないので、うまく逆変換できていることが確認できます。

3. ２次元DFTの変形実装

次に、このベーシックな二次元DFT実装を変形した実装を考えます。

3.1. 多段DFT実装

先の2次元DFTの式

$$
R(u, v) = \sum_{J=0}^{N} \sum_{I=0}^{M} c(I, J) \exp(T i(\frac{u I}{M} + \frac{v J}{N}))
$$

は、$\exp(A+B) = \exp(A)\exp(B)$であるため、$\sum$を$I$と$J$で分離すると、

$$
R(u, v) = \sum_{J=0}^{N} \left(\sum_{I=0}^{M} c(I, J) \exp(T i(\frac{u I}{M}) \right) \exp(T i\frac{v J}{N}))
$$

となります。

この内側の $\sum$ 部分は、 $0\dots N-1$で$J$をそれぞれ固定することで、N個の1次元DFTとなっています。
この部分を

$$
D(u, J) = \sum_{I=0}^{M} c(I, J) \exp(T i(\frac{u I}{M})
$$

と置く(D全体は2階複素数ベクトル)と。

$$
R(u, v) = \sum_{J=0}^{N} D(u, J) \exp(T i\frac{v J}{N}))
$$

となります。
$0\dots M-1$で$u$をそれぞれ固定するとM個のDFTとなり、$R_u(v)$はM要素の1階ベクトル$D_u(J)$へののDFTとみなせます。

この２段DFTをdft2dc(t, c2d)として実装したものが次のコードです。

const dft1dc = (t, c1d) => c1d.map((_, u) => v1sum(
  c1d.map((c, i) => mul(c, expi(t * u * i / c1d.length)))));

const dft2dc = (t, c2d) => {
  const d2d = c2d.map(c1d => dft1dc(t, c1d));
  return d2d.map((_, v) => v2sum(
    d2d.map((c1d, j) => v1mul(c1d,  expi(t * v * j / d2d.length)))));
};

const [dft2d, idft2d] = makeDft2d(dft2dc);

dft1dcは、１次元DFTそのものです。
dft2dcでははじめに1次元DFTを各1階ベクトルにかけ、その結果のd2dに対し1階ベクトルの列に対してDFTをかけています。
後半のベクトル列へのDFTの式は、dft1dcでのv1sumとmulを、v2sumとv1mulに置き換えただけの式となります。

この2段DFTの計算量のオーダーは、前半のM要素DFTがN個と、後半のN要素DFTがM個なので、$MN(M+N)$になります。

この処理のテストコードは、以下の通りです。

{
  const F2d = dft2d(v2d0), f2d = idft2d(F2d);
  cmp2d("dft2d", v2d0, f2d);
  cmp2d("idft2d", F2d0, F2d);
}

誤差は範囲内でアサーションにはかかりませんが、多段化することでexpi１つが2つの積に分離されるので、ベーシックな実装での結果とは浮動小数値としては全く同じにはなりません。

この実装は二次元データのための二段実装ですが、多次元データ一般でのDFTにも適用可能な手法です。

3.2. 二次元DFTの行列転置での実装

二段DFTの後半の(M要素)1階ベクトルへのDFTというのは、通常の配列を横方向とすると、縦方向のM列に対して1次元DFTをかけるのと同じことをしています。

二次元データを行列とみなし、行列の転置を行うことで、1階ベクトル演算のかわりに前半と同様に1次元DFTを使って、２次元DFTを実装することができます。

この行列転置を用いた実装の二次元DFTをdft2dtc(t, c2d)として実装したものが以下のコードです。

const dft2dtc = (t, c2d) => transpose(
  transpose(c2d.map(c1d => dft1dc(t, c1d))).map(c1d => dft1dc(t, c1d)));

const [dft2dt, idft2dt] = makeDft2d(dft2dtc);

やってることは、

各行に対し1次元DFTをかける
その結果を転置する
その転置した各行に対して1次元DFTをかる
その結果を転置して返す

です。
1次元DFTでのメモリのシーケンシャルな走査だけになるので、環境によっては２段DFTより実行効率が良くなる可能性があります。

テストコードは以下で、結果は誤差の範囲内で一致します。

{
  const F2d = dft2dt(v2d0), f2d = idft2dt(F2d);
  cmp2d("dft2dt", v2d0, f2d);
  cmp2d("idft2dt", F2d0, F2d);
}

4. 再帰版２次元FFT実装

多次元版FFTは、多次元版DFT変形実装と同じ構成で実装できます。
つまり、多階ベクトル演算でのFFT実装の多段化であり、そして、行列転置での実装です。

Cooley-Tukey型FFTの制約として各段で要素数が２のべき乗である必要があります。

4.1. 多段での再帰FFT実装

1次元再帰FFTfft1drc(t, c1d)とそれを使った2次元再帰FFTfft2drc(t, c2d)は以下のようになります。

const fft1drc = (t, c1d) => {
  if (c1d.length === 1) return c1d;
  const l1d = fft1drc(t, c1d.filter((_, i) => i % 2 === 0));
  const r1d = fft1drc(t, c1d.filter((_, i) => i % 2 === 1));
  const re1d = r1d.map((c, i) => mul(c, expi(t * i / c1d.length)));
  return v1add(l1d, re1d).concat(v1sub(l1d, re1d));
};

const fft2drc = (t, c2d) => {
  if (c2d.length === 1) return [fft1drc(t, c2d[0])];
  const l2d = fft2drc(t, c2d.filter((_, i) => i % 2 === 0));
  const r2d = fft2drc(t, c2d.filter((_, i) => i % 2 === 1));
  const re2d = r2d.map((c1d, i) => v1mul(c1d, expi(t * i / c2d.length)));
  return v2add(l2d, re2d).concat(v2sub(l2d, re2d));
};

const [fft2dr, ifft2dr] = makeDft2d(fft2drc);

二次元版再帰FFTは、1次元版再帰FFTと同じ5行のコードになります。
二次元版と一次元版はやってることはほぼ同じでmulがv1mulに、v1addやv1subがv2addやv2subになる点などは、多段DFTでの場合と同じです。

注意点は、2次元FFTの中のどこで1次元FFTを呼び出すかですが、長さ1のときの戻り値に対して使用します。

計算量のオーダーは、$MN\log M\log N$です。

テストコードは以下で、結果は誤差の範囲内で一致します。

{
  const F2d = fft2dr(v2d0), f2d = ifft2dr(F2d);
  cmp2d("fft2dtr", v2d0, f2d);
  cmp2d("ifft2dr", F2d0, F2d);
}

4.2. 行列転置による2次元再帰FFT実装

行列転置を使った2次元FFT実装fft2drtc(t, c2d)も、2次元DFT実装のときと全く同じもので、ただ1次元DFT関数の違うだけです。

const fft2drtc = (t, c2d) => transpose(
  transpose(c2d.map(c1d => fft1drc(t, c1d))).map(c1d => fft1drc(t, c1d)));

const [fft2drt, ifft2drt] = makeDft2d(fft2drtc);

テストコードは以下で、結果は誤差の範囲内で一致します。

{
  const F2d = fft2drt(v2d0), f2d = ifft2drt(F2d);
  cmp2d("fft2dtrt", v2d0, f2d);
  cmp2d("ifft2drt", F2d0, F2d);
}

5. ループ版2次元FFT実装

1次元でのループ版FFTの動作の詳細については、「高速フーリエ変換FFTを理解する」にあります。
この1次元FFT実装を元に、2次元FFTのループ版を作ります。

5.1. 多段でのループFFT実装

1次元ループFFTfft1dlc(t, c1d)とそれを使った2次元ループFFTfft2dlc(t, c2d)は以下のようになります。

const revbit = (k, n0) => {
  const s1 = ((n0 & 0xaaaaaaaa) >>> 1) | ((n0 & 0x55555555) << 1);
  const s2 = ((s1 & 0xcccccccc) >>> 2) | ((s1 & 0x33333333) << 2);
  const s3 = ((s2 & 0xf0f0f0f0) >>> 4) | ((s2 & 0x0f0f0f0f) << 4);
  const s4 = ((s3 & 0xff00ff00) >>> 8) | ((s3 & 0x00ff00ff) << 8);
  const s5 = ((s4 & 0xffff0000) >>> 16) | ((s4 & 0x0000ffff) << 16);
  return s5 >>> (32 - k);
};

const fft1dlc = (t, c1d) =>  {
  const n = c1d.length, k = Math.log2(n);
  const r1d = c1d.map((_, i) => c1d[revbit(k, i)]);
  for (let nh = 1; nh < n; nh *= 2) {
    t /= 2;
    for (let s = 0; s < n; s += nh * 2) {
      for (let i = 0; i < nh; i++) {
        const l = r1d[s + i], re = mul(r1d[s + i + nh], expi(t * i));
        r1d[s + i] = add(l, re), r1d[s + i + nh] = sub(l, re);
      }
    }
  }
  return r1d;
};

const fft2dlc = (t, c2d) =>  {
  const n = c2d.length, k = Math.log2(n);
  const r2d = c2d.map((_, i) => fft1dlc(t, c2d[revbit(k, i)]));
  for (let nh = 1; nh < n; nh *= 2) {
    t /= 2;
    for (let s = 0; s < n; s += nh * 2) {
      for (let i = 0; i < nh; i++) {
        const l1d = r2d[s + i], re1d = v1mul(r2d[s + i + nh], expi(t * i));
        r2d[s + i] = v1add(l1d, re1d), r2d[s + i + nh] = v1sub(l1d, re1d);
      }
    }
  }
  return r2d;
};

const [fft2dl, ifft2dl] = makeDft2d(fft2dlc);

ループ版２次元FFTでも、1次元FFTと2次元FFTの関係は、再帰版と同様の関係です。

ビット反転によるインデックスの入れ替え方も一緒ですが、2次元FFTでは、このビット反転インデックスへの入れ替え時が、1次元FFTを適用する位置になります。

注意点としては、C言語などでの場合、ベクトルl1d、re1d相当をローカル変数としてスタックメモリに置くと、要素数次第でスタックリミットに達する可能性があるので、ローカル変数は使わずに、リザルトメモリ上で掛け算や足し算を行うようにする必要があるでしょう。

テストコードは以下で、結果は誤差の範囲内で一致します。

{
  const F2d = fft2dl(v2d0), f2d = ifft2dl(F2d);
  cmp2d("fft2dtl", v2d0, f2d);
  cmp2d("ifft2dl", F2d0, F2d);
}

5.2. 行列転置による2次元ループFFT実装

行列転置を使った2次元FFT実装fft2dltc(t, c2d)も、2次元DFT実装のときと全く同じもので、ただ1次元DFT関数の違うだけです。

const fft2dltc = (t, c2d) => transpose(
  transpose(c2d.map(c1d => fft1dlc(t, c1d))).map(c1d => fft1dlc(t, c1d)));

const [fft2dlt, ifft2dlt] = makeDft2d(fft2dltc);

(正方行列でないと、インプレースで転置するのはメモリアクセス的に非効率で、転置行列のワーキングメモリを用意するべきかも)

テストコードは以下で、結果は誤差の範囲内で一致します。

{
  const F2d = fft2dlt(v2d0), f2d = ifft2dlt(F2d);
  cmp2d("fft2dtlt", v2d0, f2d);
  cmp2d("ifft2dlt", F2d0, F2d);
}

6. まとめ

要素が複素数のベクトルデータへの１次元版のDFT/FFTを理解していれば、要素が複素数ベクトルな2階のベクトルとみなすことで、2次元でのDFT/FFTを作ることができました。

同様に、要素がn階複素数ベクトルなn+1階ベクトルとして扱うことで、n次元版DFT/FFTを用いてn+1次元版のDFT/FFTを実装できます。3次元版DFT/FFT実装については、

https://gist.github.com/bellbind/b7583ca138a6a1a14763756cd75ccc86

のfft3d.jsにおいておきました。