LoginSignup
44
32

More than 3 years have passed since last update.

@bellbind

高速フーリエ変換FFTを理解する

Last updated at Posted at 2019-04-17

この文書は、離散フーリエ変換DFTの理解をもとに、高速フーリエ変換FFTの振る舞いを理解することを目標とするものです。

はじめにDFTが何をするものかを簡単に説明し、そのDFT処理をJavaScriptコードに書き下し、FFT実装の検証対象となる具体例でのDFT実行結果を示します。
そして、再帰版FFT実装を天下り的に提示し、具体例を用いてFFTの中で何が行われるかを解析します。
最後に、この解析結果にもとづいて、ループ版FFTを導出します。

1. 離散フーリエ級数(DFT)

離散フーリエ変換(Discrete Fourier Transform, DFT) は、一定期間の時系列の数値列を、その期間を周期的に繰り返す合成波とみなし、その周波数成分の複素数値の列(低周波数から高周波数の順)を算出する仕組みです。

各周波数における複素数値のうち、実数成分が偶関数であるcos波の係数であり、虚数成分が奇関数であるsin波の係数となります。
この係数の量が、時系列データの分析のために用いることができます。たとえば、一番大きい係数の周波数を取り出して音階を算出する、などが行えます。

離散フーリエ変換では、時系列データがN個な場合、区間Nで1/2周期とする周波数を基本に、0倍(定数成分)からN-1倍までのN個の周波数成分を算出します。
時系列データを$f(t), t = 0 \dots N-1$、周波数成分を$F(k), k=0 \dots N-1$とすると、この関係は数式で以下のようになります。

\begin{align}
F(k) = & \sum^{N-1}_{n=0} f(n) \exp(\frac{-2 \pi i}{N} n k) \\
f(t) = & \frac{1}{N} \sum^{N-1}_{n=0} F(n) \exp(\frac{2 \pi i}{N} n t) \\
\end{align}

$f$も$F$もN要素の複素数値配列として、前者の右辺の計算で$F(k)$を計算することを「離散フーリエ変換(DFT)」、
後者の右辺の計算で$f(t)$を計算することを、「逆離散フーリエ変換(Inverse DFT,IDFT)」と呼びます。

右辺にある$\exp(i\theta)$の部分は複素数であり、オイラーの公式$\exp(i\theta) = cos{\theta} + i sin{\theta}$となります。
離散フーリエ変換のプログラミングでは、計算でsinやcosを使い、実数成分と虚数成分のペアとして、各成分ごとにコーディングすることも多いです。
この場合、複素数の足し算と掛け算を、以下のように、各成分ごとに行います。

\begin{align}
a + b = & (a_x + i a_y) + (b_x + i b_y) \\
      = & (a_x + b_x) + i(a_y + b_y) \\
a \times b = & (a_x + i a_y) \times (b_x + i b_y) \\
           = & (a_x \times b_x - a_y \times b_y) + i(a_x \times b_y + a_y \times b_x) \\
\end{align}

この実数虚数成分を分けての計算の場合には、さらに「時系列データは実数値のみ」として、以下のように虚数部の計算を省くこともあります(ただし、ここではこの実数限定な実装は行いません)。

  • DFT: $f(t)$の虚数部が0であることから、$exp(i)$の複素数との積は、実数部虚数部ともに$f(t)$をかけるだけにする
  • IDFT: $f(t)$の実部側($a_x \times b_x - a_y \times b_y$)のみ計算する

2. JavaScriptでのDFT実装コード

まず、数値２要素Arrayを複素数として、複素数計算のヘルパーを用意します。
$\exp(i\theta)$は、実数thetaを引数に取り、複素数を返す関数expi(theta)として用意します。
$\sum$に対応するものとして、複素数配列csの総和を取る関数isum(cs)も用意します。

[1]
function expi(theta) {return [Math.cos(theta), Math.sin(theta)];}
function iadd([ax, ay], [bx, by]) {return [ax + bx, ay + by];}
function isub([ax, ay], [bx, by]) {return [ax - bx, ay - by];}
function imul([ax, ay], [bx, by]) {return [ax * bx - ay * by, ax * by + ay * bx];}
function isum(cs) {return cs.reduce((s, c) => iadd(s, c), [0, 0]);}

複素数の時系列データfを受け取り、複素数の周波数成分配列を返すDFT関数dft(f)は、以下のようになります:

[2]
function dft(f) {
    const N = f.length, T = -2 * Math.PI / N;
    return [...Array(N).keys()].map(k => isum(
        f.map((fn, n) => imul(fn, expi(T * n * k)))
    ));
}

複素数の周波数成分配列Fを受け取り、複素数の時系列データを返すIDFT関数idft(F)は、以下のようになります:

[3]
function idft(F) {
    const N = F.length, T = 2 * Math.PI / N;
    return [...Array(N).keys()].map(t => isum(
        F.map((Fn, n) => imul(Fn, expi(T * n * t)))
    )).map(([r, i]) => [r / N, i / N]);
}

以下、実数値の時系列データfr0を用意し、複素数化したf0を用いて、DFT，IDFTを行う実行例です:

[4]
{
    const fr0 = [1,3,4,2, 5,6,2,4, 0,1,3,4, 5,62,2,3];
    const f0 = fr0.map(r => [r, 0]);

    const F = dft(f0);
    const f1 = idft(F);
    const fr1 = f1.map(([r]) => r);

    console.log("fr0:", fr0);
    console.log("F:", F);
    console.log("f1:", f1);
    console.log("fr1:", fr1.map(Math.round));
}

fr0: [ 1, 3, 4, 2, 5, 6, 2, 4, 0, 1, 3, 4, 5, 62, 2, 3 ]
F: [ [ 107, 0 ],
  [ 23.295891661412693, 51.72985580737281 ],
  [ -53.54772721475247, 42.96194077712561 ],
  [ -49.21391810443097, -25.674384455895552 ],
  [ -7.869471018244225e-14, -59 ],
  [ 49.79970454205781, -24.2601708935226 ],
  [ 35.54772721475254, 48.96194077712551 ],
  [ -19.8816780990393, 53.14406936974603 ],
  [ -63, 9.583546756334772e-14 ],
  [ -19.881678099039497, -53.14406936974596 ],
  [ 35.54772721475236, -48.961940777125726 ],
  [ 49.799704542057995, 24.260170893522123 ],
  [ 2.3608413054732673e-13, 58.99999999999999 ],
  [ -49.21391810443097, 25.67438445589565 ],
  [ -53.54772721475294, -42.96194077712513 ],
  [ 23.29589166141229, -51.72985580737299 ] ]
