mp3 encoder uzura3 の作製と反省

# 前置き

昔 Fortran 90/95 で mp3 encoder を書いていたのですが、うまくいかない所があって 15 年くらい放置していました。最近 fpm (Fortran Package Manager) の練習としてファイルを整理していた時に、ようやくおかしい点に気づきました。（Fortran のいい所は、15 年ぶりにいじってもバージョン違いの問題もなく、最新の環境でそのまま作業できるところです。）

本記事では、mp3 encoder の仕組みの概説とともに、なぜ中々誤りに気付けなかった反省します。何らかの教訓となれば幸いです。

# プログラム

プログラムは fpm 対応にして github に上げてあります。

https://github.com/cure-honey/uzura3_fpm

パラメータ調整はまだ途中ですが、正しい方向に進んでいるのではないかと思います。

## 動機

私が mp3 encoder を書き始めたのは気分転換と Fortran 90 の腕試しのためで、対角化や線型方程式を解くことに飽きて、何か線形代数と違ったものをやりたくなったからです。

オーディオの知識は全くなく、デジベル(dB)の定義も常用対数の何かだったこと以外全く忘れていました。CD のデータ形式についても初めて知りました。心理音響解析の記事を読んで、人間の主観に聞こえている音は、実際の音からわずかに遅れていて大きな音が鳴ると先回りして音を聞こえなくするなど主観的には因果律が破られていることも初めて知りました。中国共産党の海外テレビニュースのライブ検閲のようです。

休みの日にちょこちょこプログラムをして、3 を基数とする FFT を書いたり、名前しか知らなかった Huffman code や CRC のルーチンを書いたり、最初の頃はすべてが新鮮で本当に楽しかったです。　


# mp3

mp3 とは MPEG I audio layer III の略称です。名前の示すようには動画の音声トラックとして定義されていますが、64kbps ISDN ２回線での 2ch stereo 音楽ストリーミングを一つのターゲットとして、128kbps で CD 並みの音質を出すことを目的としたようです。これは CD に比して 10 倍以上のデータ圧縮率となります。

mp3 は、ISDN での音楽ストリーミングには利用されませんでしたが、コンピュータ上で音楽を聴くという目的で、1990 年代後半から 2000 年代前半にかけて爆発的な人気を博しました。丁度インターネット勃興期でもあり、CD から rip した音楽データを mp3 にエンコードしてネットに流すことが広く行われて、ネットでのファイル共有や著作権、デジタルデータの所有権など様々な社会問題を提起することになりました。

MPEG I の音声フォーマットとしては　layer III の他に Layer I, II があります。これらは圧縮率は低いのですが、polyphase filter band と呼ばれるバンドパス・フィルターを用いた素直で簡潔な実装しやすいアルゴリズムになっています。一方 Layer III は、polyphase filter band をかけたあとさらに離散コサイン変換 (MDCT: Modified Cosine Transformation) をかけるという、木に竹を接いだような構成をしており、その他もかなり複雑になっています。

# mp3 encoder と Fortran

さて mp3 というとマルチメディア関連のイメージで Fortran が扱うような主題ではないイメージがあります。Fortran は実数（複素数）のような連続量の近似的計算に主に用いられ、情報系の離散数学的な離散量の計算には関心が薄い気がします。

しかし形式面から見ると mp3 encoder は Fortran の得意とする主題に近いことが分かります。

### 1. 時系列解析

CD の音声データは、44.1kHz のサンプリングレートで等間隔に並んだ 16bit 2ch の振幅データになっており、mp3 encoder はこれを周波数領域に変換して処理を行いますが、これはまさに時系列解析そのもので、Fortran の得意とするものになっています。周波数領域への変換は、行列積と FFT で実現されます。

### 2. 心理音響解析

心理音響解析とは、人間の耳の分解能や感受性を考慮して、ある音響下での音の検知閾値を求めるものです。

静音下で聞こえる音でも、同時に他の音が存在すると聞こえないことがあります。これを masking と言います。会話中に近くで掃除機をかけられると、相手が掃除機に負けない大声を出さないかぎり声は聞こえなくなり、話を続けようと、話を止めようと掃除機の音しかきこえません。これを応用すると、聞く側の認知には変化を与えずに、音を省いたり付け足せることが可能と分かります。どの程度の大きさの音まで認知されないかを求めるのが心理音響解析となります。

