第 2 章実装編 84 wav ファイルの読み込み 2.6.2 公開 API 変換され、 Ox7FFF に丸められて書き込まれる。 の書き込みの際、 0X7FFF0000 と書き込んだデータは 16bit 整数に ( 右 16bit シフトして ) 45 / * ファイルから / ファイルハンドルを作成 * / リスト 2.34 wav ファイル読み込み av. 五 ) wav ファイルの読み込みを行う関数をリスト 2.34 に示す。 49 struct WAVFi1e* WAV—Create (const struct WAVFi1eFormat* format) ; 48 / * フォーマットを指定して新規にルアイルハンドルを作成 * / リスト 2.35 WAVFi1e を波形フォーマット指定で新規作成 (wav. h) ともできる。 wav ファイルの新規作成時に使用する。 なお、 WAVFi1e 構造体は波形フォーマット情報 WAVFi1eFormat だけから新しく作るこ ンタが返り値として得られる ( 読み込みに失敗した場合は NULL が返る ) 。 引数としてファイル名を指定し、読み込みに成功した場合は WAVFi1e 構造体へのポイ 46 struct WAVFi1e* WAV—CreateFromFi1e (const char* filename) ; 55 WAVApiResu1t WAV—WriteToFi1e ( 54 / * ファイル書き出し * / リスト 2.36 av ファイル書き込み ( wav. h ) り値で得られる。 WAVFi1e をファイル名指定で書き込む関数をリスト 2.36 に示す。書き込みの成否は返 av ファイルの書き込み プログラム wavbitchange. c をリスト 2.37 に示す。 本モジュールの簡単な使用例として、 wav ファイルの量子化 bit 幅を変換して書き出す 2.6.3 使用例 56 const char* filename , const struct WAVFi1e* wavfile) ; 1 2 3 4 5 リスト 2.37 wav ファイルの bit 幅変換 (wavbitchange. c) #include "wav . h" #include く stdio . h> #include く stdlib . h> int main (int argc , char** argv)

2.7 main. c ALA バージョン 1.0.0 の段階では、 ALA-FORMAT-VERSION マクロはリスト 2.46 で定義し ALA-FORMAT-VERSION マクロ定数と不一致だった場合は、無条件でデコード失敗とする。 デコードするサンプル数は、リスト 2.47 に示す式で決定している。 デコードするサンカレ数の決定 23 #define ALA-FORMAT—VERS 工 ON 1 22 / * フォーマット / ヾーション * / リスト 2.46 ALA_FORMAT-VERSION マクロ (main. c ている。 91 361 362 リスト 2.47 デコードするサンプル数の決定 ( ma 洫・ c ) / * デコードサンプル数の確定 * / num—decode—samples = ALAUT 工 L 工 TY—MIN (num—block—samples , dec—offset—sample) ; num-samples ヘッダから取得したプロックサンプル数Ⅱー block ー s ples を基に、エンコードと全 く同じ処理でデコードサンプル数を決定する。 2.7.3 引数処理 コード (do-decode 関数 ) 処理を呼び分ける。 main 関数の処理をリスト 2.48 に示す。 main 関数はコマンドラインの引数を処理し、エンコード (do-encode 関数 ) 処理とデ 441 442 443 444 445 446 447 448 449 450 451 452 453 454 455 456 457 458 459 460 リスト 2.48 main 関数 (main ・ c) / * メインエントリ * / int main (int argc , char** argv) const char* option ; const char* input—file ; const char* output—file ; / * 引数が足らない * / if (argc く 4 ) { print—usage (argv) ; return 1 ; / * 引数文字列の取得 * / option = argvC1] ; argv [ 2 ] ; input—file output—file = argv [ 3 ] ; / * 工ンコード / デコード呼び分け * / if (strcmp(option,