f1: [ [ 0.9999999999999822, -7.993605777301127e-15 ],
  [ 2.9999999999999187, 3.68594044175552e-14 ],
  [ 4.000000000000027, 6.17284001691587e-14 ],
  [ 2.0000000000000044, 1.9095836023552692e-14 ],
  [ 5.000000000000019, -8.659739592076221e-15 ],
  [ 5.999999999999996, -1.509903313490213e-14 ],
  [ 1.9999999999999831, 2.2426505097428162e-14 ],
  [ 3.999999999999982, 3.907985046680551e-14 ],
  [ -1.2434497875801753e-14, 8.704148513061227e-14 ],
  [ 1.000000000000095, 1.9539925233402755e-14 ],
  [ 3.0000000000000027, -1.865174681370263e-14 ],
  [ 4.000000000000051, -1.9539925233402755e-14 ],
  [ 5.000000000000107, -1.2212453270876722e-14 ],
  [ 61.99999999999999, -2.816635813474022e-13 ],
  [ 1.9999999999999196, -3.6415315207705135e-14 ],
  [ 3.000000000000007, -1.3766765505351941e-14 ] ]
fr1: [ 1, 3, 4, 2, 5, 6, 2, 4, -0, 1, 3, 4, 5, 62, 2, 3 ]

DFTしたものをIDFTし誤差を丸めたfr1は、もとの値fr0の値にもどっているのが確認できました(ただし、0の符号を除く)。

3. DFTとIDFTの計算部分の共通化

DFTとIDFTの処理の違いは、expiに渡すTの符号が違う部分と、IDFTの最後のNで割る部分だけです。

この共通部分のコードをdftc(c, T)としてまとめると、以下のようになります。

[5]
function dftc(c, T) {
    return [...Array(c.length).keys()].map(i => isum(
        c.map((cn, n) => imul(cn, expi(T * n * i)))
    ));
}

function dft(f) {
    const N = f.length, T = -2 * Math.PI / N;
    return dftc(f, T);
}
function idft(F) {
    const N = F.length, T = 2 * Math.PI / N;
    return dftc(F, T).map(([r, i]) => [r / N, i / N]);
}

以降のFFTでも逆変換との計算部分の共通化をしたものとして、プログラムを考えていきます。

4. 高速フーリエ変換(FFT)

高速フーリエ変換は、要素数に制限を加えることによって可能となった、離散フーリエ変換の計算量の少ない実装です。
時系列データから周波数成分列データへの変換を「高速フーリエ変換(FFT)」、その逆変換を「逆高速フーリエ変換(Inverse FFT, IFFT)」と呼びます。

前述の離散フーリエ変換DFTやその逆変換IDFTは、ともに結果の各成分ごとに、N回の複素数の積を取りその総和を計算するため、その計算量はO(N^2)となります。
一方、高速フーリエ変換(Fast Fourier Transform, FFT)では、Nが2のべき乗のときに使える手法で、変換と逆変換の計算量をともにO(NlogN)に抑える事ができます。
たとえばN=8のとき、DFTでは8x8=64回計算するところを、FFTでは8 x log2(8) = 8x3 = 24回の計算でできるようになります。

ここでは、先に再帰関数によるFFT関数の実装を与え、その実行結果を分析することによって、FFTではどういう計算方法になっているのかについて理解していきます。

5. JavaScriptでの再帰版FFT実装コード

FFTとIFFTの共通部分を再帰関数fftrec(c, T, N, s, w)で実装します。
再帰呼び出しでは、データ列を二分させて呼び出します。つまり、log2(N)段再帰され、各段のトータルはN回の計算になることで、O(NlogN)の計算になります。

  • cとTはDFTのときと同じ、データ列とexpiに渡す定数です。Tは再帰呼び出しで2倍にしていきます。
  • Nは、その再帰ステップで処理するデータ列の要素数(2のべき乗の数)で、再帰呼出しで半減させていきます。
  • sは、その最期ステップで扱うcでの先頭インデックスです。
  • wは、2の再帰段数乗(1 => 2 => 4 => 8 => ...)です。再帰の終端(N=1)のときに、インデックスsとs+wを隣り合わせにし、計算させます(たすき掛けの相手)。

sとwは、データ列cでのインデックスでの、演算で使う２つの複素数を取り出すためのデータであり、コードではわかりにくいものとなります。
このあとで説明するように、具体的なデータを使って実行した結果を見れば、これが何をするためのコードだったかについては、すぐわかるかと思います。

まずは、再帰版FFT実装コードを提示します。FFT関数をfft0(f)、IFFT関数をifft0(F)にします。

[6]
function fftrec(c, T, N, s = 0, w = 1) {
    if (N === 1) return [c[s]];
    const Nh = N / 2, Td = T * 2, wd = w * 2;
    const rec = fftrec(c, Td, Nh, s, wd).concat(fftrec(c, Td, Nh, s + w, wd));
    for (let i = 0; i < Nh; i++) {
        const l = rec[i], re = imul(rec[i + Nh], expi(T * i));
        [rec[i], rec[i + Nh]] = [iadd(l, re), isub(l, re)];
    }
    return rec;
}

function fft0(f) {
    const N = f.length, T = -2 * Math.PI / N;
    return fftrec(f, T, N);
}
function ifft0(F) {
    const N = F.length, T = 2 * Math.PI / N;
    return fftrec(F, T, N).map(([r, i]) => [r / N, i / N]);
}

複素数計算は、forループ部分でのみ行っています。
各再帰関数内ではNh回計算します。
Nhは最初はNで、再帰ごとに半減します。一方再帰呼び出しでは2回再帰呼び出しをするので、最初の1つから始まって再帰ごとに2倍の数になっていきます。
この結果、各再帰段でトータルすると、結局はどの段でもN回複素数の計算していることになります。

再帰段数はlog2(N)なので、計算量がO(NlogN)なコードとなっていることがわかるでしょう。

このfft0とifft0を、DFTと同じ、16要素の時系列データfr0で実行してみます。

[7]
{
    const fr0 = [1,3,4,2, 5,6,2,4, 0,1,3,4, 5,62,2,3];
    const f0 = fr0.map(r => [r, 0]);

    const F = fft0(f0);
    const f1 = ifft0(F);
    const fr1 = f1.map(([r]) => r);

    console.log("fr0:", fr0);
    console.log("F:", F);
    console.log("f1:", f1);
    console.log("fr1:", fr1.map(Math.round));
}

