はじめに
AIきりたんに衝撃を受け、その仕組みについて調べたことのアウトプットです。
波形画像を用意しようと思って1年以上温めてしまったので、おそらくもう自分は画像を用意する気がないのだなと思い、文字ばかりですが思い切って公開します。またその都合で参考文献などを忘れてしまい、引用できていませんが、ご了承ください。
NEUTRINOを支える技術のうち、最も骨子となっているWOLRDの用いている統計的音声分析の一つである、ケプストラム分析に主眼を置いて解説していきます。
いくつか数学（極限まで）やシステム（伝達関数あたり）の知識が本文中で必要となりますが、そういったことを求めている場合は役に立つかと思います。
また、微妙に長いので、である調で本文は記述します。
スペクトル分析手法
スペクトル分析の手法には次のような種類がある。
パラメトリック手法
- 線形回帰分析（LPC）
ノンパラメトリック手法
- ケプストラム分析
- メルケプストラム分析
- 一般化ケプストラム分析
- メル一般化ケプストラム分析
パラメトリック手法
本記事の主眼はノンパラメトリック手法であるケプストラム分析なので、軽くだけ説明する。
全体を俯瞰するものの詳しく説明しないので、わからないところがあれば適宜調べるか、ケプストラム分析までジャンプするとよい。
線形回帰分析
線形回帰分析（LPC; Linear Predictive Coefficient、線形予測係数またはLinear Predictive Coding、線形予測符号化）は、音声波形を確率過程としてなんらかの数式を用いてモデル化する。
まず、音声波形は音波の振幅の時系列データで表せる。
また、こと周期的な波形に関しては、ある時刻$t$での振幅$A_t$は、その前の時刻$t-1$での振幅$A_{t-1}$を用いて予測できそうである。
1秒前は雨だったのだから、今も雨と考えるようなものと考えると良い。
つまり、ある時刻での振幅を用いることで、次がどのようであるかを確率的に論じる。
これを線形回帰分析では、過去の情報の線形結合で表す。
つまり、式に表すと次のようになる。
$$
A_t = c_0 + c_1A_{t-1}
$$
1時刻前の情報だけでは十分に予測ができないが、もっといくつか過去の情報を何らかの係数を掛けて足し合わせていけばより高い精度で予測ができそうである。
すなわち、
$$
A_t = c_0 + \sum_{n=0}^{p}c_pA_{t-p}
$$
と、過去$p$時刻分の情報の線形結合で、現在の振幅を表現する。
この係数$c_p$は、実際に得られたデータとの差
$$
\epsilon_t = c_0 + \sum_{n=0}^{p}c_pA_{t-p} - A_t
$$
を小さくするように回帰分析すれば得られる。
ところで、この誤差$\epsilon_t$は、モデルの予測出力との差異であるので、係数$c_p$が決まっている状態で誤差$\epsilon_t$が分かれば、モデル出力と足し合わせることで波形を復元可能である。
よって、モデルの係数$c_p$を決定した後、予測誤差$\epsilon_t$だけで音声波形を表現できると言い換えることができる。
これで現時刻の振幅$A_t$を入力することで予測誤差$\epsilon_t$が得られる線形モデル（音声波形を入力し予測誤差を出力するフィルタ）が得られた。
このフィルタの伝達関数は、ただの線形結合の式で与えられ、解が全て零点となる全零フィルタである。
また、予測誤差から音声波形を得るフィルタはこの逆フィルタで与えられ、伝達関数は逆数となり、こちらは解がすべて極となる全極フィルタとなる。
このように、自らの次の状態を自らの過去の状態で予測する分析手法は、自己回帰（AR; Auto Regression）分析と呼ばれる。
予測誤差を全極フィルタに通すことで音声波形が得られる、ということはすなわち声道特性フィルタに駆動音源を入力したとみることができる。
このため、波形を入力して駆動音源$\epsilon$を得る解析フィルタと、駆動音源に声道特性の色付けを行う復元フィルタが得られていると考えられる。
声道を音響管の多数の連結であると考え、これを自己回帰過程としてモデル化したのが統計的音声解析におけるLPCである。
すなわち、解析用のフィルタ（全零フィルタ）の係数$c_p$を回帰により決定することで、入力音声波形を駆動音源へと変換し、音声波形の復元には単純に逆フィルタ（全極フィルタ）を通して行う。
この中間で得られる駆動音源は、この音響管多数接続のモデルに対しての声道特性を示すことになり、これを用いて音声の解析を行うのである。
ノンパラメトリック手法
ノンパラメトリック手法では、あらかじめ決められたモデルの定数を決定するパラメトリック手法とは異なり、データに特定の処理を施して得られるデータを用いる手法である。
ケプストラム分析
ケプストラムとは
ノンパラメトリック手法によるスペクトル分析においては、ケプストラム分析が主流である。
ケプストラム（cepstrum）とは、スペクトル（スペクトラム、spectrum）の頭4文字を逆に並べ変えたアナグラムである。
spectrum = "spectrum"
cepstrum = spectrum.substring(0, 4).split("").reverse().join("") + spectrum.substring(4)
その名の通り、時系列波形データのフーリエ変換によって得られるスペクトラムを、さらに逆フーリエ変換するためにこのように呼ばれる（ただし、現段階でこの説明はある程度の誤解を含んでいる）。
ケプストラム分析の概要
まず、音声波形の成り立ちを考える。
声は、喉の声門の素早い開閉により発生するブザー音が、声帯、口腔などの状態によって特定の周波数が共鳴することで発声される。
管楽器経験者なら、マウスピースだけで音を鳴らした時のブザー音が、楽器内部を通ることでその楽器特有の音色へ変化する体験をしたことがあるだろうが、それがいい例である。なおこの場合声門に対応するのは唇やリードの振動である。
耳に聞こえる音の高さ（ピッチや音高という、音程ではない）は、声門の開閉間隔によって決まる。声門の1秒間の開閉回数は周波数の次元であり、これを基本周波数と呼ぶ。これを用いてこのブザー音音源は、近似的に基本周波数の逆数の間隔で並ぶパルス列で表される。
また、この声門の振動によるブザー音を音源として、声帯などを通るうちに色付けされるため、色付けするこの要素は、声道特性となど呼ばれる。
すなわち、パルス列音源が声道特性によってフィルタリングされたのが声というわけである。
これを入出力システムととらえれば、入力パルス列音源を$X(s)$、声道特性フィルタの伝達関数を$H(s)$として、声自体の周波数領域表現、すなわち声の波形スペクトル$Y(s)$は
$$
Y(s)=H(s)X(s)
$$
と表される。周期 $T$ のパルス列の周波数軸表現では、周期 $f$ の、やはりパルス列となる。
声の波形をフーリエ変換することでスペクトルを求めれば、声道特性とパルス列の積が得られることがわかる。声の特徴を持つのは声道特性フィルタ $H(s)$ であるのだから、これを分離すれば、その声の特徴を捉えることができる。そこで、積の対数は、それぞれの対数同士の和であるという対数関数の性質
$$
\log Y(s) = \log H(s)X(s) = \log H(s) + \log X(s)
$$
を利用する。
すなわち、得られたスペクトルの対数を取れば、声道特性フィルタの対数と基本周波数パルス列の対数の、それぞれの和となる。これを対数スペクトルなどと呼ぶ。
ところで、声の特徴は、一周期の間にどのように波形が変化するかによってみられるので、聞こえる声の音高の周波数と、声の特徴の周波数は明確に異なり、声の特徴の方が高い周波数を持っているということが考えられる。
対数スペクトルでは、高い周波数を持つ声道特性フィルタが、これに比べて低い基本周波数のパルス列の影響で全体的に揺らいでいるような波形が見られる。
すなわち、対数スペクトルで見たとき、その対数スペクトル波形の低周波成分が声道特性、高周波成分がパルス列による成分であると考えられる。
周波数成分を見るにはその波形のスペクトルを見ればよいのは、声の波形を分析する過程と同じ考えである。
よって、対数スペクトル波形をもう一度フーリエ変換し、さらにスペクトル波形を得る。
ところでフーリエ変換は時間軸と周波数軸の変換であり、逆に周波数軸から時間軸へ戻すのが逆フーリエ変換であったが、これらの演算は実質的には同一であるので、対数スペクトル波形を逆フーリエ変換することでも、やはり周波数成分を分析することができる。このため、音声分析においてはここで逆フーリエ変換で演算する。
ここで得られる波形は一度対数を取っているから厳密ではないものの、時間軸へと戻っている。
これを表すため、横軸はケフレンシー（quefrency、やはり周波数frequencyのアナグラムである）と呼ばれ、個々で得られる波形は__ケプストラム__と呼ばれる。
ケプストラムを見れば、低いケフレンシーの領域でいくつかの波形が見られ、ギャップが空いたのち基本周波数成分のケフレンシー領域で再びピークが見られるような波形が得られる。
これらは声道特性の周波数と、基本周波数の差が十分あるため容易に識別可能で、基本周波数より少し小さい値で分けてやれば声道特性と音源特性が分離できる。
このときにフィルタを通すことでこれらを分離するのだが、これも通常の時間軸のフィルタとは異なるため、リフター（lifter、やはりフィルターfilterのアナグラムである）と呼ばれ、その処理はリフタリングと呼ばれる。
リフタリングにより分離した声道特性と音源特性は、それぞれ変換の逆過程を経ることで時間軸へ戻すことができる。
実際には音声データはマイクによって離散情報として得られるため、フーリエ変換は離散フーリエ変換が用いられ、アルゴリズムとしてはFFTが用いられる。
また、音声を声道特性と音源特性に分離する際、音高を示す音源特性は音高を変えるだけで変化し、声質や音色を示す声道特性は声に乗せる勘定や母音の変化によって変化してしまうので、離散フーリエ変換時のウィンドウサイズは音高による周波数より小さい必要がある。
音声においては200Hz程度のような低い音は出せないだろうということで、その周期5msが採用されることが多いようである。音声を5msに区切って分析していくため、この5msの音声データをフレームと呼ぶことがある。
一般化ケプストラム分析
一般化対数関数
一般化ケプストラム分析では、ケプストラム分析において対数スペクトル波形を算出する際の対数の演算を、一般化対数関数に置き換える。
この一般化対数関数ついて、少し解説する。解説はするが、なぜ一般化対数関数が登場したかのお気持ちは想像であるので承知いただきたい。
まず、自然対数の底やネイピア数として知られる $e$ について、次のような極限で求まることが知られている。
e = \lim_{h\to\infty}\left(1 + \frac{1}{h}\right)^h = \lim_{h\to 0}\left(1 + h\right)^\frac{1}{h}
また、指数関数$\exp(x)=e^x$については、次の極限で求まる。
$$\exp(x)=e^x=\lim_{h\to 0}\left(1+hx\right)^\frac{1}{h}$$
さて、ここで極限を取らずに、次の関数$s_\gamma^{-1}(x)$
$$s_\gamma^{-1}(x)=(1+\gamma x)^{1/\gamma}$$
を考える（$h\to\gamma$としている）。
この関数では、$\gamma=0$（厳密には$\gamma\to0$）のとき指数関数を表す。このため、指数関数の拡張と考えられ、パラメータ$\gamma$を調整することで、似た性質を持ちながら異なる関数として機能する。
ここで、わざと$s_\gamma^{-1}(x)$と表現しているのは、一般化対数関数を$s_\gamma(x)$と表現したいからであり、先ほどの一般化指数関数の逆関数でもって一般化対数関数を考えることができる。すなわち一般化対数関数$s_\gamma(x)$は
$$s_\gamma(x)=\frac{x^\gamma - 1}{\gamma}$$
で定義される関数である。この関数でも、やはり$\gamma\to 0$の極限により通常の対数関数$\log x$と一致する。
しかし、厳密にはこれでは$\gamma=0$が議論できないので
$$
s_\gamma(x)=\begin{cases}
\log x & (x=0) \\ \
\cfrac{x^\gamma - 1}{\gamma} & (\mathrm{otherwise})
\end{cases}
$$
となり、さらに $\gamma$ はもともと $0 $に収束させることで対数として表現できることから、$|\gamma|<1$ の条件が付される。
すなわち、一般化ケプストラム分析には、$\gamma$ というパラメータが存在し、$\gamma=0$ のときに通常のケプストラム分析となるような、ケプストラム分析の拡張と考えられる。この $\gamma$ を調整することで通常のケプストラム分析より柔軟に分析することが可能になる。
メルケプストラム分析（メル周波数ケプストラム分析）
メルケプストラム分析は、先の一般化ケプストラム分析とは枝が異なり、通常のケプストラム分析にメル尺度やメル周波数などという、人間の周波数による聴覚特性の考え方を導入したものである。
メル尺度
人間の聴覚では、低い周波数においては敏感に、高い周波数では鈍感に反応する。すなわち、周波数の対数の尺度で音高を識別する。
音楽のオクターブと、周波数のオクターブ（2倍周波数）というのは偶然の一致ではなく、実際に周波数が倍なのである。
そのあたりはフェヒナーの法則でも示されており、簡単に確認できるだろう。
しかしそれほど人間の聴覚は単純ではないようで、フェヒナーの法則は一定の領域ではうまく近似できるものの、可聴域全体にわたってはうまく近似できないことも知られている。そこで実際に、ある周波数に対してどのような音高を感じるかを表すための尺度がメル尺度である。
メル尺度では、実際に周波数の対数スケールで表され、具体的には
$$
m=m_0\log\left(\frac{f}{f_0}+1\right)
$$
と表される。これは周波数がゼロのときに$m=0$として、ある周波数$f_0$を基準として、そのときのメル尺度を$m_0$としただけの数式である（$(f, m) = (0, 0), (f_0, m_0)$の2点を通る自然対数関数のもっとも単純なものである）。周波数 $f_0$ におけるメル尺度 $m_0$ 自体は、 $1000[\mathrm{Hz}]$ のときに $m=1000$となるように定められる。具体的な基準周波数 $f_0$ としては $700[\mathrm{Hz}]$ や $1000[\mathrm{Hz}]$ などが選ばれ、よく見るのは $f_0 = 700[\mathrm{Hz}], m_0=1127.01048[\mathrm{mel}]$ である。
基準周波数 $f_0$ の取り方にもよるが、このメル尺度においては、1000Hzの聴覚的オクターブ上が4000Hzあたりとなる場合がある。このことはフェヒナーの法則と反するが、実際に人間の聴覚特性は、メル尺度の方がよく従うことが知られている。すなわち、周波数を倍にしてもオクターブ上に聞こえないということが起こりうる。このことは、耳の構造による共振によるものであると考えられている。
この尺度を周波数の代わりに使用することで、より人間の聴覚特性に近い分析ができる。
また、高周波領域は人間にとって差を感じ難いという性質をよく反映することから、メル周波数領域では高周波成分の密度が減り、結果として元の音声の情報を保ったまま高音部の情報を減らすことができるのが利点である。
メルケプストラム分析の概要
通常のケプストラム分析における、周波数軸上のデータを、メルケプストラム分析ではメル尺度のスケールへ変換してから分析を行う。
このとき、周波数が高くなるほど密度が低くなるような三角波の重ね合わせを乗算することによってスケールの変換を行い、その三角波をメルフィルタバンクと呼ぶことが多い。
この結果得られる値は、メルフィルタバンクのチャネル数と一致し、それぞれメルケプストラム係数、あるいはメル周波数ケプストラム係数（MFCC; Mel-Frequency Cepstrum Coefficients）と呼ばれる。
反対に言えば、MFCCの次数は、たかだかメルフィルタバンクのチャネル数程度であり、それだけで声の特徴を捉えることができるという点において非常に優れていると言える。
メル一般化ケプストラム分析
ケプストラムに一般化対数の考えを導入した、一般化ケプストラム分析と、メル尺度の考えを導入したメルケプストラムは別の流れである。メル一般化ケプストラム分析では、これらをつなぎあわせたものである。
すなわち、スペクトル波形を得るときにメルフィルタバンクを適用し、そこから対数スペクトル波形を得る過程で一般化対数関数を用いるのである。
この結果得られるベクトルは、メル一般化ケプストラム係数（MGC; Mel-Generalized Cepstrum）と呼ばれる。
この方式では60次元程度で元の音声を高品質に復元できるといわれている。
この手法では、$\alpha$ と $\gamma$ の二つのパラメータが存在し、これらをうまく調節することで高品質な音響特徴量の抽出が可能である。
また、特定の $\alpha$ と $\gamma$ の組み合わせにおいて、他の分析法と等価となるので、これらの手法を統合した手法であると言える。
AIきりたんで有名なNEUTRINOが用いているWORLD法に登場するスペクトル包絡とはこのことであり、実際にNEUTRINOが出力するファイル「*.mgc」はこのメル一般化ケプストラム係数である（このように声を基本周波数、スペクトル包絡、非周期性指標に分解するような方式はVocoder（ボコーダー）と呼ばれる）。
おわりに
本記事を読んでより統計的音声分析に興味を持ってもらえたら幸いです。
実際にこれらの分析を行うためのシステムがWORLDになっています。オープンソースでありソースコードはGitHub（https://github.com/mmorise/World）に公開されています。各分析に対応するメソッドはジョジョにインスパイアされた名前で公開されており、それぞれ論文も見ることができます。
WORLDのようなアルゴリズムを用いることによって、音声波形を直接扱うのではなくこれらの分析結果を特徴量として用いて機械学習させるのが近年の潮流であり、NEUTRINOがそれにあたります。
音声認識や声質変換などの技術も、この統計的分析手法が基盤になっているものが多くあります。
また、何か間違いや疑問などありましたら、私もよりこの分野について詳しくなりたいのでコメントなどぜひ残していただければと思います。
余談
ところで、*.mgcファイルはCの順番でdouble型の表現で記述されており、5ms間のフレーム（WORLDのデフォルト値）におけるmgcの値（メルフィルタバンク60次元分）が並んでいます。
PythonならばNumPyにCの順番でバイナリファイルを読むnumpy.fromfileと、Cにおけるdoubleと対応する型numpy.float64を用いてnp.fromfile("path/to/mgc", dtype=np.float64)とすれば読むことができ、さらにリストの添え字をフレームと一致させるためreshape(-1, 60)を呼べば容易に扱うことができます。60次元もあるので可視化はしにくいですが、とりわけmgcの0次はケフレンシー領域の直流成分になりますから、Pythonで適当なコードを書けば音量を見ることができ、変化させることもできます。