COURSE 01 · AUDIO SIGNAL PROCESSING
「音とは何か」から MFCC まで
関連する予備知識がまったくない人に向けた講座です。全体を通して 1 つのプロジェクトに取り組みます：1 秒間の音楽をプログラムに聴かせて、クラシック・ジャズ・ロックのどれかを判定させます。
23 講247 枚の図8 種類の実験23 本のスクリプト
講座の本文は現在、中国語版のみです。各講のタイトルと導入文、メニューやボタンは日本語で表示しています。
00始める前に
01—05音とデジタル録音を知る
- 01プログラムがクラシック・ジャズ・ロックを聞き分けるには、まずどの一歩を越えなければならないのか？プログラムに猫の写真を渡すと、受け取るのはきちんと並んだ色の格子です。車が通り過ぎる音を渡すと、受け取るのは時間とともに次々に出てくる長い数字の列です。どちらも「分類」と呼びますが、音のほうには工程が 1 つ多くあります。まず誰かがその数字の列を、互いに比べられる証拠に整理しておかないと、プログラムには比べるものがありません。この工程を音声信号処理と呼び、この 23 講で取り組むことです。この講では具体的なアルゴリズムは扱わず、3 つのことだけを決めます。プログラムにどんな課題をこなさせるのか、23 講をどんな順序で進むのか、どんな予備知識が必要なのか。
- 02音はどうやって生まれ、どうやって録音の中のあの線になるのか？前の講で課題が決まりました。プログラムにクラシック・ジャズ・ロックを聞き分けさせることです。でもプログラムが受け取るのは長い数字の列です。この数字はどこから来るのでしょうか。この講はいちばん最初から話します。物体が振動すると周りの空気を押し、その押し合いが輪のように広がってマイクに届きます。マイクは一瞬ごとの押す強さを 1 つの数として記録します。その数を時間順に描いたものが、録音ソフトに出てくる震える線で、波形と呼びます。この線には 3 つのことが同時に隠れています。どれだけ速く振動するか、どれだけ大きく振動するか、そしてどんな音色に聞こえるか。この講では前の 2 つを説明し、音色は第 03 講に回します。
- 03音量バーは動いていないのに、曲を変えるとはっきりうるさくなるのはなぜ？前の講で波形の周波数と振幅を説明しましたが、まだ 2 つ残っています。音の強さはそもそもどう測るのか、そして音色とは何か。この講では音源から出発して、一段ずつ耳へ近づいていきます。音源が 1 秒間に出すエネルギー（音響パワー）→ あなたのいる位置まで届いたときの 1 平方メートルあたりの量（音の強さ）→ やっと聞こえる音から痛みを感じる音まで、10 兆倍もの差がある → だからデシベルを使って読みやすい数に縮める → それでもデシベルは、あなたが感じる音の大きさ（ラウドネス）とは同じではない。最後に「音量は同じなのに違う」という問いに答えます。高さと大きさをそろえても残る違いを音色と呼び、それは 3 つの要素で決まります。
- 04連続した音は、どうやって数字の列になり、また音に戻るのか？ここまでの 2 講で扱ったのは、どれも連続して変化するものでした。空気の圧力、周波数、音の大きさ。でもコンピューターはどれもそのままでは保存できません。無限に多くの時刻も、無限に多くの値も記録できないからです。この講では変換の流れを最後までたどります。マイクが連続して変化する電圧を出す → ADC という部品がそれを数字の列に変える → あとは自由に計算する → 最後に DAC という部品が数字を連続した電圧に戻し、スピーカーを動かす。ADC がすることは 2 つだけです。どれくらいの間隔で測るか（サンプリング）、1 回ごとにどれくらい細かく測るか（量子化）。この 2 つにはそれぞれ代償があり、その性質はまったく違います。一方はあとから補えますが、もう一方は二度と取り戻せません。
- 05録音から何を計算してプログラムに渡すべきか？——音声特徴量の 5 つの軸前の講までで、30 秒の音楽は 60 万個あまりの数字の列になりました（1 秒に 2 万 2 千回あまり、それを 30 秒分）。この数字の列をそのまま分類プログラムに渡してもうまくいきません。長すぎるし、細かすぎるし、1 つの数だけでは何もわからないからです。そこでまず、いくつかの要約を計算します。この要約を音声特徴量と呼びます。「どの音声特徴量を使えばいいか」と検索すると長い用語リストが出てきますが、それぞれの名前が答えている問いは同じではありません。この講ではリストを暗記するのではなく、分類の地図を示します。どんな音声特徴量も、5 つの軸で位置を決められます。抽象度、時間の範囲、音楽のどの性質か、信号の領域、そしてどんな方法で計算するか。
06—10録音を計算できる小さな区間に分ける
- 06録音全体から、比べられる数字の列になるまでに、どんな段階があるのか？前の講では「何を計算するか」を 5 つの軸に分けました。そのうち信号の領域の軸は、時間領域・周波数領域・時間周波数領域の 3 段階です。この講では前の 2 つについて、処理の流れを最初から最後まで説明します。ADC から出てきた数字の列から始まり、最後にプログラムに渡す数個の数字までのあいだに、どんな段階があり、なぜどれも欠かせないのか。2 つの流れは前半が同じで（どちらもまず区切る）、分かれ道は周波数領域のほうにあります。そちらには「窓をかける」段階が 1 つ多いのです。そして窓をかけることで新しい問題が生まれ、「フレームは重ならなければならない」という結論が導かれます。
- 07周波数の分析をせず、あの曲線だけを見て何がわかるのか？前の講で処理の流れはできましたが、「フレームごとに計算する」の欄がまだ空いています。この講では時間領域の 3 つの答えでそこを埋めます。振幅包絡（このフレームの最大値はいくつか）、RMS＝二乗平均平方根（このフレーム全体の強さはどれくらいか）、ゼロ交差率（このフレームで中心線を何回横切るか）。どれも周波数の分析はいらず、曲線そのものを見るだけです。この 3 つは入れ替えのきく名前ではなく、3 つの違う問いです。それぞれに独自の公式、弱点、使いどころがあります。この講では 3 つの公式を、8 個の数からなる小さなフレームでそれぞれ手計算します。ペンがあれば自分で確かめられます。
- 08振幅包絡の公式を、本当に動くコードにする前の講で 3 つの公式を説明しましたが、実装はまだ 1 行も書いていません。この講では最初の 1 つ、振幅包絡を一から作ります。3 つの音楽を読み込み、基本情報を確認し、波形を描き、包絡を自分で書き、フレーム番号を秒に換算し、最後に包絡を波形に重ねて 3 つのジャンルを比べます。途中で、見落としやすい abs（絶対値）に気づきます。曲全体の最大値ではほとんど影響が見えないのに、フレームごとに調べると、3 つの音楽でフレームの 18% から 36% が低く見積もられてしまいます。
- 09まずライブラリで計算し、次に自分で書き、両方を一致させるこの講では残り 2 つの時間領域の特徴量を実装します。やり方は前の講と違います。まず librosa で 1 行で計算し、次に一から自分で書き、2 つの結果をフレームごとに突き合わせます。これは形だけの作業ではありません。値が合わなければ、「1 フレームがどのサンプルを含むか」についてあなたとライブラリの理解が違うということで、その先の比較はすべて成り立ちません。RMS は両者がぴったり一致します（フレームごとの最大差 8×10⁻⁹）。ゼロ交差率はわざと一致させず、比はちょうど 1024/1023 になります。これこそ第 07 講で残した「K で割るか、K−1 で割るか」の問題です。
- 10コンピューターは 1 本の曲線から、どんな高さの成分でできているかをどうやって見つけるのか？ここまでの 4 講で、「曲線だけを見る」方法で問えることはすべて問いました。でも、それでは決して答えられないこともあります。ピアノの鍵盤を 3 つ同時に押しても、マイクが記録するのは起伏する 1 本の曲線だけで、波形からは 3 つの音が入っていることがまったくわかりません。この講では視点を変えます。時間とともにどう起伏するかではなく、どんな高さの成分でできているかを問います。やり方は意外なほど素朴です。周波数のわかっている波を用意して、1 点ずつ掛け合わせ、平均をとります。合っていれば平均はゼロにならず、合っていなければプラスとマイナスが打ち消し合ってゼロになります。すべての周波数で試して得られる曲線が、スペクトルです。
11—15波形から「時間と周波数」の表し方へ
- 111 つの数に、「どれくらい強いか」と「どこから始まるか」を同時に入れるには？前の講には宿題が残りました。探りの波 1 本だと、始まる位置のずれで結果が消えてしまうので、2 本必要になる、ということです。そのため周波数ごとに得られるのは 1 つではなく 2 つの数です。強さと、どこから始まるか。ところが「強さ」は普通の数で、普通の数 1 つに 2 つの情報は入りません。この講で導入する複素数は、まさにこのためのものです。複素数は平面上の 1 つの点で、点はもともと 2 つの情報を持っています。原点からどれだけ離れているか、どの方向にあるか。前者が強さ、後者が始まる位置に対応します。この講では抽象的な数学は扱わず、前の講で測った 2 つの読み取り値を、実際に 1 つの複素数にまとめます。
- 12「どれくらい強いか」と「どこから始まるか」を、1 つの係数に書き込む前の講で複素数の準備ができました。絶対値が強さを、角度が始まる位置を表します。この講ではそれを正式にフーリエ変換に組み込みます。周波数ごとに 1 つの複素数の係数を計算し、その絶対値をとれば振幅スペクトル、角度をとれば位相スペクトルになります。見た目はこわい公式も、分解すればやっていることは第 10 講と同じです。波を掛けて、平均をとる。ただオイラーの公式によって、2 本の探りの波が e の指数 1 つにまとめられているだけです。最後に、変換して元に戻す往復を一通りやります。1 点ごとの最大誤差は 2.5×10⁻¹³ です。ところが位相を捨ててから戻すと、誤差は 2.05 になります。8 兆倍の違いです。
- 13有限個の数字しか持たないコンピューターで、どうやってフーリエ変換をするのか？前の講の公式は、時間が連続していて、無限に多くの周波数を調べられることを前提にしていました。実際のコンピューターが持っているのは有限個のサンプルだけです。公式を本当に動かすには、時間と周波数をそれぞれ 1 回ずつ制限する必要があります。N 個のサンプルを残し、計算する周波数の区切り（ビン）も N 個だけにします。こうして得られる計算を離散フーリエ変換、略して DFT と呼びます。
- 14Python で実際の音のスペクトルを読む：4 つの音はなぜ違って見えるのか？同じ音を弾いても、バイオリン、ピアノ、サックスはやはり違って聞こえます。ノイズはまた別の姿をしています。この講では 4 つの録音を Python に渡し、周波数ごとの強さの分布を比べ、横軸を本物の Hz として読めるようにします。
- 15曲全体のスペクトルでは、音がいつ鳴ったかがわからない：短時間フーリエ変換2 つの音がどちらも低い音と高い音を順に鳴らし、順番だけが逆だったとします。曲全体の周波数の結果は、まったく同じに見えることがあります。「どちらの音が先に出たか」をプログラムに取り戻させるには、録音を互いに重なる小さな区間に分け、位置を移すたびに計算し直す必要があります。
16—20スペクトルをモデルの入力に整える
- 16同じ数字なのに、描き方を変えないと見えない：Python でスペクトログラムを描く前の講で、録音を「どの高さの成分が、どの時間に、どれくらい強いか」という大きな表にしました。この講ではそれを実際に描きます。ところが最初に描いた図はほとんど真っ黒です。計算を間違えたのではなく、色のつけ方がよくないのです。この表のマス目の 99.70% は、いちばん強いマスの 100 分の 1 にも届きません。比例のまま色を割り当てると、それらはすべて同じ黒に押し込まれます。デシベルで色をつけ、縦軸を倍数の目盛りに変えると、同じ数字がすべて見えるようになります。
- 17聞こえ方に合わせて、周波数の物差しを刻み直すには：メル尺度前の講で縦軸を倍数の目盛りに変えたのは、音楽では 1 オクターブが周波数 2 倍にあたるからでした。でも人の耳は、完全に倍数どおりに聞いているわけではありません。これを示す古い実験があります。C2 から C4 に上がるときと、G6 から A6 に上がるときでは、振動の速さの差は 196 と 192 で、ほとんど同じです。ところが前者は 2 オクターブ上がったように聞こえ、後者は全音 1 つ分にしか聞こえません。同じ差なのに、聞こえ方は 12 倍違います。この講では周波数の軸を、聴覚の実験から刻んだ物差しに取り替え、さらに 10 個の三角形で 1025 行を 10 行に圧縮します。
- 18形は合っているのに、数字がまったく違う：Python でメルスペクトログラムを取り出す前の講では、メルフィルターの一式（10 個の三角形で、それぞれが一定の周波数帯だけを担当し、真ん中がいちばん敏感で、両端に向かって次第に効かなくなる）を 5 つの手順で自作しました。この講では既製のライブラリにも作らせ、2 つの結果を並べます。形はまったく同じなのに、数字は 1 つも合いません。違いを追いかけていくと、3 つの初期設定が見つかります。ライブラリが既定で使う最高周波数は 8000 ではないこと、既定で各三角形を面積が等しくなるように縮めていること、そもそも同じメルの公式を使っていないこと。この 3 か所を直すと、最大差は 0.9999 から 0.0147 に下がります。残ったわずかな差にも、きちんと理由があります。
- 19「誰が声を出しているか」と「どんな音を出しているか」を分ける：MFCC とは何か人の声には 2 つのものが同時に入っています。声帯がどれくらい速く振動するか（音の高さを決める）と、口がどんな形になっているか（「あ」なのか「い」なのかを決める）。スペクトル上では両者は掛け算の関係にあり、絡み合っていてどちらも取り出せません。この講でやることは 1 つだけです。対数をとると、掛け算が足し算になります。足し算になった 2 つは、もう一度変換すれば 1 本の軸の両端に分かれます。この講ではその道のりの各段階を測ります。包絡の 98.3% が左側に、励振の 99.0% が右側にあり、1 か所で切れば分けられます。最後の DCT は、隣り合うメル帯どうしの相関 0.95 を 0.15 まで下げ、最初の 13 個の係数に元のスペクトルの変化の 98.2% が収まります。
- 201 行のコードが代わりにやっている 4 つの手順：Python で MFCC を取り出す前の講で MFCC の流れ全体を手で書きました。この講ではそれをライブラリに任せます。librosa.feature.mfcc は 1 行で 13 行の係数を返し、librosa.feature.delta がさらに 1 次と 2 次の変化率を出します。つなげれば、よく言われる 1 フレームあたり 39 次元になります。いつものように、1 行で結果を得たら最初にやるのは、それを前の講で手書きした 4 つの手順に分解し直して、桁まで突き合わせることです。4 つの初期設定をすべて正しく直して、ようやく最大差が 0.0000 になります。途中で 2 つのことにぶつかります。delta は「次のフレームから前のフレームを引く」ではないこと。そして 13 個の係数のどれ 1 つをとっても、単独では 3 つのジャンルを分けられないことです。
21—23少ない数値でスペクトルをまとめる
- 21スペクトルの 1 列を 3 つの数にまとめるには：帯域エネルギー比、スペクトル重心、帯域幅スペクトログラムの 1 列には、千を超える周波数のマスがあります。情報はそろっていますが、3 つの簡単な問いにそのまま答えるのには向いていません。高い周波数に比べて低い周波数がどれくらいあるか？ 周波数の重心はどこか？ 周波数は重心のまわりにどれくらい広がっているか？ この講では 3 つの問いに 1 つずつ数を割り当てます。帯域エネルギー比、スペクトル重心、スペクトル帯域幅です。さらに答えのわかっているスペクトルを 3 列作り、この 3 つの数がそれぞれ別の側面を見ていて、どれも他の代わりにはならないことを確かめます。
- 22Python で帯域エネルギー比を実装する：2000 Hz はどのマスで切るべきかBER（帯域エネルギー比）の公式は割り算 1 つだけですが、コードにはエラーを出さない落とし穴が 2 つあります。「時間フレームの数」を「周波数のマスの数」と取り違えること、そして近似したマスの間隔を使って 2000 Hz を間違った側で切ってしまうことです。この講では 30 秒の音楽 3 つについて、フレームごとに BER を計算します。正しい境界は第 186 マスです。1292 個の時間フレームを周波数のマスの数と取り違えると、境界は第 234 マス、実際には約 2519 Hz にずれます。ドビュッシーの曲の BER の中央値は、そのせいで 4.37 dB 大きくなります。
- 23スペクトル重心と帯域幅：周波数の中心と広がりをどう測るかスペクトル重心と帯域幅は、ライブラリの関数 2 行で済まされてしまうことがよくあります。この講では逆の順序で進めます。まずライブラリの既定の設定で 3 つの音楽の時間変化を求め、次に 2 行の関数を公式に分解して、自分でフレームごとに計算し直します。重心の最大差はわずか 1.917e-05 Hz、帯域幅の最大差はわずか 1.156e-05 Hz です。突き合わせる過程で、取り違えやすい定義が 1 つ明らかになります。平均絶対距離は p = 1 ですが、ライブラリ関数が既定で計算するのは二乗平均平方根距離 p = 2 で、この 3 つの音楽では両者に 1.30〜1.54 倍の差が出ます。