fr0: [ 1, 3, 4, 2, 5, 6, 2, 4, 0, 1, 3, 4, 5, 62, 2, 3 ]
F: [ [ 107, 0 ],
  [ 23.29589166141268, 51.729855807372815 ],
  [ -53.5477272147525, 42.961940777125584 ],
  [ -49.21391810443094, -25.67438445589562 ],
  [ 3.6127080574846916e-15, -59 ],
  [ 49.799704542057846, -24.260170893522517 ],
  [ 35.54772721475249, 48.96194077712559 ],
  [ -19.8816780990396, 53.14406936974591 ],
  [ -63, 0 ],
  [ -19.881678099039586, -53.14406936974591 ],
  [ 35.5477272147525, -48.961940777125584 ],
  [ 49.799704542057846, 24.260170893522528 ],
  [ -3.6127080574846916e-15, 59 ],
  [ -49.21391810443094, 25.67438445589561 ],
  [ -53.54772721475249, -42.96194077712559 ],
  [ 23.295891661412693, -51.729855807372815 ] ]
f1: [ [ 1, 0 ],
  [ 3.000000000000001, -1.3877787807814457e-15 ],
  [ 3.999999999999999, -1.0143540619928357e-16 ],
  [ 2, -3.0531133177191805e-16 ],
  [ 5, 0 ],
  [ 6, -1.7763568394002505e-15 ],
  [ 2, 4.592425496802574e-17 ],
  [ 3.9999999999999996, -2.498001805406602e-16 ],
  [ 0, 0 ],
  [ 0.9999999999999991, 1.3877787807814457e-15 ],
  [ 3.000000000000001, 9.586896263232085e-18 ],
  [ 4.000000000000001, 3.608224830031759e-16 ],
  [ 5, 0 ],
  [ 62, 1.7763568394002505e-15 ],
  [ 2, 4.592425496802574e-17 ],
  [ 2.9999999999999996, 1.942890293094024e-16 ] ]
fr1: [ 1, 3, 4, 2, 5, 6, 2, 4, 0, 1, 3, 4, 5, 62, 2, 3 ]

DFTのときとほぼ同様の結果となりました。

6. 再帰版FFTの実行結果を分析する

この再帰実装FFTのfftrecは、再帰呼び出しの前後で、別々の処理を行っています:

  • 前半: データの並べ替えをする
  • 後半: 複素数の計算を行う

まず、前半のデータの並べ替えを解析します。

6.1. FFT前半部: 並べ替え部分の実行結果を分析する

fftrecの並び替えの部分だけを切り出したものをreorderとして切り出したものは、以下のようになります:

[8]
function reorder(c, N = c.length, s = 0, w = 1) {
    if (N === 1) return [c[s]];
    const Nh = N / 2, wd = w * 2;
    const rec = reorder(c, Nh, s, wd).concat(reorder(c, Nh, s + w, wd));
    return rec;
}
console.log(reorder([0,1,2,3]));
console.log(reorder([0,1,2,3,4,5,6,7]));
console.log(reorder([0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15]));

[ 0, 2, 1, 3 ]
[ 0, 4, 2, 6, 1, 5, 3, 7 ]
[ 0, 8, 4, 12, 2, 10, 6, 14, 1, 9, 5, 13, 3, 11, 7, 15 ]

この結果の[0,2,1,3]というのは、rreorderによって[c[0], c[2], c[1], c[3]]になった、ということを意味しています。

この並びは一見わかりにくい結果ですが、数値を２進数で表現することで、その意味がすぐわかるようになります。
整数をk桁２進数文字列にする関数toBin(k)(n)は、以下の通りになります。

[9]
function toBin(k) {
    return n => n.toString(2).padStart(k, "0");
}

これを使ってreorderするものと、した結果をそれぞれ２進数で表示してみます。

[10]
console.log([0,1,2,3].map(toBin(2)));
console.log(reorder([0,1,2,3]).map(toBin(2)));

console.log([0,1,2,3,4,5,6,7].map(toBin(3)));
console.log(reorder([0,1,2,3,4,5,6,7]).map(toBin(3)));

console.log([0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15].map(toBin(4)));
console.log(reorder([0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15]).map(toBin(4)));

[ '00', '01', '10', '11' ]
[ '00', '10', '01', '11' ]
[ '000', '001', '010', '011', '100', '101', '110', '111' ]
[ '000', '100', '010', '110', '001', '101', '011', '111' ]
[ '0000',
  '0001',
  '0010',
  '0011',
  '0100',
  '0101',
  '0110',
  '0111',
  '1000',
  '1001',
  '1010',
  '1011',
  '1100',
  '1101',
  '1110',
  '1111' ]
[ '0000',
  '1000',
  '0100',
  '1100',
  '0010',
  '1010',
  '0110',
  '1110',
  '0001',
  '1001',
  '0101',
  '1101',
  '0011',
  '1011',
  '0111',
  '1111' ]

このように、reorderの再帰処理でやってることは、インデックスを２進数での前後を逆転させた位置に移したものだったのです。

次に、再帰実装ではなく、k桁２進数の順序逆転関数revBit(k, n)を実装し、それによって、reorderを再帰なしで実装してみます。

[11]
function revBit(k, n) {
    let r = 0;
    for (let i = 0; i < k; i++) r = (r << 1) | ((n >>> i) & 1);
    return r;
}
[12]
function reorder(c, N = c.length) {
    const k = Math.log2(N);
    console.assert(Number.isInteger(k), "c.length should be power of 2");
    return c.map((_, i) => c[revBit(k, i)]);
}

このコードで、再帰関数のときと同様の結果が、map(ループ)によって得られます。

整数が32bitであることを前提にすると、ビット操作でrevBitをループなしで以下のように実装できます。

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

これは、"ab cd ef gh"を1つづつ入れ替えると"ba dc fe hg"になり、次に2つづつで入れ替えると"dcba hgfe"となり、
最後に4つで入れ替えると"hgfe dcba"となるのと同じやり方です。そして、返す直前にk桁に切り詰めています。

6.2. FFT後半部: 複素数計算部分の実行結果を分析する

ここでは、要素数N=8での計算について考えていきます。

まず、DFTでの計算について、N=8のときで確認します:

