5 目次 第 1 章基本理論編 線形予測 1.1.1 予測とデータ圧縮 1.1.2 線形予測の理論 線形予測 Levinson-Durbin 再帰 (Levinson-Durbin recursion) の導出 Levinson-Durbin 再帰アルゴリズム . PARCOR 格子型フィルター PARCOR 係数による信号合成 PARCOR 係数の値域 1.1.3 線形予測の実装 実験 1.2 工ントロピー符号 1.2.1 情報とその符号化 一意に復号可能な符号 . McMillan ( マクミラン ) の不等式 1.2.2 情報量とエントロピー 情報量の定義 . 情報量の導出 . 工ントロピー 平均符号長とエントロピー 様々なファイルのエントロピー計測 . 1.2.3 工ントロピー符号の例 。 , 7 符号とその確率分布 . GoIomb ( ゴロム ) 符号 GoIomb ( ゴロム ) 符号のパラメータ設定 . 0 0 っ 4 っ -4 ー -8 一 0 1 00 ー冖ー -8 1 -1 ワ 1 -4 1 ・ 1 つ「ー -1 1 -1 1 -1 1 1 ワ 1 っ 4 っっっ 1 っっ 1 っ 4 っつひりっ 0 っっ -4 一 4 -4 ・ 4 1.1

14 より、行列形式で 1 ) ② ( 0 ) 月 ( 1 ) R(M) R(M ー 1 ) と表せられる。以下、 ( 1 ) 町 2 ) R(M) M = 月 ( 0 ) ( 1 ) R(M ー 1 ) 1 ) 0 ) R(M ー 2 ) 1 ) ( 0 ) R(M ー 2 ) R(M ー 1 ) R(M ー 2 ) 0 ) R(M ー 1 ) R(M ー 2 ) 月 ( 0 ) 第 1 章 1 aM(1) ( 2 ) aM(M) CLM 基本理論編 1 佖 M ( 1 ) M ② ( 1.7 ) aM(M) として、 M M = 0 を解くことを考える。 Levinson-Durbin 再帰 (Levinson-Durbin recursion) の導出 前節で求めた連立方程式 M M = びをもう少し整理していく。数値解法的には、 M は 正方行列にしておくのが望ましい。そこで、 M の一番上の行に ' 国 ( 0 津 ( 1 ).. 田 ( M ) ] を追 加すると、 R(O) 月 ( 1 ) Mäk 町 1 ) 月 ( 0 ) aM(M) 1 R(M) aM(1) aM ② R(M ー 1 ) ( 0 ) 0 0 0 R(M) R(M ー 1 ) と変形できる。よって、 を解くことに帰着できる。 1 ) 0 ) 1 ) 月 ( 0 ) R(M) R(M ー 1 ) 係数。 M ( 1 ) , ・ . ,aM(M) を求める問題は、次の連立方程式 1.8 aM(M) 1 R(M) 佖 M ( 1 ) 佖 M ( 2 ) R(M ー 1 ) 町 0 ) ー 0 aM(m)R(m) 0 0 0 ( 1 ・ 8 ) 連立方程式 1.8 を高速に解くアルゴリズムが、 Levinson-Durbin 再帰法 (Levinson- Durbin recursion) である。以下、 NM を式 1.9 で定義する。 0 ) 町 1 ) R(M) = を . = 0 。 M ( m 津 ( m ) とし、また自己分散行列 1 ) 0 ) R(M ー 1 ) R(M) R(M ー 1 ) R(O) ( 1.9 )

