コンピューターの中の音は、1秒あたり44100個ほど並んだ数の列 です。波形として眺めれば、時間に沿って上下する一本の線にすぎません。ところが音楽プレーヤーのスペクトラム表示には、低い音から高い音まで棒が並びます。同じ音なのに、どうやって一本の線を「どの周波数がどれだけ含まれるか」という成分表へ変えているのでしょうか。
その計算が離散フーリエ変換（DFT）です。「フーリエ変換」と聞くと難しい数式が先に立ちますが、中心でしていることは、調べたい周波数の波を掛けて、全部足すことです。
この記事では、この「掛けて足す」だけに焦点を絞ります。なぜ一致する周波数だけが残るのか、なぜ cos と sin の2本が要るのか、そして周波数の目盛りがなぜ飛び飛びになるのかを、図と操作できるデモで確かめていきます。
掛けて足すと、なぜ成分が出るのか
まず、測りたい波形の各サンプルに、調べたい周波数の cos 波を掛けます。この照らし合わせ用の波を、ここでは「探り波」と呼ぶことにします。
測りたい波形と探り波の周波数が一致していれば、山と山が同じ場所に来るので、正×正で積は正になります。谷と谷でも、負×負なので積は正です。各サンプルの積を足していくと、正の値が同じ側へ積み上がります。
周波数が違えば、山と山が重なる場所もあれば、山と谷が重なる場所もあります。積は正と負を行き来し、全部足すころには互いに打ち消し合ってゼロに近づきます。
これがDFTの核です。入力の N 個のサンプルを x[n]、調べる周波数の目盛りを k とすると、ひとつの k について x[n] と探り波を N 回掛け、N 個の積を足します。k を変えながら同じ照らし合わせを繰り返せば、倍音の記事
で足し算の側から見た成分表を、今度は分解の側から作れることになります。
cos だけでは取り逃す ── sin と対で測る
ただし、cos の探り波1本だけでは困る場合があります。測りたい成分が同じ周波数でも、cos からちょうど90度ずれていたらどうなるでしょう。波は確かに含まれているのに、cos と掛けた総和はゼロになります。周波数は合っていても、山と谷の位置、つまり位相がずれているからです。
そこで、cos と90度ずれた sin でも同時に測ります。cos との総和を実部、sin との総和を虚部として、次のように2方向を合わせます。
$$|X[k]| = \sqrt{\operatorname{Re}(X[k])^2 + \operatorname{Im}(X[k])^2}$$位相0度では cos 側がすべてを受け持ち、90度では sin 側がすべてを受け持ちます。45度では両者が約0.707ずつです。受け持ち方は変わっても、2乗和の平方根を取った大きさは常に1になります。
cos と sin を別々の波と見る代わりに、複素平面上で回る一本の矢印の横成分と縦成分だと考えることもできます。DFTでは、この回る矢印を回転因子 W として次の形で使います。
つまり W の正体は、cos と sin の探り波をひと組にしたものです。サンプル番号 n が進むごとに矢印が回り、目盛り k に対応する周波数で入力と照らし合わせます。変換後の値 X[k] が複素数になるのは、成分の大きさだけでなく位相も取り逃さないためです。
測れる周波数は飛び飛びになる
DFTに渡す N 点は、連続する音の一部分を四角く切り出した窓です。この窓に1周期、2周期、3周期と整数個ぴったり収まる探り波を並べると、周波数の目盛り k ができます。サンプリング周波数を fs とすれば、隣り合う目盛りの間隔は fs/N、目盛り k の周波数は k·fs/N です。
たとえば fs=44100、N=1024 なら、目盛りの間隔は 44100/1024 = 43.06640625 Hz です。DFTが「離散」と呼ばれるのは、入力が飛び飛びのサンプルであるだけでなく、調べる周波数もこのような飛び飛びの目盛りになるからです。
窓に2.5周期のような半端な回数しか入らない成分は、どの目盛りとも完全には一致しません。そのため値がひとつの X[k] に集まらず、両隣を含む複数の目盛りへ漏れます。ここでは、有限の長さで音を切り出すとこうした事情が生まれる、とだけ覚えておけば十分です。
直交性 ── だから成分が混ざらない
窓に整数個ぴったり入る探り波には、もうひとつ大切な性質があります。異なる目盛りの波どうしを掛けて N 点ぶん足すと、正と負が正確に打ち消し合うのです。この関係を直交と呼びます。
$$\sum_{n=0}^{N-1} \cos\left(\frac{2\pi k_1n}{N}\right) \cos\left(\frac{2\pi k_2n}{N}\right) = 0 \qquad (k_1 \ne k_2,\ \ 0 < k_1, k_2 < N/2)$$実際に N=64、k₁=3、k₂=5 で計算すると、総和は 10⁻¹⁵ 程度の値になります。これは計算機が小数を扱うときの丸め誤差で、ゼロと見てよい大きさです。一方、同じ目盛りどうし（k₁=k₂=3）なら 32 が残ります。理論値の N/2 そのものです。
異なる目盛りは互いの測定へ入り込まず、同じ目盛りだけが値を残します。だから X[k] はそれぞれ独立した測定結果になり、スペクトル、つまり音の成分表として意味を持ちます。
式に N/2 までという但し書きを付けたのには理由があります。cos は k と N−k で見分けがつきません。N=64 なら k=3 と k=61 はまったく同じ波になり、掛けて足すと 32 が残ってしまいます。つまり目盛りの後半は前半の鏡像です。実際の音のような実数の入力では、成分表の後半は前半の折り返しでしかなく、意味を持つ目盛りは半分だけになります。スペクトラム表示が標本化周波数の半分までしか描かれないのは、エイリアシングの記事
で見た事情に加えて、この折り返しがあるためです。
探り波を動かしてみる
下のデモでは、既定の「倍音入り」に 220、440、660、1100 Hz の成分が入っています。まず再生し、パートAの探り周波数をゆっくり動かしてください。含まれる周波数へ近づくと「大きさ」のバーが跳ね上がり、外れると正負の積が打ち消し合います。「探り波も小さく鳴らす」を使うと、耳でも2つの周波数の関係を確かめられます。
パートBの「1ステップ進める」は、目盛り k ごとの掛け算と足し算を一回ずつ進めます。最後まで埋めれば、選んだ音の倍音構成と同じ成分表になります。パートCでは位相を回し、cos と sin の値が変わっても「大きさ」だけが動かないことを確認できます。
うまく表示されないときは、デモを別タブで開く 。
どう実現するか
cos と sin の照らし合わせを回転因子 W にまとめると、DFTの定義は一行で書けます。
プログラムへ直しても、していることは式のままです。
for (let k = 0; k < N; k++) {
let re = 0, im = 0;
for (let n = 0; n < N; n++) {
const angle = 2 * Math.PI * n * k / N;
re += x[n] * Math.cos(angle);
im -= x[n] * Math.sin(angle);
}
X[k] = { re, im, magnitude: Math.hypot(re, im) };
}
この定義では、振幅1の余弦波が目盛り k₀ にぴったり乗ると、|X[k₀]| は N/2 になります。N=64、k₀=5 で実際に計算すると 32 でした。元の振幅へ戻して表示したければ 2/N を掛ければよく、2×32/64 でちょうど1に戻ります。ただし折り返しの相手を持たない両端、つまり直流（k=0）といちばん上の目盛り（k=N/2）だけは、2倍せず 1/N を掛けます。この2本は鏡像の片割れではなく、それ自体で完結しているからです。
同じ波の位相を0度、45度、90度へ回すと、実部と虚部の受け持ちは (32, 0)、(22.63, 22.63)、(0, 32) と移り変わります。それでも |X[5]| は3通りとも32のままです。cos と sin を対にしておけば、位相がどこにあっても成分の大きさを取り逃さない。数値の上でもそのとおりになります。
ここまでが、音を成分表へ変える計算のすべてです。ただしこの手順どおりに計算すると、N が大きくなったとき手数が耐えがたいほど増えます。同じ答えを桁違いに少ない手数で出す段取りが次回の主題です。
連載「サウンドプログラミング」
← 前回 ノイズを止めた瞬間、鳴っていない音が残る