F(k) = f(0) \exp(\frac{-2\pi i}{8} 0 k) + f(1) \exp(\frac{-2\pi i}{8} 1 k) + 
       f(2) \exp(\frac{-2\pi i}{8} 2 k) + f(3) \exp(\frac{-2\pi i}{8} 3 k) +
       f(4) \exp(\frac{-2\pi i}{8} 4 k) + f(5) \exp(\frac{-2\pi i}{8} 5 k) +
       f(6) \exp(\frac{-2\pi i}{8} 6 k) + f(7) \exp(\frac{-2\pi i}{8} 7 k)

ここで、$E(n) = \exp(\frac{-2\pi i}{8} n)$とすると、$\exp(\frac{-2\pi i}{8} t k)$は短く、E(t * k)とかけます。
このとき、$F(k)$は以下のようになります:

F(k) = E(0 * k)f(0) + E(1 * k)f(1) + E(2 * k)f(2) + E(3 * k)f(3) + E(4 * k)f(4) + E(5 * k)f(5) + E(6 * k)f(6) + E(7 * k)f(7)

まず、この係数E(n * k)に着目し、表にすると以下のようになっています:

F\f f(0) f(1) f(2) f(3) f(4) f(5) f(6) f(7)
F(0) E(0) E(0) E(0) E(0) E(0) E(0) E(0) E(0)
F(1) E(0) E(1) E(2) E(3) E(4) E(5) E(6) E(7)
F(2) E(0) E(2) E(4) E(6) E(8) E(10) E(12) E(14)
F(3) E(0) E(3) E(6) E(9) E(12) E(15) E(18) E(21)
F(4) E(0) E(4) E(8) E(12) E(16) E(20) E(24) E(28)
F(5) E(0) E(5) E(10) E(15) E(20) E(25) E(30) E(35)
F(6) E(0) E(6) E(12) E(18) E(24) E(30) E(36) E(42)
F(7) E(0) E(7) E(14) E(21) E(28) E(35) E(42) E(49)

E(n)はべき乗でもあるから、E(n)E(m) = E(n + m)です。

また、$\exp(-2\pi i) = \exp(\frac{-2\pi i}{8} 8) = 1$、$\exp(-\pi i) = \exp(\frac{-2\pi}{8} 4) = -1$といった関係があります。
これはE(0) = E(8) = 1、E(4) = -1になるので、E(n + 8) = E(n)E(8) = E(n)、E(n + 4) = E(n)E(4) = -E(n)という関係を持ちます。

E(n + 8) = E(n)から、先の表のE(n)のnを8で割った余りにすると以下のようになります:

F\f f(0) f(1) f(2) f(3) f(4) f(5) f(6) f(7)
F(0) E(0) E(0) E(0) E(0) E(0) E(0) E(0) E(0)
F(1) E(0) E(1) E(2) E(3) E(4) E(5) E(6) E(7)
F(2) E(0) E(2) E(4) E(6) E(0) E(2) E(4) E(6)
F(3) E(0) E(3) E(6) E(1) E(4) E(7) E(2) E(5)
F(4) E(0) E(4) E(0) E(4) E(0) E(4) E(0) E(4)
F(5) E(0) E(5) E(2) E(7) E(4) E(1) E(6) E(3)
F(6) E(0) E(6) E(4) E(2) E(0) E(6) E(4) E(2)
F(7) E(0) E(7) E(6) E(5) E(4) E(3) E(2) E(1)

この表に、先に調べた、前半でのビット前後反転した順でのデータ列の入れ替え[ 0, 4, 2, 6, 1, 5, 3, 7 ]を行うと、以下のようになります:

F\f f(0) f(4) f(2) f(6) f(1) f(5) f(3) f(7)
F(0) E(0) E(0) E(0) E(0) E(0) E(0) E(0) E(0)
F(1) E(0) E(4) E(2) E(6) E(1) E(5) E(3) E(7)
F(2) E(0) E(0) E(4) E(4) E(2) E(2) E(6) E(6)
F(3) E(0) E(4) E(6) E(2) E(3) E(7) E(1) E(5)
F(4) E(0) E(0) E(0) E(0) E(4) E(4) E(4) E(4)
F(5) E(0) E(4) E(2) E(6) E(5) E(1) E(7) E(3)
F(6) E(0) E(0) E(4) E(4) E(6) E(6) E(2) E(2)
F(7) E(0) E(4) E(6) E(2) E(7) E(3) E(5) E(1)

ここで、E(n)のnだけ並べてみます。

F\f 0 4 2 6 1 5 3 7
0 0 0 0 0 0 0 0 0
1 0 4 2 6 1 5 3 7
2 0 0 4 4 2 2 6 6
3 0 4 6 2 3 7 1 5
4 0 0 0 0 4 4 4 4
5 0 4 2 6 5 1 7 3
6 0 0 4 4 6 6 2 2
7 0 4 6 2 7 3 5 1

この並び替えによって、F(０)とF(4)のあいだなどで、共通部分のあるパターンがなんとなく見えてきたかもしれません。
次は、このパターンについて見ていきます。

先に示したE(n + 8) = E(n)、E(n + 4) = -E(n)より、E(n)には:

  • -E(0) = E(4)
  • -E(1) = E(5)
  • -E(2) = E(6)
  • -E(3) = E(7)

かつ:

  • -E(4) = E(0)
  • -E(5) = E(1)
  • -E(6) = E(2)
  • -E(7) = E(3)

の関係があります。つまりE(n)のnには、0と4、1と5、2と6、3と7は、互いにマイナスをかけた関係にあります。

ここからたとえば、F(0)とF(4)は、F(0)の前側の0426部分をL0、後ろ側の1537部分をR0とすると、

  • F(0) = L0 + R0
  • F(4) = L0 - R0 = L0 + E(4)R0

という関係になっています。同様に、

  • F(1) = L1 + R1
  • F(5) = L1 - R1
  • F(2) = L2 + R2
  • F(6) = L2 - R2
  • F(3) = L3 + R3
  • F(7) = L3 - R3

ともなっています。

最後に、このLRの計算を、再帰の一番奥の段から、足し合わせていくときの計算ステップを具体的に並べていきます。

再帰版FFTの後半の計算コードは、同一段数の再帰処理でまとめることで、各再帰段ごとで全要素を1回づつ更新しています。
(実際の動作では、葉がN要素な二分木での深さ優先の帰りがけ順で、各段としてはバラバラに計算されます。
しかし、二分木のそれぞれの枝での計算範囲は完全に分割されているので、幅優先で処理しても同じ結果になります。)
N=8では、Nを半減させてN=1になるまで、3段の(2回づつの)再帰呼び出しが、再帰関数中で行われることになります。