1.1 線形予測 1 ュ .3 線形予測の実装 Levinson-Durbin 法の C 言語による実装をリスト 1.1 に示す。 #include く float . > 17 { 23 Levinson-Durbin 再帰法の実装例 1 2 3 4 5 6 7 9 10 11 12 13 14 15 16 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 リスト 1.1 く stdio . h> く math . > #include く stdlib . h> く string. h> く stdint . 五 > #include #include #include #include 8 / * ( 標本 ) 自己相関の計算 * / static int32_t calc—auto—correlation(double *auto—corr, num—sample , uint32—t max—order) ; / * 乙 e ⅵれ s 。れ一加再帰計算 * / static int32_t const double *data, uint32_t levinson—durbin—recursion(double *lpc—coef , uint32-t max-order) ; const double *auto_corr, int main (void) int32—t num—sample int32—t max—delay int32—t smpl , i ; = 300 ; / * サンプル数 */ 10 ; / * c 係数の数 * / (double *)malloc (sizeof (double) * num—sample) ; double *data (double *)malloc (sizeof (double) * num—sample) ; double *predict (double *)malloc (sizeof (double) * (max—delay + 1 ) ) ; double *auto_cor sin(smpl * 0.01 ) + 0.5 * cos ( 4. Of * sin(smpl * 0.05 ) ) ; smpl く num—sample ; smpl + + ) { (double *)malloc(sizeof (double) * (max—delay + 1 ) ) ; data [smpl] for (smpl = 0 ; / * 波形の生成 * / double e て ro て ; double *coff levinson—durbin—recursion(coff , auto—cor, max—delay) ; calc—auto—correlation(auto—cor, data, num—sample , max—delay + 1 ) ; / * 自己相関・ e ⅵれ s 。れ初再帰計算 * / / * 以降は予測 * / } else { data [smpl] ; predict [smpl] / * 最初の m 。 de 乙ステップ分は元信号を単純コピー */ if (smpl く max—delay) { f0 て (smpl 0 ; smpl く num-sample ; smpl + + ) { / * 予測テスト * /