人間の耳は鼓膜で音を周波数分解していますが、これはフーリエ変換で表わせます。また masking は鼓膜の振動子が独立ではなく、隣接する周波数の振動子と連結した連結振動子になっていると考えると素朴には理解できます。ある周波数の入力が他の周波数を mask するのは、周波数空間での畳み込みで表わせせることが分かっています。

結局、心理音響解析も主として FFT と行列畳み込みに還元されます。 

### 3. １変数方程式の求解

周波数変換された音声データは、振幅を量子化した上で Huffman code で圧縮します。量子化間隔が荒ければ圧縮に必要な bit 数は小さくなりますが誤差は増えます。量子化間隔が細かければ必要な bit 数が増えますが誤差は減ります。出力の bit rate から許される bit 数の下で、最も量子間隔を細かくすると、量子化誤差を最小にできます。

これは量子化間隔 q に対して入力データ x を Huffman code 化した bit 数 $G(x; q)$ が許容 bit 数 M に等しくなるよう $ G(x; q) = M $ をみたす q を定める 1 変数方程式の問題と見なせます。（実際は q は整数値しか取らないので $G(x; q)\le M \lt G(x;q+1)$ の不等式になります。）

この時、$G(x; q)$ は微分を求められないので Newton 法は使えず、数値計算の授業で一番最初によく出てくるはさみうち法での求解になります。

### 4. 極値問題

3．の量子化による誤差は、2．の心理音響解析で求められた masking の範囲内であれば、ずれていても主観的には問題無いことになりますが、masking の範囲を超えている場合、scale factor というパラメータを変化させることにより量子化誤差を減らす仕組みが用意されています。