以下は、この再帰段数を一度に巻き戻すことの順に、計算される内容と結果を列挙したものです。
入力値の時系列データの値は、それぞれf(0) = f0、f(1) = f1、 ...、f(7) = f7と記述します。

再帰段3

前半の入れ替え直後の初期値です。

  • F0(0) = f0
  • F0(1) = f4
  • F0(2) = f2
  • F0(3) = f6
  • F0(4) = f1
  • F0(5) = f5
  • F0(6) = f3
  • F0(7) = f7

再帰段2

Nhが1、つまりループはi = 0のみで、ループ中の複素数計算は:

  • L1 = L0 + E(i*4)R0 = E(0)L0 + E(i*4+0)R0
  • R1 = L0 - E(i*4)R0 = E(0)L0 + E(i*4+4)R0

より、以下の計算が行われます:

  • s = 0, i = 0
    • F1(0) = F0(0) + F0(1) = E(0)f0 + E(0)f4
    • F1(1) = F0(0) - F0(1) = E(0)f0 + E(4)f4
  • s = 2, i = 0
    • F1(2) = F0(2) + F0(3) = E(0)f2 + E(0)f6
    • F1(3) = F0(2) - F0(3) = E(0)f2 + E(4)f6
  • s = 4, i = 0
    • F1(4) = F0(4) + F0(5) = E(0)f1 + E(0)f5
    • F1(5) = F0(4) - F0(5) = E(0)f1 + E(4)f5
  • s = 6, i = 0
    • F1(6) = F0(6) + F0(7) = E(0)f3 + E(0)f7
    • F1(7) = F0(6) - F0(7) = E(0)f3 + E(4)f7

この結果の値を順に並べると:

  • F1(0) = E(0)f0 + E(0)f4
  • F1(1) = E(0)f0 + E(4)f4
  • F1(2) = E(0)f2 + E(0)f6
  • F1(3) = E(0)f2 + E(4)f6
  • F1(4) = E(0)f1 + E(0)f5
  • F1(5) = E(0)f1 + E(4)f5
  • F1(6) = E(0)f3 + E(0)f7
  • F1(7) = E(0)f3 + E(4)f7

となります。

再帰段1

Nhが２，つまりループはi = 0,1で、ループ中の複素数計算は:

  • L2 = L1 + E(i*2)R1 = E(0)L1 + E(i*2)R1
  • R2 = L1 - E(i*2)R1 = E(0)L1 + E(i*2+4)R1

より、以下の計算が行われます:

  • s = 0, i = 0
    • F2(0) = F1(0) + E(0)F1(2) = E(0)f0 + E(0)f4 + E(0)f2 + E(0)f6
    • F2(2) = F1(0) - E(0)F1(2) = E(0)f0 + E(0)f4 + E(4)f2 + E(4)f6
  • s = 0, i = 1
    • F2(1) = F1(1) + E(2)F1(3) = E(0)f0 + E(4)f4 + E(2)f2 + E(6)f6
    • F2(3) = F1(1) - E(2)F1(3) = E(0)f0 + E(4)f4 + E(6)f2 + E(2)f6
  • s = 4, i = 0
    • F2(4) = F1(4) + E(0)F1(6) = E(0)f1 + E(0)f5 + E(0)f3 + E(0)f7
    • F2(6) = F1(4) - E(0)F1(6) = E(0)f1 + E(0)f5 + E(4)f3 + E(4)f7
  • s = 4, i = 1
    • F2(5) = F1(5) + E(2)F1(7) = E(0)f1 + E(4)f5 + E(2)f3 + E(6)f7
    • F2(7) = F1(5) - E(2)F1(7) = E(0)f1 + E(4)f5 + E(6)f3 + E(2)f7

この結果の値を順に並べると:

  • F2(0) = E(0)f0 + E(0)f4 + E(0)f2 + E(0)f6
  • F2(1) = E(0)f0 + E(4)f4 + E(2)f2 + E(6)f6
  • F2(2) = E(0)f0 + E(0)f4 + E(4)f2 + E(4)f6
  • F2(3) = E(0)f0 + E(4)f4 + E(6)f2 + E(2)f6
  • F2(4) = E(0)f1 + E(0)f5 + E(0)f3 + E(0)f7
  • F2(5) = E(0)f1 + E(4)f5 + E(2)f3 + E(6)f7
  • F2(6) = E(0)f1 + E(0)f5 + E(4)f3 + E(4)f7
  • F2(7) = E(0)f1 + E(4)f5 + E(6)f3 + E(2)f7

となります。

再帰段0

Nhが4、つまりループはi = 0,1,2,3で、ループ中の複素数計算は:

  • L3 = L2 + E(i*1)R2 = E(0)L2 + E(i*1)R2
  • R3 = L2 - E(i*1)R2 = E(0)L2 + E(i*1+4)R2

より、以下の計算が行われます:

  • s = 0, i = 0
    • F3(0) = F2(0) + E(0)F2(4) = E(0)f0 + E(0)f4 + E(0)f2 + E(0)f6 + E(0)f1 + E(0)f5 + E(0)f3 + E(0)f7
    • F3(4) = F2(0) - E(0)F2(4) = E(0)f0 + E(0)f4 + E(0)f2 + E(0)f6 + E(4)f1 + E(4)f5 + E(4)f3 + E(4)f7
  • s = 0, i = 1
    • F3(1) = F2(1) + E(1)F2(5) = E(0)f0 + E(4)f4 + E(2)f2 + E(6)f6 + E(1)f1 + E(5)f5 + E(3)f3 + E(7)f7
    • F3(5) = F2(1) - E(1)F2(5) = E(0)f0 + E(4)f4 + E(2)f2 + E(6)f6 + E(5)f1 + E(1)f5 + E(7)f3 + E(3)f7
  • s = 0, i = 2
    • F3(2) = F2(2) + E(2)F2(6) = E(0)f0 + E(0)f4 + E(4)f2 + E(4)f6 + E(2)f1 + E(2)f5 + E(6)f3 + E(6)f7
    • F3(6) = F2(2) - E(2)F2(6) = E(0)f0 + E(0)f4 + E(4)f2 + E(4)f6 + E(6)f1 + E(6)f5 + E(2)f3 + E(2)f7
  • s = 0, i = 3
    • F3(3) = F2(3) + E(3)F2(7) = E(0)f0 + E(4)f4 + E(6)f2 + E(2)f6 + E(3)f1 + E(7)f5 + E(1)f3 + E(5)f7
    • F3(7) = F2(3) - E(3)F2(7) = E(0)f0 + E(4)f4 + E(6)f2 + E(2)f6 + E(7)f1 + E(3)f5 + E(5)f3 + E(1)f7