1.1 線形予測 Levinson-Durbin 再帰法は、数学的帰納法によく似ており、次のステップにより係数を 求めていく。 * 6 1. 初期化ステップ : 1 次の係数の ( 1 ) を求める。 2. 再帰ステップ : ん次の係数 ( 1 ) , ... , ( ん ) から、ん十 1 次の係数 + 1 ( 1 ) , ・ 1 ) を求める。 こでは、 1. および 2. の場合の解をそれぞれ見ていく。 ■初期化ステップ M = 1 のとき、 , ル 1 は次の様に定義される。 0 月 1 1 15 ・ , + 1 ( ん十 佖 1 ( 1 ) ' 実際にを計算してみると、 ル 1 1 0 0 十月 1 ( 1 ) 1 十召 0 佖 1 ( 1 ) 従って、 el = 十月 1 の ( 1 ) 、及び十召 0 の ( 1 ) = 0 からの ( 1 ) = と求められ る * 7 ・再帰ステップ仮定として、 ( 0 ) 月 ( 1 ) 召 ( ん ) ( 1 ) ( 0 ) ん一 1 ) ん ) 月 ( ん一 1 ) ( 0 ) が成立していたとする。 M = ん十 1 の時、行列ルん + 1 は ( ん ) ん + 1 ) ん一 1 ) ( ん ) 月 ( 0 ) ( 1 ) ( 1 ) 月 ( 0 ) 咒 ( ん + 1 ) 月 ( ん ) 月 ( 1 ) 町 1 ) 0 ) 1 ( 1 ) ( 2 ) 0 0 ( ん ) 0 0 ) 1 ) 1 ) ( 0 ) 町ん ) 月 ( ん一 1 ) ん + 1 ) 月 ( ん ) + 1 町ん + 1 ) ん ) 、 6 参考資料 [ 19 ] で筆者は、「 Levinson-Durbin 帰納法と言ったほうがいいんじゃないか」と書いてあった。 Ro= を > 0 より、確率論的に至る所ゼロ除算の心配はない・・・が、デジタル信号においては が全て 0 のときが往々にして起こる。この場合は、係数を全て 0 にする等の特殊処理が必要である。 * 7

1.1 線形予測 十 1 ん十 1 = ルん + 1 ( + 1 十ん + 1 + 1 ) を計算することで導出できる。 十 1 ( ん十 1 十ん十 1 「ん十 1 ) = △十 1 ん十 1 十ん十十 1 「ん十 1 鉄十ん + 1 を 0 ん ( の月 ( ん十 1 ーの 17 0 0 0 津 ( ん + 1 = 0 0 ( ん + 1 ーの 十 1 ん十 1 = 十 1 ( 江ん十 1 十ん十 1 ん十 1 ) とすれば、 こで /\ ん + 1 0 0 ek ー を 0 津 ( ん + 1 ーの一を 0 滝 ( ん + 1 ーの十ん + 1 鉄 ( 1 ーん + 1 ) 鉄 0 0 となって鉄 + 1 を求めることができる上、 M = ん十 1 の時でも再帰が成立することが分 かる。 Levinson-Durbin 再帰アルゴリズム 前節までの導出結果をまとめると、 1. ん = 1 の時 : ( 1 ) の ( 1 ) = 2. んが求まった時、ん十 1 は : /\ ん十 1 e ん十 1 0 ん十 1 町 0 ) el = 月 ( 0 ) + 町 1 ) 佖 1 ( 1 ) = 0 ( の町ん + 1 ーの 社ん十 1 十ん十 1 「ん十 1 ( 1.11 ) ( 1 ・ 12 ) ( 1 ・ 13 ) ( 1 ・ 14 ) ( 1.15 )

104 索引 。符号 , 41 一意に復号可能な符号 , 28 後ろ向き誤差 , 18 MS 処理 , 59 工ンコード , 27 工ントロピー , 34 工ントロピー符号 , 9 工ントロピー符号化 , 27 符号 , 41 格子型フィルター , 20 固定小数 , 52 Golomb 符号 , 43 残差 , 9 自己相関 , 13 情報 , 27 情報源 , 27 情報量 , 31 線形予測 , 12 デコード , 28 PARCOR 係数 , 18 標本自己相関 , 18 復号 , 28 符号 , 27 符号化 , 27 プリエンファシス , 62 前向き誤差 , 18 マクミランの不等式 , 29 窓関数 , 87 無記憶情報源 , 37 Rice 符号 , 43 Levinson-Durbin 再帰法 , 14

102 [ 19 ] [ 20 ] [ 21 ] [ 22 ] [ 23 ] [ 24 ] [ 25 ] [ 26 ] [ 27 ] [ 28 ] [ 29 ] [ 30 ] [ 31 ] [ 32 ] [ 33 ] [ 34 ] [ 35 ] [ 36 ] [ 37 ] [ 39 ] [ 40 ] 参考文献 dvorakjp/hinshutu. htm Linear Prediction and Levinson-Durbin Algorithmhttp : //www.emptyloop ・ com/tech110tes/A\%20tutoria1\%200n\%201inear\%20prediction\%20and\ %20Levinson—Durbin. pdf The Canterbury Corpus, http : //corpus. canterbury. ac. nz/descriptions/ 7 符号、 5 符号、ゴロム符号による圧縮効果 , https://naoya-2.hatenadiary ・ org/e Ⅱ try / 20090804 / 1249380645 Zhu Li, Lec 03 Entropy and Coding Ⅱ Hoffman and Golomb C0ding http ://1.web.umkc.edu/1izhu/teaching/2016sp.video.com/unication/ notes/Iec03. pdf 浮動小数点を利用する際に知っておきたいこと , https : //blogs ・ msdn ・ microsoft.com/jpvsblog/2014/10/28/93/ 固定小数 , http : //www ・ sage-p ・ com/compone/toda/fixdec.htm 第 10 回「固定小数点の演算」 ( 201302 ) , http://www.hdlab ・ co ・ jp/web/ a0600nepoint / 201302. php 7.11 Wavpack, http : //www.wavpack. com/WavPack.pdf WavPack Audio Compression http : / / w .wavpack. com FLAC - Free Lossless Audio Codec, https ://xiph.org/flac/index.html FLAC - format, https : //xiph. org/flac/format . html Lossless AudiO Homepage, http : //www. lossless—audio . com OptimFROG, http : //losslessaudio. org TAK, http : //thbeck. de/Tak/Tak. html TTA, http ://tausoft.org/wiki/True—Audi0—C0dec—Overview MPEG-4 Audio Lossless Coding (ALS) , https : //www.nue.tu-berlin. de/ menue/research/research_topic/compression_and_transmission/mpeg—4— audi0-10ss1ess-c0ding-a1s/parameter/en/#c230252 Monkey's Audio - a fast and powerful lossless audio compressor https : / / . monkeysaudio. com/index.html ビットを数える・探すアルゴリズム , http://www.nminoru ・ jp/-nminoru/ programming/bitcount. html nlz. c. txt, https : //www.hackersdelight.org/hdcodetxt/nlz.c.txt [ 38 ] 数値解析 ( 帝京大学 ) 2016 年度前期浮動小数点数補足 http://y-m ・ jp/wp-content/ up10ads/2016/04/na—2nd—f10atingpoint . pdf Vorbis I specification, https ://xiph.org/vorbis/doc/Vorbis- 工 -spec. html Burg' s Method, Algorithm and Recursion, http : //www.emptyloop.com/ techn0tes/A%20tutoria1%200n%20Burg' s%20meth0d, %20a1gorithm%20and%

24 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 } 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 第 1 章 predict [smpl] fo て (i 1 ; i く = max-delay; i + + ) { predict [smpl] (coff [i] * data[smpl-i] ) ; = 0 . 0f ; / * 誤差計算・結果表示 */ error = 0 . 0f ; for (smpl 0 ; smpl く num-sample ; smpl + + ) { e てて 0 て + = pow(predict [smpl] ー data[smpl] , 2 ) ; printf ( "NO : %d 凵 Data : 凵 %f 凵 Predict : 凵 % f 凵 \ Ⅱ " , smpl , data [smpl] , 基本理論編 predict [ smpl] ) ; printf ( " Er て 0 て凵 : 凵 % f 凵 \ Ⅱ” , sqrt (error / num_sample) ) ; free (data) ; free (predict) ; free(auto-cor) ; free(coff) ; return 0 ; static 土Ⅱ t32 ー t levinson_durbin_recursion(double *lpc_coef , *auto—corr, uint32—t max_order) 土Ⅱ t32 ー t k, 土 ; double lambda ; double *u_vec , *V_vec, *a_vec , *e_vec ; if (lpc-coef = = NULL Ⅱ auto—corr = = NULL) { const double fprintf (stderr , "Data 凵 0て凵 result 凵 pointer 凵 point 凵 to 凵 NULL. 凵 \ Ⅱ” ) ; return ー 1 ; * 0 次自己相関 ( 信号のニ乗和 ) が 0 に近い場合、入力信号は無音と判定 * = > 予測誤差 , p び系数は全て 0 として無音出力システムを予測 . return 0 ; / * 初期化 * / 0 . 0f ; lpc—coef [i] 0 ; i く max—order + 1 ; i + + ) { for (i if (fabs (auto—corr [ 0 ] ) く FLT_EPS 工 LON) { (double *)malloc (sizeof (double) a_vec 0 ーの 0 ーん + プを含めると ma 0 der + 2 * / (double *)malloc(sizeof (double) e—O, e ーん + プを含めると ma 0 de + 2 * / (double *)malloc (sizeof (double) e_vec u_vec * (max- * (max- * (max- 0 て de て order order

2.9 性能向上の指針 Monkey ' s Audi0 #include く math . > 筆者は [ 40 ] を参考に、リスト 2.49 に示すように C 言語版の簡易実装を行った。 き誤差と後ろ向き誤差の二乗和の最小化により定式化される。 べースの手法よりもより精度の高い解析ができる [ 41 ] 事が知られている。 Burg 法は前向 していたが、他にも線形予測手法の別の定式化が存在する。特に Burg 法は自己相関関数 ALA では自己相関を計算することによる Levinson-Durbin 法べースの線形予測を使用 Burg 法 予測についての取り組みを本節で述べる。 2.9.2 予測手法 したソースが公開されている [ 42 ] が、筆者はこの実装をまだ理解しきれていない。 り、マルチプラットフォーム対応ができていない。 ffmpeg がリバースエンジニアリング が公式に配布している SDK には W ⅲ dows 向け exe とライプラリが含まれるだけ * 14 であ ると言える。しかし、 TAK は実装がクローズの上、商用利用を禁止している。更に、 TAK 能なコーデックである。事実上、 2019 年の段階では最強のロスレス音声コーデックであ TAK[32] は Monkey's Audio に並ぶ圧縮率と FLAC 並のデコード速度を持つ非常に高性 TAK (Tom' s verlustfreier Audiokompressor) ンであるものの、筆者としては難解であり読み解くことができなかった。 い。高圧縮率の秘訣はプロックに分けないことで達成していると思われる。実装はオープ ロックに分けずに符号化するため、ストリーミングエンコード / デコードが事実上できな れ、しばしば他コーデックの圧縮率と比較される。 Monkey's Audio は音声データをプ Monkey's Audio[35] はロスレス音声コーデックの中でも圧縮率が高い、 13 ことが知ら 95 1 2 3 4 5 6 7 8 リスト 2.49 Burg 法の C 言語実装 く stdint . 五 > く stdio . h> #include く stdlib . h> く string. h> く float . > く assert . 五> #include #include # inc lude #include #include 、 13 圧縮率の高いロスレス音声コーデックとしては他に La [ 30 ] と OptimFROG[31] が挙げられるが、エンコー ド / デコード速度が異常に遅いため、 こでは解説しない。 wine を使用すれば Linux 環境で動かすことは可能。 * 14

1.1 線形予測 これらの式 1.18 , 1.19 を変形することで、線形予測の別の定式化を行うことを考える。 Levinson-Durbin 法の係数更新式 1.13 に注目すると、äM = äM 十 /\MÜM = は行毎に次の計算をしていることが観察できる。 aM m) = 0M-1(m) ー 7MaM-1(M ー m) ( この式 1.20 を前向き誤差の式 1.18 に代入すると、 ( 1.20 ) 19 社ィーイむ M ん ( れ ) = を 0M(m)y(n ー m) を。 M ー 1 ( m ル ( れ一 m) rn=0 m=0 m=0 m=0 rn=0 rn=0 を M ー 1 ( m ル ( れ一 m) をー 1 ( M ー rn ル ( れ一 m) ( ・ . ・ m → m + 1 ) m=l ・ aM-1(M) = 0 ) ー 7M aM ー 1 ( M ー 1 ー m ル ( 〃ー 1 ー m ) 佖 M ー 1 ( m ル ( れ一 m) -7M aM(M ー m)y(n-m) m=0 = ん一 1 ( れ ) ーー 1 ( れ一 1 ) 示している。後ろ向き誤差についても、前向き誤差のときと同じく、後ろ向き誤差の式 式 1.21 は、 PARCOR 係数 7 。 , ( m = 1 , ... , M ) から、前向き誤差を計算できることを M—I m=0 m=0 rn=0 aM-1(M ー m ル ( れ一 m) 。 M ー 1 ( m ル ( れ一 m) を aM-1(M ー m ル ( れ一 m) ー加 aM-1(m)y(n ー m) を {aM-1(M ー m) ー し M ( れ ) = aM(M ー m ル ( れ一 m) 1.19 に式 1.20 を代入すれば、 m=0 M—I m=0 ( 1.21 ) 。 M ー 1 ( M ー 1 ー m ル ( れ一 1 ー m ) -7M を aM-1(m)y(n-m) ( ・ . ・ m → m + 1 ) = し M ー 1 ( 〃ー 1 ) ーん一 1 ( れ ) rn=0 m=0 すフィルターによって計算できる。このフィルターは、乗算器を交差させた特有の構造に となり、後ろ向き誤差も PARCOR 係数から計算できる。式 1.21 , 1.22 は、図 1.3 に示 ( 1 ・ 22 )