(long block の場合）scale factor は 21 個あるので、21次元のパラメータ空間のベクトル x を動かして周波数空間での音声データ $F(f)$ とそれを量子化間隔 q で量子化した $F'(f,x;q)$ が心理音響解析で求めた許容誤差 $\delta(f)$ 以下になるように $ |F(f) - F'(f,  q; x)| \le \delta(f) $ 誤差を小さくしてゆく極値探索問題と見なせます。

mp3 encoder の実際の処理はこの 4 つに還元されます。プログラムとしては、この前後に wav ファイルと mp3 ファイルの I/O がくっつく形になります。


こうしてみると、mp3 encoder は Fortran による数値計算でよく見かけるものの組み合わせのアナロジーで理解できることが分かります。そうして、このアナロジーが mp3 encoder の制作を容易にするとともに、先入観により誤りを見つけるのを困難にしていました。

# 先入観による誤り

さて長く解決できなかった問題は、4 番目の残差ベクトルの最小化の極値探索の部分です。残差の部分を関数空間の重み付きの距離問題 $\int |f(x)-g(x)|^p w(x) dx$ のイメージで捕えて、なんとなく $f(x),g(x)$ も滑らかな連続関数と思い込んでしまいました。このため二重に誤りを招きました。

ここにきてデジタル信号処理の基礎知識の欠落が問題になりました。

### 離散値の問題

ひとつは残差を減らすように 21次元の scale factor を変えるところで、連続関数でグラジエント ($\nabla $) を取るイメージで、残差が許容値を超えている場合対応する scale factor 全てを一斉に変化させましたが、量子化されているために最小変化量が 1 で変化が大きくなりすぎて、極値を通り過ぎて探索に失敗していました。連続変数なら $(1,1,1)/\sqrt3$ と規格化される所が、最小変化量が１なので $(1,1,1)$ で変化して行き過ぎるようなものです。離散値を取るデジタルの場合、いちどにパラメータ１個変化させるのが基本のようです。 

### dB の問題

信号処理ではよく縦軸が対数を取った dB になっていますが、dB のセンスが無かったために、残差の距離の計算は対数を取る前の真値での差の（重み付きの）相加平均で良いと思い込んでいました。確かにダイナミックレンジが大きいので、対数を取るのはもっともだとは思ったのですが、対数値で相加平均を取ると真値では相乗平均になってしまうので、関数間の距離が積で定義されることになって ($\prod |f(x)/g(x)|dx$) 何か納得がいかない気がした覚えがあります。

普通ならばとりあえず分野の常識に従ってよく考えてみるところですが、「dB よくワカラン」で済ませたのが失敗でした。よくよく考えると、周波数空間での波形は正弦波のフーリエ変換がデルタ関数になることから分かるように、ギザギザした不連続でいたるところ微分不能な形をイメージする方が適切で、差よりも比で距離を定義する方が合理的な気もします。実際、対数で残差を定めた方が良い結果が得られます。


# 心理音響解析 psychoacoustic analysis

上記の修正によって、残差を小さくする極値探索の部分も期待した動きをするようになったので、残差の許容値(masking 閾値と同じもの)を与える心理音響解析ルーチンも新たに作りました。心理音響解析というと、難しそうで長年ビビッていましたが、背後の理屈はともかく形式的には FFT を何本かづつ束ねたあと耳に固有の spreading function という単純な関数と畳み込みを取るだけのもので、Fortran では数行で書けてしまいます。

モデルとしては、昔にネットの掲示板で紹介されていた F. Baumgarte, C. Ferekidis and H. Fuchs の "A Nonlinear Psychoacoustic Model Applied to the ISO MPEG Layer 3 Coder" を用いました。

一般に心理音響解析は、経験的なパラメータが多い複雑なモデルが多いのですが、この論文では静音時の耳の感受性での規格化と、畳み込みで非線形和を取ることにより最小限のパラメータで実験事実を説明するもので、本物くさい香りを漂わせています。

この論文中ですでに指摘されているのですが、mp3 固有の事情で encode に使えるビットを、フレームと呼ばれる１単位当たりで、左右のチャンネルと granule と呼ばれるさらに時間的に細かく２分割した部分とで４つに分けて利用するのですが、このビット分配が音質に非常に重要になっています。

皮肉なことに scale factor 部分でつまづきましたが、一旦出来てしまうと scale factor で極値探索して残差を減らすより、ビット分配を工夫する方がよほど大事な感じを受けました。いずれにせよ、これらのためには心理音響解析が適切な許容誤差の見積もりを与えてくれることが重要になります。

# 残る謎　polyphase filterbank + MDCT vs FFT 

mp3 では音声の時系列データを周波数領域に変換するのに、Layer I/II と共通の polyphase filterbank 　というバンドパスフィルターを掛けて 32 個のバンドに分解した後で、さらに MDCT を掛けてより細かな周波数領域に分解しています。この木に竹を接いだようなスペクトルを、同じ時系列データに直接 FFT を掛けて同じ周波数間隔にしたスペクトルと比べると、結構な違いがあります。

正弦波を入力として与えると、FFT の方は期待される周波数位置にピークが立ってかつ大きさは時間的に変化せず定常で一定値を取ります。これに対して polyphase filterbank + MDCT ではピーク位置は広がりながらも期待する位置に出ますが、強度が時間的に 2 桁位大きく変動しています。

またフーリエ変換はユニタリー変換なので、入力と出力の二乗強度和が一致するはずで、計算値もそうなっていますが、polyphase filterbank + MDCT の方は二乗強度和が一致せず大きさも比率も変動しています。

これが私の誤りや勘違いでないとすると、FFT によっていくら正確に心理音響解析をしても、これを mp3 に適用するところでブレ幅が大きくて、何をやっているのかよく分からないことになります。少なくとも真値での大小比較は意味がなく、相対比率での比較ぐらいしか意味がないのかもしれません。

この辺は、未だ謎のままです。 


# 反省

仕様上の致命的なバグがまだ残っているかもしれませんが、個々のルーチン的は昔の時点で正しく機能しているのではなかと思います。

最後に mp3 encoder 作製に長期間かかった点を反省したいと思います。 

1．遊びのプログラムなので締め切りもなく気が向いた時にいじる形で楽しくやれてはいたが、よく分からなくても専門外だからといって追及しない、原因をとことん探らない、きちんと記録ノートなどを取らないなど、やってはいけないとされることをやってしまって見事にドツボにはまってしまった。

2. 基礎知識が足りず、分野の常識に欠け、よく分からない所をうやむやにしたので、文献値や先行研究の再現を怠った。結果的に最初に確認すべき規格化や値の原点の確認などを怠った。

3．連続関数や Fortran の数値計算のアナロジーに頼り過ぎた。アナロジーは助けになるとともに、また足かせにもなる。

以上、いくつか反省を書いてみました。今回プログラムを修正するにあたっては、よく文献を読まず、ノートも殴り書きの記録しか取らず、規格化因子等をうやむやにしてしまったので、口先はともかく行動では反省が足りないようです。