この結果の値を順に並べると:

  • F3(0) = E(0)f0 + E(0)f4 + E(0)f2 + E(0)f6 + E(0)f1 + E(0)f5 + E(0)f3 + E(0)f7
  • F3(1) = E(0)f0 + E(4)f4 + E(2)f2 + E(6)f6 + E(1)f1 + E(5)f5 + E(3)f3 + E(7)f7
  • F3(2) = E(0)f0 + E(0)f4 + E(4)f2 + E(4)f6 + E(2)f1 + E(2)f5 + E(6)f3 + E(6)f7
  • F3(3) = E(0)f0 + E(4)f4 + E(6)f2 + E(2)f6 + E(3)f1 + E(7)f5 + E(1)f3 + E(5)f7
  • F3(4) = E(0)f0 + E(0)f4 + E(0)f2 + E(0)f6 + E(4)f1 + E(4)f5 + E(4)f3 + E(4)f7
  • F3(5) = E(0)f0 + E(4)f4 + E(2)f2 + E(6)f6 + E(5)f1 + E(1)f5 + E(7)f3 + E(3)f7
  • F3(6) = E(0)f0 + E(0)f4 + E(4)f2 + E(4)f6 + E(6)f1 + E(6)f5 + E(2)f3 + E(2)f7
  • F3(7) = E(0)f0 + E(4)f4 + E(6)f2 + E(2)f6 + E(7)f1 + E(3)f5 + E(5)f3 + E(1)f7

となります。

このE(n)を表にすると、

F3\ft f0 f4 f2 f6 f1 f5 f3 f7
F3(0) E(0) E(0) E(0) E(0) E(0) E(0) E(0) E(0)
F3(1) E(0) E(4) E(2) E(6) E(1) E(5) E(3) E(7)
F3(2) E(0) E(0) E(4) E(4) E(2) E(2) E(6) E(6)
F3(3) E(0) E(4) E(6) E(2) E(3) E(7) E(1) E(5)
F3(4) E(0) E(0) E(0) E(0) E(4) E(4) E(4) E(4)
F3(5) E(0) E(4) E(2) E(6) E(5) E(1) E(7) E(3)
F3(6) E(0) E(0) E(4) E(4) E(6) E(6) E(2) E(2)
F3(7) E(0) E(4) E(6) E(2) E(7) E(3) E(5) E(1)

となり、先の表と同じ値となっていることがわかります。

7. JavaScriptによるループ版FFT実装コード

上述の再帰段数ごとの計算を、要素数が任意の２のべき乗数Nとして、ループによって直接的に実装したものが、ループ版FFTとなります。

FFTとIFFTの共通処理部分をfftin1(c, T, N)としました。ただし、引数Tの値は再帰版と違い、Nで割っていません。

最内ループ内のコードでの変数(i, Nh, l, re)は、再帰版FFTでのループ処理での変数と合わせてあります。

[14]
function fftin1(c, T, N) {
    const k = Math.log2(N);
    const rec = c.map((_, i) => c[revBit(k, i)]);
    for (let Nh = 1; Nh < N; Nh *= 2) {
        T /= 2;
        for (let s = 0; s < N; s += Nh * 2) {
            for (let i = 0; i < Nh; i++) {
                const l = rec[s + i], re = imul(rec[s + i + Nh], expi(T * i));
                [rec[s + i], rec[s + i + Nh]] = [iadd(l, re), isub(l, re)];
            }
        }
    }
    return rec;
}

function fft1(f) {
    const N = f.length, T = -2 * Math.PI;
    return fftin1(f, T, N);
}
function ifft1(F) {
    const N = F.length, T = 2 * Math.PI;
    return fftin1(F, T, N).map(([r, i]) => [r / N, i / N]);
}

最外ループでは、Nhを倍々にしていくことから、log2(N)回実行されます。
その中のループは２つ合わせて、中の計算部分をN/2回実行するものとなります。

この実装で、再帰版FFTのときと同じfr0を使って計算してみましょう。

[15]
{
    const fr0 = [1,3,4,2, 5,6,2,4, 0,1,3,4, 5,62,2,3];
    const f0 = fr0.map(r => [r, 0]);

    const F = fft1(f0);
    const f1 = ifft1(F);
    const fr1 = f1.map(([r]) => r);

    console.log("fr0:", fr0);
    console.log("F:", F);
    console.log("f1:", f1);
    console.log("fr1:", fr1.map(Math.round));
}

fr0: [ 1, 3, 4, 2, 5, 6, 2, 4, 0, 1, 3, 4, 5, 62, 2, 3 ]
F: [ [ 107, 0 ],
  [ 23.29589166141268, 51.729855807372815 ],
  [ -53.5477272147525, 42.961940777125584 ],
  [ -49.21391810443094, -25.67438445589562 ],
  [ 3.6127080574846916e-15, -59 ],
  [ 49.799704542057846, -24.260170893522517 ],
  [ 35.54772721475249, 48.96194077712559 ],
  [ -19.8816780990396, 53.14406936974591 ],
  [ -63, 0 ],
  [ -19.881678099039586, -53.14406936974591 ],
  [ 35.5477272147525, -48.961940777125584 ],
  [ 49.799704542057846, 24.260170893522528 ],
  [ -3.6127080574846916e-15, 59 ],
  [ -49.21391810443094, 25.67438445589561 ],
  [ -53.54772721475249, -42.96194077712559 ],
  [ 23.295891661412693, -51.729855807372815 ] ]
f1: [ [ 1, 0 ],
  [ 3.000000000000001, -1.3877787807814457e-15 ],
  [ 3.999999999999999, -1.0143540619928357e-16 ],
  [ 2, -3.0531133177191805e-16 ],
  [ 5, 0 ],
  [ 6, -1.7763568394002505e-15 ],
  [ 2, 4.592425496802574e-17 ],
  [ 3.9999999999999996, -2.498001805406602e-16 ],
  [ 0, 0 ],
  [ 0.9999999999999991, 1.3877787807814457e-15 ],
  [ 3.000000000000001, 9.586896263232085e-18 ],
  [ 4.000000000000001, 3.608224830031759e-16 ],
  [ 5, 0 ],
  [ 62, 1.7763568394002505e-15 ],
  [ 2, 4.592425496802574e-17 ],
  [ 2.9999999999999996, 1.942890293094024e-16 ] ]