56 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 第 2 章 int32—t* residual) 実装編 const int32—t* data, uint32—t num—samples , const int32—t* parcor—coef , uint32—t order, uint32—t samp , 0 て d ; int32—t* forward—residual ; int32_t* backward_residual ; int32—t mul—temp ; / * 丸め誤差軽減のための加算定数 = 0.5 * / (IUL くく 14 ) ; const int32_t half / * 引数チェック * / if (lpc = = NULL Ⅱ data = = NULL Ⅱ parcor—coef = = NULL Ⅱ residual = = NULL) { return ALAPRED 工 CTOR—APIRESULT_INVAL 工 D—ARGUMENT ; / * 次数チェック * / if ( 0 て de て > lpc—>max-order) { return ALAPRED 工 CTOR_APIRESULT-EXCEED_MAX—ORDER; / * オート変数にポインタをコピー * / lpc—>forward—residual ; forward_residual lpc—>backward—residual ; backward_residual / * 誤差計算 * / fo て (samp = 0 ; samp く num—samples ; samp + + ) { / * 格子型フィルタにデータ入力 * / forward—residual [ 0 ] data Csamp] ; / * 前向き誤差計算 * / for (ord = 1 ; 0 て d く = order; 0 て d + + ) { mul—temp return ALAPRED 工 CTOR-APIRESULT_OK ; residual Csamp] forward—residual [ 0 て de て ] ; / * 残差信号 * / backward—residual [ 0 ] = data Csamp] ; / * 後ろ向き誤差計算部にデータ入力 */ backward—residual [ 0 て d ] backward—residual [ord ー 1 ] ー mul—temp ; forward—residual [ 0 て d ー 1 ] + half , 15 ) ; (int32—t) ALAUTILITY—SHIFT_RIGHT_ARITHMETIC (parcor—coef Cord] * mul—temp for ( 0 て d = 0 て de て ; ord > = 1 ; ord——) { / * 後ろ向き誤差計算 * / forward—residual [ 0 て d ] forward—residual [ 0 て d ー 1 ] ー mul—temp ; backward—residual [ 0 て d ー 1 ] + half , 15 ) ; (int32—t) ALAUTIL 工 TY-SHIFT-R 工 GHT—AR 工 THMETIC (parcor-coef [ 0 て d ] *

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 + + ) { / * 予測テスト * /

2.9 性能向上の指針 / * f-vec, とー ec の更新 * / fO て ( 土 = 0 ; 土く num—samples ー k ー 1 ; i + + ) { f-vecCi + k + 1 ] ; tmp2 = b—vec[i] ; tmpl f—vec[i + k + 1 ] tmpl + mu * tmp2; b—vec は ] = mu * tmpl + tmp2; 97 8 一 9 0 1 2 3 4 5 6 5 5 6 6 6 6 6 6 6 のんの更新 * / Dk = ( 1 . 0f ー mu * (u) * Dk ー f—vec[k + 1 ] * f—vec[k + 1 ] b—vec [ k ー 2 ] * b—vec [num—samples ー k ー 2 ] ; num—samples 67 68 69 / * 係数コピー * / sizeof (double) * (coef-order + 1 ) ) ; 70 memcpy(coef , a—vec, 71 free(b—vec) ; 72 free(f—vec) ; 73 free(a—vec) ; 74 free (parcor—coef) ; 75 76 } 77 78 int main(void) 79 32 * 8192 ; const Ⅱ土Ⅱ t32 ー t num—salnples 80 const uint32_t coef—order 256 ; 81 82 uint32-t i , j ; 83 double* data; 84 double* predict ; 85 double* coef ; 86 double err ; 87 88 data = malloc (sizeof (double) * num—samples) ; 89 predict = ma110c (sizeof (double) * num—samples) ; coef = ma110c (sizeof (double) * (coef—order + 1 ) ) ; 90 91 / * データ作成 * / 92 fo て (i = 0 ; i く num—samples; i + + ) { 93 sin(i * 0 .01f ) + cos ( 4.0f * sin(i * 0.05f ) ) ; data[i] 94 95 96 97 / * B 社四アルゴリズム * / 98 burg-algorithm(data, num—samples , coef , coef—order) ; 99 / * 予測 */ 100 0 . 0f ; 101 coef—order ; i く num—samples ; i 十十 ) { for (i 102 predict [ 幻 = 0 . 0f ; 103 1 ; j く = coef-order; j + + ) { for (j 104 predict Ci] coef Cj] * data[i 105 106 err + = pow(data [i] 107 predict [i] , 2 ) ;

2.3 384 385 386 387 388 389 390 391 392 393 394 395 396 397 398 399 400 401 402 403 404 405 406 ala—predictor. c リスト 2.11 プリエンファシスフィルター処理 ( ala ー p て edicto て . c ) / * プリエンファシス ( t32 , マ ce ) */ ALAPredictorApiResu1t ALAEmphasisFi1ter—PreEmphasisInt32 ( int32—t* data, uint32—t num—samples , int32—t coef—shift) uint32—t smpl ; int32-t prev—int32, tmp-int32; const int32—t coef—numer = (int32-t) ( ( 1 くく coef-shift) / * 引数チェック * / if (data = = NULL) { return ALAPRED 工 CTOR—AP 工 RESULT—INVAL 工 D-ARGUMENT ; / * フィルタ適用 * / prev—int32 0 ; smpl く num-samples ; smpl + + ) { for (smpl tmp—int32 = data Csmp1] ; (int32—t) ALAUTIL 工 TY-SHIFT—RIGHT—AR 工 T}NET 工 C (prev—int32 * data [smpl] coef—numer, coef—shift) ; prev—int32 tmp—int32 ; return ALAPRED 工 CTOR—APIRESULT—OK ; 63 入力データ data をその場で ( ⅲ - place で ) 書き換える関数である。フィルター計算は固 定小数乗算を使用しており、式 2.7 の″に該当する係数との乗算を、 ( 1 くく coef-shift) 1 との乗算処理と、 coef-shift の右シフト処理に分けて実現している。 ALA バージョ ン 1.0.0 では、 coef-shift は 5 に設定されているため、フィルター係数″の値は ″ = ( 25 ー 1 ) / 25 = 31 / 32 = 0.96875 に相当する。 デ工ンファシスフィルター処理 デ工ンファシスフィルター処理の実装をリスト 2.12 に示す。 408 409 410 411 412 413 414 415 416 417 418 419 リスト 2.12 デ工ンファシスフィルター処理 (ala_predictor. c) / * デ工ンファシス ( をれ t32 , をれマ乙 ace ) * / ALAPredictorApiResu1t ALAEmphasisFilter—DeEmphasisInt32 ( int32—t* data, uint32—t num—samples , int32—t coef—shift) uint32—t smpl ; const int32-t coef—numer = (int32-t) ( ( 1 くく coef-shift) / * 引数チェック */ if (data = = NULL) { return ALAPRED 工 CTOR_AP 工 RESULT_ 工 NVAL 工 D—ARGUMENT ;

2.5 ala_coder . c 77 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 } / * 商部分を取得 * / quot BitStream—GetBit (strm, &bit) ; while (bit quot 十十 ; BitStream—GetBit (strm, &bit) ; } else { rest / * 1 の剰余は 0 */ if (rice—parameter / * 剰余部分を取得 * / return (rice—parameter * quot + rest) ; rest (uint32—t) bitsbuf ; BitStream—GetBits (strm, ALAUti1ity—Log2Cei1(rice—parameter) , uint64_t bitsbuf ; / * ライス符号の剰余部分取得 * / &bitsbuf 符号化処理 ( リスト 2.19 ) の逆の操作を行っているだけである。まず商 quot を住符号 の定義 ( 1 が出現するまでの bit 数 ) に従って取得し、次に Rice 符号のパラメータを元に 剰余 rest を取得する。 2.5.2 残差の符号化と復号 残差の符号化と復号処理は、 Rice 符号のパラメータを入力データに適応しながら変更 し、残差を出力、あるいは取得する処理に要約される。 Rice 符号のパラメータの更新処理 は式 1.54 で見たように平均値の更新処理に帰着される。平均値は有限のサンプル数では 厳密に計算できないので、替わりに推定平均値の計算を行う。符号化と復号処理で全く同 じ様に推定平均値を更新することで、パラメータを適応的に変更しながらでも、可逆な符 号化が実現できる。 残差の符号化 残差の符号化処理をリスト 2.22 に示す。 110 111 112 113 リスト 2.22 残差の符号化処理 (ala-coder ・ c) / * 符号付き整数配列の符号化 * / ALACoderApiResu1t ALACoder-PutDataArray ( struct ALACoder* coder, struct BitStream* strm, const 土Ⅱ t32 ー t * * data, uint32_t num—channels , uint32_t num—samples)

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

96 9 10 11 uint32 t i, k; 12 13 14 15 16 double Dk , mu ; 17 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 44 45 46 47 48 49 50 51 52 53 54 55 56 57 double* coef , uint32—t coef—order) static void burg—algorithm (const double* data , double double double double double tmpl , tmp2 ; *parcor_coef ; *b_vec ; *f_vec; *a_vec ; / * べクトル領域割り当て * / 第 2 章実装編 uint32-t num-samples , a—vec = malloc (sizeof (double) * (coef—order + 1 ) ) ; parcor—coef = malloc(sizeof (double) * (coef—order + 1 ) ) ; f—vec = malloc (sizeof (double) * num—samples) ; b—vec = malloc (sizeof (double) * num—samples) ; / * 各べクトル初期化 */ fO て (k = 0 ; k く coef—order + 1 ; k + + ) { a-vec [k] = 0 . 0f ; a-vec [ 0 ] 1 . 0f ; memcpy(f-vec , data, memcpy (b-vec , data , / * りんの初期化 * / Dk = 0 .0f ; sizeof (double) * num—samples) ; sizeof (double) * num—samples) ; for (i 0 ; i く num-samples; i + + ) { Dk + = 2.0f * f—vec[i] * f—vec[i] ; Dk ー f—vec [ 0 ] * f—vec [ 0 ] + b—vec [num—samples 1 ] * b—vec [num-samples / * 四再帰アルゴリズム * / fO て (k = 0 ; k く coef—order; k + + ) { / * 反射 ( 月 C ) 係数の計算 * / mu = 0 . 0f ; for (i 0 ; i く num—samples ー k ー 1 ; mu + = f—vec[i + k + 1 ] * b—vec[i] ; ー 2 . 0f / Dk ; mu * = parcor—coef [k] assert(fabs(parcor-coef[k]) く 1 . Of ) ; / * 。ー ec の更新 * / for (i tmpl = a—vec[i] ; tmp2 = a-vec[k + 1 a-vec Ci] tmpl + mu * tmp2 ; a—vec Ck + 1 i] = mu * tmpl + tmp2; 幻 ;

1.1 線形予測 25 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 V_vec for (i (double *)malloc (sizeof (double) * (max—order + 2 ) ) ; 0 ; i く max—order + 2 ; i + + ) { u-vec Ci] / * 最初のステップの係数をセット * / v-vec Ci] a-vec [i] = 0 . 0f ; a-vec [ 0 ] e-vec [ 0 ] a-vecC1] e-vec [ 1 ] u-vec [ 0 ] v-vec [ 0 ] 1 . 0f ; auto—corr [ 0 ] ; auto-corr [ 1 ] / auto-corr [ 0 ] ; ー co てて [ 0 ] + auto-corr[l] * a-vec[l] ; = 0 . 0f ; 1 . 0f ; auto u-vec[l] = 0 . 0f ; v-vec[l] 1 . 0f ; for (k = 1 ; k く max—order; k + + ) { 0 ; i く k + 1 ; i + + ) { lambda + = a_vec [i] * auto—corr [k + 1 / * 再帰処理 * / lambda = 0 .0f ; fo て ( 土 lambda / = (—e e—vec [k + 1 ] ( 1 ー lambda * lambda) * e-vec [k] ; / * u-vec, セ ec の更新 */ 0 ; i く k ; i + + ) { for ( 土 u—vec [i + 1 ] = v—vec [k u-vec [ 0 ] 1 .0f ; u—vec[k + v-vec [ 0 ] 0 . 0f ; v—vec[k + 1 ] / * es t の更新 * / 1 ] 0 ; i く k + 2 ; i + + ) { for (i a-vecCi + 1 ] ; = 0 . Of ; = 1 . Of ; / * 結果の取得 * / memcpy(lpc-coef , a—vec , sizeof (double) * (max-order + 1 ) ) ; a-vec [i] = u—vec [i] + lambda * v—vec [i] ; free (u—vec) ; free (v—vec) ; free(a—vec) ; free(e—vec) ; return 0 ; static int32—t calc—auto—correlation(double *auto—corr, data, uint32—t num—sample , uint32—t max—order) int32—t smpl , i ; if (max—order > num-sample) { const double * fprintf (stderr , "Max 凵 order (%d) 凵 is 凵 larger 凵 than 凵 number 凵 0f 凵 samples (%d) ・凵