離散フーリエ変換の高速化（高速フーリエ変換）をソースの変形だけで導いてみます。 「あまり得意ではない」といっても、DFTの性質や、複素数の極形式の演算とオイラーの公式くらいはわかっている必要があります。 最終的には実離散高速フーリエ変換を目指しますが、書いていたら思いのほか大きくなってしまったので、今回は \(N^2\) 回の計算をする複素DFT関数から、ループを用いた複素FFT関数を導くところまでです。
以下のサイトが色々と参考になります。
つか、まぁ、「FFTの概略と設計法」に書いてある「機械的な変形」をていねいに説明しただけなんだが・・・
function formatNumeric(x) {
return x.toFixed(6).padStart(11);
}
class Complex {
constructor(re = 0, im = 0) {
this.re = re;
this.im = im;
}
toString() {
return `${formatNumeric(this.re)} ${formatNumeric(this.im)}`;
}
}
const src = Array(16).fill(0).map(_ => new Complex(Math.random() - 0.5));
src.map(z => LOG(z));
STATUS("finished");
DFTは複素数の方が楽なので実数化は後で考えることにして、とりあえずComplexクラスを作ってテストデータを入れる。
LOGとかSTATUSというのはこのスクリプトを実行している外側のHTMLファイルで定義されていて、なんかログ文字列を出力する関数だと思えばいい。
興味がある人はjs/fft.htmを見ればすぐに分かる。
ログ出力の時にconsole.logと同じように文字列化されるので、Complex.toStringも定義しておく。
function formatNumeric(x) {
return x.toFixed(6).padStart(11);
}
class Complex {
constructor(re = 0, im = 0) {
this.re = re;
this.im = im;
}
toString() {
return `${formatNumeric(this.re)} ${formatNumeric(this.im)}`;
}
add(z) {
this.re += z.re;
this.im += z.im;
return this;
}
mul(z) {
const re = this.re;
this.re = re * z.re - this.im * z.im;
this.im = re * z.im + this.im * z.re;
return this;
}
}
function W(x) {
return new Complex(
Math.cos(2 * Math.PI * x),
Math.sin(2 * Math.PI * x)
);
}
async function DFT(src) {
const n = src.length;
const dst = [];
for (let k = 0; k < n; k++) {
await STATUS(k);
const sum = new Complex();
for (let i = 0; i < n; i++)
sum.add(W(-i * k / n).mul(src[i]));
dst.push(sum);
}
return dst;
}
const src = Array(16).fill(0).map(_ => new Complex(Math.random() - 0.5));
DFT(src).then(dft => {
dft.forEach(z => LOG(z));
STATUS("finished");
});
複素数の演算が増えている。
Wというのはいわゆる回転因子（twiddle factor）というやつで、オイラーの公式の引数の
\(2\pi\)
を省略した感じになっている。
つまり、
Wをシンプルに定義して、使う側で符号を付けることにする。
回転因子は（複素）指数関数の一種なので、指数法則が使えることに注意。
DFT関数は定義通りに書いただけで、引数srcにComplexの配列を与えると、複素DFTして結果をComplexの配列として返す。
ここでのポイントは、実DFTなので、
ことを確認することである。 以後、同じように石橋を叩きがならコードを確認していく。
一番よく使われる基数2のCooley-Tukey型の変形をするつもりなので、要素数は \(N=2^n\) 個にしておく必要がある。 16点くらいがご機嫌である。 多ければ結果の確認が面倒になる。 少ないと基数を大きくしたり、FFTを流用して高速DCTなんかを計算するときに、書き換えごとにどんどん要素数が半分になっていって、最後 \(N=1\) になってしまいました、ということが起きる。
DFT関数がasyncになっているが、これはDFTの演算が長いとJavaScriptのUIスレッドがブロックしてしまい途中経過が表示されなくなるからで、STATUSでsetTimeoutを実行して処理をJSエンジン側に譲り渡すようにしてある。
トップレベルではawaitは使えないからasync関数も使えないと勘違いしている人がいるが、実際にはasync/awaitはプロミスのシンタックスシュガーで、async関数はプロミスを返しているだけなので、awaitは使えなくても普通にPromise.thenが使える。
逆に、awaitにプロミスを指定してもちゃんとプロミスを待って普通に動く。
念の為に言っておくと、本当に普通のプロミスなので、thenに指定する関数をasyncにする必要はない。
というか、アロー関数にasync付けるのが面倒だったので、最後のSTATUSにawaitを付けてなかったりする。
まぁfinishedの文字通り、一番最後の呼び出しになるので。
function formatNumeric(x) {
return x.toFixed(6).padStart(11);
}
class Complex {
constructor(re = 0, im = 0) {
this.re = re;
this.im = im;
}
toString() {
return `${formatNumeric(this.re)} ${formatNumeric(this.im)}`;
}
add(z) {
this.re += z.re;
this.im += z.im;
return this;
}
mul(z) {
const re = this.re;
this.re = re * z.re - this.im * z.im;
this.im = re * z.im + this.im * z.re;
return this;
}
scale(r) {
this.re /= r;
this.im /= r;
return this;
}
}
function W(x) {
return new Complex(
Math.cos(2 * Math.PI * x),
Math.sin(2 * Math.PI * x)
);
}
async function DFT(src, inv = -1) {
const n = src.length;
const dst = [];
for (let k = 0; k < n; k++) {
await STATUS(inv > 0 ? "inv" : "fwd", k);
const sum = new Complex();
for (let i = 0; i < n; i++)
sum.add(W(inv * i * k / n).mul(src[i]));
dst.push(sum);
}
return dst;
}
const src = Array(16).fill(0).map(_ => new Complex(Math.random() - 0.5));
DFT(src)
.then(dft => DFT(dft, 1))
.then(idft => {
idft.map(z => z.scale(src.length))
.forEach((z, i) => LOG(`${src[i]} -> ${z}`));
STATUS("finished");
});
結果が合っているかどうか確認するため、逆DFTする。 ここでのポイントはもちろん値がちゃんと元に戻っているか確認することである。
DFTの行きと帰りの違いは原理的には回転因子Wの回転方向だけなので、invという引数を設けて、正変換では-1、逆変換では+1を指定する。
実際には結果は
\(N\)
倍になっているので、これを元に戻すためにComplex.scaleという関数を作ってある。
プログラム上のポイントはasync関数DFTの結果をさらにDFTに渡すので、thenを数珠繋ぎにしてるのと、ComplexにtoStringが定義してあるので、テンプレート文字列の中にComplexがそのまま書けることくらい。
function formatNumeric(x) {
return x.toFixed(6).padStart(11);
}
class Complex {
constructor(re = 0, im = 0) {
this.re = re;
this.im = im;
}
toString() {
return `${formatNumeric(this.re)} ${formatNumeric(this.im)}`;
}
add(z) {
this.re += z.re;
this.im += z.im;
return this;
}
sub(z) {
this.re -= z.re;
this.im -= z.im;
return this;
}
mul(z) {
const re = this.re;
this.re = re * z.re - this.im * z.im;
this.im = re * z.im + this.im * z.re;
return this;
}
scale(r) {
this.re /= r;
this.im /= r;
return this;
}
abs() {
return Math.hypot(this.re, this.im);
}
}
function W(x) {
return new Complex(
Math.cos(2 * Math.PI * x),
Math.sin(2 * Math.PI * x)
);
}
async function DFT(src, inv = -1) {
const n = src.length;
const dst = [];
for (let k = 0; k < n; k++) {
await STATUS(inv > 0 ? "inv" : "fwd", k);
const sum = new Complex();
for (let i = 0; i < n; i++)
sum.add(W(inv * i * k / n).mul(src[i]));
dst.push(sum);
}
return dst;
}
const src = Array(16).fill(0).map(_ => new Complex(Math.random() - 0.5));
DFT(src)
.then(dft => DFT(dft, 1))
.then(idft => {
idft.map(z => z.scale(src.length))
.forEach((z, i) => LOG(z.sub(src[i]).abs().toFixed(8)));
STATUS("finished");
});
せっかくスクリプトを使っているのだから、結果の誤差もスクリプトに計算してもらう。
入力と逆DFTの差をとって、複素数だから絶対値を表示すればいいだろう。
絶対値の計算は真面目に実部虚部を2乗して足して平方根を計算してもいいが、実はMath.hypotという便利な関数がある。
マイナーだけど。
プログラム上のポイントは、差の計算をする時にsrc[i].sub(z)ではなくz.sub(src[i])と書いているところで、sub関数はthisオブジェクトの値を更新してしまうから、こう書いておかないと後でsrcを使い回す時に困ってしまう。
function formatNumeric(x) {
return x.toFixed(6).padStart(11);
}
class Complex {
constructor(re = 0, im = 0) {
this.re = re;
this.im = im;
}
toString() {
return `${formatNumeric(this.re)} ${formatNumeric(this.im)}`;
}
add(z) {
this.re += z.re;
this.im += z.im;
return this;
}
sub(z) {
this.re -= z.re;
this.im -= z.im;
return this;
}
mul(z) {
const re = this.re;
this.re = re * z.re - this.im * z.im;
this.im = re * z.im + this.im * z.re;
return this;
}
scale(r) {
this.re /= r;
this.im /= r;
return this;
}
abs() {
return Math.hypot(this.re, this.im);
}
}
function W(x) {
return new Complex(
Math.cos(2 * Math.PI * x),
Math.sin(2 * Math.PI * x)
);
}
async function DFT(src, inv = -1) {
const n = src.length;
const m = n / 2;
const dst = [];
for (let k = 0; k < n; k++) {
await STATUS(inv > 0 ? "inv" : "fwd", k);
const sum = new Complex();
for (let i = 0; i < m; i++) {
sum.add(W(inv * (i * 2) * k / n).mul(src[i * 2]));
sum.add(W(inv * (i * 2 + 1) * k / n).mul(src[i * 2 + 1]));
}
dst.push(sum);
}
return dst;
}
const src = Array(16).fill(0).map(_ => new Complex(Math.random() - 0.5));
DFT(src)
.then(dft => DFT(dft, 1))
.then(idft => {
idft.map(z => z.scale(src.length))
.forEach((z, i) => LOG(z.sub(src[i]).abs().toFixed(8)));
STATUS("finished");
});
見ての通り分けただけ。 Cooley-Tukey型FFTの分解には時間間引きと周波数間引きがあるが、最終的には実FFTを計算したいので、時間間引きで計算する。
function formatNumeric(x) {
return x.toFixed(6).padStart(11);
}
class Complex {
constructor(re = 0, im = 0) {
this.re = re;
this.im = im;
}
toString() {
return `${formatNumeric(this.re)} ${formatNumeric(this.im)}`;
}
add(z) {
this.re += z.re;
this.im += z.im;
return this;
}
sub(z) {
this.re -= z.re;
this.im -= z.im;
return this;
}
mul(z) {
const re = this.re;
this.re = re * z.re - this.im * z.im;
this.im = re * z.im + this.im * z.re;
return this;
}
scale(r) {
this.re /= r;
this.im /= r;
return this;
}
abs() {
return Math.hypot(this.re, this.im);
}
}
function W(x) {
return new Complex(
Math.cos(2 * Math.PI * x),
Math.sin(2 * Math.PI * x)
);
}
async function DFT(src, inv = -1) {
const n = src.length;
const m = n / 2;
const dst = [];
for (let k = 0; k < n; k++) {
await STATUS(inv > 0 ? "inv" : "fwd", k);
const sum = new Complex();
for (let i = 0; i < m; i++) {
sum.add(W(inv * i * 2 * k / n).mul(src[i * 2]));
sum.add(W(inv * i * 2 * k / n).mul(W(inv * k / n)).mul(src[i * 2 + 1]));
}
dst.push(sum);
}
return dst;
}
const src = Array(16).fill(0).map(_ => new Complex(Math.random() - 0.5));
DFT(src)
.then(dft => DFT(dft, 1))
.then(idft => {
idft.map(z => z.scale(src.length))
.forEach((z, i) => LOG(z.sub(src[i]).abs().toFixed(8)));
STATUS("finished");
});
Wは指数関数の一種なので指数法則が使える。
これを利用して展開する。 この段階ではとりあえず展開しただけである。
function formatNumeric(x) {
return x.toFixed(6).padStart(11);
}
class Complex {
constructor(re = 0, im = 0) {
this.re = re;
this.im = im;
}
toString() {
return `${formatNumeric(this.re)} ${formatNumeric(this.im)}`;
}
add(z) {
this.re += z.re;
this.im += z.im;
return this;
}
sub(z) {
this.re -= z.re;
this.im -= z.im;
return this;
}
mul(z) {
const re = this.re;
this.re = re * z.re - this.im * z.im;
this.im = re * z.im + this.im * z.re;
return this;
}
scale(r) {
this.re /= r;
this.im /= r;
return this;
}
abs() {
return Math.hypot(this.re, this.im);
}
}
function W(x) {
return new Complex(
Math.cos(2 * Math.PI * x),
Math.sin(2 * Math.PI * x)
);
}
async function DFT(src, inv = -1) {
const n = src.length;
const m = n / 2;
const dst = [];
for (let k = 0; k < n; k++) {
await STATUS(inv > 0 ? "inv" : "fwd", k);
const sum = new Complex();
for (let i = 0; i < m; i++) {
sum.add(W(inv * i * k / m).mul(src[i * 2]));
sum.add(W(inv * i * k / m).mul(W(inv * k / n)).mul(src[i * 2 + 1]));
}
dst.push(sum);
}
return dst;
}
const src = Array(16).fill(0).map(_ => new Complex(Math.random() - 0.5));
DFT(src)
.then(dft => DFT(dft, 1))
.then(idft => {
idft.map(z => z.scale(src.length))
.forEach((z, i) => LOG(z.sub(src[i]).abs().toFixed(8)));
STATUS("finished");
});
m=n/2、つまりn=m*2を使って書き換える。
2*k/nは2*k/(m*2)なので、2が約分できてk/mになる。
function formatNumeric(x) {
return x.toFixed(6).padStart(11);
}
class Complex {
constructor(re = 0, im = 0) {
this.re = re;
this.im = im;
}
toString() {
return `${formatNumeric(this.re)} ${formatNumeric(this.im)}`;
}
add(z) {
this.re += z.re;
this.im += z.im;
return this;
}
sub(z) {
this.re -= z.re;
this.im -= z.im;
return this;
}
mul(z) {
const re = this.re;
this.re = re * z.re - this.im * z.im;
this.im = re * z.im + this.im * z.re;
return this;
}
scale(r) {
this.re /= r;
this.im /= r;
return this;
}
abs() {
return Math.hypot(this.re, this.im);
}
}
function W(x) {
return new Complex(
Math.cos(2 * Math.PI * x),
Math.sin(2 * Math.PI * x)
);
}
async function DFT(src, inv = -1) {
const n = src.length;
const m = n / 2;
const dst = [];
for (let k = 0; k < n; k++) {
await STATUS(inv > 0 ? "inv" : "fwd", k);
const even = new Complex(), odd = new Complex();
for (let i = 0; i < m; i++) {
even.add(W(inv * i * k / m).mul(src[i * 2]));
odd.add(W(inv * i * k / m).mul(src[i * 2 + 1]));
}
dst.push(even.add(odd.mul(W(inv * k / n))));
}
return dst;
}
const src = Array(16).fill(0).map(_ => new Complex(Math.random() - 0.5));
DFT(src)
.then(dft => DFT(dft, 1))
.then(idft => {
idft.map(z => z.scale(src.length))
.forEach((z, i) => LOG(z.sub(src[i]).abs().toFixed(8)));
STATUS("finished");
});
W(inv*k/n)の項にはiが含まれていないから、ループ内では定数である。
つまり、ループの外に出せる。
今は一種の総和を計算しているので、乗算してから総和を取っていたものを、総和を取ってから乗算するように書き換える。
乗算が必要なのは奇数項だけなので、偶奇別々に和をとって、奇数項に回転因子を乗じてから両者を足す。
function formatNumeric(x) {
return x.toFixed(6).padStart(11);
}
class Complex {
constructor(re = 0, im = 0) {
this.re = re;
this.im = im;
}
toString() {
return `${formatNumeric(this.re)} ${formatNumeric(this.im)}`;
}
add(z) {
this.re += z.re;
this.im += z.im;
return this;
}
sub(z) {
this.re -= z.re;
this.im -= z.im;
return this;
}
mul(z) {
const re = this.re;
this.re = re * z.re - this.im * z.im;
this.im = re * z.im + this.im * z.re;
return this;
}
scale(r) {
this.re /= r;
this.im /= r;
return this;
}
abs() {
return Math.hypot(this.re, this.im);
}
}
function W(x) {
return new Complex(
Math.cos(2 * Math.PI * x),
Math.sin(2 * Math.PI * x)
);
}
async function DFT(src, inv = -1) {
const n = src.length;
const m = n / 2;
const dst = Array(n);
for (let k = 0; k < m; k++) {
await STATUS(inv > 0 ? "inv" : "fwd", k);
let even = new Complex(), odd = new Complex();
for (let i = 0; i < m; i++) {
even.add(W(inv * i * k / m).mul(src[i * 2]));
odd.add(W(inv * i * k / m).mul(src[i * 2 + 1]));
}
dst[k] = even.add(odd.mul(W(inv * k / n)));
even = new Complex();
odd = new Complex();
for (let i = 0; i < m; i++) {
even.add(W(inv * i * (k + m) / m).mul(src[i * 2]));
odd.add(W(inv * i * (k + m) / m).mul(src[i * 2 + 1]));
}
dst[k + m] = even.add(odd.mul(W(inv * (k + m) / n)));
}
return dst;
}
const src = Array(16).fill(0).map(_ => new Complex(Math.random() - 0.5));
DFT(src)
.then(dft => DFT(dft, 1))
.then(idft => {
idft.map(z => z.scale(src.length))
.forEach((z, i) => LOG(z.sub(src[i]).abs().toFixed(8)));
STATUS("finished");
});
最終的には長さが半分のDFTふたつに分けたいので、kのループ回数も半分にする必要がある。
ループ1回ごとにiのループを2回実行するようにして、kのループ回数を半分にする。
今は時間間引きなので、周波数側であるkのループはkとk+n/2、つまり前後に分けて計算する。
先ほどと同様に回転因子を展開して、
function formatNumeric(x) {
return x.toFixed(6).padStart(11);
}
class Complex {
constructor(re = 0, im = 0) {
this.re = re;
this.im = im;
}
toString() {
return `${formatNumeric(this.re)} ${formatNumeric(this.im)}`;
}
add(z) {
this.re += z.re;
this.im += z.im;
return this;
}
sub(z) {
this.re -= z.re;
this.im -= z.im;
return this;
}
mul(z) {
const re = this.re;
this.re = re * z.re - this.im * z.im;
this.im = re * z.im + this.im * z.re;
return this;
}
scale(r) {
this.re /= r;
this.im /= r;
return this;
}
abs() {
return Math.hypot(this.re, this.im);
}
}
function W(x) {
return new Complex(
Math.cos(2 * Math.PI * x),
Math.sin(2 * Math.PI * x)
);
}
async function DFT(src, inv = -1) {
const n = src.length;
const m = n / 2;
const dst = Array(n);
for (let k = 0; k < m; k++) {
await STATUS(inv > 0 ? "inv" : "fwd", k);
let even = new Complex(), odd = new Complex();
for (let i = 0; i < m; i++) {
even.add(W(inv * i * k / m).mul(src[i * 2]));
odd.add(W(inv * i * k / m).mul(src[i * 2 + 1]));
}
dst[k] = even.add(odd.mul(W(inv * k / n)));
even = new Complex();
odd = new Complex();
for (let i = 0; i < m; i++) {
even.add(W(inv * i * k / m).mul(W(inv * i * m / m)).mul(src[i * 2]));
odd.add(W(inv * i * k / m).mul(W(inv * i * m / m)).mul(src[i * 2 + 1]));
}
dst[k + m] = even.add(odd.mul(W(inv * k / n).mul(W(inv * m / n))));
}
return dst;
}
const src = Array(16).fill(0).map(_ => new Complex(Math.random() - 0.5));
DFT(src)
.then(dft => DFT(dft, 1))
.then(idft => {
idft.map(z => z.scale(src.length))
.forEach((z, i) => LOG(z.sub(src[i]).abs().toFixed(8)));
STATUS("finished");
});
約分する。
function formatNumeric(x) {
return x.toFixed(6).padStart(11);
}
class Complex {
constructor(re = 0, im = 0) {
this.re = re;
this.im = im;
}
toString() {
return `${formatNumeric(this.re)} ${formatNumeric(this.im)}`;
}
add(z) {
this.re += z.re;
this.im += z.im;
return this;
}
sub(z) {
this.re -= z.re;
this.im -= z.im;
return this;
}
mul(z) {
const re = this.re;
this.re = re * z.re - this.im * z.im;
this.im = re * z.im + this.im * z.re;
return this;
}
scale(r) {
this.re /= r;
this.im /= r;
return this;
}
abs() {
return Math.hypot(this.re, this.im);
}
}
function W(x) {
return new Complex(
Math.cos(2 * Math.PI * x),
Math.sin(2 * Math.PI * x)
);
}
async function DFT(src, inv = -1) {
const n = src.length;
const m = n / 2;
const dst = Array(n);
for (let k = 0; k < m; k++) {
await STATUS(inv > 0 ? "inv" : "fwd", k);
let even = new Complex(), odd = new Complex();
for (let i = 0; i < m; i++) {
even.add(W(inv * i * k / m).mul(src[i * 2]));
odd.add(W(inv * i * k / m).mul(src[i * 2 + 1]));
}
dst[k] = even.add(odd.mul(W(inv * k / n)));
even = new Complex();
odd = new Complex();
for (let i = 0; i < m; i++) {
even.add(W(inv * i * k / m).mul(W(inv * i)).mul(src[i * 2]));
odd.add(W(inv * i * k / m).mul(W(inv * i)).mul(src[i * 2 + 1]));
}
dst[k + m] = even.add(odd.mul(W(inv * k / n).mul(W(inv / 2))));
}
return dst;
}
const src = Array(16).fill(0).map(_ => new Complex(Math.random() - 0.5));
DFT(src)
.then(dft => DFT(dft, 1))
.then(idft => {
idft.map(z => z.scale(src.length))
.forEach((z, i) => LOG(z.sub(src[i]).abs().toFixed(8)));
STATUS("finished");
});
W(inv*i)という項が出てくるが、回転因子のパラメータが整数ならば、その回転因子の値は1なので、乗算が丸ごと省略できる。
function formatNumeric(x) {
return x.toFixed(6).padStart(11);
}
class Complex {
constructor(re = 0, im = 0) {
this.re = re;
this.im = im;
}
toString() {
return `${formatNumeric(this.re)} ${formatNumeric(this.im)}`;
}
add(z) {
this.re += z.re;
this.im += z.im;
return this;
}
sub(z) {
this.re -= z.re;
this.im -= z.im;
return this;
}
mul(z) {
const re = this.re;
this.re = re * z.re - this.im * z.im;
this.im = re * z.im + this.im * z.re;
return this;
}
scale(r) {
this.re /= r;
this.im /= r;
return this;
}
abs() {
return Math.hypot(this.re, this.im);
}
}
function W(x) {
return new Complex(
Math.cos(2 * Math.PI * x),
Math.sin(2 * Math.PI * x)
);
}
async function DFT(src, inv = -1) {
const n = src.length;
const m = n / 2;
const dst = Array(n);
for (let k = 0; k < m; k++) {
await STATUS(inv > 0 ? "inv" : "fwd", k);
let even = new Complex(), odd = new Complex();
for (let i = 0; i < m; i++) {
even.add(W(inv * i * k / m).mul(src[i * 2]));
odd.add(W(inv * i * k / m).mul(src[i * 2 + 1]));
}
dst[k] = even.add(odd.mul(W(inv * k / n)));
even = new Complex();
odd = new Complex();
for (let i = 0; i < m; i++) {
even.add(W(inv * i * k / m).mul(src[i * 2]));
odd.add(W(inv * i * k / m).mul(src[i * 2 + 1]));
}
dst[k + m] = even.add(odd.mul(W(inv * k / n).mul(W(inv / 2))));
}
return dst;
}
const src = Array(16).fill(0).map(_ => new Complex(Math.random() - 0.5));
DFT(src)
.then(dft => DFT(dft, 1))
.then(idft => {
idft.map(z => z.scale(src.length))
.forEach((z, i) => LOG(z.sub(src[i]).abs().toFixed(8)));
STATUS("finished");
});
さらに、W(inv/2)という項は正変換か逆変換かに応じてW(-0.5)かW(0.5)になるが、
\(\cos(\mp \pi)+j\sin(\mp \pi)\)
のことなので、どっちにしろ
\(-1\)
である。
したがって、乗算して加算していたものが、単なる減算になる。
function formatNumeric(x) {
return x.toFixed(6).padStart(11);
}
class Complex {
constructor(re = 0, im = 0) {
this.re = re;
this.im = im;
}
toString() {
return `${formatNumeric(this.re)} ${formatNumeric(this.im)}`;
}
add(z) {
this.re += z.re;
this.im += z.im;
return this;
}
sub(z) {
this.re -= z.re;
this.im -= z.im;
return this;
}
mul(z) {
const re = this.re;
this.re = re * z.re - this.im * z.im;
this.im = re * z.im + this.im * z.re;
return this;
}
scale(r) {
this.re /= r;
this.im /= r;
return this;
}
abs() {
return Math.hypot(this.re, this.im);
}
}
function W(x) {
return new Complex(
Math.cos(2 * Math.PI * x),
Math.sin(2 * Math.PI * x)
);
}
async function DFT(src, inv = -1) {
const n = src.length;
const m = n / 2;
const dst = Array(n);
for (let k = 0; k < m; k++) {
await STATUS(inv > 0 ? "inv" : "fwd", k);
let even = new Complex(), odd = new Complex();
for (let i = 0; i < m; i++) {
even.add(W(inv * i * k / m).mul(src[i * 2]));
odd.add(W(inv * i * k / m).mul(src[i * 2 + 1]));
}
dst[k] = even.add(odd.mul(W(inv * k / n)));
even = new Complex();
odd = new Complex();
for (let i = 0; i < m; i++) {
even.add(W(inv * i * k / m).mul(src[i * 2]));
odd.add(W(inv * i * k / m).mul(src[i * 2 + 1]));
}
dst[k + m] = even.sub(odd.mul(W(inv * k / n)));
}
return dst;
}
const src = Array(16).fill(0).map(_ => new Complex(Math.random() - 0.5));
DFT(src)
.then(dft => DFT(dft, 1))
.then(idft => {
idft.map(z => z.scale(src.length))
.forEach((z, i) => LOG(z.sub(src[i]).abs().toFixed(8)));
STATUS("finished");
});
function formatNumeric(x) {
return x.toFixed(6).padStart(11);
}
class Complex {
constructor(re = 0, im = 0) {
if (re instanceof Complex) {
this.re = re.re;
this.im = re.im;
return;
}
this.re = re;
this.im = im;
}
toString() {
return `${formatNumeric(this.re)} ${formatNumeric(this.im)}`;
}
add(z) {
this.re += z.re;
this.im += z.im;
return this;
}
sub(z) {
this.re -= z.re;
this.im -= z.im;
return this;
}
mul(z) {
const re = this.re;
this.re = re * z.re - this.im * z.im;
this.im = re * z.im + this.im * z.re;
return this;
}
scale(r) {
this.re /= r;
this.im /= r;
return this;
}
abs() {
return Math.hypot(this.re, this.im);
}
}
function W(x) {
return new Complex(
Math.cos(2 * Math.PI * x),
Math.sin(2 * Math.PI * x)
);
}
async function DFT(src, inv = -1) {
const n = src.length;
const m = n / 2;
const dst = Array(n);
for (let k = 0; k < m; k++) {
await STATUS(inv > 0 ? "inv" : "fwd", k);
let even = new Complex(), odd = new Complex();
for (let i = 0; i < m; i++) {
even.add(W(inv * i * k / m).mul(src[i * 2]));
odd.add(W(inv * i * k / m).mul(src[i * 2 + 1]));
}
odd.mul(W(inv * k / n));
dst[k] = new Complex(even).add(odd);
dst[k + m] = even.sub(odd);
}
return dst;
}
const src = Array(16).fill(0).map(_ => new Complex(Math.random() - 0.5));
DFT(src)
.then(dft => DFT(dft, 1))
.then(idft => {
idft.map(z => z.scale(src.length))
.forEach((z, i) => LOG(z.sub(src[i]).abs().toFixed(8)));
STATUS("finished");
});
ここでよく見るとふたつのiループは全く同じになっているので、片方を丸ごと削除できる。
また、oddに乗じられている回転因子も全く同じなので、この乗算も共通にできる。
JavaScriptの仕様から、代入はオブジェクトの複製ではなく共有なので、evenかoddのどちらかは複製する必要がある点に注意。
function formatNumeric(x) {
return x.toFixed(6).padStart(11);
}
class Complex {
constructor(re = 0, im = 0) {
if (re instanceof Complex) {
this.re = re.re;
this.im = re.im;
return;
}
this.re = re;
this.im = im;
}
toString() {
return `${formatNumeric(this.re)} ${formatNumeric(this.im)}`;
}
add(z) {
this.re += z.re;
this.im += z.im;
return this;
}
sub(z) {
this.re -= z.re;
this.im -= z.im;
return this;
}
mul(z) {
const re = this.re;
this.re = re * z.re - this.im * z.im;
this.im = re * z.im + this.im * z.re;
return this;
}
scale(r) {
this.re /= r;
this.im /= r;
return this;
}
abs() {
return Math.hypot(this.re, this.im);
}
}
function W(x) {
return new Complex(
Math.cos(2 * Math.PI * x),
Math.sin(2 * Math.PI * x)
);
}
async function DFT(src, inv = -1) {
const n = src.length;
const m = n / 2;
const dst = Array(n);
const tmp = [];
for (let k = 0; k < m; k++) {
await STATUS(inv > 0 ? "inv" : "fwd", k);
let even = new Complex(), odd = new Complex();
for (let i = 0; i < m; i++) {
even.add(W(inv * i * k / m).mul(src[i * 2]));
odd.add(W(inv * i * k / m).mul(src[i * 2 + 1]));
}
tmp.push(even, odd);
}
for (let k = 0; k < m; k++) {
const even = tmp[k * 2];
const odd = tmp[k * 2 + 1].mul(W(inv * k / n));
dst[k] = new Complex(even).add(odd);
dst[k + m] = even.sub(odd);
}
return dst;
}
const src = Array(16).fill(0).map(_ => new Complex(Math.random() - 0.5));
DFT(src)
.then(dft => DFT(dft, 1))
.then(idft => {
idft.map(z => z.scale(src.length))
.forEach((z, i) => LOG(z.sub(src[i]).abs().toFixed(8)));
STATUS("finished");
});
DFTっぽい部分とそうでない部分に分けておく。
DFTっぽくない部分を便宜上「後処理」と呼ぶことにする。
evenとoddは一度配列tmpにとっておき、後処理部分でそれを読み出して処理する。
その名の通り、evenはtmpの添字が偶数、oddは奇数になる。
function formatNumeric(x) {
return x.toFixed(6).padStart(11);
}
class Complex {
constructor(re = 0, im = 0) {
if (re instanceof Complex) {
this.re = re.re;
this.im = re.im;
return;
}
this.re = re;
this.im = im;
}
toString() {
return `${formatNumeric(this.re)} ${formatNumeric(this.im)}`;
}
add(z) {
this.re += z.re;
this.im += z.im;
return this;
}
sub(z) {
this.re -= z.re;
this.im -= z.im;
return this;
}
mul(z) {
const re = this.re;
this.re = re * z.re - this.im * z.im;
this.im = re * z.im + this.im * z.re;
return this;
}
scale(r) {
this.re /= r;
this.im /= r;
return this;
}
abs() {
return Math.hypot(this.re, this.im);
}
}
function W(x) {
return new Complex(
Math.cos(2 * Math.PI * x),
Math.sin(2 * Math.PI * x)
);
}
async function DFT(src, inv = -1) {
const n = src.length;
const m = n / 2;
const dst = Array(n);
const tmp = [];
for (let k = 0; k < m; k++) {
await STATUS(inv > 0 ? "inv" : "fwd", k);
let even = new Complex(), odd = new Complex();
for (let i = 0; i < m; i++)
even.add(W(inv * i * k / m).mul(src[i * 2]));
for (let i = 0; i < m; i++)
odd.add(W(inv * i * k / m).mul(src[i * 2 + 1]));
tmp.push(even, odd);
}
for (let k = 0; k < m; k++) {
const even = tmp[k * 2];
const odd = tmp[k * 2 + 1].mul(W(inv * k / n));
dst[k] = new Complex(even).add(odd);
dst[k + m] = even.sub(odd);
}
return dst;
}
const src = Array(16).fill(0).map(_ => new Complex(Math.random() - 0.5));
DFT(src)
.then(dft => DFT(dft, 1))
.then(idft => {
idft.map(z => z.scale(src.length))
.forEach((z, i) => LOG(z.sub(src[i]).abs().toFixed(8)));
STATUS("finished");
});
なんとなく見えてきただろうか。
function formatNumeric(x) {
return x.toFixed(6).padStart(11);
}
class Complex {
constructor(re = 0, im = 0) {
if (re instanceof Complex) {
this.re = re.re;
this.im = re.im;
return;
}
this.re = re;
this.im = im;
}
toString() {
return `${formatNumeric(this.re)} ${formatNumeric(this.im)}`;
}
add(z) {
this.re += z.re;
this.im += z.im;
return this;
}
sub(z) {
this.re -= z.re;
this.im -= z.im;
return this;
}
mul(z) {
const re = this.re;
this.re = re * z.re - this.im * z.im;
this.im = re * z.im + this.im * z.re;
return this;
}
scale(r) {
this.re /= r;
this.im /= r;
return this;
}
abs() {
return Math.hypot(this.re, this.im);
}
}
function W(x) {
return new Complex(
Math.cos(2 * Math.PI * x),
Math.sin(2 * Math.PI * x)
);
}
async function DFT(src, inv = -1) {
const n = src.length;
const m = n / 2;
const dst = Array(n);
const tmp = [];
for (let k = 0; k < m; k++) {
await STATUS(inv > 0 ? "inv" : "fwd", k);
let even = new Complex();
for (let i = 0; i < m; i++)
even.add(W(inv * i * k / m).mul(src[i * 2]));
tmp.push(even);
}
for (let k = 0; k < m; k++) {
await STATUS(inv > 0 ? "inv" : "fwd", k);
let odd = new Complex();
for (let i = 0; i < m; i++)
odd.add(W(inv * i * k / m).mul(src[i * 2 + 1]));
tmp.push(odd);
}
for (let k = 0; k < m; k++) {
const even = tmp[k];
const odd = tmp[k + m].mul(W(inv * k / n));
dst[k] = new Complex(even).add(odd);
dst[k + m] = even.sub(odd);
}
return dst;
}
const src = Array(16).fill(0).map(_ => new Complex(Math.random() - 0.5));
DFT(src)
.then(dft => DFT(dft, 1))
.then(idft => {
idft.map(z => z.scale(src.length))
.forEach((z, i) => LOG(z.sub(src[i]).abs().toFixed(8)));
STATUS("finished");
});
3つあるkのループのうち、前のふたつをよく見ると、長さがn/2のDFTを計算になっている。
片方は入力の偶数項、もう一方は奇数項のみの計算をしている。
先ほどまでは偶数項・奇数項をひとつずつ交互にtmpにpushしていたが、今度は偶数項・奇数項をそれぞれまとめてtmpにずるずるとつなげているので、tmpの前半が偶数項のDFT、後半が奇数項のDFTである。
それを後処理することで最終的な結果が求まる。
それをはっきりさせるためにこの部分をDFT関数の呼び出しに書き換えたいのだが、そうすると結果の検証部分にも影響が出る。
そこで、DFT関数は素朴な
\(N^2\)
ループのDFTに戻し、今いぢっている関数をFFTと命名して、正逆DFTを行って元に戻るのを比較するのではなく、DFT関数とFFT関数の結果を比較するように変更しておく。
function formatNumeric(x) {
return x.toFixed(6).padStart(11);
}
class Complex {
constructor(re = 0, im = 0) {
if (re instanceof Complex) {
this.re = re.re;
this.im = re.im;
return;
}
this.re = re;
this.im = im;
}
toString() {
return `${formatNumeric(this.re)} ${formatNumeric(this.im)}`;
}
add(z) {
this.re += z.re;
this.im += z.im;
return this;
}
sub(z) {
this.re -= z.re;
this.im -= z.im;
return this;
}
mul(z) {
const re = this.re;
this.re = re * z.re - this.im * z.im;
this.im = re * z.im + this.im * z.re;
return this;
}
scale(r) {
this.re /= r;
this.im /= r;
return this;
}
abs() {
return Math.hypot(this.re, this.im);
}
}
function W(x) {
return new Complex(
Math.cos(2 * Math.PI * x),
Math.sin(2 * Math.PI * x)
);
}
async function DFT(src, inv = -1) {
const n = src.length;
const dst = [];
for (let k = 0; k < n; k++) {
await STATUS(inv > 0 ? "inv" : "fwd", k);
const sum = new Complex();
for (let i = 0; i < n; i++)
sum.add(W(inv * i * k / n).mul(src[i]));
dst.push(sum);
}
return dst;
}
async function FFT(src, inv = -1) {
const n = src.length;
const m = n / 2;
const dst = Array(n);
const tmp = [];
for (let k = 0; k < m; k++) {
await STATUS(inv > 0 ? "inv" : "fwd", k);
let even = new Complex();
for (let i = 0; i < m; i++)
even.add(W(inv * i * k / m).mul(src[i * 2]));
tmp.push(even);
}
for (let k = 0; k < m; k++) {
await STATUS(inv > 0 ? "inv" : "fwd", k);
let odd = new Complex();
for (let i = 0; i < m; i++)
odd.add(W(inv * i * k / m).mul(src[i * 2 + 1]));
tmp.push(odd);
}
for (let k = 0; k < m; k++) {
const even = tmp[k];
const odd = tmp[k + m].mul(W(inv * k / n));
dst[k] = new Complex(even).add(odd);
dst[k + m] = even.sub(odd);
}
return dst;
}
const src = Array(16).fill(0).map(_ => new Complex(Math.random() - 0.5));
DFT(src).then(async dft => {
const fft = await FFT(src);
fft.forEach((z, i) => LOG(z.sub(dft[i]).abs().toFixed(8)));
STATUS("finished");
});
その上でふたつ出てきたちびっ子DFTをDFT関数の呼び出しに置き換える。
入力はあらかじめ偶数項・奇数項に振り分けておく必要がある。
function formatNumeric(x) {
return x.toFixed(6).padStart(11);
}
class Complex {
constructor(re = 0, im = 0) {
if (re instanceof Complex) {
this.re = re.re;
this.im = re.im;
return;
}
this.re = re;
this.im = im;
}
toString() {
return `${formatNumeric(this.re)} ${formatNumeric(this.im)}`;
}
add(z) {
this.re += z.re;
this.im += z.im;
return this;
}
sub(z) {
this.re -= z.re;
this.im -= z.im;
return this;
}
mul(z) {
const re = this.re;
this.re = re * z.re - this.im * z.im;
this.im = re * z.im + this.im * z.re;
return this;
}
scale(r) {
this.re /= r;
this.im /= r;
return this;
}
abs() {
return Math.hypot(this.re, this.im);
}
}
function W(x) {
return new Complex(
Math.cos(2 * Math.PI * x),
Math.sin(2 * Math.PI * x)
);
}
async function DFT(src, inv = -1) {
const n = src.length;
const dst = [];
for (let k = 0; k < n; k++) {
await STATUS(inv > 0 ? "inv" : "fwd", k);
const sum = new Complex();
for (let i = 0; i < n; i++)
sum.add(W(inv * i * k / n).mul(src[i]));
dst.push(sum);
}
return dst;
}
async function FFT(src, inv = -1) {
const n = src.length;
const m = n / 2;
const dst = Array(n);
const even = [], odd = []
for (let k = 0; k < m; k++) {
even.push(src[k * 2]);
odd.push(src[k * 2 + 1]);
}
let tmp = await DFT(even, inv);
tmp = tmp.concat(await DFT(odd, inv));
for (let k = 0; k < m; k++) {
const even = tmp[k];
const odd = tmp[k + m].mul(W(inv * k / n));
dst[k] = new Complex(even).add(odd);
dst[k + m] = even.sub(odd);
}
return dst;
}
const src = Array(16).fill(0).map(_ => new Complex(Math.random() - 0.5));
DFT(src).then(async dft => {
const fft = await FFT(src);
fft.forEach((z, i) => LOG(z.sub(dft[i]).abs().toFixed(8)));
STATUS("finished");
});
これで偶数項・奇数項に分けてDFTしていることがはっきり分かるだろう。
function formatNumeric(x) {
return x.toFixed(6).padStart(11);
}
class Complex {
constructor(re = 0, im = 0) {
if (re instanceof Complex) {
this.re = re.re;
this.im = re.im;
return;
}
this.re = re;
this.im = im;
}
toString() {
return `${formatNumeric(this.re)} ${formatNumeric(this.im)}`;
}
add(z) {
this.re += z.re;
this.im += z.im;
return this;
}
sub(z) {
this.re -= z.re;
this.im -= z.im;
return this;
}
mul(z) {
const re = this.re;
this.re = re * z.re - this.im * z.im;
this.im = re * z.im + this.im * z.re;
return this;
}
scale(r) {
this.re /= r;
this.im /= r;
return this;
}
abs() {
return Math.hypot(this.re, this.im);
}
}
function W(x) {
return new Complex(
Math.cos(2 * Math.PI * x),
Math.sin(2 * Math.PI * x)
);
}
async function DFT(src, inv = -1) {
const n = src.length;
const dst = [];
for (let k = 0; k < n; k++) {
await STATUS(inv > 0 ? "inv" : "fwd", k);
const sum = new Complex();
for (let i = 0; i < n; i++)
sum.add(W(inv * i * k / n).mul(src[i]));
dst.push(sum);
}
return dst;
}
async function FFT(src, inv = -1) {
const n = src.length;
const m = n / 2;
const dst = Array(n);
const even = [], odd = []
for (let k = 0; k < m; k++) {
even.push(src[k * 2]);
odd.push(src[k * 2 + 1]);
}
let tmp = await (async (src, inv = -1) => {
const n = src.length;
const m = n / 2;
const dst = Array(n);
const even = [], odd = []
for (let k = 0; k < m; k++) {
even.push(src[k * 2]);
odd.push(src[k * 2 + 1]);
}
let tmp = await DFT(even, inv);
tmp = tmp.concat(await DFT(odd, inv));
for (let k = 0; k < m; k++) {
const even = tmp[k];
const odd = tmp[k + m].mul(W(inv * k / n));
dst[k] = new Complex(even).add(odd);
dst[k + m] = even.sub(odd);
}
return dst;
})(even, inv);
tmp = tmp.concat(await (async (src, inv = -1) => {
const n = src.length;
const m = n / 2;
const dst = Array(n);
const even = [], odd = []
for (let k = 0; k < m; k++) {
even.push(src[k * 2]);
odd.push(src[k * 2 + 1]);
}
let tmp = await DFT(even, inv);
tmp = tmp.concat(await DFT(odd, inv));
for (let k = 0; k < m; k++) {
const even = tmp[k];
const odd = tmp[k + m].mul(W(inv * k / n));
dst[k] = new Complex(even).add(odd);
dst[k + m] = even.sub(odd);
}
return dst;
})(odd, inv));
for (let k = 0; k < m; k++) {
const even = tmp[k];
const odd = tmp[k + m].mul(W(inv * k / n));
dst[k] = new Complex(even).add(odd);
dst[k + m] = even.sub(odd);
}
return dst;
}
const src = Array(16).fill(0).map(_ => new Complex(Math.random() - 0.5));
DFT(src).then(async dft => {
const fft = await FFT(src);
fft.forEach((z, i) => LOG(z.sub(dft[i]).abs().toFixed(8)));
STATUS("finished");
});
このままFFT関数を再帰呼び出ししても結果を求めることができるが、ここではループに展開することを目標にする。
そこで、もう一段分解したらどうなるかを調べるため、DFT関数の呼び出しをFFT関数の中身でそっくり置き換える。
アロー関数って便利だね。
JavaScriptはこういう検証にも案外向いているのかもしれない。
inv引数は別にアロー関数の引数として渡す必要がないので削除してしまおう。
function formatNumeric(x) {
return x.toFixed(6).padStart(11);
}
class Complex {
constructor(re = 0, im = 0) {
if (re instanceof Complex) {
this.re = re.re;
this.im = re.im;
return;
}
this.re = re;
this.im = im;
}
toString() {
return `${formatNumeric(this.re)} ${formatNumeric(this.im)}`;
}
add(z) {
this.re += z.re;
this.im += z.im;
return this;
}
sub(z) {
this.re -= z.re;
this.im -= z.im;
return this;
}
mul(z) {
const re = this.re;
this.re = re * z.re - this.im * z.im;
this.im = re * z.im + this.im * z.re;
return this;
}
scale(r) {
this.re /= r;
this.im /= r;
return this;
}
abs() {
return Math.hypot(this.re, this.im);
}
}
function W(x) {
return new Complex(
Math.cos(2 * Math.PI * x),
Math.sin(2 * Math.PI * x)
);
}
async function DFT(src, inv = -1) {
const n = src.length;
const dst = [];
for (let k = 0; k < n; k++) {
await STATUS(inv > 0 ? "inv" : "fwd", k);
const sum = new Complex();
for (let i = 0; i < n; i++)
sum.add(W(inv * i * k / n).mul(src[i]));
dst.push(sum);
}
return dst;
}
async function FFT(src, inv = -1) {
const n = src.length;
const m = n / 2;
const dst = Array(n);
const even = [], odd = []
for (let k = 0; k < m; k++) {
even.push(src[k * 2]);
odd.push(src[k * 2 + 1]);
}
let tmp = await (async (src) => {
const n = src.length;
const m = n / 2;
const dst = Array(n);
const even = [], odd = []
for (let k = 0; k < m; k++) {
even.push(src[k * 2]);
odd.push(src[k * 2 + 1]);
}
let tmp = await DFT(even, inv);
tmp = tmp.concat(await DFT(odd, inv));
for (let k = 0; k < m; k++) {
const even = tmp[k];
const odd = tmp[k + m].mul(W(inv * k / n));
dst[k] = new Complex(even).add(odd);
dst[k + m] = even.sub(odd);
}
return dst;
})(even);
tmp = tmp.concat(await (async (src) => {
const n = src.length;
const m = n / 2;
const dst = Array(n);
const even = [], odd = []
for (let k = 0; k < m; k++) {
even.push(src[k * 2]);
odd.push(src[k * 2 + 1]);
}
let tmp = await DFT(even, inv);
tmp = tmp.concat(await DFT(odd, inv));
for (let k = 0; k < m; k++) {
const even = tmp[k];
const odd = tmp[k + m].mul(W(inv * k / n));
dst[k] = new Complex(even).add(odd);
dst[k + m] = even.sub(odd);
}
return dst;
})(odd));
for (let k = 0; k < m; k++) {
const even = tmp[k];
const odd = tmp[k + m].mul(W(inv * k / n));
dst[k] = new Complex(even).add(odd);
dst[k + m] = even.sub(odd);
}
return dst;
}
const src = Array(16).fill(0).map(_ => new Complex(Math.random() - 0.5));
DFT(src).then(async dft => {
const fft = await FFT(src);
fft.forEach((z, i) => LOG(z.sub(dft[i]).abs().toFixed(8)));
STATUS("finished");
});
src引数は外側ではevenとoddだが、これも内側の関数で展開しよう。
function formatNumeric(x) {
return x.toFixed(6).padStart(11);
}
class Complex {
constructor(re = 0, im = 0) {
if (re instanceof Complex) {
this.re = re.re;
this.im = re.im;
return;
}
this.re = re;
this.im = im;
}
toString() {
return `${formatNumeric(this.re)} ${formatNumeric(this.im)}`;
}
add(z) {
this.re += z.re;
this.im += z.im;
return this;
}
sub(z) {
this.re -= z.re;
this.im -= z.im;
return this;
}
mul(z) {
const re = this.re;
this.re = re * z.re - this.im * z.im;
this.im = re * z.im + this.im * z.re;
return this;
}
scale(r) {
this.re /= r;
this.im /= r;
return this;
}
abs() {
return Math.hypot(this.re, this.im);
}
}
function W(x) {
return new Complex(
Math.cos(2 * Math.PI * x),
Math.sin(2 * Math.PI * x)
);
}
async function DFT(src, inv = -1) {
const n = src.length;
const dst = [];
for (let k = 0; k < n; k++) {
await STATUS(inv > 0 ? "inv" : "fwd", k);
const sum = new Complex();
for (let i = 0; i < n; i++)
sum.add(W(inv * i * k / n).mul(src[i]));
dst.push(sum);
}
return dst;
}
async function FFT(src, inv = -1) {
const n = src.length;
const m = n / 2;
const dst = Array(n);
let tmp = await (async () => {
const n = src.length / 2;
const m = n / 2;
const dst = Array(n);
const even = [], odd = []
for (let k = 0; k < m; k++) {
even.push(src[(k * 2) * 2]);
odd.push(src[(k * 2 + 1) * 2]);
}
let tmp = await DFT(even, inv);
tmp = tmp.concat(await DFT(odd, inv));
for (let k = 0; k < m; k++) {
const even = tmp[k];
const odd = tmp[k + m].mul(W(inv * k / n));
dst[k] = new Complex(even).add(odd);
dst[k + m] = even.sub(odd);
}
return dst;
})();
tmp = tmp.concat(await (async () => {
const n = src.length / 2;
const m = n / 2;
const dst = Array(n);
const even = [], odd = []
for (let k = 0; k < m; k++) {
even.push(src[(k * 2) * 2 + 1]);
odd.push(src[(k * 2 + 1) * 2 + 1]);
}
let tmp = await DFT(even, inv);
tmp = tmp.concat(await DFT(odd, inv));
for (let k = 0; k < m; k++) {
const even = tmp[k];
const odd = tmp[k + m].mul(W(inv * k / n));
dst[k] = new Complex(even).add(odd);
dst[k + m] = even.sub(odd);
}
return dst;
})());
for (let k = 0; k < m; k++) {
const even = tmp[k];
const odd = tmp[k + m].mul(W(inv * k / n));
dst[k] = new Complex(even).add(odd);
dst[k + m] = even.sub(odd);
}
return dst;
}
const src = Array(16).fill(0).map(_ => new Complex(Math.random() - 0.5));
DFT(src).then(async dft => {
const fft = await FFT(src);
fft.forEach((z, i) => LOG(z.sub(dft[i]).abs().toFixed(8)));
STATUS("finished");
});
すると外側のevenとoddが不要になる。
ここからが佳境である。
ループの中のevenとoddはsrcを偶奇に振り分けているだけである。
そんな無駄なコピーはしないで、実際の計算の時に初めて読み込めばいいではないか。
そのためにDFT関数にも協力してもらい、srcを読み込むときに添え字の変換関数を通すようにして、偶数項の計算をしたければ変換関数としてk=>2*kを渡すようにする。
奇数項の場合はk=>2*k+1である。
srcを飛び飛びに読んでしまうと、DFT点数をsrcの要素数から判断できなくなるため、nも別途渡す必要がある。
function formatNumeric(x) {
return x.toFixed(6).padStart(11);
}
class Complex {
constructor(re = 0, im = 0) {
if (re instanceof Complex) {
this.re = re.re;
this.im = re.im;
return;
}
this.re = re;
this.im = im;
}
toString() {
return `${formatNumeric(this.re)} ${formatNumeric(this.im)}`;
}
add(z) {
this.re += z.re;
this.im += z.im;
return this;
}
sub(z) {
this.re -= z.re;
this.im -= z.im;
return this;
}
mul(z) {
const re = this.re;
this.re = re * z.re - this.im * z.im;
this.im = re * z.im + this.im * z.re;
return this;
}
scale(r) {
this.re /= r;
this.im /= r;
return this;
}
abs() {
return Math.hypot(this.re, this.im);
}
}
function W(x) {
return new Complex(
Math.cos(2 * Math.PI * x),
Math.sin(2 * Math.PI * x)
);
}
async function DFT(src, inv = -1, n = null, idx = k => k) {
n ??= src.length;
const dst = [];
for (let k = 0; k < n; k++) {
await STATUS(inv > 0 ? "inv" : "fwd", k);
const sum = new Complex();
for (let i = 0; i < n; i++)
sum.add(W(inv * i * k / n).mul(src[idx(i)]));
dst.push(sum);
}
return dst;
}
async function FFT(src, inv = -1) {
const n = src.length;
const m = n / 2;
const dst = Array(n);
let tmp = await (async () => {
const n = src.length / 2;
const m = n / 2;
const dst = Array(n);
let tmp = await DFT(src, inv, m, k => (k * 2) * 2);
tmp = tmp.concat(await DFT(src, inv, m, k => (k * 2 + 1) * 2));
for (let k = 0; k < m; k++) {
const even = tmp[k];
const odd = tmp[k + m].mul(W(inv * k / n));
dst[k] = new Complex(even).add(odd);
dst[k + m] = even.sub(odd);
}
return dst;
})();
tmp = tmp.concat(await (async () => {
const n = src.length / 2;
const m = n / 2;
const dst = Array(n);
let tmp = await DFT(src, inv, m, k => (k * 2) * 2 + 1);
tmp = tmp.concat(await DFT(src, inv, m, k => (k * 2 + 1) * 2 + 1));
for (let k = 0; k < m; k++) {
const even = tmp[k];
const odd = tmp[k + m].mul(W(inv * k / n));
dst[k] = new Complex(even).add(odd);
dst[k + m] = even.sub(odd);
}
return dst;
})());
for (let k = 0; k < m; k++) {
const even = tmp[k];
const odd = tmp[k + m].mul(W(inv * k / n));
dst[k] = new Complex(even).add(odd);
dst[k + m] = even.sub(odd);
}
return dst;
}
const src = Array(16).fill(0).map(_ => new Complex(Math.random() - 0.5));
DFT(src).then(async dft => {
const fft = await FFT(src);
fft.forEach((z, i) => LOG(z.sub(dft[i]).abs().toFixed(8)));
STATUS("finished");
});
DFT関数は結果の比較対象としても使っているから、互換性を保たなければならない。
nは省略したらsrc.lengthである。
変換関数は恒等変換k=>kでよい。
アロー関数の引数がなくなったので、今度は戻り値に手をつける。
function formatNumeric(x) {
return x.toFixed(6).padStart(11);
}
class Complex {
constructor(re = 0, im = 0) {
if (re instanceof Complex) {
this.re = re.re;
this.im = re.im;
return;
}
this.re = re;
this.im = im;
}
toString() {
return `${formatNumeric(this.re)} ${formatNumeric(this.im)}`;
}
add(z) {
this.re += z.re;
this.im += z.im;
return this;
}
sub(z) {
this.re -= z.re;
this.im -= z.im;
return this;
}
mul(z) {
const re = this.re;
this.re = re * z.re - this.im * z.im;
this.im = re * z.im + this.im * z.re;
return this;
}
scale(r) {
this.re /= r;
this.im /= r;
return this;
}
abs() {
return Math.hypot(this.re, this.im);
}
}
function W(x) {
return new Complex(
Math.cos(2 * Math.PI * x),
Math.sin(2 * Math.PI * x)
);
}
async function DFT(src, inv = -1, n = null, idx = k => k) {
n ??= src.length;
const dst = [];
for (let k = 0; k < n; k++) {
await STATUS(inv > 0 ? "inv" : "fwd", k);
const sum = new Complex();
for (let i = 0; i < n; i++)
sum.add(W(inv * i * k / n).mul(src[idx(i)]));
dst.push(sum);
}
return dst;
}
async function FFT(src, inv = -1) {
const n = src.length;
const m = n / 2;
const dst = Array(n);
await (async () => {
const n = src.length / 2;
const m = n / 2;
let tmp = await DFT(src, inv, m, k => (k * 2) * 2);
tmp = tmp.concat(await DFT(src, inv, m, k => (k * 2 + 1) * 2));
for (let k = 0; k < m; k++) {
const even = tmp[k];
const odd = tmp[k + m].mul(W(inv * k / n));
dst[k] = new Complex(even).add(odd);
dst[k + m] = even.sub(odd);
}
})();
await (async () => {
const n = src.length / 2;
const m = n / 2;
let tmp = await DFT(src, inv, m, k => (k * 2) * 2 + 1);
tmp = tmp.concat(await DFT(src, inv, m, k => (k * 2 + 1) * 2 + 1));
for (let k = 0; k < m; k++) {
const even = tmp[k];
const odd = tmp[k + m].mul(W(inv * k / n));
dst[n + k] = new Complex(even).add(odd);
dst[n + k + m] = even.sub(odd);
}
})();
for (let k = 0; k < m; k++) {
const even = dst[k];
const odd = dst[k + m].mul(W(inv * k / n));
dst[k] = new Complex(even).add(odd);
dst[k + m] = even.sub(odd);
}
return dst;
}
const src = Array(16).fill(0).map(_ => new Complex(Math.random() - 0.5));
DFT(src).then(async dft => {
const fft = await FFT(src);
fft.forEach((z, i) => LOG(z.sub(dft[i]).abs().toFixed(8)));
STATUS("finished");
});
外側のdstに直接出力してしまえばよい。
これでtmpが消滅し、awaitもスッキリ。
function formatNumeric(x) {
return x.toFixed(6).padStart(11);
}
class Complex {
constructor(re = 0, im = 0) {
if (re instanceof Complex) {
this.re = re.re;
this.im = re.im;
return;
}
this.re = re;
this.im = im;
}
toString() {
return `${formatNumeric(this.re)} ${formatNumeric(this.im)}`;
}
add(z) {
this.re += z.re;
this.im += z.im;
return this;
}
sub(z) {
this.re -= z.re;
this.im -= z.im;
return this;
}
mul(z) {
const re = this.re;
this.re = re * z.re - this.im * z.im;
this.im = re * z.im + this.im * z.re;
return this;
}
scale(r) {
this.re /= r;
this.im /= r;
return this;
}
abs() {
return Math.hypot(this.re, this.im);
}
}
function W(x) {
return new Complex(
Math.cos(2 * Math.PI * x),
Math.sin(2 * Math.PI * x)
);
}
async function DFT(src, inv = -1, n = null, idx = k => k) {
n ??= src.length;
const dst = [];
for (let k = 0; k < n; k++) {
await STATUS(inv > 0 ? "inv" : "fwd", k);
const sum = new Complex();
for (let i = 0; i < n; i++)
sum.add(W(inv * i * k / n).mul(src[idx(i)]));
dst.push(sum);
}
return dst;
}
async function FFT(src, inv = -1) {
const n = src.length;
const m = n / 2;
const dst = Array(n);
await (async () => {
const n = src.length / 2;
const m = n / 2;
let tmp = await DFT(src, inv, m, k => (k * 2) * 2);
tmp = tmp.concat(await DFT(src, inv, m, k => (k * 2 + 1) * 2));
for (let k = 0; k < m; k++) {
const even = tmp[k];
const odd = tmp[k + m].mul(W(inv * k / n));
dst[k] = new Complex(even).add(odd);
dst[k + m] = even.sub(odd);
}
tmp = await DFT(src, inv, m, k => (k * 2) * 2 + 1);
tmp = tmp.concat(await DFT(src, inv, m, k => (k * 2 + 1) * 2 + 1));
for (let k = 0; k < m; k++) {
const even = tmp[k];
const odd = tmp[k + m].mul(W(inv * k / n));
dst[n + k] = new Complex(even).add(odd);
dst[n + k + m] = even.sub(odd);
}
})();
for (let k = 0; k < m; k++) {
const even = dst[k];
const odd = dst[k + m].mul(W(inv * k / n));
dst[k] = new Complex(even).add(odd);
dst[k + m] = even.sub(odd);
}
return dst;
}
const src = Array(16).fill(0).map(_ => new Complex(Math.random() - 0.5));
DFT(src).then(async dft => {
const fft = await FFT(src);
fft.forEach((z, i) => LOG(z.sub(dft[i]).abs().toFixed(8)));
STATUS("finished");
});
ふたつあるアロー関数のnとmの値は同じで、tmpはアロー関数の外側では使わないので、変数を使い回すことができる。
function formatNumeric(x) {
return x.toFixed(6).padStart(11);
}
class Complex {
constructor(re = 0, im = 0) {
if (re instanceof Complex) {
this.re = re.re;
this.im = re.im;
return;
}
this.re = re;
this.im = im;
}
toString() {
return `${formatNumeric(this.re)} ${formatNumeric(this.im)}`;
}
add(z) {
this.re += z.re;
this.im += z.im;
return this;
}
sub(z) {
this.re -= z.re;
this.im -= z.im;
return this;
}
mul(z) {
const re = this.re;
this.re = re * z.re - this.im * z.im;
this.im = re * z.im + this.im * z.re;
return this;
}
scale(r) {
this.re /= r;
this.im /= r;
return this;
}
abs() {
return Math.hypot(this.re, this.im);
}
}
function W(x) {
return new Complex(
Math.cos(2 * Math.PI * x),
Math.sin(2 * Math.PI * x)
);
}
async function DFT(src, inv = -1, n = null, idx = k => k) {
n ??= src.length;
const dst = [];
for (let k = 0; k < n; k++) {
await STATUS(inv > 0 ? "inv" : "fwd", k);
const sum = new Complex();
for (let i = 0; i < n; i++)
sum.add(W(inv * i * k / n).mul(src[idx(i)]));
dst.push(sum);
}
return dst;
}
async function FFT(src, inv = -1) {
const n = src.length;
const m = n / 2;
const dst = Array(n);
await (async () => {
const n = src.length / 2;
const m = n / 2;
let tmp = await DFT(src, inv, m, k => (k * 2) * 2);
tmp = tmp.concat(await DFT(src, inv, m, k => (k * 2 + 1) * 2));
tmp = tmp.concat(await DFT(src, inv, m, k => (k * 2) * 2 + 1));
tmp = tmp.concat(await DFT(src, inv, m, k => (k * 2 + 1) * 2 + 1));
for (let k = 0; k < m; k++) {
const even = tmp[k];
const odd = tmp[k + m].mul(W(inv * k / n));
dst[k] = new Complex(even).add(odd);
dst[k + m] = even.sub(odd);
}
for (let k = 0; k < m; k++) {
const even = tmp[n + k];
const odd = tmp[n + k + m].mul(W(inv * k / n));
dst[n + k] = new Complex(even).add(odd);
dst[n + k + m] = even.sub(odd);
}
})();
for (let k = 0; k < m; k++) {
const even = dst[k];
const odd = dst[k + m].mul(W(inv * k / n));
dst[k] = new Complex(even).add(odd);
dst[k + m] = even.sub(odd);
}
return dst;
}
const src = Array(16).fill(0).map(_ => new Complex(Math.random() - 0.5));
DFT(src).then(async dft => {
const fft = await FFT(src);
fft.forEach((z, i) => LOG(z.sub(dft[i]).abs().toFixed(8)));
STATUS("finished");
});
tmpに全部のDFT結果を押し込んでから後処理を行えば、DFTはDFTで、後処理は後処理で固めることができる。
ふたつ目のアロー関数から来たtmpの添え字が変わることになる。
ここでまたDFT関数に協力してもらい、出力先の配列と、出力する位置を指定できるようにする。
互換性を保つため、出力先が指定されていなければ内部で確保した配列を返す。
tmpはconcatをやめてまとめて確保する。
function formatNumeric(x) {
return x.toFixed(6).padStart(11);
}
class Complex {
constructor(re = 0, im = 0) {
if (re instanceof Complex) {
this.re = re.re;
this.im = re.im;
return;
}
this.re = re;
this.im = im;
}
toString() {
return `${formatNumeric(this.re)} ${formatNumeric(this.im)}`;
}
add(z) {
this.re += z.re;
this.im += z.im;
return this;
}
sub(z) {
this.re -= z.re;
this.im -= z.im;
return this;
}
mul(z) {
const re = this.re;
this.re = re * z.re - this.im * z.im;
this.im = re * z.im + this.im * z.re;
return this;
}
scale(r) {
this.re /= r;
this.im /= r;
return this;
}
abs() {
return Math.hypot(this.re, this.im);
}
}
function W(x) {
return new Complex(
Math.cos(2 * Math.PI * x),
Math.sin(2 * Math.PI * x)
);
}
async function DFT(src, inv = -1, n = null, idx = k => k, dst = null, offset = 0) {
n ??= src.length;
dst ??= Array(n);
for (let k = 0; k < n; k++) {
await STATUS(inv > 0 ? "inv" : "fwd", k);
const sum = new Complex();
for (let i = 0; i < n; i++)
sum.add(W(inv * i * k / n).mul(src[idx(i)]));
dst[k + offset] = sum;
}
return dst;
}
async function FFT(src, inv = -1) {
const n = src.length;
const m = n / 2;
const dst = Array(n);
await (async () => {
const n = src.length / 2;
const m = n / 2;
const tmp = Array(src.length);
await DFT(src, inv, m, k => (k * 2) * 2, tmp, 0);
await DFT(src, inv, m, k => (k * 2 + 1) * 2, tmp, m);
await DFT(src, inv, m, k => (k * 2) * 2 + 1, tmp, m * 2);
await DFT(src, inv, m, k => (k * 2 + 1) * 2 + 1, tmp, m * 3);
for (let k = 0; k < m; k++) {
const even = tmp[k];
const odd = tmp[k + m].mul(W(inv * k / n));
dst[k] = new Complex(even).add(odd);
dst[k + m] = even.sub(odd);
}
for (let k = 0; k < m; k++) {
const even = tmp[n + k];
const odd = tmp[n + k + m].mul(W(inv * k / n));
dst[n + k] = new Complex(even).add(odd);
dst[n + k + m] = even.sub(odd);
}
})();
for (let k = 0; k < m; k++) {
const even = dst[k];
const odd = dst[k + m].mul(W(inv * k / n));
dst[k] = new Complex(even).add(odd);
dst[k + m] = even.sub(odd);
}
return dst;
}
const src = Array(16).fill(0).map(_ => new Complex(Math.random() - 0.5));
DFT(src).then(async dft => {
const fft = await FFT(src);
fft.forEach((z, i) => LOG(z.sub(dft[i]).abs().toFixed(8)));
STATUS("finished");
});
そうするとtmpの代わりにdstが使えるようになり、tmpを削除できる。
function formatNumeric(x) {
return x.toFixed(6).padStart(11);
}
class Complex {
constructor(re = 0, im = 0) {
if (re instanceof Complex) {
this.re = re.re;
this.im = re.im;
return;
}
this.re = re;
this.im = im;
}
toString() {
return `${formatNumeric(this.re)} ${formatNumeric(this.im)}`;
}
add(z) {
this.re += z.re;
this.im += z.im;
return this;
}
sub(z) {
this.re -= z.re;
this.im -= z.im;
return this;
}
mul(z) {
const re = this.re;
this.re = re * z.re - this.im * z.im;
this.im = re * z.im + this.im * z.re;
return this;
}
scale(r) {
this.re /= r;
this.im /= r;
return this;
}
abs() {
return Math.hypot(this.re, this.im);
}
}
function W(x) {
return new Complex(
Math.cos(2 * Math.PI * x),
Math.sin(2 * Math.PI * x)
);
}
async function DFT(src, inv = -1, n = null, idx = k => k, dst = null, offset = 0) {
n ??= src.length;
dst ??= Array(n);
for (let k = 0; k < n; k++) {
await STATUS(inv > 0 ? "inv" : "fwd", k);
const sum = new Complex();
for (let i = 0; i < n; i++)
sum.add(W(inv * i * k / n).mul(src[idx(i)]));
dst[k + offset] = sum;
}
return dst;
}
async function FFT(src, inv = -1) {
const n = src.length;
const m = n / 2;
const dst = Array(n);
await (async () => {
const n = src.length / 2;
const m = n / 2;
await DFT(src, inv, m, k => (k * 2) * 2, dst, 0);
await DFT(src, inv, m, k => (k * 2 + 1) * 2, dst, m);
await DFT(src, inv, m, k => (k * 2) * 2 + 1, dst, m * 2);
await DFT(src, inv, m, k => (k * 2 + 1) * 2 + 1, dst, m * 3);
for (let k = 0; k < m; k++) {
const even = dst[k];
const odd = dst[k + m].mul(W(inv * k / n));
dst[k] = new Complex(even).add(odd);
dst[k + m] = even.sub(odd);
}
for (let k = 0; k < m; k++) {
const even = dst[n + k];
const odd = dst[n + k + m].mul(W(inv * k / n));
dst[n + k] = new Complex(even).add(odd);
dst[n + k + m] = even.sub(odd);
}
})();
for (let k = 0; k < m; k++) {
const even = dst[k];
const odd = dst[k + m].mul(W(inv * k / n));
dst[k] = new Complex(even).add(odd);
dst[k + m] = even.sub(odd);
}
return dst;
}
const src = Array(16).fill(0).map(_ => new Complex(Math.random() - 0.5));
DFT(src).then(async dft => {
const fft = await FFT(src);
fft.forEach((z, i) => LOG(z.sub(dft[i]).abs().toFixed(8)));
STATUS("finished");
});
function formatNumeric(x) {
return x.toFixed(6).padStart(11);
}
class Complex {
constructor(re = 0, im = 0) {
if (re instanceof Complex) {
this.re = re.re;
this.im = re.im;
return;
}
this.re = re;
this.im = im;
}
toString() {
return `${formatNumeric(this.re)} ${formatNumeric(this.im)}`;
}
add(z) {
this.re += z.re;
this.im += z.im;
return this;
}
sub(z) {
this.re -= z.re;
this.im -= z.im;
return this;
}
mul(z) {
const re = this.re;
this.re = re * z.re - this.im * z.im;
this.im = re * z.im + this.im * z.re;
return this;
}
scale(r) {
this.re /= r;
this.im /= r;
return this;
}
abs() {
return Math.hypot(this.re, this.im);
}
}
function W(x) {
return new Complex(
Math.cos(2 * Math.PI * x),
Math.sin(2 * Math.PI * x)
);
}
async function DFT(src, inv = -1, n = null, idx = k => k, dst = null, offset = 0) {
n ??= src.length;
dst ??= Array(n);
for (let k = 0; k < n; k++) {
await STATUS(inv > 0 ? "inv" : "fwd", k);
const sum = new Complex();
for (let i = 0; i < n; i++)
sum.add(W(inv * i * k / n).mul(src[idx(i)]));
dst[k + offset] = sum;
}
return dst;
}
async function FFT(src, inv = -1) {
const n = src.length;
const m = n / 2;
const dst = Array(n);
await (async () => {
const n = src.length / 2;
const m = n / 2;
await DFT(src, inv, m, k => (k * 2) * 2, dst, 0);
await DFT(src, inv, m, k => (k * 2 + 1) * 2, dst, m);
await DFT(src, inv, m, k => (k * 2) * 2 + 1, dst, m * 2);
await DFT(src, inv, m, k => (k * 2 + 1) * 2 + 1, dst, m * 3);
for (let j = 0; j < 2; j++) {
for (let k = 0; k < m; k++) {
const even = dst[n * j + k];
const odd = dst[n * j + k + m].mul(W(inv * k / n));
dst[n * j + k] = new Complex(even).add(odd);
dst[n * j + k + m] = even.sub(odd);
}
}
})();
for (let k = 0; k < m; k++) {
const even = dst[k];
const odd = dst[k + m].mul(W(inv * k / n));
dst[k] = new Complex(even).add(odd);
dst[k + m] = even.sub(odd);
}
return dst;
}
const src = Array(16).fill(0).map(_ => new Complex(Math.random() - 0.5));
DFT(src).then(async dft => {
const fft = await FFT(src);
fft.forEach((z, i) => LOG(z.sub(dft[i]).abs().toFixed(8)));
STATUS("finished");
});
ふたつある後処理をループにまとめる。
後処理のkのループはなんとなく元からあった一番外側の後処理と共通にできそうに見える。
そこで、変数名をかち合わないように変更してから、
function formatNumeric(x) {
return x.toFixed(6).padStart(11);
}
class Complex {
constructor(re = 0, im = 0) {
if (re instanceof Complex) {
this.re = re.re;
this.im = re.im;
return;
}
this.re = re;
this.im = im;
}
toString() {
return `${formatNumeric(this.re)} ${formatNumeric(this.im)}`;
}
add(z) {
this.re += z.re;
this.im += z.im;
return this;
}
sub(z) {
this.re -= z.re;
this.im -= z.im;
return this;
}
mul(z) {
const re = this.re;
this.re = re * z.re - this.im * z.im;
this.im = re * z.im + this.im * z.re;
return this;
}
scale(r) {
this.re /= r;
this.im /= r;
return this;
}
abs() {
return Math.hypot(this.re, this.im);
}
}
function W(x) {
return new Complex(
Math.cos(2 * Math.PI * x),
Math.sin(2 * Math.PI * x)
);
}
async function DFT(src, inv = -1, n = null, idx = k => k, dst = null, offset = 0) {
n ??= src.length;
dst ??= Array(n);
for (let k = 0; k < n; k++) {
await STATUS(inv > 0 ? "inv" : "fwd", k);
const sum = new Complex();
for (let i = 0; i < n; i++)
sum.add(W(inv * i * k / n).mul(src[idx(i)]));
dst[k + offset] = sum;
}
return dst;
}
async function FFT(src, inv = -1) {
const n = src.length;
const m = n / 2;
const dst = Array(n);
await (async () => {
const h = src.length / 2;
const q = h / 2;
await DFT(src, inv, q, k => (k * 2) * 2, dst, 0);
await DFT(src, inv, q, k => (k * 2 + 1) * 2, dst, q);
await DFT(src, inv, q, k => (k * 2) * 2 + 1, dst, q * 2);
await DFT(src, inv, q, k => (k * 2 + 1) * 2 + 1, dst, q * 3);
for (let j = 0; j < 2; j++) {
for (let k = 0; k < q; k++) {
const even = dst[h * j + k];
const odd = dst[h * j + k + q].mul(W(inv * k / h));
dst[h * j + k] = new Complex(even).add(odd);
dst[h * j + k + q] = even.sub(odd);
}
}
})();
for (let k = 0; k < m; k++) {
const even = dst[k];
const odd = dst[k + m].mul(W(inv * k / n));
dst[k] = new Complex(even).add(odd);
dst[k + m] = even.sub(odd);
}
return dst;
}
const src = Array(16).fill(0).map(_ => new Complex(Math.random() - 0.5));
DFT(src).then(async dft => {
const fft = await FFT(src);
fft.forEach((z, i) => LOG(z.sub(dft[i]).abs().toFixed(8)));
STATUS("finished");
});
アロー関数の外に出してみる。
function formatNumeric(x) {
return x.toFixed(6).padStart(11);
}
class Complex {
constructor(re = 0, im = 0) {
if (re instanceof Complex) {
this.re = re.re;
this.im = re.im;
return;
}
this.re = re;
this.im = im;
}
toString() {
return `${formatNumeric(this.re)} ${formatNumeric(this.im)}`;
}
add(z) {
this.re += z.re;
this.im += z.im;
return this;
}
sub(z) {
this.re -= z.re;
this.im -= z.im;
return this;
}
mul(z) {
const re = this.re;
this.re = re * z.re - this.im * z.im;
this.im = re * z.im + this.im * z.re;
return this;
}
scale(r) {
this.re /= r;
this.im /= r;
return this;
}
abs() {
return Math.hypot(this.re, this.im);
}
}
function W(x) {
return new Complex(
Math.cos(2 * Math.PI * x),
Math.sin(2 * Math.PI * x)
);
}
async function DFT(src, inv = -1, n = null, idx = k => k, dst = null, offset = 0) {
n ??= src.length;
dst ??= Array(n);
for (let k = 0; k < n; k++) {
await STATUS(inv > 0 ? "inv" : "fwd", k);
const sum = new Complex();
for (let i = 0; i < n; i++)
sum.add(W(inv * i * k / n).mul(src[idx(i)]));
dst[k + offset] = sum;
}
return dst;
}
async function FFT(src, inv = -1) {
const n = src.length;
const m = n / 2;
const dst = Array(n);
const q = m / 2;
await (async () => {
await DFT(src, inv, q, k => (k * 2) * 2, dst, 0);
await DFT(src, inv, q, k => (k * 2 + 1) * 2, dst, q);
await DFT(src, inv, q, k => (k * 2) * 2 + 1, dst, q * 2);
await DFT(src, inv, q, k => (k * 2 + 1) * 2 + 1, dst, q * 3);
})();
const h = src.length / 2;
for (let j = 0; j < 2; j++) {
for (let k = 0; k < q; k++) {
const even = dst[h * j + k];
const odd = dst[h * j + k + q].mul(W(inv * k / h));
dst[h * j + k] = new Complex(even).add(odd);
dst[h * j + k + q] = even.sub(odd);
}
}
for (let k = 0; k < m; k++) {
const even = dst[k];
const odd = dst[k + m].mul(W(inv * k / n));
dst[k] = new Complex(even).add(odd);
dst[k + m] = even.sub(odd);
}
return dst;
}
const src = Array(16).fill(0).map(_ => new Complex(Math.random() - 0.5));
DFT(src).then(async dft => {
const fft = await FFT(src);
fft.forEach((z, i) => LOG(z.sub(dft[i]).abs().toFixed(8)));
STATUS("finished");
});
そうすると、アロー関数の中身はきれいに添え字変換つきのDFTだけになる。 これで2回目の分解も終わりである。
DFTが倍に、後処理が1段増える。
function formatNumeric(x) {
return x.toFixed(6).padStart(11);
}
class Complex {
constructor(re = 0, im = 0) {
if (re instanceof Complex) {
this.re = re.re;
this.im = re.im;
return;
}
this.re = re;
this.im = im;
}
toString() {
return `${formatNumeric(this.re)} ${formatNumeric(this.im)}`;
}
add(z) {
this.re += z.re;
this.im += z.im;
return this;
}
sub(z) {
this.re -= z.re;
this.im -= z.im;
return this;
}
mul(z) {
const re = this.re;
this.re = re * z.re - this.im * z.im;
this.im = re * z.im + this.im * z.re;
return this;
}
scale(r) {
this.re /= r;
this.im /= r;
return this;
}
abs() {
return Math.hypot(this.re, this.im);
}
}
function W(x) {
return new Complex(
Math.cos(2 * Math.PI * x),
Math.sin(2 * Math.PI * x)
);
}
async function DFT(src, inv = -1, n = null, idx = k => k, dst = null, offset = 0) {
n ??= src.length;
dst ??= Array(n);
for (let k = 0; k < n; k++) {
await STATUS(inv > 0 ? "inv" : "fwd", k);
const sum = new Complex();
for (let i = 0; i < n; i++)
sum.add(W(inv * i * k / n).mul(src[idx(i)]));
dst[k + offset] = sum;
}
return dst;
}
async function FFT(src, inv = -1) {
const n = src.length;
const m = n / 2;
const dst = Array(n);
await (async () => {
const qq = n / 8;
await DFT(src, inv, qq, k => ((k * 2) * 2) * 2, dst, 0);
await DFT(src, inv, qq, k => ((k * 2 + 1) * 2) * 2, dst, qq);
await DFT(src, inv, qq, k => ((k * 2) * 2 + 1) * 2, dst, qq * 2);
await DFT(src, inv, qq, k => ((k * 2 + 1) * 2 + 1) * 2, dst, qq * 3);
await DFT(src, inv, qq, k => ((k * 2) * 2) * 2 + 1, dst, qq * 4);
await DFT(src, inv, qq, k => ((k * 2 + 1) * 2) * 2 + 1, dst, qq * 5);
await DFT(src, inv, qq, k => ((k * 2) * 2 + 1) * 2 + 1, dst, qq * 6);
await DFT(src, inv, qq, k => ((k * 2 + 1) * 2 + 1) * 2 + 1, dst, qq * 7);
})();
let h = src.length / 4;
let q = h / 2;
for (let j = 0; j < 4; j++) {
for (let k = 0; k < q; k++) {
const even = dst[h * j + k];
const odd = dst[h * j + k + q].mul(W(inv * k / h));
dst[h * j + k] = new Complex(even).add(odd);
dst[h * j + k + q] = even.sub(odd);
}
}
h *= 2;
q *= 2;
for (let j = 0; j < 2; j++) {
for (let k = 0; k < q; k++) {
const even = dst[h * j + k];
const odd = dst[h * j + k + q].mul(W(inv * k / h));
dst[h * j + k] = new Complex(even).add(odd);
dst[h * j + k + q] = even.sub(odd);
}
}
for (let k = 0; k < m; k++) {
const even = dst[k];
const odd = dst[k + m].mul(W(inv * k / n));
dst[k] = new Complex(even).add(odd);
dst[k + m] = even.sub(odd);
}
return dst;
}
const src = Array(16).fill(0).map(_ => new Complex(Math.random() - 0.5));
DFT(src).then(async dft => {
const fft = await FFT(src);
fft.forEach((z, i) => LOG(z.sub(dft[i]).abs().toFixed(8)));
STATUS("finished");
});
添え字変換関数がややこしいことになっている。
元々の関数がk=>F(k)だったとしたら、k=>F(k*2)とk=>F(k)*2の2通りが考えられる。
これは添え字変換関数を導入したあたりでのDFTの偶奇分解を見れば、k=>(k*2)*2+1とk=>(k*2+1)*2+1に分ける操作をしているので、置き換えるのは関数の外側ではなく、引数の方である。
つまり、前者のk=>F(k*2)が正しく、DFT関数の呼び出しを2行に複製して、kを(k*2)と(k*2+1)で置き換えればよい。
点数は半分になり、出力オフセットは単純に増分が半分になるだけである。
後処理は既にアロー関数の外に出してある。 処理点数を半分に、ループ回数を倍にすればいい。 あとでまとめてループにする予定なので、点数を示す変数は新たに作るのではなく、後処理が1段終わったら倍にしてそのまま同じ変数を使い回す。
そうすると、後処理の最後の段もループ回数1回のループで囲めば、後処理3段分をひとつのループにまとめられそうに見える。 とりあえず変数を揃えて、
function formatNumeric(x) {
return x.toFixed(6).padStart(11);
}
class Complex {
constructor(re = 0, im = 0) {
if (re instanceof Complex) {
this.re = re.re;
this.im = re.im;
return;
}
this.re = re;
this.im = im;
}
toString() {
return `${formatNumeric(this.re)} ${formatNumeric(this.im)}`;
}
add(z) {
this.re += z.re;
this.im += z.im;
return this;
}
sub(z) {
this.re -= z.re;
this.im -= z.im;
return this;
}
mul(z) {
const re = this.re;
this.re = re * z.re - this.im * z.im;
this.im = re * z.im + this.im * z.re;
return this;
}
scale(r) {
this.re /= r;
this.im /= r;
return this;
}
abs() {
return Math.hypot(this.re, this.im);
}
}
function W(x) {
return new Complex(
Math.cos(2 * Math.PI * x),
Math.sin(2 * Math.PI * x)
);
}
async function DFT(src, inv = -1, n = null, idx = k => k, dst = null, offset = 0) {
n ??= src.length;
dst ??= Array(n);
for (let k = 0; k < n; k++) {
await STATUS(inv > 0 ? "inv" : "fwd", k);
const sum = new Complex();
for (let i = 0; i < n; i++)
sum.add(W(inv * i * k / n).mul(src[idx(i)]));
dst[k + offset] = sum;
}
return dst;
}
async function FFT(src, inv = -1) {
const n = src.length;
const dst = Array(n);
await (async () => {
const qq = n / 8;
await DFT(src, inv, qq, k => ((k * 2) * 2) * 2, dst, 0);
await DFT(src, inv, qq, k => ((k * 2 + 1) * 2) * 2, dst, qq);
await DFT(src, inv, qq, k => ((k * 2) * 2 + 1) * 2, dst, qq * 2);
await DFT(src, inv, qq, k => ((k * 2 + 1) * 2 + 1) * 2, dst, qq * 3);
await DFT(src, inv, qq, k => ((k * 2) * 2) * 2 + 1, dst, qq * 4);
await DFT(src, inv, qq, k => ((k * 2 + 1) * 2) * 2 + 1, dst, qq * 5);
await DFT(src, inv, qq, k => ((k * 2) * 2 + 1) * 2 + 1, dst, qq * 6);
await DFT(src, inv, qq, k => ((k * 2 + 1) * 2 + 1) * 2 + 1, dst, qq * 7);
})();
let h = src.length / 4;
let q = h / 2;
for (let j = 0; j < 4; j++) {
for (let k = 0; k < q; k++) {
const even = dst[h * j + k];
const odd = dst[h * j + k + q].mul(W(inv * k / h));
dst[h * j + k] = new Complex(even).add(odd);
dst[h * j + k + q] = even.sub(odd);
}
}
h *= 2;
q *= 2;
for (let j = 0; j < 2; j++) {
for (let k = 0; k < q; k++) {
const even = dst[h * j + k];
const odd = dst[h * j + k + q].mul(W(inv * k / h));
dst[h * j + k] = new Complex(even).add(odd);
dst[h * j + k + q] = even.sub(odd);
}
}
h *= 2;
q *= 2;
for (let k = 0; k < q; k++) {
const even = dst[k];
const odd = dst[k + q].mul(W(inv * k / h));
dst[k] = new Complex(even).add(odd);
dst[k + q] = even.sub(odd);
}
return dst;
}
const src = Array(16).fill(0).map(_ => new Complex(Math.random() - 0.5));
DFT(src).then(async dft => {
const fft = await FFT(src);
fft.forEach((z, i) => LOG(z.sub(dft[i]).abs().toFixed(8)));
STATUS("finished");
});
えいやっと書き換える。
function formatNumeric(x) {
return x.toFixed(6).padStart(11);
}
class Complex {
constructor(re = 0, im = 0) {
if (re instanceof Complex) {
this.re = re.re;
this.im = re.im;
return;
}
this.re = re;
this.im = im;
}
toString() {
return `${formatNumeric(this.re)} ${formatNumeric(this.im)}`;
}
add(z) {
this.re += z.re;
this.im += z.im;
return this;
}
sub(z) {
this.re -= z.re;
this.im -= z.im;
return this;
}
mul(z) {
const re = this.re;
this.re = re * z.re - this.im * z.im;
this.im = re * z.im + this.im * z.re;
return this;
}
scale(r) {
this.re /= r;
this.im /= r;
return this;
}
abs() {
return Math.hypot(this.re, this.im);
}
}
function W(x) {
return new Complex(
Math.cos(2 * Math.PI * x),
Math.sin(2 * Math.PI * x)
);
}
async function DFT(src, inv = -1, n = null, idx = k => k, dst = null, offset = 0) {
n ??= src.length;
dst ??= Array(n);
for (let k = 0; k < n; k++) {
await STATUS(inv > 0 ? "inv" : "fwd", k);
const sum = new Complex();
for (let i = 0; i < n; i++)
sum.add(W(inv * i * k / n).mul(src[idx(i)]));
dst[k + offset] = sum;
}
return dst;
}
async function FFT(src, inv = -1) {
const n = src.length;
const dst = Array(n);
await (async () => {
const qq = n / 8;
await DFT(src, inv, qq, k => ((k * 2) * 2) * 2, dst, 0);
await DFT(src, inv, qq, k => ((k * 2 + 1) * 2) * 2, dst, qq);
await DFT(src, inv, qq, k => ((k * 2) * 2 + 1) * 2, dst, qq * 2);
await DFT(src, inv, qq, k => ((k * 2 + 1) * 2 + 1) * 2, dst, qq * 3);
await DFT(src, inv, qq, k => ((k * 2) * 2) * 2 + 1, dst, qq * 4);
await DFT(src, inv, qq, k => ((k * 2 + 1) * 2) * 2 + 1, dst, qq * 5);
await DFT(src, inv, qq, k => ((k * 2) * 2 + 1) * 2 + 1, dst, qq * 6);
await DFT(src, inv, qq, k => ((k * 2 + 1) * 2 + 1) * 2 + 1, dst, qq * 7);
})();
let h = src.length / 4;
let q = h / 2;
for (let m = 4; m >= 1; m /= 2) {
for (let j = 0; j < m; j++) {
for (let k = 0; k < q; k++) {
const even = dst[h * j + k];
const odd = dst[h * j + k + q].mul(W(inv * k / h));
dst[h * j + k] = new Complex(even).add(odd);
dst[h * j + k + q] = even.sub(odd);
}
}
h *= 2;
q *= 2;
}
return dst;
}
const src = Array(16).fill(0).map(_ => new Complex(Math.random() - 0.5));
DFT(src).then(async dft => {
const fft = await FFT(src);
fft.forEach((z, i) => LOG(z.sub(dft[i]).abs().toFixed(8)));
STATUS("finished");
});
3つバラバラだった後処理が、無事にひとつのループになった。 これでだいぶ見慣れたFFTのプログラムに近くなってきた。
jのループの中にある乗算をなくす。
このループで全体を演算する必要があるから、jの最終値はnになるはずである。
function formatNumeric(x) {
return x.toFixed(6).padStart(11);
}
class Complex {
constructor(re = 0, im = 0) {
if (re instanceof Complex) {
this.re = re.re;
this.im = re.im;
return;
}
this.re = re;
this.im = im;
}
toString() {
return `${formatNumeric(this.re)} ${formatNumeric(this.im)}`;
}
add(z) {
this.re += z.re;
this.im += z.im;
return this;
}
sub(z) {
this.re -= z.re;
this.im -= z.im;
return this;
}
mul(z) {
const re = this.re;
this.re = re * z.re - this.im * z.im;
this.im = re * z.im + this.im * z.re;
return this;
}
scale(r) {
this.re /= r;
this.im /= r;
return this;
}
abs() {
return Math.hypot(this.re, this.im);
}
}
function W(x) {
return new Complex(
Math.cos(2 * Math.PI * x),
Math.sin(2 * Math.PI * x)
);
}
async function DFT(src, inv = -1, n = null, idx = k => k, dst = null, offset = 0) {
n ??= src.length;
dst ??= Array(n);
for (let k = 0; k < n; k++) {
await STATUS(inv > 0 ? "inv" : "fwd", k);
const sum = new Complex();
for (let i = 0; i < n; i++)
sum.add(W(inv * i * k / n).mul(src[idx(i)]));
dst[k + offset] = sum;
}
return dst;
}
async function FFT(src, inv = -1) {
const n = src.length;
const dst = Array(n);
await (async () => {
const qq = n / 8;
await DFT(src, inv, qq, k => ((k * 2) * 2) * 2, dst, 0);
await DFT(src, inv, qq, k => ((k * 2 + 1) * 2) * 2, dst, qq);
await DFT(src, inv, qq, k => ((k * 2) * 2 + 1) * 2, dst, qq * 2);
await DFT(src, inv, qq, k => ((k * 2 + 1) * 2 + 1) * 2, dst, qq * 3);
await DFT(src, inv, qq, k => ((k * 2) * 2) * 2 + 1, dst, qq * 4);
await DFT(src, inv, qq, k => ((k * 2 + 1) * 2) * 2 + 1, dst, qq * 5);
await DFT(src, inv, qq, k => ((k * 2) * 2 + 1) * 2 + 1, dst, qq * 6);
await DFT(src, inv, qq, k => ((k * 2 + 1) * 2 + 1) * 2 + 1, dst, qq * 7);
})();
let h = src.length / 4;
let q = h / 2;
for (let m = 4; m >= 1; m /= 2) {
for (let j = 0; j < n; j += h) {
for (let k = 0; k < q; k++) {
const even = dst[j + k];
const odd = dst[j + k + q].mul(W(inv * k / h));
dst[j + k] = new Complex(even).add(odd);
dst[j + k + q] = even.sub(odd);
}
}
h *= 2;
q *= 2;
}
return dst;
}
const src = Array(16).fill(0).map(_ => new Complex(Math.random() - 0.5));
DFT(src).then(async dft => {
const fft = await FFT(src);
fft.forEach((z, i) => LOG(z.sub(dft[i]).abs().toFixed(8)));
STATUS("finished");
});
mのループはqかhを見ていれば終了が分かる。
hはどこまで行くかというと、jのループは最後の後処理で1度しか回らないわけだから、nが正解である。
h（Half、のつもり）なのにnまで行くのは気に入らない分かりづらいので、mで統一する。
「各段のnのようなもの」と捉えられるから、この方が分かりやすい。
function formatNumeric(x) {
return x.toFixed(6).padStart(11);
}
class Complex {
constructor(re = 0, im = 0) {
if (re instanceof Complex) {
this.re = re.re;
this.im = re.im;
return;
}
this.re = re;
this.im = im;
}
toString() {
return `${formatNumeric(this.re)} ${formatNumeric(this.im)}`;
}
add(z) {
this.re += z.re;
this.im += z.im;
return this;
}
sub(z) {
this.re -= z.re;
this.im -= z.im;
return this;
}
mul(z) {
const re = this.re;
this.re = re * z.re - this.im * z.im;
this.im = re * z.im + this.im * z.re;
return this;
}
scale(r) {
this.re /= r;
this.im /= r;
return this;
}
abs() {
return Math.hypot(this.re, this.im);
}
}
function W(x) {
return new Complex(
Math.cos(2 * Math.PI * x),
Math.sin(2 * Math.PI * x)
);
}
async function DFT(src, inv = -1, n = null, idx = k => k, dst = null, offset = 0) {
n ??= src.length;
dst ??= Array(n);
for (let k = 0; k < n; k++) {
await STATUS(inv > 0 ? "inv" : "fwd", k);
const sum = new Complex();
for (let i = 0; i < n; i++)
sum.add(W(inv * i * k / n).mul(src[idx(i)]));
dst[k + offset] = sum;
}
return dst;
}
async function FFT(src, inv = -1) {
const n = src.length;
const dst = Array(n);
await (async () => {
const qq = n / 8;
await DFT(src, inv, qq, k => ((k * 2) * 2) * 2, dst, 0);
await DFT(src, inv, qq, k => ((k * 2 + 1) * 2) * 2, dst, qq);
await DFT(src, inv, qq, k => ((k * 2) * 2 + 1) * 2, dst, qq * 2);
await DFT(src, inv, qq, k => ((k * 2 + 1) * 2 + 1) * 2, dst, qq * 3);
await DFT(src, inv, qq, k => ((k * 2) * 2) * 2 + 1, dst, qq * 4);
await DFT(src, inv, qq, k => ((k * 2 + 1) * 2) * 2 + 1, dst, qq * 5);
await DFT(src, inv, qq, k => ((k * 2) * 2 + 1) * 2 + 1, dst, qq * 6);
await DFT(src, inv, qq, k => ((k * 2 + 1) * 2 + 1) * 2 + 1, dst, qq * 7);
})();
for (let m = n/4; m <= n; m *= 2) {
const q = m / 2;
for (let j = 0; j < n; j += m) {
for (let k = 0; k < q; k++) {
const even = dst[j + k];
const odd = dst[j + k + q].mul(W(inv * k / m));
dst[j + k] = new Complex(even).add(odd);
dst[j + k + q] = even.sub(odd);
}
}
}
return dst;
}
const src = Array(16).fill(0).map(_ => new Complex(Math.random() - 0.5));
DFT(src).then(async dft => {
const fft = await FFT(src);
fft.forEach((z, i) => LOG(z.sub(dft[i]).abs().toFixed(8)));
STATUS("finished");
});
そうすると、q（Quarter）の方が実はHalfなので、qもhに書き換えよう。
function formatNumeric(x) {
return x.toFixed(6).padStart(11);
}
class Complex {
constructor(re = 0, im = 0) {
if (re instanceof Complex) {
this.re = re.re;
this.im = re.im;
return;
}
this.re = re;
this.im = im;
}
toString() {
return `${formatNumeric(this.re)} ${formatNumeric(this.im)}`;
}
add(z) {
this.re += z.re;
this.im += z.im;
return this;
}
sub(z) {
this.re -= z.re;
this.im -= z.im;
return this;
}
mul(z) {
const re = this.re;
this.re = re * z.re - this.im * z.im;
this.im = re * z.im + this.im * z.re;
return this;
}
scale(r) {
this.re /= r;
this.im /= r;
return this;
}
abs() {
return Math.hypot(this.re, this.im);
}
}
function W(x) {
return new Complex(
Math.cos(2 * Math.PI * x),
Math.sin(2 * Math.PI * x)
);
}
async function DFT(src, inv = -1, n = null, idx = k => k, dst = null, offset = 0) {
n ??= src.length;
dst ??= Array(n);
for (let k = 0; k < n; k++) {
await STATUS(inv > 0 ? "inv" : "fwd", k);
const sum = new Complex();
for (let i = 0; i < n; i++)
sum.add(W(inv * i * k / n).mul(src[idx(i)]));
dst[k + offset] = sum;
}
return dst;
}
async function FFT(src, inv = -1) {
const n = src.length;
const dst = Array(n);
await (async () => {
const qq = n / 8;
await DFT(src, inv, qq, k => ((k * 2) * 2) * 2, dst, 0);
await DFT(src, inv, qq, k => ((k * 2 + 1) * 2) * 2, dst, qq);
await DFT(src, inv, qq, k => ((k * 2) * 2 + 1) * 2, dst, qq * 2);
await DFT(src, inv, qq, k => ((k * 2 + 1) * 2 + 1) * 2, dst, qq * 3);
await DFT(src, inv, qq, k => ((k * 2) * 2) * 2 + 1, dst, qq * 4);
await DFT(src, inv, qq, k => ((k * 2 + 1) * 2) * 2 + 1, dst, qq * 5);
await DFT(src, inv, qq, k => ((k * 2) * 2 + 1) * 2 + 1, dst, qq * 6);
await DFT(src, inv, qq, k => ((k * 2 + 1) * 2 + 1) * 2 + 1, dst, qq * 7);
})();
for (let m = n/4; m <= n; m *= 2) {
const h = m / 2;
for (let j = 0; j < n; j += m) {
for (let k = 0; k < h; k++) {
const even = dst[j + k];
const odd = dst[j + k + h].mul(W(inv * k / m));
dst[j + k] = new Complex(even).add(odd);
dst[j + k + h] = even.sub(odd);
}
}
}
return dst;
}
const src = Array(16).fill(0).map(_ => new Complex(Math.random() - 0.5));
DFT(src).then(async dft => {
const fft = await FFT(src);
fft.forEach((z, i) => LOG(z.sub(dft[i]).abs().toFixed(8)));
STATUS("finished");
});
あとは8個あるDFT関数をループに直すだけなのだが、そのためには変換関数をどうにかしなければならない。
変換関数は先ほど検討した通り、一段ごとにkを(k*2)と(k*2+1)で置き換えればいいので、とりあえず文字列で最終形を確認してしまおう。
function formatNumeric(x) {
return x.toFixed(6).padStart(11);
}
class Complex {
constructor(re = 0, im = 0) {
if (re instanceof Complex) {
this.re = re.re;
this.im = re.im;
return;
}
this.re = re;
this.im = im;
}
toString() {
return `${formatNumeric(this.re)} ${formatNumeric(this.im)}`;
}
add(z) {
this.re += z.re;
this.im += z.im;
return this;
}
sub(z) {
this.re -= z.re;
this.im -= z.im;
return this;
}
mul(z) {
const re = this.re;
this.re = re * z.re - this.im * z.im;
this.im = re * z.im + this.im * z.re;
return this;
}
scale(r) {
this.re /= r;
this.im /= r;
return this;
}
abs() {
return Math.hypot(this.re, this.im);
}
}
function W(x) {
return new Complex(
Math.cos(2 * Math.PI * x),
Math.sin(2 * Math.PI * x)
);
}
async function DFT(src, inv = -1, n = null, idx = k => k, dst = null, offset = 0) {
n ??= src.length;
dst ??= Array(n);
for (let k = 0; k < n; k++) {
await STATUS(inv > 0 ? "inv" : "fwd", k);
const sum = new Complex();
for (let i = 0; i < n; i++)
sum.add(W(inv * i * k / n).mul(src[idx(i)]));
dst[k + offset] = sum;
}
return dst;
}
async function FFT(src, inv = -1) {
const n = src.length;
const dst = Array(n);
let f = [ "k * 2", "k * 2 + 1" ];
f = f.flatMap(k => [ k.replace("k", "(k * 2)"), k.replace("k", "(k * 2 + 1)") ]);
f = f.flatMap(k => [ k.replace("k", "(k * 2)"), k.replace("k", "(k * 2 + 1)") ]);
console.log(f);
await (async () => {
const qq = n / 8;
await DFT(src, inv, qq, k => ((k * 2) * 2) * 2, dst, 0);
await DFT(src, inv, qq, k => ((k * 2 + 1) * 2) * 2, dst, qq);
await DFT(src, inv, qq, k => ((k * 2) * 2 + 1) * 2, dst, qq * 2);
await DFT(src, inv, qq, k => ((k * 2 + 1) * 2 + 1) * 2, dst, qq * 3);
await DFT(src, inv, qq, k => ((k * 2) * 2) * 2 + 1, dst, qq * 4);
await DFT(src, inv, qq, k => ((k * 2 + 1) * 2) * 2 + 1, dst, qq * 5);
await DFT(src, inv, qq, k => ((k * 2) * 2 + 1) * 2 + 1, dst, qq * 6);
await DFT(src, inv, qq, k => ((k * 2 + 1) * 2 + 1) * 2 + 1, dst, qq * 7);
})();
for (let m = n/4; m <= n; m *= 2) {
const h = m / 2;
for (let j = 0; j < n; j += m) {
for (let k = 0; k < h; k++) {
const even = dst[j + k];
const odd = dst[j + k + h].mul(W(inv * k / m));
dst[j + k] = new Complex(even).add(odd);
dst[j + k + h] = even.sub(odd);
}
}
}
return dst;
}
const src = Array(16).fill(0).map(_ => new Complex(Math.random() - 0.5));
DFT(src).then(async dft => {
const fft = await FFT(src);
fft.forEach((z, i) => LOG(z.sub(dft[i]).abs().toFixed(8)));
STATUS("finished");
});
実行結果にはFFTの誤差しか出ていないが、ブラウザのログの方に関数の展開結果が出ている。
通常、配列を何か変換したいときはArrayを使う。
要素を削りたい場合はArrayだ。
要素を増やすためにはどうすればいいかというと、実はArrayが使える。
これはa.map()（ただしflatは1段分のみ）と同じなのだが、なぜ別の関数になっているかといえば、mapみたいな処理したいんだけど、要素数が変わる場合に便利だからなんだと思う。
要素を増やしたければ配列を返せばそれをバラして追加してくれるし、空の配列を返すと要素を削除することもできる。
実際の展開の方は、初期値を["k"]にしてflatMapを3回実行してもいいのだが、そうするとkの周りに余計なカッコがつくので、ソースとの比較のためにこうしてある。
ソースと同じ結果になっているのを確認したら、これを実際の関数に直そう。
function formatNumeric(x) {
return x.toFixed(6).padStart(11);
}
class Complex {
constructor(re = 0, im = 0) {
if (re instanceof Complex) {
this.re = re.re;
this.im = re.im;
return;
}
this.re = re;
this.im = im;
}
toString() {
return `${formatNumeric(this.re)} ${formatNumeric(this.im)}`;
}
add(z) {
this.re += z.re;
this.im += z.im;
return this;
}
sub(z) {
this.re -= z.re;
this.im -= z.im;
return this;
}
mul(z) {
const re = this.re;
this.re = re * z.re - this.im * z.im;
this.im = re * z.im + this.im * z.re;
return this;
}
scale(r) {
this.re /= r;
this.im /= r;
return this;
}
abs() {
return Math.hypot(this.re, this.im);
}
}
function W(x) {
return new Complex(
Math.cos(2 * Math.PI * x),
Math.sin(2 * Math.PI * x)
);
}
async function DFT(src, inv = -1, n = null, idx = k => k, dst = null, offset = 0) {
n ??= src.length;
dst ??= Array(n);
for (let k = 0; k < n; k++) {
await STATUS(inv > 0 ? "inv" : "fwd", k);
const sum = new Complex();
for (let i = 0; i < n; i++)
sum.add(W(inv * i * k / n).mul(src[idx(i)]));
dst[k + offset] = sum;
}
return dst;
}
async function FFT(src, inv = -1) {
const n = src.length;
const dst = Array(n);
let f = [ k => k * 2, k => k * 2 + 1 ];
f = f.flatMap(v => [ k => v(k * 2), k => v(k * 2 + 1) ]);
f = f.flatMap(v => [ k => v(k * 2), k => v(k * 2 + 1) ]);
f.forEach(v => console.log(v(0)));
await (async () => {
const qq = n / 8;
await DFT(src, inv, qq, k => ((k * 2) * 2) * 2, dst, 0);
await DFT(src, inv, qq, k => ((k * 2 + 1) * 2) * 2, dst, qq);
await DFT(src, inv, qq, k => ((k * 2) * 2 + 1) * 2, dst, qq * 2);
await DFT(src, inv, qq, k => ((k * 2 + 1) * 2 + 1) * 2, dst, qq * 3);
await DFT(src, inv, qq, k => ((k * 2) * 2) * 2 + 1, dst, qq * 4);
await DFT(src, inv, qq, k => ((k * 2 + 1) * 2) * 2 + 1, dst, qq * 5);
await DFT(src, inv, qq, k => ((k * 2) * 2 + 1) * 2 + 1, dst, qq * 6);
await DFT(src, inv, qq, k => ((k * 2 + 1) * 2 + 1) * 2 + 1, dst, qq * 7);
})();
for (let m = n/4; m <= n; m *= 2) {
const h = m / 2;
for (let j = 0; j < n; j += m) {
for (let k = 0; k < h; k++) {
const even = dst[j + k];
const odd = dst[j + k + h].mul(W(inv * k / m));
dst[j + k] = new Complex(even).add(odd);
dst[j + k + h] = even.sub(odd);
}
}
}
return dst;
}
const src = Array(16).fill(0).map(_ => new Complex(Math.random() - 0.5));
DFT(src).then(async dft => {
const fft = await FFT(src);
fft.forEach((z, i) => LOG(z.sub(dft[i]).abs().toFixed(8)));
STATUS("finished");
});
念のために書いておくと、fは関数の配列なので、flatMapのコールバック関数の引数vは関数になる。
ここから新しいアロー関数をふたつ作って配列にして返すと、flatMapはそれを展開して、元のfの倍の関数を含む配列を返してくる。
これを繰り返す。
実際にforEachで回して、関数だから適当に引数を指定して結果を確認すると、いわゆるビット逆順変換をやっていることが分かる。
ビット逆順変換はみな色々と苦労するところで、C言語だと再帰を使うのが一番簡単なのだが、JavaScriptだとflatMapでできることに気づいたのは、ある意味、この記事を書いた収穫である。
ということは、mapとflatがある処理系ならば同様に計算できるということでもある。
年をとるとこういうことにふと気づくことがある。
やっぱり経験を重ねてきたのが大きいのかな。
function formatNumeric(x) {
return x.toFixed(6).padStart(11);
}
class Complex {
constructor(re = 0, im = 0) {
if (re instanceof Complex) {
this.re = re.re;
this.im = re.im;
return;
}
this.re = re;
this.im = im;
}
toString() {
return `${formatNumeric(this.re)} ${formatNumeric(this.im)}`;
}
add(z) {
this.re += z.re;
this.im += z.im;
return this;
}
sub(z) {
this.re -= z.re;
this.im -= z.im;
return this;
}
mul(z) {
const re = this.re;
this.re = re * z.re - this.im * z.im;
this.im = re * z.im + this.im * z.re;
return this;
}
scale(r) {
this.re /= r;
this.im /= r;
return this;
}
abs() {
return Math.hypot(this.re, this.im);
}
}
function W(x) {
return new Complex(
Math.cos(2 * Math.PI * x),
Math.sin(2 * Math.PI * x)
);
}
async function DFT(src, inv = -1, n = null, idx = k => k, dst = null, offset = 0) {
n ??= src.length;
dst ??= Array(n);
for (let k = 0; k < n; k++) {
await STATUS(inv > 0 ? "inv" : "fwd", k);
const sum = new Complex();
for (let i = 0; i < n; i++)
sum.add(W(inv * i * k / n).mul(src[idx(i)]));
dst[k + offset] = sum;
}
return dst;
}
async function FFT(src, inv = -1) {
const n = src.length;
const dst = Array(n);
let f = [ k => k * 2, k => k * 2 + 1 ];
f = f.flatMap(v => [ k => v(k * 2), k => v(k * 2 + 1) ]);
f = f.flatMap(v => [ k => v(k * 2), k => v(k * 2 + 1) ]);
const qq = n / 8;
await Promise.all(f.map((f, i) => DFT(src, inv, qq, f, dst, qq * i)));
for (let m = n/4; m <= n; m *= 2) {
const h = m / 2;
for (let j = 0; j < n; j += m) {
for (let k = 0; k < h; k++) {
const even = dst[j + k];
const odd = dst[j + k + h].mul(W(inv * k / m));
dst[j + k] = new Complex(even).add(odd);
dst[j + k + h] = even.sub(odd);
}
}
}
return dst;
}
const src = Array(16).fill(0).map(_ => new Complex(Math.random() - 0.5));
DFT(src).then(async dft => {
const fft = await FFT(src);
fft.forEach((z, i) => LOG(z.sub(dft[i]).abs().toFixed(8)));
STATUS("finished");
});
ここまで来ればDFTの呼び出しをループにまとめるのは比較的簡単である。
ただ、DFTはasync関数で、forEachで回すとawaitできないという罠がある。
for ofで回せばよいのだが、qqの計算にiが欲しい。
そこで、iが取れるmapで回して、それをPromise.allで待つという、これまた大変にいい加減なことをしている。
DFT関数はasync関数だから、実際の戻り値はプロミスであり、mapの結果はプロミスの配列になるので、これで動く。
これもasync・awaitとプロミスの関係をきちんと理解していないと思いつかないだろう。
ただ、JavaScriptはシングルスレッドとは言え、実行順がどうなるか分からないので、STATUSによる表示は使い物にならない。
それぞれのブロックは重なっていないはずなので結果は大丈夫なはずである。
これはマルチスレッド化の可能性も示していることになる。
function formatNumeric(x) {
return x.toFixed(6).padStart(11);
}
class Complex {
constructor(re = 0, im = 0) {
if (re instanceof Complex) {
this.re = re.re;
this.im = re.im;
return;
}
this.re = re;
this.im = im;
}
toString() {
return `${formatNumeric(this.re)} ${formatNumeric(this.im)}`;
}
add(z) {
this.re += z.re;
this.im += z.im;
return this;
}
sub(z) {
this.re -= z.re;
this.im -= z.im;
return this;
}
mul(z) {
const re = this.re;
this.re = re * z.re - this.im * z.im;
this.im = re * z.im + this.im * z.re;
return this;
}
scale(r) {
this.re /= r;
this.im /= r;
return this;
}
abs() {
return Math.hypot(this.re, this.im);
}
}
function W(x) {
return new Complex(
Math.cos(2 * Math.PI * x),
Math.sin(2 * Math.PI * x)
);
}
async function DFT(src, inv = -1, n = null, idx = k => k, dst = null, offset = 0) {
n ??= src.length;
dst ??= Array(n);
for (let k = 0; k < n; k++) {
await STATUS(inv > 0 ? "inv" : "fwd", k);
const sum = new Complex();
for (let i = 0; i < n; i++)
sum.add(W(inv * i * k / n).mul(src[idx(i)]));
dst[k + offset] = sum;
}
return dst;
}
async function FFT(src, inv = -1) {
const n = src.length;
const dst = Array(n);
let f = [ k => k * 2, k => k * 2 + 1 ];
while (f.length < n)
f = f.flatMap(v => [ k => v(k * 2), k => v(k * 2 + 1) ]);
await Promise.all(f.map((f, i) => DFT(src, inv, 1, f, dst, i)));
for (let m = 2; m <= n; m *= 2) {
const h = m / 2;
for (let j = 0; j < n; j += m) {
for (let k = 0; k < h; k++) {
const even = dst[j + k];
const odd = dst[j + k + h].mul(W(inv * k / m));
dst[j + k] = new Complex(even).add(odd);
dst[j + k + h] = even.sub(odd);
}
}
}
return dst;
}
const src = Array(16).fill(0).map(_ => new Complex(Math.random() - 0.5));
DFT(src).then(async dft => {
const fft = await FFT(src);
fft.forEach((z, i) => LOG(z.sub(dft[i]).abs().toFixed(8)));
STATUS("finished");
});
DFT関数呼び出しの要素数が1になるようにループ回数を調整する。
要素数が1なのだから、DFT関数の呼び出しはn回になり、変換関数もn個必要になるので、n個できるまでループすればいい。
mの初期値はh=m/2という式があるので、最低が1ということはなく、2が正解である。
ここで、1点DFTというのはどういう演算なのか？ と考えると、単純に入力をそのまま返す演算だと分かる。
したがって、DFT関数は呼び出す必要すらなく、入力を並び替えてdstに設定すればよい。
function formatNumeric(x) {
return x.toFixed(6).padStart(11);
}
class Complex {
constructor(re = 0, im = 0) {
if (re instanceof Complex) {
this.re = re.re;
this.im = re.im;
return;
}
this.re = re;
this.im = im;
}
toString() {
return `${formatNumeric(this.re)} ${formatNumeric(this.im)}`;
}
add(z) {
this.re += z.re;
this.im += z.im;
return this;
}
sub(z) {
this.re -= z.re;
this.im -= z.im;
return this;
}
mul(z) {
const re = this.re;
this.re = re * z.re - this.im * z.im;
this.im = re * z.im + this.im * z.re;
return this;
}
scale(r) {
this.re /= r;
this.im /= r;
return this;
}
abs() {
return Math.hypot(this.re, this.im);
}
}
function W(x) {
return new Complex(
Math.cos(2 * Math.PI * x),
Math.sin(2 * Math.PI * x)
);
}
async function DFT(src, inv = -1, n = null, idx = k => k, dst = null, offset = 0) {
n ??= src.length;
dst ??= Array(n);
for (let k = 0; k < n; k++) {
await STATUS(inv > 0 ? "inv" : "fwd", k);
const sum = new Complex();
for (let i = 0; i < n; i++)
sum.add(W(inv * i * k / n).mul(src[idx(i)]));
dst[k + offset] = sum;
}
return dst;
}
async function FFT(src, inv = -1) {
const n = src.length;
const dst = Array(n);
let f = [ k => k * 2, k => k * 2 + 1 ];
while (f.length < n)
f = f.flatMap(v => [ k => v(k * 2), k => v(k * 2 + 1) ]);
f.forEach((v, i) => dst[i] = new Complex(src[v(0)]));
for (let m = 2; m <= n; m *= 2) {
const h = m / 2;
for (let j = 0; j < n; j += m) {
for (let k = 0; k < h; k++) {
const even = dst[j + k];
const odd = dst[j + k + h].mul(W(inv * k / m));
dst[j + k] = new Complex(even).add(odd);
dst[j + k + h] = even.sub(odd);
}
}
}
return dst;
}
const src = Array(16).fill(0).map(_ => new Complex(Math.random() - 0.5));
DFT(src).then(async dft => {
const fft = await FFT(src);
fft.forEach((z, i) => LOG(z.sub(dft[i]).abs().toFixed(8)));
STATUS("finished");
});
変換関数は点数で決まるため、同じ点数のFFTを繰り返す場合は最初の1度だけ計算すればよく、値もあらかじめ計算しておくのが普通である。
function formatNumeric(x) {
return x.toFixed(6).padStart(11);
}
class Complex {
constructor(re = 0, im = 0) {
if (re instanceof Complex) {
this.re = re.re;
this.im = re.im;
return;
}
this.re = re;
this.im = im;
}
toString() {
return `${formatNumeric(this.re)} ${formatNumeric(this.im)}`;
}
add(z) {
this.re += z.re;
this.im += z.im;
return this;
}
sub(z) {
this.re -= z.re;
this.im -= z.im;
return this;
}
mul(z) {
const re = this.re;
this.re = re * z.re - this.im * z.im;
this.im = re * z.im + this.im * z.re;
return this;
}
scale(r) {
this.re /= r;
this.im /= r;
return this;
}
abs() {
return Math.hypot(this.re, this.im);
}
}
function W(x) {
return new Complex(
Math.cos(2 * Math.PI * x),
Math.sin(2 * Math.PI * x)
);
}
async function DFT(src, inv = -1) {
const n = src.length;
const dst = [];
for (let k = 0; k < n; k++) {
await STATUS(inv > 0 ? "inv" : "fwd", k);
const sum = new Complex();
for (let i = 0; i < n; i++)
sum.add(W(inv * i * k / n).mul(src[i]));
dst.push(sum);
}
return dst;
}
async function FFT(src, inv = -1) {
const n = src.length;
const dst = Array(n);
let f = [ k => k * 2, k => k * 2 + 1 ];
while (f.length < n)
f = f.flatMap(v => [ k => v(k * 2), k => v(k * 2 + 1) ]);
const idx = f.map(v => v(0));
console.log(idx);
idx.forEach((v, i) => dst[i] = new Complex(src[v]));
for (let m = 2; m <= n; m *= 2) {
const h = m / 2;
for (let j = 0; j < n; j += m) {
for (let k = 0; k < h; k++) {
const even = dst[j + k];
const odd = dst[j + k + h].mul(W(inv * k / m));
dst[j + k] = new Complex(even).add(odd);
dst[j + k + h] = even.sub(odd);
}
}
}
return dst;
}
const src = Array(16).fill(0).map(_ => new Complex(Math.random() - 0.5));
DFT(src).then(async dft => {
const fft = await FFT(src);
fft.forEach((z, i) => LOG(z.sub(dft[i]).abs().toFixed(8)));
STATUS("finished");
});
FFT関数からDFT関数を呼び出すこともなくなったので、DFT関数も元の素朴な形に戻してある。
変換関数は要するにビット逆順の添え字テーブルなのだから、以下のように生成でき、関数をわざわざ作る必要もない。
function formatNumeric(x) {
return x.toFixed(6).padStart(11);
}
class Complex {
constructor(re = 0, im = 0) {
if (re instanceof Complex) {
this.re = re.re;
this.im = re.im;
return;
}
this.re = re;
this.im = im;
}
toString() {
return `${formatNumeric(this.re)} ${formatNumeric(this.im)}`;
}
add(z) {
this.re += z.re;
this.im += z.im;
return this;
}
sub(z) {
this.re -= z.re;
this.im -= z.im;
return this;
}
mul(z) {
const re = this.re;
this.re = re * z.re - this.im * z.im;
this.im = re * z.im + this.im * z.re;
return this;
}
scale(r) {
this.re /= r;
this.im /= r;
return this;
}
abs() {
return Math.hypot(this.re, this.im);
}
}
function W(x) {
return new Complex(
Math.cos(2 * Math.PI * x),
Math.sin(2 * Math.PI * x)
);
}
async function DFT(src, inv = -1) {
const n = src.length;
const dst = [];
for (let k = 0; k < n; k++) {
await STATUS(inv > 0 ? "inv" : "fwd", k);
const sum = new Complex();
for (let i = 0; i < n; i++)
sum.add(W(inv * i * k / n).mul(src[i]));
dst.push(sum);
}
return dst;
}
async function FFT(src, inv = -1) {
await STATUS("FFT: doing", src.length);
const n = src.length;
const dst = Array(n);
let idx = [ 0 ];
while (idx.length < n)
idx = idx.flatMap(v => [ v, v + idx.length]);
console.log(idx);
idx.forEach((v, i) => dst[i] = new Complex(src[v]));
for (let m = 2; m <= n; m *= 2) {
const h = m / 2;
for (let j = 0; j < n; j += m) {
for (let k = 0; k < h; k++) {
const even = dst[j + k];
const odd = dst[j + k + h].mul(W(inv * k / m));
dst[j + k] = new Complex(even).add(odd);
dst[j + k + h] = even.sub(odd);
}
}
}
await STATUS("FFT: done");
return dst;
}
const src = Array(16).fill(0).map(_ => new Complex(Math.random() - 0.5));
DFT(src).then(async dft => {
const fft = await FFT(src);
fft.forEach((z, i) => LOG(z.sub(dft[i]).abs().toFixed(8)));
STATUS("finished");
});
これが一般的なFFTルーチンになる。 そして、ここからさらに無駄な演算を省く最適化を行うのが普通である。 本当はなんとなく最適化した実FFTまで一気に書くつもりだったが、ていねいにdiffを載せてたらすごく長くなってしまったので、続きは次回。
16 Feb 2026: 新規作成
ご意見・ご要望の送り先は あかもず仮店舗 の末尾をご覧ください。
Copyright (C) 2026 akamoz.jp