fr1: [ 1, 3, 4, 2, 5, 6, 2, 4, 0, 1, 3, 4, 5, 62, 2, 3 ]

結果は、再帰版FFTのときと同じ値になりました。

8. まとめ

DFTの実装の理解を前提に、再帰版FFTの実装を提示し、再帰版FFTの振る舞いを具体例で解析することで、FFTの仕組みを解明し、最後にその仕組をループで直接実行したループ版FFTを作成しました。

再帰版FFTの前半部の解析では、インデックスの2進数表現を用いました。

内部で２度自身を呼び出す再帰関数は、O(NlogN)な処理として、葉がN要素の完全二分木のデータ構造への処理を行っているとみなせます。
この隠れた二分木を表出させるのが、再帰で扱う数値の２進数表現になります。

たとえば、ハノイの塔は、この内部で二度自身を呼び出す再帰関数になります。
これも２進数表現化を応用することで、固定段のループによって実装させることができます。

付録: 任意長データのためのFFT

２のべき乗のFFT/IFFTを用いて、任意長のデータを扱うFFTを作ることができます。
この任意長FFTは、Bluestein's algorithmと呼ばれています。

２つの数列(関数)をかけ合わせて総和(積分)を取ることを、畳み込み(convolution)といいます。
フーリエ変換やラプラス変換は、べき乗タイプ($\exp(C x)$)の数列(関数)との畳み込みを行っています。
このタイプの畳み込みによる変換では、２つの数列(関数)に対して、変換した後同士の積をとり、その逆変換をすると、もとの2数列(関数)の畳み込みになっている、という性質があります。

Bluestein's algorithmでは、DFTの$\exp(\frac{-2\pi i}{N}n k)$を分解して、その畳み込みに変形します。
この変形後の畳み込みは、2のべき乗長にデータを拡大するのに適した形になります。

まずDFTの畳み込みの式が、リンク先での$a_n$と$b_n$との畳込みの式へと変換されます。
そこから、この$a_n$と$b_n$を、2のべき乗長に拡大します。
$a_n$は後ろを単に０で埋めるだけですが、$b_n$のほうは後ろ側に$b_n$を反転したものを入れます。
これは元からの$a_n$との畳み込みがされる範囲を埋めています。
(先頭の畳込みは$ab_0 = a_0 b_0 + a_1 b_{-1} + a_2 b_{-2} + \dots + a_{N-1} b_{1-N}$であり、この$b_{-1} \dots b_{1-N}$が拡大した$b_n$の後ろ側に必要)。

この２のべき乗長に拡大された$a_n$と$b_n$の畳み込みを、2のべき乗長のFFTとIFFTを使って行い、その結果から、先頭からもとの長さだけ切り出し、$b$の複素共役とかけることで、DFTでの結果とほぼ同じものになる、という仕組みです。

このBluestein's algorithmのコードは以下のようになります(fft2(f)、ifft2(F)):

[16]
function conj([x, y]) {return [x, -y];}

function fftin2(c, T, N) {
    const Nd = N * 2;
    const bc = c.map((_, n) => expi(T * n * n / Nd));
    const b = bc.map(conj);
    const a = c.map((cn, n) => imul(cn, bc[n]));

    //const N2 = 1 << Math.ceil(Math.log2(Nd - 1));
    const N2 = 1 << (32 - Math.clz32(Nd - 1));
    const a2 = a.concat(Array(N2 - N).fill([0, 0]));
    const b2 = b.concat(Array(N2 - Nd + 1).fill([0, 0]), b.slice(1).reverse());
    const A2 = fft1(a2);
    const B2 = fft1(b2);
    const AB2 = A2.map((A2n, n) => imul(A2n, B2[n]));
    const ab2 = ifft1(AB2);
    return bc.map((bcn, n) => imul(bcn, ab2[n]));
}

function fft2(f) {
    const N = f.length, T = -2 * Math.PI; 
    return fftin2(f, T, N);
}
function ifft2(F) {
    const N = F.length, T = 2 * Math.PI;
    return fftin2(F, T, N).map(([x, y]) => [x / N, y / N]);
}
[17]
{
    // 15 elements example
    const fr0 = [1,3,4,2, 5,6,2,4, 0,1,3,4, 5,62,2];
    const f0 = fr0.map(r => [r, 0]);

    const F = fft2(f0);
    const f1 = ifft2(F);
    const fr1 = f1.map(([r]) => r);

    console.log("fr0:", fr0);
    console.log("F:", F);
    console.log("f1:", f1);
    console.log("fr1:", fr1.map(Math.round));
}

fr0: [ 1, 3, 4, 2, 5, 6, 2, 4, 0, 1, 3, 4, 5, 62, 2 ]
F: [ [ 104.00000000000001, 1.4210854715202004e-14 ],
  [ 40.113068713053366, 40.536802661213386 ],
  [ -16.938440097404285, 64.08647185352932 ],
  [ -47.041019662496865, 29.02599135062092 ],
  [ -57.94588444305084, -15.35101781191083 ],
  [ -35.500000000000064, -52.82754963085073 ],
  [ 20.04101966249682, -49.09166758334534 ],
  [ 52.77125582740197, -26.917243683390566 ],
  [ 52.771255827401916, 26.91724368339068 ],
  [ 20.041019662496826, 49.091667583345306 ],
  [ -35.49999999999999, 52.82754963085077 ],
  [ -57.94588444305086, 15.351017811910662 ],
  [ -47.04101966249681, -29.02599135062106 ],
  [ -16.938440097404428, -64.08647185352933 ],
  [ 40.11306871305353, -40.53680266121326 ] ]
f1: [ [ 1.0000000000000246, -8.763360407707903e-15 ],
  [ 3.0000000000000235, -1.4092430925908654e-14 ],
  [ 4.000000000000012, -6.158037043254202e-15 ],
  [ 2.000000000000004, 6.2764608325475514e-15 ],
  [ 4.999999999999981, 5.921189464667502e-15 ],
  [ 5.999999999999992, -2.0368891758456206e-14 ],
  [ 2.0000000000000098, -8.052817671947802e-15 ],
  [ 4.000000000000021, 9.473903143468002e-16 ],
  [ 1.1817018103590317e-14, 9.119920035791959e-15 ],
  [ 0.9999999999999932, 1.0835776720341527e-14 ],
  [ 2.9999999999999916, -1.4210854715202004e-14 ],
  [ 3.999999999999985, 5.8027656753741514e-15 ],
  [ 4.9999999999999964, 8.763360407707903e-15 ],
  [ 62.00000000000001, 7.579122514774402e-15 ],
  [ 2.000000000000005, 1.4802973661668755e-14 ] ]
