前回の続きです。 実数入力の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) {
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");
});
ソースを少し保守しておく。
"use strict";
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]);
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));
FFT(src).then(async fft => {
const dft = await DFT(src);
fft.forEach((z, i) => LOG(z.sub(dft[i]).abs().toFixed(8)));
STATUS("finished");
});
strictモードにして、FFTを先に計算するようにする。
"use strict";が効いているかどうかは、グローバルスコープでx=1;と書いてみれば分かる。
単なる文字列だからよく間違えるんだな、これが。
DFTはコードとしては安定しているが、点数を増やすと計算に時間がかかる。 FFTはゴリゴリ書き換えているのでよくエラーが混ざり込む。 エラーの発生をDFTの計算完了まで待たされるのは時間の無駄なので、FFTを先に計算する。
まず、jとkのループは機械的に入れ替えが可能だ。
"use strict";
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]);
idx.forEach((v, i) => dst[i] = new Complex(src[v]));
for (let m = 2; m <= n; m *= 2) {
const h = m / 2;
for (let k = 0; k < h; k++) {
for (let j = 0; j < n; j += m) {
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));
FFT(src).then(async fft => {
const dft = await DFT(src);
fft.forEach((z, i) => LOG(z.sub(dft[i]).abs().toFixed(8)));
STATUS("finished");
});
すると、回転因子はjには依存しないから、jのループの中では不変定数になり、ループの外に出せる。
"use strict";
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]);
idx.forEach((v, i) => dst[i] = new Complex(src[v]));
for (let m = 2; m <= n; m *= 2) {
const h = m / 2;
for (let k = 0; k < h; k++) {
const w = W(inv * k / m);
for (let j = 0; j < n; j += m) {
const even = dst[j + k];
const odd = dst[j + k + h].mul(w);
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));
FFT(src).then(async fft => {
const dft = await DFT(src);
fft.forEach((z, i) => LOG(z.sub(dft[i]).abs().toFixed(8)));
STATUS("finished");
});
これだけで回転因子の計算回数は激減するし、実用上は回転因子は表引きにすることが多いが、表引きも楽になる。
実FFTは \(F(k)\) と \(F(N-k)\) が複素共役の関係になる。 \(k=0\) と \(k=N/2\) の時は実数になり、複素共役になる相手がいないことに注意。 この性質が色々なところで面倒を引き起こす。
各ループではちびっ子FFTの点数はm/2、つまりh点なので、その半分のqについて複素共役対称になる。
しかし、k=0とk=qの点は相手が実数になり、複素共役になる相手がいない。
そこで、kのループをk=0・k=qとそれ以外に分ける。
"use strict";
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]);
idx.forEach((v, i) => dst[i] = new Complex(src[v]));
for (let m = 2; m <= n; m *= 2) {
const h = m / 2;
const q = h / 2;
let k = 0;
for (; k < 1; k++) {
const w = W(inv * k / m);
for (let j = 0; j < n; j += m) {
const even = dst[j + k];
const odd = dst[j + k + h].mul(w);
dst[j + k] = new Complex(even).add(odd);
dst[j + k + h] = even.sub(odd);
}
}
for (; k < q; k++) {
const w = W(inv * k / m);
for (let j = 0; j < n; j += m) {
const even = dst[j + k];
const odd = dst[j + k + h].mul(w);
dst[j + k] = new Complex(even).add(odd);
dst[j + k + h] = even.sub(odd);
}
}
for (; k <= q; k++) {
const w = W(inv * k / m);
for (let j = 0; j < n; j += m) {
const even = dst[j + k];
const odd = dst[j + k + h].mul(w);
dst[j + k] = new Complex(even).add(odd);
dst[j + k + h] = even.sub(odd);
}
}
for (; k < h; k++) {
const w = W(inv * k / m);
for (let j = 0; j < n; j += m) {
const even = dst[j + k];
const odd = dst[j + k + h].mul(w);
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));
FFT(src).then(async fft => {
const dft = await DFT(src);
fft.forEach((z, i) => LOG(z.sub(dft[i]).abs().toFixed(8)));
STATUS("finished");
});
3つ目のループはループの形はしているが、実際にはk=qの1回しか実行されない。
k=0の場合を整理する。
"use strict";
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]);
idx.forEach((v, i) => dst[i] = new Complex(src[v]));
for (let m = 2; m <= n; m *= 2) {
const h = m / 2;
const q = h / 2;
for (let j = 0; j < n; j += m) {
const even = dst[j];
const odd = dst[j + h];
dst[j] = new Complex(even).add(odd);
dst[j + h] = even.sub(odd);
}
let k = 1;
for (; k < q; k++) {
const w = W(inv * k / m);
for (let j = 0; j < n; j += m) {
const even = dst[j + k];
const odd = dst[j + k + h].mul(w);
dst[j + k] = new Complex(even).add(odd);
dst[j + k + h] = even.sub(odd);
}
}
for (; k <= q; k++) {
const w = W(inv * k / m);
for (let j = 0; j < n; j += m) {
const even = dst[j + k];
const odd = dst[j + k + h].mul(w);
dst[j + k] = new Complex(even).add(odd);
dst[j + k + h] = even.sub(odd);
}
}
for (; k < h; k++) {
const w = W(inv * k / m);
for (let j = 0; j < n; j += m) {
const even = dst[j + k];
const odd = dst[j + k + h].mul(w);
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));
FFT(src).then(async fft => {
const dft = await DFT(src);
fft.forEach((z, i) => LOG(z.sub(dft[i]).abs().toFixed(8)));
STATUS("finished");
});
回転因子が1になるので乗算が省略できる。
k=qの場合も整理する。
"use strict";
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]);
idx.forEach((v, i) => dst[i] = new Complex(src[v]));
for (let m = 2; m <= n; m *= 2) {
const h = m / 2;
const q = h / 2;
for (let j = 0; j < n; j += m) {
const even = dst[j];
const odd = dst[j + h];
dst[j] = new Complex(even).add(odd);
dst[j + h] = even.sub(odd);
}
if (h < 2)
continue;
let k = 1;
for (; k < q; k++) {
const w = W(inv * k / m);
for (let j = 0; j < n; j += m) {
const even = dst[j + k];
const odd = dst[j + k + h].mul(w);
dst[j + k] = new Complex(even).add(odd);
dst[j + k + h] = even.sub(odd);
}
}
const w = W(inv * q / m);
for (let j = 0; j < n; j += m) {
const even = dst[j + q];
const odd = dst[j + q + h].mul(w);
dst[j + q] = new Complex(even).add(odd);
dst[j + q + h] = even.sub(odd);
}
for (k++; k < h; k++) {
const w = W(inv * k / m);
for (let j = 0; j < n; j += m) {
const even = dst[j + k];
const odd = dst[j + k + h].mul(w);
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));
FFT(src).then(async fft => {
const dft = await DFT(src);
fft.forEach((z, i) => LOG(z.sub(dft[i]).abs().toFixed(8)));
STATUS("finished");
});
if文が増えており、hが1の場合はそれより下のループを実行しないようにしている。
このif文がないとq=0.5について計算しようとしてしまい、おかしな結果になってしまう。
回転因子を約分して、
"use strict";
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]);
idx.forEach((v, i) => dst[i] = new Complex(src[v]));
for (let m = 2; m <= n; m *= 2) {
const h = m / 2;
const q = h / 2;
for (let j = 0; j < n; j += m) {
const even = dst[j];
const odd = dst[j + h];
dst[j] = new Complex(even).add(odd);
dst[j + h] = even.sub(odd);
}
if (h < 2)
continue;
let k = 1;
for (; k < q; k++) {
const w = W(inv * k / m);
for (let j = 0; j < n; j += m) {
const even = dst[j + k];
const odd = dst[j + k + h].mul(w);
dst[j + k] = new Complex(even).add(odd);
dst[j + k + h] = even.sub(odd);
}
}
const w = W(inv / 4);
for (let j = 0; j < n; j += m) {
const even = dst[j + q];
const odd = dst[j + q + h].mul(w);
dst[j + q] = new Complex(even).add(odd);
dst[j + q + h] = even.sub(odd);
}
for (k++; k < h; k++) {
const w = W(inv * k / m);
for (let j = 0; j < n; j += m) {
const even = dst[j + k];
const odd = dst[j + k + h].mul(w);
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));
FFT(src).then(async fft => {
const dft = await DFT(src);
fft.forEach((z, i) => LOG(z.sub(dft[i]).abs().toFixed(8)));
STATUS("finished");
});
"use strict";
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]);
idx.forEach((v, i) => dst[i] = new Complex(src[v]));
for (let m = 2; m <= n; m *= 2) {
const h = m / 2;
const q = h / 2;
for (let j = 0; j < n; j += m) {
const even = dst[j];
const odd = dst[j + h];
dst[j] = new Complex(even).add(odd);
dst[j + h] = even.sub(odd);
}
if (h < 2)
continue;
let k = 1;
for (; k < q; k++) {
const w = W(inv * k / m);
for (let j = 0; j < n; j += m) {
const even = dst[j + k];
const odd = dst[j + k + h].mul(w);
dst[j + k] = new Complex(even).add(odd);
dst[j + k + h] = even.sub(odd);
}
}
const w = new Complex(0, inv);
for (let j = 0; j < n; j += m) {
const even = dst[j + q];
const odd = dst[j + q + h].mul(w);
dst[j + q] = new Complex(even).add(odd);
dst[j + q + h] = even.sub(odd);
}
for (k++; k < h; k++) {
const w = W(inv * k / m);
for (let j = 0; j < n; j += m) {
const even = dst[j + k];
const odd = dst[j + k + h].mul(w);
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));
FFT(src).then(async fft => {
const dft = await DFT(src);
fft.forEach((z, i) => LOG(z.sub(dft[i]).abs().toFixed(8)));
STATUS("finished");
});
この時点ではeven・oddとも実数で、odd側に虚数単位がかかって加減算されることで、共役な複素数が生まれる。
"use strict";
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]);
idx.forEach((v, i) => dst[i] = new Complex(src[v]));
for (let m = 2; m <= n; m *= 2) {
const h = m / 2;
const q = h / 2;
for (let j = 0; j < n; j += m) {
const even = dst[j];
const odd = dst[j + h];
dst[j] = new Complex(even).add(odd);
dst[j + h] = even.sub(odd);
}
if (h < 2)
continue;
let k = 1;
for (; k < q; k++) {
const w = W(inv * k / m);
for (let j = 0; j < n; j += m) {
const even = dst[j + k];
const odd = dst[j + k + h].mul(w);
dst[j + k] = new Complex(even).add(odd);
dst[j + k + h] = even.sub(odd);
}
}
for (let j = 0; j < n; j += m) {
const even = new Complex(dst[j + q].re, 0);
const odd = new Complex(0, inv * dst[j + q + h].re);
dst[j + q] = new Complex(even).add(odd);
dst[j + q + h] = even.sub(odd);
}
for (k++; k < h; k++) {
const w = W(inv * k / m);
for (let j = 0; j < n; j += m) {
const even = dst[j + k];
const odd = dst[j + k + h].mul(w);
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));
FFT(src).then(async fft => {
const dft = await DFT(src);
fft.forEach((z, i) => LOG(z.sub(dft[i]).abs().toFixed(8)));
STATUS("finished");
});
実部・虚部が簡単に得られ、複素共役になると分かっているのだから、結果を書き込む時に複素共役の複素数を作ればよい。
"use strict";
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]);
idx.forEach((v, i) => dst[i] = new Complex(src[v]));
for (let m = 2; m <= n; m *= 2) {
const h = m / 2;
const q = h / 2;
for (let j = 0; j < n; j += m) {
const even = dst[j];
const odd = dst[j + h];
dst[j] = new Complex(even).add(odd);
dst[j + h] = even.sub(odd);
}
if (h < 2)
continue;
let k = 1;
for (; k < q; k++) {
const w = W(inv * k / m);
for (let j = 0; j < n; j += m) {
const even = dst[j + k];
const odd = dst[j + k + h].mul(w);
dst[j + k] = new Complex(even).add(odd);
dst[j + k + h] = even.sub(odd);
}
}
for (let j = 0; j < n; j += m) {
const even = dst[j + q].re;
const odd = inv * dst[j + q + h].re;
dst[j + q] = new Complex(even, odd);
dst[j + q + h] = new Complex(even, -odd);
}
for (k++; k < h; k++) {
const w = W(inv * k / m);
for (let j = 0; j < n; j += m) {
const even = dst[j + k];
const odd = dst[j + k + h].mul(w);
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));
FFT(src).then(async fft => {
const dft = await DFT(src);
fft.forEach((z, i) => LOG(z.sub(dft[i]).abs().toFixed(8)));
STATUS("finished");
});
これで複素共役であることもはっきりする。
残ったループをひとつにまとめたいので、ループ変数の範囲を変更する。
このとき、同時に計算する相手を複素共役にしないと意味がないので、計算するデータはd[q+k]ではなく、d[h-k]である。
そこで、とりあえずkを機械的にh-kで置き換えてしまう。
"use strict";
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]);
idx.forEach((v, i) => dst[i] = new Complex(src[v]));
for (let m = 2; m <= n; m *= 2) {
const h = m / 2;
const q = h / 2;
for (let j = 0; j < n; j += m) {
const even = dst[j];
const odd = dst[j + h];
dst[j] = new Complex(even).add(odd);
dst[j + h] = even.sub(odd);
}
if (h < 2)
continue;
for (let k = 1; k < q; k++) {
const w = W(inv * k / m);
for (let j = 0; j < n; j += m) {
const even = dst[j + k];
const odd = dst[j + k + h].mul(w);
dst[j + k] = new Complex(even).add(odd);
dst[j + k + h] = even.sub(odd);
}
}
for (let j = 0; j < n; j += m) {
const even = dst[j + q].re;
const odd = inv * dst[j + q + h].re;
dst[j + q] = new Complex(even, odd);
dst[j + q + h] = new Complex(even, -odd);
}
for (let k = 1; k < q; k++) {
const w = W(inv * (h - k) / m);
for (let j = 0; j < n; j += m) {
const even = dst[j + (h - k)];
const odd = dst[j + (h - k) + h].mul(w);
dst[j + (h - k)] = new Complex(even).add(odd);
dst[j + (h - k) + h] = even.sub(odd);
}
}
}
await STATUS("FFT: done");
return dst;
}
const src = Array(16).fill(0).map(_ => new Complex(Math.random() - 0.5));
FFT(src).then(async fft => {
const dft = await DFT(src);
fft.forEach((z, i) => LOG(z.sub(dft[i]).abs().toFixed(8)));
STATUS("finished");
});
かっこを外して整理する。 回転因子も指数法則で展開しておく。
"use strict";
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]);
idx.forEach((v, i) => dst[i] = new Complex(src[v]));
for (let m = 2; m <= n; m *= 2) {
const h = m / 2;
const q = h / 2;
for (let j = 0; j < n; j += m) {
const even = dst[j];
const odd = dst[j + h];
dst[j] = new Complex(even).add(odd);
dst[j + h] = even.sub(odd);
}
if (h < 2)
continue;
for (let k = 1; k < q; k++) {
const w = W(inv * k / m);
for (let j = 0; j < n; j += m) {
const even = dst[j + k];
const odd = dst[j + k + h].mul(w);
dst[j + k] = new Complex(even).add(odd);
dst[j + k + h] = even.sub(odd);
}
}
for (let j = 0; j < n; j += m) {
const even = dst[j + q].re;
const odd = inv * dst[j + q + h].re;
dst[j + q] = new Complex(even, odd);
dst[j + q + h] = new Complex(even, -odd);
}
for (let k = 1; k < q; k++) {
const w = W(inv * h / m).mul(W(-inv * k / m));
for (let j = 0; j < n; j += m) {
const even = dst[j + h - k];
const odd = dst[j + m - k].mul(w);
dst[j + h - k] = new Complex(even).add(odd);
dst[j + m - k] = even.sub(odd);
}
}
}
await STATUS("FFT: done");
return dst;
}
const src = Array(16).fill(0).map(_ => new Complex(Math.random() - 0.5));
FFT(src).then(async fft => {
const dft = await DFT(src);
fft.forEach((z, i) => LOG(z.sub(dft[i]).abs().toFixed(8)));
STATUS("finished");
});
回転因子を約分すると、
"use strict";
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]);
idx.forEach((v, i) => dst[i] = new Complex(src[v]));
for (let m = 2; m <= n; m *= 2) {
const h = m / 2;
const q = h / 2;
for (let j = 0; j < n; j += m) {
const even = dst[j];
const odd = dst[j + h];
dst[j] = new Complex(even).add(odd);
dst[j + h] = even.sub(odd);
}
if (h < 2)
continue;
for (let k = 1; k < q; k++) {
const w = W(inv * k / m);
for (let j = 0; j < n; j += m) {
const even = dst[j + k];
const odd = dst[j + k + h].mul(w);
dst[j + k] = new Complex(even).add(odd);
dst[j + k + h] = even.sub(odd);
}
}
for (let j = 0; j < n; j += m) {
const even = dst[j + q].re;
const odd = inv * dst[j + q + h].re;
dst[j + q] = new Complex(even, odd);
dst[j + q + h] = new Complex(even, -odd);
}
for (let k = 1; k < q; k++) {
const w = W(inv / 2).mul(W(-inv * k / m));
for (let j = 0; j < n; j += m) {
const even = dst[j + h - k];
const odd = dst[j + m - k].mul(w);
dst[j + h - k] = new Complex(even).add(odd);
dst[j + m - k] = even.sub(odd);
}
}
}
await STATUS("FFT: done");
return dst;
}
const src = Array(16).fill(0).map(_ => new Complex(Math.random() - 0.5));
FFT(src).then(async fft => {
const dft = await DFT(src);
fft.forEach((z, i) => LOG(z.sub(dft[i]).abs().toFixed(8)));
STATUS("finished");
});
\(e^{\mp j\pi}\) は常に \(-1\) なので、
"use strict";
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]);
idx.forEach((v, i) => dst[i] = new Complex(src[v]));
for (let m = 2; m <= n; m *= 2) {
const h = m / 2;
const q = h / 2;
for (let j = 0; j < n; j += m) {
const even = dst[j];
const odd = dst[j + h];
dst[j] = new Complex(even).add(odd);
dst[j + h] = even.sub(odd);
}
if (h < 2)
continue;
for (let k = 1; k < q; k++) {
const w = W(inv * k / m);
for (let j = 0; j < n; j += m) {
const even = dst[j + k];
const odd = dst[j + k + h].mul(w);
dst[j + k] = new Complex(even).add(odd);
dst[j + k + h] = even.sub(odd);
}
}
for (let j = 0; j < n; j += m) {
const even = dst[j + q].re;
const odd = inv * dst[j + q + h].re;
dst[j + q] = new Complex(even, odd);
dst[j + q + h] = new Complex(even, -odd);
}
for (let k = 1; k < q; k++) {
const w = W(-inv * k / m);
for (let j = 0; j < n; j += m) {
const even = dst[j + h - k];
const odd = dst[j + m - k].mul(w);
dst[j + h - k] = new Complex(even).sub(odd);
dst[j + m - k] = even.add(odd);
}
}
}
await STATUS("FFT: done");
return dst;
}
const src = Array(16).fill(0).map(_ => new Complex(Math.random() - 0.5));
FFT(src).then(async fft => {
const dft = await DFT(src);
fft.forEach((z, i) => LOG(z.sub(dft[i]).abs().toFixed(8)));
STATUS("finished");
});
乗算を省略でき、加減算が入れ替わる。 ループをまとめるために変数名を変更しておいて、
"use strict";
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]);
idx.forEach((v, i) => dst[i] = new Complex(src[v]));
for (let m = 2; m <= n; m *= 2) {
const h = m / 2;
const q = h / 2;
for (let j = 0; j < n; j += m) {
const even = dst[j];
const odd = dst[j + h];
dst[j] = new Complex(even).add(odd);
dst[j + h] = even.sub(odd);
}
if (h < 2)
continue;
for (let k = 1; k < q; k++) {
const w0 = W(inv * k / m);
for (let j = 0; j < n; j += m) {
const z0 = dst[j + k];
const z1 = dst[j + k + h].mul(w0);
dst[j + k] = new Complex(z0).add(z1);
dst[j + k + h] = z0.sub(z1);
}
}
for (let j = 0; j < n; j += m) {
const even = dst[j + q].re;
const odd = inv * dst[j + q + h].re;
dst[j + q] = new Complex(even, odd);
dst[j + q + h] = new Complex(even, -odd);
}
for (let k = 1; k < q; k++) {
const w1 = W(-inv * k / m);
for (let j = 0; j < n; j += m) {
const z2 = dst[j + h - k];
const z3 = dst[j + m - k].mul(w1);
dst[j + h - k] = new Complex(z2).sub(z3);
dst[j + m - k] = z2.add(z3);
}
}
}
await STATUS("FFT: done");
return dst;
}
const src = Array(16).fill(0).map(_ => new Complex(Math.random() - 0.5));
FFT(src).then(async fft => {
const dft = await DFT(src);
fft.forEach((z, i) => LOG(z.sub(dft[i]).abs().toFixed(8)));
STATUS("finished");
});
えいやっとまとめる。
"use strict";
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]);
idx.forEach((v, i) => dst[i] = new Complex(src[v]));
for (let m = 2; m <= n; m *= 2) {
const h = m / 2;
const q = h / 2;
for (let j = 0; j < n; j += m) {
const even = dst[j];
const odd = dst[j + h];
dst[j] = new Complex(even).add(odd);
dst[j + h] = even.sub(odd);
}
if (h < 2)
continue;
for (let k = 1; k < q; k++) {
const w0 = W(inv * k / m);
const w1 = W(-inv * k / m);
for (let j = 0; j < n; j += m) {
const z0 = dst[j + k];
const z1 = dst[j + k + h].mul(w0);
const z2 = dst[j + h - k];
const z3 = dst[j + m - k].mul(w1);
dst[j + k] = new Complex(z0).add(z1);
dst[j + k + h] = z0.sub(z1);
dst[j + h - k] = new Complex(z2).sub(z3);
dst[j + m - k] = z2.add(z3);
}
}
for (let j = 0; j < n; j += m) {
const even = dst[j + q].re;
const odd = inv * dst[j + q + h].re;
dst[j + q] = new Complex(even, odd);
dst[j + q + h] = new Complex(even, -odd);
}
}
await STATUS("FFT: done");
return dst;
}
const src = Array(16).fill(0).map(_ => new Complex(Math.random() - 0.5));
FFT(src).then(async fft => {
const dft = await DFT(src);
fft.forEach((z, i) => LOG(z.sub(dft[i]).abs().toFixed(8)));
STATUS("finished");
});
回転因子は符号が逆なだけなので、w0とw1は複素共役だということが分かる。
"use strict";
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]);
idx.forEach((v, i) => dst[i] = new Complex(src[v]));
for (let m = 2; m <= n; m *= 2) {
const h = m / 2;
const q = h / 2;
for (let j = 0; j < n; j += m) {
const even = dst[j];
const odd = dst[j + h];
dst[j] = new Complex(even).add(odd);
dst[j + h] = even.sub(odd);
}
if (h < 2)
continue;
for (let k = 1; k < q; k++) {
const w0 = W(inv * k / m);
const w1 = new Complex(w0.re, -w0.im);
for (let j = 0; j < n; j += m) {
const z0 = dst[j + k];
const z1 = dst[j + k + h].mul(w0);
const z2 = dst[j + h - k];
const z3 = dst[j + m - k].mul(w1);
dst[j + k] = new Complex(z0).add(z1);
dst[j + k + h] = z0.sub(z1);
dst[j + h - k] = new Complex(z2).sub(z3);
dst[j + m - k] = z2.add(z3);
}
}
for (let j = 0; j < n; j += m) {
const even = dst[j + q].re;
const odd = inv * dst[j + q + h].re;
dst[j + q] = new Complex(even, odd);
dst[j + q + h] = new Complex(even, -odd);
}
}
await STATUS("FFT: done");
return dst;
}
const src = Array(16).fill(0).map(_ => new Complex(Math.random() - 0.5));
FFT(src).then(async fft => {
const dft = await DFT(src);
fft.forEach((z, i) => LOG(z.sub(dft[i]).abs().toFixed(8)));
STATUS("finished");
});
ここで、もし結果の後半を省略したらどうなるかを考える。 結果の後半を出力しないということは、今計算しているDFTの前段のちびっ子FFT達も結果の後半を出力してくれない、ということだ。 つまり、入力がそもそも半分しかない。 ここから全ての計算を行い、後半のデータを省略して出力することになる。 状況をはっきりさせるため、実際にデータの有無をコメントで書くとこうなる。
"use strict";
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]);
idx.forEach((v, i) => dst[i] = new Complex(src[v]));
for (let m = 2; m <= n; m *= 2) {
const h = m / 2;
const q = h / 2;
for (let j = 0; j < n; j += m) {
const even = dst[j];
const odd = dst[j + h];
dst[j] = new Complex(even).add(odd);
dst[j + h] = even.sub(odd);
}
if (h < 2)
continue;
for (let k = 1; k < q; k++) {
const w0 = W(inv * k / m);
const w1 = new Complex(w0.re, -w0.im);
for (let j = 0; j < n; j += m) {
const z0 = dst[j + k]; // exists
const z1 = dst[j + k + h].mul(w0); // exists
const z2 = dst[j + h - k]; // doesn't exist
const z3 = dst[j + m - k].mul(w1); // doesn't exist
dst[j + k] = new Complex(z0).add(z1); // needed
dst[j + k + h] = z0.sub(z1); // not needed
dst[j + h - k] = new Complex(z2).sub(z3); // needed
dst[j + m - k] = z2.add(z3); // not needed
}
}
for (let j = 0; j < n; j += m) {
const even = dst[j + q].re; // exists
const odd = inv * dst[j + q + h].re; // exists
dst[j + q] = new Complex(even, odd); // needed
dst[j + q + h] = new Complex(even, -odd); // not needed
}
}
await STATUS("FFT: done");
return dst;
}
const src = Array(16).fill(0).map(_ => new Complex(Math.random() - 0.5));
FFT(src).then(async fft => {
const dft = await DFT(src);
fft.forEach((z, i) => LOG(z.sub(dft[i]).abs().toFixed(8)));
STATUS("finished");
});
入力は4つに分けた真ん中より前の1/4と、最後の1/4、つまり、h-kとm-kと書かれている部分が存在せず、出力は後半、つまりh+kとm-kと書かれた部分を出力しなくてよい。
実際にえいやっと書き換えるとこうなる。
"use strict";
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;
}
conj() {
this.im = -this.im;
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]);
idx.forEach((v, i) => dst[i] = new Complex(src[v]));
for (let m = 2; m <= n; m *= 2) {
const h = m / 2;
const q = h / 2;
for (let j = 0; j < n; j += m) {
const even = dst[j];
const odd = dst[j + h];
dst[j] = new Complex(even).add(odd);
dst[j + h] = even.sub(odd);
}
if (h < 2)
continue;
for (let k = 1; k < q; k++) {
const w0 = W(inv * k / m);
// const w1 = new Complex(w0.re, -w0.im);
for (let j = 0; j < n; j += m) {
const z0 = dst[j + k];
const z1 = dst[j + k + h].mul(w0);
const z2 = new Complex(z0).conj(); // dst[j + h - k];
const z3 = new Complex(z1).conj(); // dst[j + m - k].mul(w1);
dst[j + k] = new Complex(z0).add(z1);
dst[j + k + h] = null;
dst[j + h - k] = new Complex(z2).sub(z3);
dst[j + m - k] = null;
}
}
for (let j = 0; j < n; j += m) {
const even = dst[j + q].re;
const odd = inv * dst[j + q + h].re;
dst[j + q] = new Complex(even, odd);
dst[j + q + h] = null;
}
}
await STATUS("FFT: done");
return dst;
}
function expandFFTResult(raw) {
const n = raw.length;
return raw.map((v, i) => v ?? new Complex(raw[n - i]).conj());
}
const src = Array(16).fill(0).map(_ => new Complex(Math.random() - 0.5));
FFT(src).then(async fft => {
const dft = await DFT(src);
fft = expandFFTResult(fft);
fft.forEach((z, i) => LOG(z.sub(dft[i]).abs().toFixed(8)));
STATUS("finished");
});
出力側は複素共役になるものをnullにしてしまい、入力側はデータが存在しない部分は存在する部分から複素共役の性質を使ってでっち上げる。
Complexクラスにconjメッソドが増えている。
これはconjugate、つまり複素共役の意味で、自身を複素共役の値に変更してthisを返す。
conjは自身を変更してしまうので、z0とz1は複素共役を取る前に複製している。
z3はdst[j + k + h]の共役複素数にw1、つまりw0の共役複素数を乗じたものだが、
dst[j + k + h] * w0の複素共役に等しい。
つまり、z1の複素共役になっている。
最終的に求めた結果も後半部分がない。
このままでは検証ができないので、expandFFTResultで後半部分を復元している。
今まで、入力はComplexを想定していたが、これを実数入力に変更する。
あとで分かりやすいように一時変数を導入してから、
"use strict";
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;
}
conj() {
this.im = -this.im;
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]);
idx.forEach((v, i) => dst[i] = new Complex(src[v]));
for (let m = 2; m <= n; m *= 2) {
const h = m / 2;
const q = h / 2;
for (let j = 0; j < n; j += m) {
const even = dst[j];
const odd = dst[j + h];
dst[j] = new Complex(even).add(odd);
dst[j + h] = even.sub(odd);
}
if (h < 2)
continue;
for (let k = 1; k < q; k++) {
const w0 = W(inv * k / m);
for (let j = 0; j < n; j += m) {
const z0 = dst[j + k];
const z1 = dst[j + k + h].mul(w0);
const z2 = new Complex(z0).conj(); // dst[j + h - k];
const z3 = new Complex(z1).conj(); // dst[j + m - k].mul(w1);
const y0 = z0.add(z1);
const y1 = z2.sub(z3);
dst[j + k] = y0;
dst[j + k + h] = null;
dst[j + h - k] = y1;
dst[j + m - k] = null;
}
}
for (let j = 0; j < n; j += m) {
const even = dst[j + q].re;
const odd = inv * dst[j + q + h].re;
const y = new Complex(even, odd);
dst[j + q] = y;
dst[j + q + h] = null;
}
}
await STATUS("FFT: done");
return dst;
}
function expandFFTResult(raw) {
const n = raw.length;
return raw.map((v, i) => v ?? new Complex(raw[n - i]).conj());
}
const src = Array(16).fill(0).map(_ => new Complex(Math.random() - 0.5));
FFT(src).then(async fft => {
const dft = await DFT(src);
fft = expandFFTResult(fft);
fft.forEach((z, i) => LOG(z.sub(dft[i]).abs().toFixed(8)));
STATUS("finished");
});
入力データを変更する。
"use strict";
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;
}
conj() {
this.im = -this.im;
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]);
idx.forEach((v, i) => dst[i] = new Complex(src[v]));
for (let m = 2; m <= n; m *= 2) {
const h = m / 2;
const q = h / 2;
for (let j = 0; j < n; j += m) {
const even = dst[j];
const odd = dst[j + h];
dst[j] = new Complex(even).add(odd);
dst[j + h] = even.sub(odd);
}
if (h < 2)
continue;
for (let k = 1; k < q; k++) {
const w0 = W(inv * k / m);
for (let j = 0; j < n; j += m) {
const z0 = dst[j + k];
const z1 = dst[j + k + h].mul(w0);
const z2 = new Complex(z0).conj(); // dst[j + h - k];
const z3 = new Complex(z1).conj(); // dst[j + m - k].mul(w1);
const y0 = z0.add(z1);
const y1 = z2.sub(z3);
dst[j + k] = y0;
dst[j + k + h] = null;
dst[j + h - k] = y1;
dst[j + m - k] = null;
}
}
for (let j = 0; j < n; j += m) {
const even = dst[j + q].re;
const odd = inv * dst[j + q + h].re;
const y = new Complex(even, odd);
dst[j + q] = y;
dst[j + q + h] = null;
}
}
await STATUS("FFT: done");
return dst;
}
function expandFFTResult(raw) {
const n = raw.length;
return raw.map((v, i) => v ?? new Complex(raw[n - i]).conj());
}
const src = Array(16).fill(0).map(_ => Math.random() - 0.5);
FFT(src).then(async fft => {
const dft = await DFT(src.map(v => new Complex(v)));
fft = expandFFTResult(fft);
fft.forEach((z, i) => LOG(z.sub(dft[i]).abs().toFixed(8)));
STATUS("finished");
});
DFT関数の入力は複素数なので、複素数に変換して渡す。
FFT関数では入力の入れ替え操作の時に複製が作成されているので、今のところ問題なく動く。
続けてdstを実数のまま使い続けるように変更する。
"use strict";
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;
}
conj() {
this.im = -this.im;
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]);
idx.forEach((v, i) => dst[i] = src[v]);
for (let m = 2; m <= n; m *= 2) {
const h = m / 2;
const q = h / 2;
for (let j = 0; j < n; j += m) {
const even = dst[j];
const odd = dst[j + h];
dst[j] = even + odd;
dst[j + h] = even - odd;
}
if (h < 2)
continue;
for (let k = 1; k < q; k++) {
const w0 = W(inv * k / m);
for (let j = 0; j < n; j += m) {
const z0 = new Complex(dst[j + k], dst[j + h - k]);
const z1 = new Complex(dst[j + k + h], dst[j + m - k]).mul(w0);
const z2 = new Complex(z0).conj(); // dst[j + h - k];
const z3 = new Complex(z1).conj(); // dst[j + m - k].mul(w1);
const y0 = z0.add(z1);
const y1 = z2.sub(z3);
dst[j + k] = y0.re;
dst[j + k + h] = y1.im;
dst[j + h - k] = y1.re;
dst[j + m - k] = y0.im;
}
}
for (let j = 0; j < n; j += m) {
const even = dst[j + q];
const odd = inv * dst[j + q + h];
const y = new Complex(even, odd);
dst[j + q] = y.re;
dst[j + q + h] = y.im;
}
}
await STATUS("FFT: done");
return dst;
}
function expandFFTResult(raw) {
const n = raw.length;
const out = Array(n);
out[0] = new Complex(raw[0]);
if (n < 2)
return out;
for (let i = 1; i < n / 2; i++) {
out[i] = new Complex(raw[i], raw[n - i]);
out[n - i] = new Complex(out[i]).conj();
}
out[n / 2] = new Complex(raw[n / 2]);
return out;
}
const src = Array(16).fill(0).map(_ => Math.random() - 0.5);
FFT(src).then(async fft => {
const dft = await DFT(src.map(v => new Complex(v)));
fft = expandFFTResult(fft);
fft.forEach((z, i) => LOG(z.sub(dft[i]).abs().toFixed(8)));
STATUS("finished");
});
今まで複素数を出力していたところに実部を、nullを出力していたところに虚部を出力する。
出力結果も仕様が変わるので、expandFFTResultで対応する。
out[k]が実部ならば、虚部はout[n-k]にある。
k[0]とk[n/2]は実数で、対応する虚部はない。
"use strict";
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;
}
conj() {
this.im = -this.im;
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]);
idx.forEach((v, i) => dst[i] = src[v]);
for (let m = 2; m <= n; m *= 2) {
const h = m / 2;
const q = h / 2;
for (let j = 0; j < n; j += m) {
const even = dst[j];
const odd = dst[j + h];
dst[j] = even + odd;
dst[j + h] = even - odd;
}
if (h < 2)
continue;
for (let k = 1; k < q; k++) {
const w0 = W(inv * k / m);
for (let j = 0; j < n; j += m) {
const z0 = new Complex(dst[j + k], dst[j + h - k]);
const z1 = new Complex(dst[j + k + h], dst[j + m - k]).mul(w0);
// const z2 = new Complex(z0).conj(); // dst[j + h - k];
// const z3 = new Complex(z1).conj(); // dst[j + m - k].mul(w1);
const y0 = new Complex(z0).add(z1);
const y1 = z0.sub(z1).conj();
dst[j + k] = y0.re;
dst[j + k + h] = y1.im;
dst[j + h - k] = y1.re;
dst[j + m - k] = y0.im;
}
}
for (let j = 0; j < n; j += m) {
const even = dst[j + q];
const odd = inv * dst[j + q + h];
const y = new Complex(even, odd);
dst[j + q] = y.re;
dst[j + q + h] = y.im;
}
}
await STATUS("FFT: done");
return dst;
}
function expandFFTResult(raw) {
const n = raw.length;
const out = Array(n);
out[0] = new Complex(raw[0]);
if (n < 2)
return out;
for (let i = 1; i < n / 2; i++) {
out[i] = new Complex(raw[i], raw[n - i]);
out[n - i] = new Complex(out[i]).conj();
}
out[n / 2] = new Complex(raw[n / 2]);
return out;
}
const src = Array(16).fill(0).map(_ => Math.random() - 0.5);
FFT(src).then(async fft => {
const dft = await DFT(src.map(v => new Complex(v)));
fft = expandFFTResult(fft);
fft.forEach((z, i) => LOG(z.sub(dft[i]).abs().toFixed(8)));
STATUS("finished");
});
z2とz3は不要。
回転因子も含め、複素演算を実数演算に展開する。 とりあえずデータを読み込む側。
"use strict";
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;
}
conj() {
this.im = -this.im;
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]);
idx.forEach((v, i) => dst[i] = src[v]);
for (let m = 2; m <= n; m *= 2) {
const h = m / 2;
const q = h / 2;
for (let j = 0; j < n; j += m) {
const even = dst[j];
const odd = dst[j + h];
dst[j] = even + odd;
dst[j + h] = even - odd;
}
if (h < 2)
continue;
for (let k = 1; k < q; k++) {
// const w0 = W(inv * k / m);
const wr = Math.cos(2 * Math.PI * inv * k / m);
const wi = Math.sin(2 * Math.PI * inv * k / m);
for (let j = 0; j < n; j += m) {
// const z0 = new Complex(dst[j + k], dst[j + h - k]);
const re0 = dst[j + k];
const im0 = dst[j + h - k];
// const z1 = new Complex(dst[j + k + h], dst[j + m - k]).mul(w0);
const re = dst[j + k + h];
const im = dst[j + m - k];
const re1 = re * wr - im * wi;
const im1 = re * wi + im * wr;
const y0 = new Complex(re0 + re1, im0 + im1);
const y1 = new Complex(re0 - re1, im0 - im1);
dst[j + k] = y0.re;
dst[j + k + h] = -y1.im;
dst[j + h - k] = y1.re;
dst[j + m - k] = y0.im;
}
}
for (let j = 0; j < n; j += m) {
const even = dst[j + q];
const odd = inv * dst[j + q + h];
const y = new Complex(even, odd);
dst[j + q] = y.re;
dst[j + q + h] = y.im;
}
}
await STATUS("FFT: done");
return dst;
}
function expandFFTResult(raw) {
const n = raw.length;
const out = Array(n);
out[0] = new Complex(raw[0]);
if (n < 2)
return out;
for (let i = 1; i < n / 2; i++) {
out[i] = new Complex(raw[i], raw[n - i]);
out[n - i] = new Complex(out[i]).conj();
}
out[n / 2] = new Complex(raw[n / 2]);
return out;
}
const src = Array(16).fill(0).map(_ => Math.random() - 0.5);
FFT(src).then(async fft => {
const dft = await DFT(src.map(v => new Complex(v)));
fft = expandFFTResult(fft);
fft.forEach((z, i) => LOG(z.sub(dft[i]).abs().toFixed(8)));
STATUS("finished");
});
y1は複素共役をdstの虚部への出力時に符号を操作することで行っていることに注意。
さらにデータを書き込む側も展開する。
ここで先ほど導入した一時変数y0・y1・yもお役御免になる。
"use strict";
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;
}
conj() {
this.im = -this.im;
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]);
idx.forEach((v, i) => dst[i] = src[v]);
for (let m = 2; m <= n; m *= 2) {
const h = m / 2;
const q = h / 2;
for (let j = 0; j < n; j += m) {
const even = dst[j];
const odd = dst[j + h];
dst[j] = even + odd;
dst[j + h] = even - odd;
}
if (h < 2)
continue;
for (let k = 1; k < q; k++) {
const wr = Math.cos(2 * Math.PI * inv * k / m);
const wi = Math.sin(2 * Math.PI * inv * k / m);
for (let j = 0; j < n; j += m) {
const re0 = dst[j + k];
const im0 = dst[j + h - k];
const re = dst[j + k + h];
const im = dst[j + m - k];
const re1 = re * wr - im * wi;
const im1 = re * wi + im * wr;
// const y0 = new Complex(re0 + re1, im0 + im1);
// const y1 = new Complex(re0 - re1, im0 - im1);
dst[j + k] = re0 + re1;
dst[j + k + h] = -(im0 - im1);
dst[j + h - k] = re0 - re1;
dst[j + m - k] = im0 + im1;
}
}
for (let j = 0; j < n; j += m) {
const even = dst[j + q];
const odd = inv * dst[j + q + h];
// const y = new Complex(even, odd);
dst[j + q] = even;
dst[j + q + h] = odd;
}
}
await STATUS("FFT: done");
return dst;
}
function expandFFTResult(raw) {
const n = raw.length;
const out = Array(n);
out[0] = new Complex(raw[0]);
if (n < 2)
return out;
for (let i = 1; i < n / 2; i++) {
out[i] = new Complex(raw[i], raw[n - i]);
out[n - i] = new Complex(out[i]).conj();
}
out[n / 2] = new Complex(raw[n / 2]);
return out;
}
const src = Array(16).fill(0).map(_ => Math.random() - 0.5);
FFT(src).then(async fft => {
const dft = await DFT(src.map(v => new Complex(v)));
fft = expandFFTResult(fft);
fft.forEach((z, i) => LOG(z.sub(dft[i]).abs().toFixed(8)));
STATUS("finished");
});
すると、最後のループのevenは読んで書いているだけなので処理を削除できる。
コメントも整理しておく。
"use strict";
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;
}
conj() {
this.im = -this.im;
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]);
idx.forEach((v, i) => dst[i] = src[v]);
for (let m = 2; m <= n; m *= 2) {
const h = m / 2;
const q = h / 2;
for (let j = 0; j < n; j += m) {
const even = dst[j];
const odd = dst[j + h];
dst[j] = even + odd;
dst[j + h] = even - odd;
}
if (h < 2)
continue;
for (let k = 1; k < q; k++) {
const wr = Math.cos(2 * Math.PI * inv * k / m);
const wi = Math.sin(2 * Math.PI * inv * k / m);
for (let j = 0; j < n; j += m) {
const re0 = dst[j + k];
const im0 = dst[j + h - k];
const re = dst[j + k + h];
const im = dst[j + m - k];
const re1 = re * wr - im * wi;
const im1 = re * wi + im * wr;
dst[j + k] = re0 + re1;
dst[j + k + h] = -(im0 - im1);
dst[j + h - k] = re0 - re1;
dst[j + m - k] = im0 + im1;
}
}
for (let j = 0; j < n; j += m) {
const odd = inv * dst[j + q + h];
dst[j + q + h] = odd;
}
}
await STATUS("FFT: done");
return dst;
}
function expandFFTResult(raw) {
const n = raw.length;
const out = Array(n);
out[0] = new Complex(raw[0]);
if (n < 2)
return out;
for (let i = 1; i < n / 2; i++) {
out[i] = new Complex(raw[i], raw[n - i]);
out[n - i] = new Complex(out[i]).conj();
}
out[n / 2] = new Complex(raw[n / 2]);
return out;
}
const src = Array(16).fill(0).map(_ => Math.random() - 0.5);
FFT(src).then(async fft => {
const dft = await DFT(src.map(v => new Complex(v)));
fft = expandFFTResult(fft);
fft.forEach((z, i) => LOG(z.sub(dft[i]).abs().toFixed(8)));
STATUS("finished");
});
この段階で、後処理や結果検証のためのDFTにはComplexクラスが残っているが、FFT関数自体はComplexクラスは不要になっている。
この状態でinvを1にすれば一応逆変換はできるが、入力、つまり周波数領域の値が実数に限るということなる。
フィルタのインパルス応答など、応用がないわけでもないが、この仕様ではあまり使い道がないので、思い切って正変換専用にする。
つまり、invを-1固定にする。
"use strict";
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;
}
conj() {
this.im = -this.im;
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) {
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]);
idx.forEach((v, i) => dst[i] = src[v]);
for (let m = 2; m <= n; m *= 2) {
const h = m / 2;
const q = h / 2;
for (let j = 0; j < n; j += m) {
const even = dst[j];
const odd = dst[j + h];
dst[j] = even + odd;
dst[j + h] = even - odd;
}
if (h < 2)
continue;
for (let k = 1; k < q; k++) {
const wr = Math.cos(2 * Math.PI * k / m);
const wi = -Math.sin(2 * Math.PI * k / m);
for (let j = 0; j < n; j += m) {
const re0 = dst[j + k];
const im0 = dst[j + h - k];
const re = dst[j + k + h];
const im = dst[j + m - k];
const re1 = re * wr - im * wi;
const im1 = re * wi + im * wr;
dst[j + k] = re0 + re1;
dst[j + k + h] = -(im0 - im1);
dst[j + h - k] = re0 - re1;
dst[j + m - k] = im0 + im1;
}
}
for (let j = 0; j < n; j += m) {
const odd = -dst[j + q + h];
dst[j + q + h] = odd;
}
}
await STATUS("FFT: done");
return dst;
}
function expandFFTResult(raw) {
const n = raw.length;
const out = Array(n);
out[0] = new Complex(raw[0]);
if (n < 2)
return out;
for (let i = 1; i < n / 2; i++) {
out[i] = new Complex(raw[i], raw[n - i]);
out[n - i] = new Complex(out[i]).conj();
}
out[n / 2] = new Complex(raw[n / 2]);
return out;
}
const src = Array(16).fill(0).map(_ => Math.random() - 0.5);
FFT(src).then(async fft => {
const dft = await DFT(src.map(v => new Complex(v)));
fft = expandFFTResult(fft);
fft.forEach((z, i) => LOG(z.sub(dft[i]).abs().toFixed(8)));
STATUS("finished");
});
回転因子の三角関数の中に現れる符号は整理して外に出してある。
最後のループに付いている負号をなくす方向で符号を調整する。
"use strict";
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;
}
conj() {
this.im = -this.im;
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) {
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]);
idx.forEach((v, i) => dst[i] = src[v]);
for (let m = 2; m <= n; m *= 2) {
const h = m / 2;
const q = h / 2;
for (let j = 0; j < n; j += m) {
const even = dst[j];
const odd = dst[j + h];
dst[j] = even + odd;
dst[j + h] = even - odd;
}
if (h < 2)
continue;
for (let k = 1; k < q; k++) {
const wr = Math.cos(2 * Math.PI * k / m);
const wi = -Math.sin(2 * Math.PI * k / m);
for (let j = 0; j < n; j += m) {
const re0 = dst[j + k];
const im0 = -dst[j + h - k];
const re = dst[j + k + h];
const im = -dst[j + m - k];
const re1 = re * wr - im * wi;
const im1 = re * wi + im * wr;
dst[j + k] = re0 + re1;
dst[j + k + h] = im0 - im1;
dst[j + h - k] = re0 - re1;
dst[j + m - k] = -(im0 + im1);
}
}
for (let j = 0; j < n; j += m) {
const odd = dst[j + q + h];
dst[j + q + h] = odd;
}
}
await STATUS("FFT: done");
return dst;
}
function expandFFTResult(raw) {
const n = raw.length;
const out = Array(n);
out[0] = new Complex(raw[0]);
if (n < 2)
return out;
for (let i = 1; i < n / 2; i++) {
out[i] = new Complex(raw[i], -raw[n - i]);
out[n - i] = new Complex(out[i]).conj();
}
out[n / 2] = new Complex(raw[n / 2]);
return out;
}
const src = Array(16).fill(0).map(_ => Math.random() - 0.5);
FFT(src).then(async fft => {
const dft = await DFT(src.map(v => new Complex(v)));
fft = expandFFTResult(fft);
fft.forEach((z, i) => LOG(z.sub(dft[i]).abs().toFixed(8)));
STATUS("finished");
});
これは虚部の符号を逆にして扱うことを意味するから、expandFFTResultでの出力の展開時も符号を調整する必要がある。
実部は前半、虚部は後半をそのまま出力した状態と考えるとこの状態になる。
すると、最後のループのoddも読んで書くだけになるので処理を削除でき、そうするとループの中が空になるので、なんと、ループを丸ごと削除できる。
"use strict";
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;
}
conj() {
this.im = -this.im;
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) {
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]);
idx.forEach((v, i) => dst[i] = src[v]);
for (let m = 2; m <= n; m *= 2) {
const h = m / 2;
const q = h / 2;
for (let j = 0; j < n; j += m) {
const even = dst[j];
const odd = dst[j + h];
dst[j] = even + odd;
dst[j + h] = even - odd;
}
if (h < 2)
continue;
for (let k = 1; k < q; k++) {
const wr = Math.cos(2 * Math.PI * k / m);
const wi = -Math.sin(2 * Math.PI * k / m);
for (let j = 0; j < n; j += m) {
const re0 = dst[j + k];
const im0 = -dst[j + h - k];
const re = dst[j + k + h];
const im = -dst[j + m - k];
const re1 = re * wr - im * wi;
const im1 = re * wi + im * wr;
dst[j + k] = re0 + re1;
dst[j + k + h] = im0 - im1;
dst[j + h - k] = re0 - re1;
dst[j + m - k] = -(im0 + im1);
}
}
}
await STATUS("FFT: done");
return dst;
}
function expandFFTResult(raw) {
const n = raw.length;
const out = Array(n);
out[0] = new Complex(raw[0]);
if (n < 2)
return out;
for (let i = 1; i < n / 2; i++) {
out[i] = new Complex(raw[i], -raw[n - i]);
out[n - i] = new Complex(out[i]).conj();
}
out[n / 2] = new Complex(raw[n / 2]);
return out;
}
const src = Array(16).fill(0).map(_ => Math.random() - 0.5);
FFT(src).then(async fft => {
const dft = await DFT(src.map(v => new Complex(v)));
fft = expandFFTResult(fft);
fft.forEach((z, i) => LOG(z.sub(dft[i]).abs().toFixed(8)));
STATUS("finished");
});
符号の調整もバカにできない。
このループをスキップするためのif文は分かりやすくqで書き換えると
"use strict";
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;
}
conj() {
this.im = -this.im;
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) {
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]);
idx.forEach((v, i) => dst[i] = src[v]);
for (let m = 2; m <= n; m *= 2) {
const h = m / 2;
const q = h / 2;
for (let j = 0; j < n; j += m) {
const even = dst[j];
const odd = dst[j + h];
dst[j] = even + odd;
dst[j + h] = even - odd;
}
if (q < 1)
continue;
for (let k = 1; k < q; k++) {
const wr = Math.cos(2 * Math.PI * k / m);
const wi = -Math.sin(2 * Math.PI * k / m);
for (let j = 0; j < n; j += m) {
const re0 = dst[j + k];
const im0 = -dst[j + h - k];
const re = dst[j + k + h];
const im = -dst[j + m - k];
const re1 = re * wr - im * wi;
const im1 = re * wi + im * wr;
dst[j + k] = re0 + re1;
dst[j + k + h] = im0 - im1;
dst[j + h - k] = re0 - re1;
dst[j + m - k] = -(im0 + im1);
}
}
}
await STATUS("FFT: done");
return dst;
}
function expandFFTResult(raw) {
const n = raw.length;
const out = Array(n);
out[0] = new Complex(raw[0]);
if (n < 2)
return out;
for (let i = 1; i < n / 2; i++) {
out[i] = new Complex(raw[i], -raw[n - i]);
out[n - i] = new Complex(out[i]).conj();
}
out[n / 2] = new Complex(raw[n / 2]);
return out;
}
const src = Array(16).fill(0).map(_ => Math.random() - 0.5);
FFT(src).then(async fft => {
const dft = await DFT(src.map(v => new Complex(v)));
fft = expandFFTResult(fft);
fft.forEach((z, i) => LOG(z.sub(dft[i]).abs().toFixed(8)));
STATUS("finished");
});
という条件になるが、この条件を満たす場合はすぐ下のループがそもそも実行されないことが分かる。
ということで、このif文も削除できる。
"use strict";
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;
}
conj() {
this.im = -this.im;
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) {
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]);
idx.forEach((v, i) => dst[i] = src[v]);
for (let m = 2; m <= n; m *= 2) {
const h = m / 2;
const q = h / 2;
for (let j = 0; j < n; j += m) {
const even = dst[j];
const odd = dst[j + h];
dst[j] = even + odd;
dst[j + h] = even - odd;
}
// if (q < 1)
// continue;
for (let k = 1; k < q; k++) {
const wr = Math.cos(2 * Math.PI * k / m);
const wi = -Math.sin(2 * Math.PI * k / m);
for (let j = 0; j < n; j += m) {
const re0 = dst[j + k];
const im0 = -dst[j + h - k];
const re = dst[j + k + h];
const im = -dst[j + m - k];
const re1 = re * wr - im * wi;
const im1 = re * wi + im * wr;
dst[j + k] = re0 + re1;
dst[j + k + h] = im0 - im1;
dst[j + h - k] = re0 - re1;
dst[j + m - k] = -(im0 + im1);
}
}
}
await STATUS("FFT: done");
return dst;
}
function expandFFTResult(raw) {
const n = raw.length;
const out = Array(n);
out[0] = new Complex(raw[0]);
if (n < 2)
return out;
for (let i = 1; i < n / 2; i++) {
out[i] = new Complex(raw[i], -raw[n - i]);
out[n - i] = new Complex(out[i]).conj();
}
out[n / 2] = new Complex(raw[n / 2]);
return out;
}
const src = Array(16).fill(0).map(_ => Math.random() - 0.5);
FFT(src).then(async fft => {
const dft = await DFT(src.map(v => new Complex(v)));
fft = expandFFTResult(fft);
fft.forEach((z, i) => LOG(z.sub(dft[i]).abs().toFixed(8)));
STATUS("finished");
});
コメントを削除して、残りの符号も調整し、自然な演算で計算できるようにする。
"use strict";
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;
}
conj() {
this.im = -this.im;
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) {
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]);
idx.forEach((v, i) => dst[i] = src[v]);
for (let m = 2; m <= n; m *= 2) {
const h = m / 2;
const q = h / 2;
for (let j = 0; j < n; j += m) {
const even = dst[j];
const odd = dst[j + h];
dst[j] = even + odd;
dst[j + h] = even - odd;
}
for (let k = 1; k < q; k++) {
const wr = Math.cos(2 * Math.PI * k / m);
const wi = Math.sin(2 * Math.PI * k / m);
for (let j = 0; j < n; j += m) {
const re0 = dst[j + k];
const im0 = dst[j + h - k];
const re = dst[j + k + h];
const im = dst[j + m - k];
const re1 = re * wr - im * wi;
const im1 = re * wi + im * wr;
dst[j + k] = re0 + re1;
dst[j + k + h] = im1 - im0;
dst[j + h - k] = re0 - re1;
dst[j + m - k] = im0 + im1;
}
}
}
await STATUS("FFT: done");
return dst;
}
function expandFFTResult(raw) {
const n = raw.length;
const out = Array(n);
out[0] = new Complex(raw[0]);
if (n < 2)
return out;
for (let i = 1; i < n / 2; i++) {
out[i] = new Complex(raw[i], -raw[n - i]);
out[n - i] = new Complex(out[i]).conj();
}
out[n / 2] = new Complex(raw[n / 2]);
return out;
}
const src = Array(16).fill(0).map(_ => Math.random() - 0.5);
FFT(src).then(async fft => {
const dft = await DFT(src.map(v => new Complex(v)));
fft = expandFFTResult(fft);
fft.forEach((z, i) => LOG(z.sub(dft[i]).abs().toFixed(8)));
STATUS("finished");
});
これが一般的な実FFTのコードになる。
あとはループを少し展開したり、基数を変更したり、k=m/8の場合は回転因子が
\((1+j)/\sqrt 2\)
になるのでもう少し複素乗算が減らせたり、回転因子の対称性からさらにループの長さを半分にして回転因子の計算回数（表引ならテーブル参照回数と容量）を半分にする、といった最適化が可能だが、この状態でもそこそこ速いので、今回はとりあえずここまでにして、その先へ進むことにする。
このコードでは関数が呼ばれてすぐにデータを添字のビット逆順で並び替えて（スクランブルして）dstに保持し、あとはdstの内容を複製することなく処理を進めている。
これを、入力された順序のまま計算し、出力の直前でスクランブルする（スクランブルを解除する）ことを考える。
これは、たとえば畳み込みの計算を実行するときなどによく使われるテクニックである。 畳み込み（高速巡回畳み込み）演算は入力をFFTし、あらかじめDFTしておいた畳み込み対象関数を乗じ（スペクトルごとに単純に乗算するだけだが、複素乗算であることには留意）、逆FFTすることで得られる。 もし、入力側でスクランブルしていると、逆FFTでは最後にスクランブルを解除することになる。 つまり、スクランブル→FFT→窓関数のFFT結果を乗算→逆FFT→スクランブル解除の手順になる。
このスクランブル位置を出力側に移すと、逆FFTでは先にスクランブルすることになるから、FFT→スクランブル解除→窓関数のFFT結果を乗算→スクランブル→逆FFTの順になる。 もし、窓関数の乗算をスクランブルしたまま実行すれば、スクランブル処理を丸々省くことができ、FFT→スクランブルしたまま窓関数のFFT結果を乗算→逆FFT、と実行できる。
まず、スクランブルを後ろへ持っていき、入力は複製だけしておく。
当然のように演算位置がずれてしまうから、とりあえずdstの添え字に全部idxによる添え字変換を書いてしまう。
"use strict";
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;
}
conj() {
this.im = -this.im;
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) {
await STATUS("FFT: doing", src.length);
const n = src.length;
let idx = [ 0 ];
while (idx.length < n)
idx = idx.flatMap(v => [ v, v + idx.length]);
const dst = Array.from(src);
for (let m = 2; m <= n; m *= 2) {
const h = m / 2;
const q = h / 2;
for (let j = 0; j < n; j += m) {
const even = dst[idx[j]];
const odd = dst[idx[j + h]];
dst[idx[j]] = even + odd;
dst[idx[j + h]] = even - odd;
}
for (let k = 1; k < q; k++) {
const wr = Math.cos(2 * Math.PI * k / m);
const wi = Math.sin(2 * Math.PI * k / m);
for (let j = 0; j < n; j += m) {
const re0 = dst[idx[j + k]];
const im0 = dst[idx[j + h - k]];
const re = dst[idx[j + k + h]];
const im = dst[idx[j + m - k]];
const re1 = re * wr - im * wi;
const im1 = re * wi + im * wr;
dst[idx[j + k]] = re0 + re1;
dst[idx[j + k + h]] = im1 - im0;
dst[idx[j + h - k]] = re0 - re1;
dst[idx[j + m - k]] = im0 + im1;
}
}
}
await STATUS("FFT: done");
return idx.map(i => dst[i]);
}
function expandFFTResult(raw) {
const n = raw.length;
const out = Array(n);
out[0] = new Complex(raw[0]);
if (n < 2)
return out;
for (let i = 1; i < n / 2; i++) {
out[i] = new Complex(raw[i], -raw[n - i]);
out[n - i] = new Complex(out[i]).conj();
}
out[n / 2] = new Complex(raw[n / 2]);
return out;
}
const src = Array(16).fill(0).map(_ => Math.random() - 0.5);
FFT(src).then(async fft => {
const dft = await DFT(src.map(v => new Complex(v)));
fft = expandFFTResult(fft);
fft.forEach((z, i) => LOG(z.sub(dft[i]).abs().toFixed(8)));
STATUS("finished");
});
元々のスクランブル処理はdstを作ってからforEachで一生懸命処理しているが、書き換え後の最後でやっているように、idxでmapすれば並び替えのついでに複製まで一発でやってくれることに気づいた。
ソースコードレビューって大事だね。
結果が正しいことを確認できたら、このidxの変換を削除していく。
まずjのループだが、元々のループでは0からmごとにnまで変化している。
これはビット逆順に直すと0からn/mまで動くことになり、元々のループの最中に最上位ビットが変化するから、ビット逆順では増分は1である。
なおかつ、回転因子のくくり出しで分かったように、jは演算内容には関係なく純粋に演算位置だけを指定している。
したがって、演算順序は気にする必要がなく、関係する場所をすべて列挙できれば問題ない。
よって、jはループ範囲と増分を書き換えることで簡単にidxの外に出せる。
また、m・h・qといった変数の加算は、たとえばmが2だったならビット逆順にすればn/4になるので、単純にn/2/mのように書き換えればidxの外に出せる。
試しに初段のループについて実施してみるとこうなる。
"use strict";
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;
}
conj() {
this.im = -this.im;
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) {
await STATUS("FFT: doing", src.length);
const n = src.length;
let idx = [ 0 ];
while (idx.length < n)
idx = idx.flatMap(v => [ v, v + idx.length]);
const dst = Array.from(src);
for (let m = 2; m <= n; m *= 2) {
const h = m / 2;
const q = h / 2;
for (let j = 0; j < n/m; j++) {
const even = dst[j];
const odd = dst[j + n/2/h];
dst[j] = even + odd;
dst[j + n/2/h] = even - odd;
}
for (let k = 1; k < q; k++) {
const wr = Math.cos(2 * Math.PI * k / m);
const wi = Math.sin(2 * Math.PI * k / m);
for (let j = 0; j < n; j += m) {
const re0 = dst[idx[j + k]];
const im0 = dst[idx[j + h - k]];
const re = dst[idx[j + k + h]];
const im = dst[idx[j + m - k]];
const re1 = re * wr - im * wi;
const im1 = re * wi + im * wr;
dst[idx[j + k]] = re0 + re1;
dst[idx[j + k + h]] = im1 - im0;
dst[idx[j + h - k]] = re0 - re1;
dst[idx[j + m - k]] = im0 + im1;
}
}
}
await STATUS("FFT: done");
return idx.map(i => dst[i]);
}
function expandFFTResult(raw) {
const n = raw.length;
const out = Array(n);
out[0] = new Complex(raw[0]);
if (n < 2)
return out;
for (let i = 1; i < n / 2; i++) {
out[i] = new Complex(raw[i], -raw[n - i]);
out[n - i] = new Complex(out[i]).conj();
}
out[n / 2] = new Complex(raw[n / 2]);
return out;
}
const src = Array(16).fill(0).map(_ => Math.random() - 0.5);
FFT(src).then(async fft => {
const dft = await DFT(src.map(v => new Complex(v)));
fft = expandFFTResult(fft);
fft.forEach((z, i) => LOG(z.sub(dft[i]).abs().toFixed(8)));
STATUS("finished");
});
大丈夫そうなので、後半のループについても同じように書き換える。
"use strict";
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;
}
conj() {
this.im = -this.im;
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) {
await STATUS("FFT: doing", src.length);
const n = src.length;
let idx = [ 0 ];
while (idx.length < n)
idx = idx.flatMap(v => [ v, v + idx.length]);
const dst = Array.from(src);
for (let m = 2; m <= n; m *= 2) {
const h = m / 2;
const q = h / 2;
for (let j = 0; j < n/m; j++) {
const even = dst[j];
const odd = dst[j + n/2/h];
dst[j] = even + odd;
dst[j + n/2/h] = even - odd;
}
for (let k = 1; k < q; k++) {
const wr = Math.cos(2 * Math.PI * k / m);
const wi = Math.sin(2 * Math.PI * k / m);
for (let j = 0; j < n/m; j++) {
const re0 = dst[j + idx[k]];
const im0 = dst[j + idx[h - k]];
const re = dst[j + idx[k + h]];
const im = dst[j + idx[m - k]];
const re1 = re * wr - im * wi;
const im1 = re * wi + im * wr;
dst[j + idx[k]] = re0 + re1;
dst[j + idx[k + h]] = im1 - im0;
dst[j + idx[h - k]] = re0 - re1;
dst[j + idx[m - k]] = im0 + im1;
}
}
}
await STATUS("FFT: done");
return idx.map(i => dst[i]);
}
function expandFFTResult(raw) {
const n = raw.length;
const out = Array(n);
out[0] = new Complex(raw[0]);
if (n < 2)
return out;
for (let i = 1; i < n / 2; i++) {
out[i] = new Complex(raw[i], -raw[n - i]);
out[n - i] = new Complex(out[i]).conj();
}
out[n / 2] = new Complex(raw[n / 2]);
return out;
}
const src = Array(16).fill(0).map(_ => Math.random() - 0.5);
FFT(src).then(async fft => {
const dft = await DFT(src.map(v => new Complex(v)));
fft = expandFFTResult(fft);
fft.forEach((z, i) => LOG(z.sub(dft[i]).abs().toFixed(8)));
STATUS("finished");
});
kはjとは逆に1ずつ増加しているので、ビット逆順にすると最上位ビットが変化することになり、ループ上限がnになる。
元々の最大値はqだったから、増分はn/qである。
これでq回繰り返せばnになる。
元々の初期値は0ではなく1なので、ビット逆順にした後も0は除外する必要がある。
とりあえず、大文字のKでループを書き換え、それをビット逆順にしたものを新たにkとして問題がないか確認してみる。
"use strict";
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;
}
conj() {
this.im = -this.im;
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) {
await STATUS("FFT: doing", src.length);
const n = src.length;
let idx = [ 0 ];
while (idx.length < n)
idx = idx.flatMap(v => [ v, v + idx.length]);
const dst = Array.from(src);
for (let m = 2; m <= n; m *= 2) {
const h = m / 2;
const q = h / 2;
for (let j = 0; j < n/m; j++) {
const even = dst[j];
const odd = dst[j + n/2/h];
dst[j] = even + odd;
dst[j + n/2/h] = even - odd;
}
for (let K = n/q; K < n; K += n/q) {
const k = idx[K];
const wr = Math.cos(2 * Math.PI * k / m);
const wi = Math.sin(2 * Math.PI * k / m);
for (let j = 0; j < n/m; j++) {
const re0 = dst[j + idx[k]];
const im0 = dst[j + idx[h - k]];
const re = dst[j + idx[k + h]];
const im = dst[j + idx[m - k]];
const re1 = re * wr - im * wi;
const im1 = re * wi + im * wr;
dst[j + idx[k]] = re0 + re1;
dst[j + idx[k + h]] = im1 - im0;
dst[j + idx[h - k]] = re0 - re1;
dst[j + idx[m - k]] = im0 + im1;
}
}
}
await STATUS("FFT: done");
return idx.map(i => dst[i]);
}
function expandFFTResult(raw) {
const n = raw.length;
const out = Array(n);
out[0] = new Complex(raw[0]);
if (n < 2)
return out;
for (let i = 1; i < n / 2; i++) {
out[i] = new Complex(raw[i], -raw[n - i]);
out[n - i] = new Complex(out[i]).conj();
}
out[n / 2] = new Complex(raw[n / 2]);
return out;
}
const src = Array(16).fill(0).map(_ => Math.random() - 0.5);
FFT(src).then(async fft => {
const dft = await DFT(src.map(v => new Complex(v)));
fft = expandFFTResult(fft);
fft.forEach((z, i) => LOG(z.sub(dft[i]).abs().toFixed(8)));
STATUS("finished");
});
続いてループ内部でjと同様に書き換えられる部分を書き換える。
"use strict";
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;
}
conj() {
this.im = -this.im;
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) {
await STATUS("FFT: doing", src.length);
const n = src.length;
let idx = [ 0 ];
while (idx.length < n)
idx = idx.flatMap(v => [ v, v + idx.length]);
const dst = Array.from(src);
for (let m = 2; m <= n; m *= 2) {
const h = m / 2;
const q = h / 2;
for (let j = 0; j < n/m; j++) {
const even = dst[j];
const odd = dst[j + n/2/h];
dst[j] = even + odd;
dst[j + n/2/h] = even - odd;
}
for (let K = n/q; K < n; K += n/q) {
const k = idx[K];
const wr = Math.cos(2 * Math.PI * k / m);
const wi = Math.sin(2 * Math.PI * k / m);
for (let j = 0; j < n/m; j++) {
const re0 = dst[j + K];
const im0 = dst[j + idx[h - k]];
const re = dst[j + K + n/2/h];
const im = dst[j + idx[m - k]];
const re1 = re * wr - im * wi;
const im1 = re * wi + im * wr;
dst[j + K] = re0 + re1;
dst[j + K + n/2/h] = im1 - im0;
dst[j + idx[h - k]] = re0 - re1;
dst[j + idx[m - k]] = im0 + im1;
}
}
}
await STATUS("FFT: done");
return idx.map(i => dst[i]);
}
function expandFFTResult(raw) {
const n = raw.length;
const out = Array(n);
out[0] = new Complex(raw[0]);
if (n < 2)
return out;
for (let i = 1; i < n / 2; i++) {
out[i] = new Complex(raw[i], -raw[n - i]);
out[n - i] = new Complex(out[i]).conj();
}
out[n / 2] = new Complex(raw[n / 2]);
return out;
}
const src = Array(16).fill(0).map(_ => Math.random() - 0.5);
FFT(src).then(async fft => {
const dft = await DFT(src.map(v => new Complex(v)));
fft = expandFFTResult(fft);
fft.forEach((z, i) => LOG(z.sub(dft[i]).abs().toFixed(8)));
STATUS("finished");
});
問題はここからである。
実は、ビット逆順の補数の計算は、1の補数は簡単にできるが、2の補数は簡単ではない。
1の補数は全ビットの1と0を反転することになるので、ビット逆順にしても1と0の反転である。
これはn-k-1のような演算に相当し、DCT4のようにオフセットのついた演算ではうまくいくことが多い。
n-kのような演算は2の補数であり、これは簡単には求められない。
データの並び替えをしていなければ1の補数と2の補数はひとつずれているだけなので、ループをまたいでインデックスを渡せばなんとかなるが、並び替えてしまっているのでそれもできない。 Split-Radix FFTのように二重ループにすれば求められそうな気がするが、定かではない。
そこで、とりあえずiという変数を導入して、h-kはidxを引いて求めてしまう。
iはkの近くで空いていたのと、imaginary、つまり虚数にも通じるのでこの名前にしてみた。
これさえ求めればn-kはh+h-kなので、idx[n-k]はidx[h+h-k] → idx[h]+i → n/2/h+iと計算できる。
コード中ではちょっと変数の順序が違うけど。
"use strict";
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;
}
conj() {
this.im = -this.im;
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) {
await STATUS("FFT: doing", src.length);
const n = src.length;
let idx = [ 0 ];
while (idx.length < n)
idx = idx.flatMap(v => [ v, v + idx.length]);
const dst = Array.from(src);
for (let m = 2; m <= n; m *= 2) {
const h = m / 2;
const q = h / 2;
for (let j = 0; j < n/m; j++) {
const even = dst[j];
const odd = dst[j + n/2/h];
dst[j] = even + odd;
dst[j + n/2/h] = even - odd;
}
for (let K = n/q; K < n; K += n/q) {
const k = idx[K];
const i = idx[h - k];
const wr = Math.cos(2 * Math.PI * k / m);
const wi = Math.sin(2 * Math.PI * k / m);
for (let j = 0; j < n/m; j++) {
const re0 = dst[j + K];
const im0 = dst[j + i];
const re = dst[j + K + n/2/h];
const im = dst[j + i + n/2/h];
const re1 = re * wr - im * wi;
const im1 = re * wi + im * wr;
dst[j + K] = re0 + re1;
dst[j + K + n/2/h] = im1 - im0;
dst[j + i] = re0 - re1;
dst[j + i + n/2/h] = im0 + im1;
}
}
}
await STATUS("FFT: done");
return idx.map(i => dst[i]);
}
function expandFFTResult(raw) {
const n = raw.length;
const out = Array(n);
out[0] = new Complex(raw[0]);
if (n < 2)
return out;
for (let i = 1; i < n / 2; i++) {
out[i] = new Complex(raw[i], -raw[n - i]);
out[n - i] = new Complex(out[i]).conj();
}
out[n / 2] = new Complex(raw[n / 2]);
return out;
}
const src = Array(16).fill(0).map(_ => Math.random() - 0.5);
FFT(src).then(async fft => {
const dft = await DFT(src.map(v => new Complex(v)));
fft = expandFFTResult(fft);
fft.forEach((z, i) => LOG(z.sub(dft[i]).abs().toFixed(8)));
STATUS("finished");
});
iと回転因子の計算に残っているkをidx[K]で書き換える。
iはidx[h - idx[K]]と二重にidxを引いてる点に注意。
"use strict";
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;
}
conj() {
this.im = -this.im;
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) {
await STATUS("FFT: doing", src.length);
const n = src.length;
let idx = [ 0 ];
while (idx.length < n)
idx = idx.flatMap(v => [ v, v + idx.length]);
const dst = Array.from(src);
for (let m = 2; m <= n; m *= 2) {
const h = m / 2;
const q = h / 2;
for (let j = 0; j < n/m; j++) {
const even = dst[j];
const odd = dst[j + n/2/h];
dst[j] = even + odd;
dst[j + n/2/h] = even - odd;
}
for (let K = n/q; K < n; K += n/q) {
// const k = idx[K];
const i = idx[h - idx[K]];
const wr = Math.cos(2 * Math.PI * idx[K] / m);
const wi = Math.sin(2 * Math.PI * idx[K] / m);
for (let j = 0; j < n/m; j++) {
const re0 = dst[j + K];
const im0 = dst[j + i];
const re = dst[j + K + n/2/h];
const im = dst[j + i + n/2/h];
const re1 = re * wr - im * wi;
const im1 = re * wi + im * wr;
dst[j + K] = re0 + re1;
dst[j + K + n/2/h] = im1 - im0;
dst[j + i] = re0 - re1;
dst[j + i + n/2/h] = im0 + im1;
}
}
}
await STATUS("FFT: done");
return idx.map(i => dst[i]);
}
function expandFFTResult(raw) {
const n = raw.length;
const out = Array(n);
out[0] = new Complex(raw[0]);
if (n < 2)
return out;
for (let i = 1; i < n / 2; i++) {
out[i] = new Complex(raw[i], -raw[n - i]);
out[n - i] = new Complex(out[i]).conj();
}
out[n / 2] = new Complex(raw[n / 2]);
return out;
}
const src = Array(16).fill(0).map(_ => Math.random() - 0.5);
FFT(src).then(async fft => {
const dft = await DFT(src.map(v => new Complex(v)));
fft = expandFFTResult(fft);
fft.forEach((z, i) => LOG(z.sub(dft[i]).abs().toFixed(8)));
STATUS("finished");
});
これで小文字のkがなくなるので、大文字のKを小文字のkに書き直す。
"use strict";
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;
}
conj() {
this.im = -this.im;
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) {
await STATUS("FFT: doing", src.length);
const n = src.length;
let idx = [ 0 ];
while (idx.length < n)
idx = idx.flatMap(v => [ v, v + idx.length]);
const dst = Array.from(src);
for (let m = 2; m <= n; m *= 2) {
const h = m / 2;
const q = h / 2;
for (let j = 0; j < n/m; j++) {
const even = dst[j];
const odd = dst[j + n/2/h];
dst[j] = even + odd;
dst[j + n/2/h] = even - odd;
}
for (let k = n/q; k < n; k += n/q) {
const i = idx[h - idx[k]];
const wr = Math.cos(2 * Math.PI * idx[k] / m);
const wi = Math.sin(2 * Math.PI * idx[k] / m);
for (let j = 0; j < n/m; j++) {
const re0 = dst[j + k];
const im0 = dst[j + i];
const re = dst[j + k + n/2/h];
const im = dst[j + i + n/2/h];
const re1 = re * wr - im * wi;
const im1 = re * wi + im * wr;
dst[j + k] = re0 + re1;
dst[j + k + n/2/h] = im1 - im0;
dst[j + i] = re0 - re1;
dst[j + i + n/2/h] = im0 + im1;
}
}
}
await STATUS("FFT: done");
return idx.map(i => dst[i]);
}
function expandFFTResult(raw) {
const n = raw.length;
const out = Array(n);
out[0] = new Complex(raw[0]);
if (n < 2)
return out;
for (let i = 1; i < n / 2; i++) {
out[i] = new Complex(raw[i], -raw[n - i]);
out[n - i] = new Complex(out[i]).conj();
}
out[n / 2] = new Complex(raw[n / 2]);
return out;
}
const src = Array(16).fill(0).map(_ => Math.random() - 0.5);
FFT(src).then(async fft => {
const dft = await DFT(src.map(v => new Complex(v)));
fft = expandFFTResult(fft);
fft.forEach((z, i) => LOG(z.sub(dft[i]).abs().toFixed(8)));
STATUS("finished");
});
idxの表引きは残ってしまっているわけだが、実用上はどうするかというと、この変換はDFT点数が決まれば決まってしまうものなので、回転因子も含めて全部テーブルにしてしまう。
三角関数も順序よく計算できるわけではないので加法定理のような公式も使えない。
テーブルにするのが一番楽である。
続いてhとqをmで置き換える。
"use strict";
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;
}
conj() {
this.im = -this.im;
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) {
await STATUS("FFT: doing", src.length);
const n = src.length;
let idx = [ 0 ];
while (idx.length < n)
idx = idx.flatMap(v => [ v, v + idx.length]);
const dst = Array.from(src);
for (let m = 2; m <= n; m *= 2) {
// const h = m / 2;
// const q = h / 2;
for (let j = 0; j < n/m; j++) {
const even = dst[j];
const odd = dst[j + n/m];
dst[j] = even + odd;
dst[j + n/m] = even - odd;
}
for (let k = 4*n/m; k < n; k += 4*n/m) {
const i = idx[m/2 - idx[k]];
const wr = Math.cos(2 * Math.PI * idx[k] / m);
const wi = Math.sin(2 * Math.PI * idx[k] / m);
for (let j = 0; j < n/m; j++) {
const re0 = dst[j + k];
const im0 = dst[j + i];
const re = dst[j + k + n/m];
const im = dst[j + i + n/m];
const re1 = re * wr - im * wi;
const im1 = re * wi + im * wr;
dst[j + k] = re0 + re1;
dst[j + k + n/m] = im1 - im0;
dst[j + i] = re0 - re1;
dst[j + i + n/m] = im0 + im1;
}
}
}
await STATUS("FFT: done");
return idx.map(i => dst[i]);
}
function expandFFTResult(raw) {
const n = raw.length;
const out = Array(n);
out[0] = new Complex(raw[0]);
if (n < 2)
return out;
for (let i = 1; i < n / 2; i++) {
out[i] = new Complex(raw[i], -raw[n - i]);
out[n - i] = new Complex(out[i]).conj();
}
out[n / 2] = new Complex(raw[n / 2]);
return out;
}
const src = Array(16).fill(0).map(_ => Math.random() - 0.5);
FFT(src).then(async fft => {
const dft = await DFT(src.map(v => new Complex(v)));
fft = expandFFTResult(fft);
fft.forEach((z, i) => LOG(z.sub(dft[i]).abs().toFixed(8)));
STATUS("finished");
});
すると、やたらとn/mという項が目につくようになる。
m・h・qはなんとなく逆数っぽいものになりそうである。
とりあえず、n/mをMと置いて書き換える。
"use strict";
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;
}
conj() {
this.im = -this.im;
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) {
await STATUS("FFT: doing", src.length);
const n = src.length;
let idx = [ 0 ];
while (idx.length < n)
idx = idx.flatMap(v => [ v, v + idx.length]);
const dst = Array.from(src);
for (let m = 2; m <= n; m *= 2) {
let M = n / m;
for (let j = 0; j < M; j++) {
const even = dst[j];
const odd = dst[j + M];
dst[j] = even + odd;
dst[j + M] = even - odd;
}
for (let k = 4*M; k < n; k += 4*M) {
const i = idx[n / (2*M) - idx[k]];
const wr = Math.cos(2 * Math.PI * idx[k] * M / n);
const wi = Math.sin(2 * Math.PI * idx[k] * M / n);
for (let j = 0; j < M; j++) {
const re0 = dst[j + k];
const im0 = dst[j + i];
const re = dst[j + k + M];
const im = dst[j + i + M];
const re1 = re * wr - im * wi;
const im1 = re * wi + im * wr;
dst[j + k] = re0 + re1;
dst[j + k + M] = im1 - im0;
dst[j + i] = re0 - re1;
dst[j + i + M] = im0 + im1;
}
}
}
await STATUS("FFT: done");
return idx.map(i => dst[i]);
}
function expandFFTResult(raw) {
const n = raw.length;
const out = Array(n);
out[0] = new Complex(raw[0]);
if (n < 2)
return out;
for (let i = 1; i < n / 2; i++) {
out[i] = new Complex(raw[i], -raw[n - i]);
out[n - i] = new Complex(out[i]).conj();
}
out[n / 2] = new Complex(raw[n / 2]);
return out;
}
const src = Array(16).fill(0).map(_ => Math.random() - 0.5);
FFT(src).then(async fft => {
const dft = await DFT(src.map(v => new Complex(v)));
fft = expandFFTResult(fft);
fft.forEach((z, i) => LOG(z.sub(dft[i]).abs().toFixed(8)));
STATUS("finished");
});
iと回転因子を求めている部分にもmが出てくるので、これも置き換えてある。
m=n/Mである。
これでループの中からmがなくなるので、ループもMで書き換える。
このループはm*M=n、つまりM=n/mとなるように書き換えればいいので、初期値がn/2、最終値はn/nだから1で、ループごとに2倍になっていたものはループごとに半分になる。
減少方向になるのでループ終了条件の不等号を書き換える必要がある（よく忘れる）。
"use strict";
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;
}
conj() {
this.im = -this.im;
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) {
await STATUS("FFT: doing", src.length);
const n = src.length;
let idx = [ 0 ];
while (idx.length < n)
idx = idx.flatMap(v => [ v, v + idx.length]);
const dst = Array.from(src);
for (let M = n / 2; M >= 1; M /= 2) {
for (let j = 0; j < M; j++) {
const even = dst[j];
const odd = dst[j + M];
dst[j] = even + odd;
dst[j + M] = even - odd;
}
for (let k = 4*M; k < n; k += 4*M) {
const i = idx[n / (2*M) - idx[k]];
const wr = Math.cos(2 * Math.PI * idx[k] * M / n);
const wi = Math.sin(2 * Math.PI * idx[k] * M / n);
for (let j = 0; j < M; j++) {
const re0 = dst[j + k];
const im0 = dst[j + i];
const re = dst[j + k + M];
const im = dst[j + i + M];
const re1 = re * wr - im * wi;
const im1 = re * wi + im * wr;
dst[j + k] = re0 + re1;
dst[j + k + M] = im1 - im0;
dst[j + i] = re0 - re1;
dst[j + i + M] = im0 + im1;
}
}
}
await STATUS("FFT: done");
return idx.map(i => dst[i]);
}
function expandFFTResult(raw) {
const n = raw.length;
const out = Array(n);
out[0] = new Complex(raw[0]);
if (n < 2)
return out;
for (let i = 1; i < n / 2; i++) {
out[i] = new Complex(raw[i], -raw[n - i]);
out[n - i] = new Complex(out[i]).conj();
}
out[n / 2] = new Complex(raw[n / 2]);
return out;
}
const src = Array(16).fill(0).map(_ => Math.random() - 0.5);
FFT(src).then(async fft => {
const dft = await DFT(src.map(v => new Complex(v)));
fft = expandFFTResult(fft);
fft.forEach((z, i) => LOG(z.sub(dft[i]).abs().toFixed(8)));
STATUS("finished");
});
2*M・4*Mという項も見えるので、D（double）とQ（quadruple）で置き換える。
"use strict";
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;
}
conj() {
this.im = -this.im;
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) {
await STATUS("FFT: doing", src.length);
const n = src.length;
let idx = [ 0 ];
while (idx.length < n)
idx = idx.flatMap(v => [ v, v + idx.length]);
const dst = Array.from(src);
for (let M = n / 2; M >= 1; M /= 2) {
const D = 2 * M;
const Q = 4 * M;
for (let j = 0; j < M; j++) {
const even = dst[j];
const odd = dst[j + M];
dst[j] = even + odd;
dst[j + M] = even - odd;
}
for (let k = Q; k < n; k += Q) {
const i = idx[n / D - idx[k]];
const wr = Math.cos(2 * Math.PI * idx[k] * M / n);
const wi = Math.sin(2 * Math.PI * idx[k] * M / n);
for (let j = 0; j < M; j++) {
const re0 = dst[j + k];
const im0 = dst[j + i];
const re = dst[j + k + M];
const im = dst[j + i + M];
const re1 = re * wr - im * wi;
const im1 = re * wi + im * wr;
dst[j + k] = re0 + re1;
dst[j + k + M] = im1 - im0;
dst[j + i] = re0 - re1;
dst[j + i + M] = im0 + im1;
}
}
}
await STATUS("FFT: done");
return idx.map(i => dst[i]);
}
function expandFFTResult(raw) {
const n = raw.length;
const out = Array(n);
out[0] = new Complex(raw[0]);
if (n < 2)
return out;
for (let i = 1; i < n / 2; i++) {
out[i] = new Complex(raw[i], -raw[n - i]);
out[n - i] = new Complex(out[i]).conj();
}
out[n / 2] = new Complex(raw[n / 2]);
return out;
}
const src = Array(16).fill(0).map(_ => Math.random() - 0.5);
FFT(src).then(async fft => {
const dft = await DFT(src.map(v => new Complex(v)));
fft = expandFFTResult(fft);
fft.forEach((z, i) => LOG(z.sub(dft[i]).abs().toFixed(8)));
STATUS("finished");
});
これで大体完成である。
FFT関数の冒頭でsrcをdstにコピーしているが、呼び出し元がsrcの内容を必要としない場合、これは単なる無駄である。
必要なら呼び出し側でコピーしてもらえばよいではないか。
また、FFT関数の最後でスクランブルを解除して返しているが、元々の趣旨を考えたら、スクランブルを解除しないで返した方がいい、ということになる。
そこでこれらをえいやっと削除したのがこれ。
"use strict";
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;
}
conj() {
this.im = -this.im;
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(dst) {
await STATUS("FFT: doing", src.length);
const n = src.length;
let idx = [ 0 ];
while (idx.length < n)
idx = idx.flatMap(v => [ v, v + idx.length]);
for (let M = n / 2; M >= 1; M /= 2) {
const D = 2 * M;
const Q = 4 * M;
for (let j = 0; j < M; j++) {
const even = dst[j];
const odd = dst[j + M];
dst[j] = even + odd;
dst[j + M] = even - odd;
}
for (let k = Q; k < n; k += Q) {
const i = idx[n / D - idx[k]];
const wr = Math.cos(2 * Math.PI * idx[k] * M / n);
const wi = Math.sin(2 * Math.PI * idx[k] * M / n);
for (let j = 0; j < M; j++) {
const re0 = dst[j + k];
const im0 = dst[j + i];
const re = dst[j + k + M];
const im = dst[j + i + M];
const re1 = re * wr - im * wi;
const im1 = re * wi + im * wr;
dst[j + k] = re0 + re1;
dst[j + k + M] = im1 - im0;
dst[j + i] = re0 - re1;
dst[j + i + M] = im0 + im1;
}
}
}
await STATUS("FFT: done");
return idx;
}
function expandFFTResult(raw) {
const n = raw.length;
const out = Array(n);
out[0] = new Complex(raw[0]);
if (n < 2)
return out;
for (let i = 1; i < n / 2; i++) {
out[i] = new Complex(raw[i], -raw[n - i]);
out[n - i] = new Complex(out[i]).conj();
}
out[n / 2] = new Complex(raw[n / 2]);
return out;
}
const src = Array(16).fill(0).map(_ => Math.random() - 0.5);
let fft = Array.from(src);
FFT(fft).then(async idx => {
const dft = await DFT(src.map(v => new Complex(v)));
fft = expandFFTResult(idx.map(i => fft[i]));
fft.forEach((z, i) => LOG(z.sub(dft[i]).abs().toFixed(8)));
STATUS("finished");
});
ビット逆順テーブルさえなんとかすれば、すぐにでもC言語に移植できそうな姿をしている。 これはインプレース演算（in-place、その場で）と呼ばれる形で、入力データの記憶領域に結果を上書きして返すことになる。 入力データは破壊されるがメモリ容量は半分で済む。 メモリアクセスの範囲も半分以下になるため、現代のCPUにおいてはメモリ容量よりもキャッシュに対する性能向上効果が期待できる。
ちなみに、ビット逆順テーブルは再帰でさらっと書くことが多い。
void bitReverse(int *buf, int n, int step, int i)
{
if (n < 2) {
*buf = i;
return;
}
n /= 2;
bitReverse(buf, n, step * 2, i);
bitReverse(buf + step, n, step * 2, i + n);
}
が、この記事を書いていて、C++ならstd::vectorで割とさらっと書けることに気づいた。
void bitReverse(std::vector<std::size_t> &v, std::size_t n)
{
v.push_back(0);
while (v.size() < n) {
std::vector<std::size_t> u;
for (auto i: v) {
u.push_back(i);
u.push_back(i + v.size());
}
v.swap(u);
}
}
size()関数を使うので、要素型をstd::size_tにしておかないと、符号が混ざってるとかなんとか怒られることがあるのがちょっとアレだけど。
std::vectorも何度も確保して捨てて、を繰り返すけど、事前にテーブルを生成する程度なら十分に実用になる。
元々の趣旨を考えると、逆FFTもできないと意味がない。 逆FFTは演算の順序を逆にすればいいだけである。 と、FFT (高速フーリエ・コサイン・サイン変換) の概略と設計法に書いてあって、このサイトを初めて読んだ当時（20年くらい前）は途方に暮れたものである。
具体的には以下の書き換えを行う。
Mのループを逆順にする。その他のループは逆にする必要はない。
a・bに対してc=a+b・d=a-bを求めている場所は、a=(c+d)/2・b=(c-d)/2とすれば元に戻る。
wr*wr+wi*wi=1になることを利用すれば、実虚それぞれにwrあるいはwiを乗じた式を適当に加減算すると、元に戻る組み合わせがある
実際に逆FFTをinvFFTという名前の関数にして超えいやっっっっっと実装するとこうなる。
"use strict";
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;
}
conj() {
this.im = -this.im;
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(dst) {
await STATUS("FFT: doing", src.length);
const n = src.length;
let idx = [ 0 ];
while (idx.length < n)
idx = idx.flatMap(v => [ v, v + idx.length]);
for (let M = n / 2; M >= 1; M /= 2) {
const D = 2 * M;
const Q = 4 * M;
for (let j = 0; j < M; j++) {
const even = dst[j];
const odd = dst[j + M];
dst[j] = even + odd;
dst[j + M] = even - odd;
}
for (let k = Q; k < n; k += Q) {
const i = idx[n / D - idx[k]];
const wr = Math.cos(2 * Math.PI * idx[k] * M / n);
const wi = Math.sin(2 * Math.PI * idx[k] * M / n);
for (let j = 0; j < M; j++) {
const re0 = dst[j + k];
const im0 = dst[j + i];
const re = dst[j + k + M];
const im = dst[j + i + M];
const re1 = re * wr - im * wi;
const im1 = re * wi + im * wr;
dst[j + k] = re0 + re1;
dst[j + k + M] = im1 - im0;
dst[j + i] = re0 - re1;
dst[j + i + M] = im0 + im1;
}
}
}
await STATUS("FFT: done");
return idx;
}
async function invFFT(dst) {
await STATUS("invFFT: doing", src.length);
const n = src.length;
let idx = [ 0 ];
while (idx.length < n)
idx = idx.flatMap(v => [ v, v + idx.length]);
for (let M = 1; M < n; M *= 2) {
const D = 2 * M;
const Q = 4 * M;
for (let j = 0; j < M; j++) {
const even = dst[j];
const odd = dst[j + M];
dst[j] = (even + odd) / 2;
dst[j + M] = (even - odd) / 2;
}
for (let k = Q; k < n; k += Q) {
const i = idx[n / D - idx[k]];
const wr = Math.cos(2 * Math.PI * idx[k] * M / n);
const wi = Math.sin(2 * Math.PI * idx[k] * M / n);
for (let j = 0; j < M; j++) {
const a = dst[j + k]; // = re0 + re1;
const b = dst[j + k + M]; // = im1 - im0;
const c = dst[j + i]; // = re0 - re1;
const d = dst[j + i + M]; // = im0 + im1;
const re0 = (a + c) / 2;
const im0 = (d - b) / 2;
const re1 = (a - c) / 2;
const im1 = (b + d) / 2;
dst[j + k] = re0;
dst[j + i] = im0;
// re1 * wr = re * wr * wr - im * wr * wi;
// re1 * wi = re * wr * wi - im * wi * wi;
// im1 * wr = re * wr * wi + im * wr * wr;
// im1 * wi = re * wi * wi + im * wr * wi;
const re = re1 * wr + im1 * wi;
const im = im1 * wr - re1 * wi;
dst[j + k + M] = re;
dst[j + i + M] = im;
}
}
}
await STATUS("invFFT: done");
return idx;
}
function expandFFTResult(raw) {
const n = raw.length;
const out = Array(n);
out[0] = new Complex(raw[0]);
if (n < 2)
return out;
for (let i = 1; i < n / 2; i++) {
out[i] = new Complex(raw[i], -raw[n - i]);
out[n - i] = new Complex(out[i]).conj();
}
out[n / 2] = new Complex(raw[n / 2]);
return out;
}
const src = Array(16).fill(0).map(_ => Math.random() - 0.5);
let raw = Array.from(src);
FFT(raw).then(async idx => {
const dft = await DFT(src.map(v => new Complex(v)));
const fft = expandFFTResult(idx.map(i => raw[i]));
fft.forEach((z, i) => LOG(z.sub(dft[i]).abs().toFixed(8)));
await invFFT(raw);
LOG(raw.map((v, i) => v - src[i]).reduce((sum, v) => sum + v * v, 0));
STATUS("finished");
});
本当はもう少し段階を踏んで、たとえばn=1から順に倍々にしていこうかと思ったが、面倒くさいその方が意外と難しかったりするので、申し訳ないがここは気合いで一気に。
実行結果は、最初の16行が正FFTの結果の誤差の絶対値、最後の1行が逆FFTのsrcに対する各要素の誤差の2乗の和である。
誤差を計算して2乗して足すところはreduce一発でも書けるのだが、一時変数を作るかかっこを付けるかしないといけないので、一度mapで誤差を求めてからreduceで2乗和を計算している。
少し変数を整理するとこうなる。
"use strict";
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;
}
conj() {
this.im = -this.im;
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(dst) {
await STATUS("FFT: doing", src.length);
const n = src.length;
let idx = [ 0 ];
while (idx.length < n)
idx = idx.flatMap(v => [ v, v + idx.length]);
for (let M = n / 2; M >= 1; M /= 2) {
const D = 2 * M;
const Q = 4 * M;
for (let j = 0; j < M; j++) {
const even = dst[j];
const odd = dst[j + M];
dst[j] = even + odd;
dst[j + M] = even - odd;
}
for (let k = Q; k < n; k += Q) {
const i = idx[n / D - idx[k]];
const wr = Math.cos(2 * Math.PI * idx[k] * M / n);
const wi = Math.sin(2 * Math.PI * idx[k] * M / n);
for (let j = 0; j < M; j++) {
const re0 = dst[j + k];
const im0 = dst[j + i];
const re = dst[j + k + M];
const im = dst[j + i + M];
const re1 = re * wr - im * wi;
const im1 = re * wi + im * wr;
dst[j + k] = re0 + re1;
dst[j + k + M] = im1 - im0;
dst[j + i] = re0 - re1;
dst[j + i + M] = im0 + im1;
}
}
}
await STATUS("FFT: done");
return idx;
}
async function invFFT(dst) {
await STATUS("invFFT: doing", src.length);
const n = src.length;
let idx = [ 0 ];
while (idx.length < n)
idx = idx.flatMap(v => [ v, v + idx.length]);
for (let M = 1; M < n; M *= 2) {
const D = 2 * M;
const Q = 4 * M;
for (let j = 0; j < M; j++) {
const even = dst[j];
const odd = dst[j + M];
dst[j] = (even + odd) / 2;
dst[j + M] = (even - odd) / 2;
}
for (let k = Q; k < n; k += Q) {
const i = idx[n / D - idx[k]];
const wr = Math.cos(2 * Math.PI * idx[k] * M / n);
const wi = Math.sin(2 * Math.PI * idx[k] * M / n);
for (let j = 0; j < M; j++) {
const a = dst[j + k];
const b = dst[j + k + M];
const c = dst[j + i];
const d = dst[j + i + M];
dst[j + k] = (a + c) / 2;
dst[j + i] = (d - b) / 2;
const re1 = (a - c) / 2;
const im1 = (d + b) / 2;
dst[j + k + M] = re1 * wr + im1 * wi;
dst[j + i + M] = im1 * wr - re1 * wi;
}
}
}
await STATUS("invFFT: done");
return idx;
}
function expandFFTResult(raw) {
const n = raw.length;
const out = Array(n);
out[0] = new Complex(raw[0]);
if (n < 2)
return out;
for (let i = 1; i < n / 2; i++) {
out[i] = new Complex(raw[i], -raw[n - i]);
out[n - i] = new Complex(out[i]).conj();
}
out[n / 2] = new Complex(raw[n / 2]);
return out;
}
const src = Array(16).fill(0).map(_ => Math.random() - 0.5);
let raw = Array.from(src);
FFT(raw).then(async idx => {
const dft = await DFT(src.map(v => new Complex(v)));
const fft = expandFFTResult(idx.map(i => raw[i]));
fft.forEach((z, i) => LOG(z.sub(dft[i]).abs().toFixed(8)));
await invFFT(raw);
LOG(raw.map((v, i) => v - src[i]).reduce((sum, v) => sum + v * v, 0));
STATUS("finished");
});
実は複素FFTの場合は/2が不要で、そうすると結果がn倍になるので最後にnで割るのが普通である。
DFT関数もn倍の値が得られるのでわざわざscaleメソッドを作ったのであった。
これに対して実FFTの場合、この辺でなくなってしまったループの演算がない分、結果が合わなくなる。
元々はちゃんと回転因子を乗じて加減算していたのだから、ここの値を2倍にしないといけない。
これを補って、/2を取り除くとこうなる。
"use strict";
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;
}
conj() {
this.im = -this.im;
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(dst) {
await STATUS("FFT: doing", src.length);
const n = src.length;
let idx = [ 0 ];
while (idx.length < n)
idx = idx.flatMap(v => [ v, v + idx.length]);
for (let M = n / 2; M >= 1; M /= 2) {
const D = 2 * M;
const Q = 4 * M;
for (let j = 0; j < M; j++) {
const even = dst[j];
const odd = dst[j + M];
dst[j] = even + odd;
dst[j + M] = even - odd;
}
for (let k = Q; k < n; k += Q) {
const i = idx[n / D - idx[k]];
const wr = Math.cos(2 * Math.PI * idx[k] * M / n);
const wi = Math.sin(2 * Math.PI * idx[k] * M / n);
for (let j = 0; j < M; j++) {
const re0 = dst[j + k];
const im0 = dst[j + i];
const re = dst[j + k + M];
const im = dst[j + i + M];
const re1 = re * wr - im * wi;
const im1 = re * wi + im * wr;
dst[j + k] = re0 + re1;
dst[j + k + M] = im1 - im0;
dst[j + i] = re0 - re1;
dst[j + i + M] = im0 + im1;
}
}
}
await STATUS("FFT: done");
return idx;
}
async function invFFT(dst) {
await STATUS("invFFT: doing", src.length);
const n = src.length;
let idx = [ 0 ];
while (idx.length < n)
idx = idx.flatMap(v => [ v, v + idx.length]);
for (let M = 1; M < n; M *= 2) {
const D = 2 * M;
const Q = 4 * M;
for (let j = 0; j < M; j++) {
const even = dst[j];
const odd = dst[j + M];
dst[j] = even + odd;
dst[j + M] = even - odd;
}
if (D < n) {
for (let j = 0; j < M; j++) {
dst[j + D] *= 2;
dst[j + D + M] *= 2;
}
}
for (let k = Q; k < n; k += Q) {
const i = idx[n / D - idx[k]];
const wr = Math.cos(2 * Math.PI * idx[k] * M / n);
const wi = Math.sin(2 * Math.PI * idx[k] * M / n);
for (let j = 0; j < M; j++) {
const a = dst[j + k];
const b = dst[j + k + M];
const c = dst[j + i];
const d = dst[j + i + M];
dst[j + k] = a + c;
dst[j + i] = d - b;
const re1 = a - c;
const im1 = d + b;
dst[j + k + M] = re1 * wr + im1 * wi;
dst[j + i + M] = im1 * wr - re1 * wi;
}
}
}
await STATUS("invFFT: done");
return idx;
}
function expandFFTResult(raw) {
const n = raw.length;
const out = Array(n);
out[0] = new Complex(raw[0]);
if (n < 2)
return out;
for (let i = 1; i < n / 2; i++) {
out[i] = new Complex(raw[i], -raw[n - i]);
out[n - i] = new Complex(out[i]).conj();
}
out[n / 2] = new Complex(raw[n / 2]);
return out;
}
const src = Array(16).fill(0).map(_ => Math.random() - 0.5);
let raw = Array.from(src);
FFT(raw).then(async idx => {
const dft = await DFT(src.map(v => new Complex(v)));
const fft = expandFFTResult(idx.map(i => raw[i]));
fft.forEach((z, i) => LOG(z.sub(dft[i]).abs().toFixed(8)));
await invFFT(raw);
LOG(raw.map((v, i) => v / raw.length - src[i]).reduce((sum, v) => sum + v * v, 0));
STATUS("finished");
});
この演算は、正FFTのmのループ初回ではやっていなかったので、逆FFTではループの最終回、つまりm=n/2の時は何もしなくてよい。
そこで、削ってしまったif文をとりあえず復活させて、変数を書き換える。
どこを2倍するかだか、kがnまで行ったすぐ後に計算するはずの要素が対象だ。
ビット逆順なのだから、kのループで使っているビットのすぐ下のビットを1にすれば、元のビット順序で桁上げしたことになる。
どこのビットを使っているかというと、Qずつ増えるのだから、Qの位置がkが変化する場合の最下位ビットになる。
そのすぐ下のビットを1にすればよいので、Q/2、つまりdst[D]について演算する。
もう一方のデータ、つまり虚部は、他のjループを見てみるとMだけ離れた要素について演算を行っているので、結局、D+jの位置とM+D+jの位置の要素を2倍すればよい。
結果はn倍になるのでnで割って比較する。
しかし、やっぱり無駄に2倍しているのが悔しい。 どれくらいの演算量があるのか検討してみる。
M=1の時dst[0]とdst[1]が書き換わる
dst[2]とdst[3]が2倍になる
dstの残りは3番目のループで処理される
M=2の時dst[0]からdst[3]が書き換わる
dst[4]からdst[7]が2倍になる
dstの残りは3番目のループで処理される
M=4の時dst[0]からdst[7]が書き換わる
dst[8]からdst[15]が2倍になる
dstの残りは3番目のループで処理される
M=n/4の時n/4回実行され、dst[0]からdst[n/2-1]が書き換わる
n/4回実行され、dst[n/2]からdst[n-1]が2倍になる
M=n/2の時（最後）n/2回実行され、dst[0]からdst[n-1]が書き換わる
大体n回実行されていることが分かるのだが、よく見るとdst[2]からdst[n-1]までが2倍の対象で、dst[0]とdst[1]は対象外だということが分かる。
それ以外の要素はどこかでふたつ目のループによって2倍され、以後は最初のループで処理されて出力されている。
これが結果が合わない原因である。
ということは、2倍のループを無くし、dst[0]とdst[1]をあらかじめ2で割っておくと全体が半分になる。
invFFTはビット逆順のデータを受け取っているので、
\(F(0)\)
と
\(F(N/2)\)
の位置、つまり、実FFTとして見たときに複素数ではなく実数になっている部分に相当する。
"use strict";
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;
}
conj() {
this.im = -this.im;
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(dst) {
await STATUS("FFT: doing", src.length);
const n = src.length;
let idx = [ 0 ];
while (idx.length < n)
idx = idx.flatMap(v => [ v, v + idx.length]);
for (let M = n / 2; M >= 1; M /= 2) {
const D = 2 * M;
const Q = 4 * M;
for (let j = 0; j < M; j++) {
const even = dst[j];
const odd = dst[j + M];
dst[j] = even + odd;
dst[j + M] = even - odd;
}
for (let k = Q; k < n; k += Q) {
const i = idx[n / D - idx[k]];
const wr = Math.cos(2 * Math.PI * idx[k] * M / n);
const wi = Math.sin(2 * Math.PI * idx[k] * M / n);
for (let j = 0; j < M; j++) {
const re0 = dst[j + k];
const im0 = dst[j + i];
const re = dst[j + k + M];
const im = dst[j + i + M];
const re1 = re * wr - im * wi;
const im1 = re * wi + im * wr;
dst[j + k] = re0 + re1;
dst[j + k + M] = im1 - im0;
dst[j + i] = re0 - re1;
dst[j + i + M] = im0 + im1;
}
}
}
await STATUS("FFT: done");
return idx;
}
async function invFFT(dst) {
await STATUS("invFFT: doing", src.length);
const n = src.length;
let idx = [ 0 ];
while (idx.length < n)
idx = idx.flatMap(v => [ v, v + idx.length]);
dst[0] /= 2;
dst[1] /= 2;
for (let M = 1; M < n; M *= 2) {
const D = 2 * M;
const Q = 4 * M;
for (let j = 0; j < M; j++) {
const even = dst[j];
const odd = dst[j + M];
dst[j] = even + odd;
dst[j + M] = even - odd;
}
for (let k = Q; k < n; k += Q) {
const i = idx[n / D - idx[k]];
const wr = Math.cos(2 * Math.PI * idx[k] * M / n);
const wi = Math.sin(2 * Math.PI * idx[k] * M / n);
for (let j = 0; j < M; j++) {
const a = dst[j + k];
const b = dst[j + k + M];
const c = dst[j + i];
const d = dst[j + i + M];
dst[j + k] = a + c;
dst[j + i] = d - b;
const re1 = a - c;
const im1 = d + b;
dst[j + k + M] = re1 * wr + im1 * wi;
dst[j + i + M] = im1 * wr - re1 * wi;
}
}
}
await STATUS("invFFT: done");
return idx;
}
function expandFFTResult(raw) {
const n = raw.length;
const out = Array(n);
out[0] = new Complex(raw[0]);
if (n < 2)
return out;
for (let i = 1; i < n / 2; i++) {
out[i] = new Complex(raw[i], -raw[n - i]);
out[n - i] = new Complex(out[i]).conj();
}
out[n / 2] = new Complex(raw[n / 2]);
return out;
}
const src = Array(16).fill(0).map(_ => Math.random() - 0.5);
let raw = Array.from(src);
FFT(raw).then(async idx => {
const dft = await DFT(src.map(v => new Complex(v)));
const fft = expandFFTResult(idx.map(i => raw[i]));
fft.forEach((z, i) => LOG(z.sub(dft[i]).abs().toFixed(8)));
await invFFT(raw);
LOG(raw.map((v, i) => v * 2 / raw.length - src[i]).reduce((sum, v) => sum + v * v, 0));
STATUS("finished");
});
全体を2で割っているので、スケーリングがnではなくn/2になる点に注意。
FFTの実装の話はとりあえずこれで終わりである。
19 Feb 2026: 新規作成
ご意見・ご要望の送り先は あかもず仮店舗 の末尾をご覧ください。
Copyright (C) 2026 akamoz.jp