fr1: [ 1, 3, 4, 2, 5, 6, 2, 4, 0, 1, 3, 4, 5, 62, 2 ]

比較対象に、同じデータのDFTでの結果も載せておきます:

[18]
{
    // 15 elements example
    const fr0 = [1,3,4,2, 5,6,2,4, 0,1,3,4, 5,62,2];
    const f0 = fr0.map(r => [r, 0]);

    const F = dft(f0);
    const f1 = idft(F);
    const fr1 = f1.map(([r]) => r);

    console.log("fr0:", fr0);
    console.log("F:", F);
    console.log("f1:", f1);
    console.log("fr1:", fr1.map(Math.round));
}

fr0: [ 1, 3, 4, 2, 5, 6, 2, 4, 0, 1, 3, 4, 5, 62, 2 ]
F: [ [ 104, 0 ],
  [ 40.1130687130533, 40.53680266121345 ],
  [ -16.93844009740454, 64.08647185352929 ],
  [ -47.04101966249686, 29.02599135062094 ],
  [ -57.945884443050794, -15.351017811910953 ],
  [ -35.49999999999983, -52.827549630850896 ],
  [ 20.041019662496918, -49.09166758334529 ],
  [ 52.771255827402115, -26.917243683390264 ],
  [ 52.77125582740182, 26.917243683390836 ],
  [ 20.041019662496755, 49.09166758334532 ],
  [ -35.50000000000034, 52.82754963085049 ],
  [ -57.94588444305084, 15.351017811910584 ],
  [ -47.04101966249679, -29.025991350621084 ],
  [ -16.938440097404428, -64.08647185352936 ],
  [ 40.113068713053934, -40.536802661212924 ] ]
f1: [ [ 1.0000000000000286, 1.0894988614988204e-14 ],
  [ 3.000000000000027, 1.8947806286936004e-14 ],
  [ 4.000000000000055, 1.5158245029548804e-14 ],
  [ 2.0000000000000298, -6.015928496102182e-14 ],
  [ 5.000000000000037, -8.052817671947802e-15 ],
  [ 6.000000000000018, -5.6369723703634616e-14 ],
  [ 1.9999999999999796, -6.987003568307652e-14 ],
  [ 3.9999999999999383, -3.789561257387201e-15 ],
  [ -1.7289873236829103e-14, -1.5158245029548804e-14 ],
  [ 0.9999999999999837, -1.8947806286936004e-14 ],
  [ 2.999999999999985, 1.8000415972589204e-14 ],
  [ 3.999999999999993, 3.600083194517841e-14 ],
  [ 5.000000000000014, -1.3026616822268503e-14 ],
  [ 62, -2.2737367544323207e-14 ],
  [ 1.999999999999967, -1.3263464400855204e-14 ] ]
fr1: [ 1, 3, 4, 2, 5, 6, 2, 4, -0, 1, 3, 4, 5, 62, 2 ]

データ列の長さが2のべき乗のときはfftin1を呼び、そうでない場合のみfftin2を呼ぶようにしたfft/ifftは以下のようになります。

[19]
function isPowerOf2(n) {return (n & (n - 1)) === 0;}
function fft(f) {
    const N = f.length, T = -2 * Math.PI;
    const fftin = isPowerOf2(N) ? fftin1 : fftin2;
    return fftin(f, T, N);
}
function ifft(F) {
    const N = F.length, T = 2 * Math.PI;
    const fftin = isPowerOf2(N) ? fftin1 : fftin2;
    return fftin(F, T, N).map(([x, y]) => [x / N, y / N]);
}

参考: このドキュメントの編集環境

このドキュメントは、jupyter-notebook形式で記述したもので、ijavascriptカーネルを用いて、nodejsによってドキュメント内のコードが実行できます。
エディタ環境として、ブラウザを用いるJupyterLabを使いました。

jupyterインストール:

$ pip3 install --upgrade jupyter

ijavascriptインストールとnodejsカーネルのインストール(macの場合):

$ brew install zeromq
$ npm install -g ijavascript --zmq-external
$ ijsinstall

注: zeromqのネイティブモジュールを含むので、nodejsを更新した場合、再度npm installする必要がでるかもしれません。

JupyterLabインストールと実行:

$ pip3 install --upgrade jupyterlab
$ jupyter-lab

成功すればブラウザタブが開き、ipynbファイルの編集を行えます。

追加: ipynbからqiitaのmarkdownにしてqiitaへ貼り付ける

この本文を記述したipynbは以下のgistへ置いてあります。

jupyter-nbconvertコマンドを使えば、ipynb形式のファイルをmarkdown形式のファイルへ変換できます:

$ jupyter-nbconvert --to markdown understanding-fft.ipynb
[NbConvertApp] Converting notebook understanding-fft.ipynb to markdown
[NbConvertApp] Writing 25680 bytes to understanding-fft.md

macosでは、このmarkdownファイルをpbcopyコマンドでクリップボードに載せ、ブラウザ上で貼り付けられます:

$ cat understanding-fft.ipynb | pbcopy

ただし、いくつかの点でqiitaとjupyter-notebookでのmarkdownには、解釈に違いがあります:

  • 改行の解釈: qiitaは改行ごとに折り返されるが、jupyter-notebookは空行で改段落になる
  • 表の解釈: qiitaでは行頭の|の前に空白を入れられない(装飾可能なテキストブロック扱いになる)が、jupyter-notebookでは行頭の|の前に空白を入れても表として処理される
  • $$な数式ブロックの解釈: qiitaでは$$なブロックでは\begin{align}環境は解釈されず数式扱いもされなくなるが、jupyter-notebookでは解釈される

この対応のために、ipynbのmarkdownで使った$$な数式ブロック部分は、mathコードブロックに書き換えるなど、の後処理が必要となるでしょう。

ここではさらに、javascriptコードブロックに、qiitaのコードファイル名表示機能を使って、実行順の:[0]をつけておきました。

44
32
3

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
44
32